We designed a deformable model of the human CNS to reproduce functional biomechanical relations between dynamically interacting CSF compartments (Figures 1A, B). An anatomically accurate analog of the spinal SASs with the transparent spinal cord including pairs of peripheral nerve roots and the translucent dural surfaces was manufactured in a multistep 3D printing and casting process with subject-specific imaging data of a 26-year-old male volunteer (Tangen et al., 2015). A mold was 3D printed using dissolvable polyvinyl alcohol (PVA). Platinum silicone casting resin of shore A8 hardness (TAP Plastics Inc.) was used for casting a flexible and transparent phantom. All parts including the dural and pial surfaces and nerve roots are deformable. Spinal CSF motion was generated by transmitting oscillatory expansion and contraction of an inflatable balloon located in the head section to the fluid. The balloon, mimicking the cerebrocranial vascular bed, in turn, was driven by a piston pump capable of generating stroke volumes up to 1 mL/beat in the frequency range of 0–127 beats per minute, as shown in detail in Supplementary Appendix SA. The connected spinal and cranial CSF spaces had no artificial openings consistent with the closed nature of the natural CNS. The mode of pulsatile flow conditions in the spinal CSF-filled spaces of the bench model reproduces pulsatile vascular bed dilation consistent with our understanding of periodic intracranial CSF displacement (Basati et al., 2012; Tangen et al., 2015). More manufacturing details of the subject-specific CNS replica (including 3D printing settings, stitching, and patching) can be found elsewhere (Tangen et al., 2016b; Lu et al., 2016; Ayansiji, 2023).

To reproduce conditions of in vivo IT procedures on the bench, infusion catheters with inner diameters of 0.2 mm, 1.0 mm, and 3.2 mm were inserted into lumbar and thoracic regions with inner diameters of 0.2 mm, 1.0 mm, and 3.2 mm. The elastic dura of the spine model has a self-sealing property, thus enabling realistic catheter insertion and placement as used in clinical practice for human therapies. A wide range of infusion parameter settings (i.e., infusion volume, flow rate, position, and duration) and systemic blood pump settings (i.e., stoke volume and frequency) enabled the implementation of a comprehensive spectrum of infusion scenarios (bolus and chronic drug pump) occurring in the physiological range of CSF oscillations. More than 50 tracer infusion experiments of trypan blue (Sigma-Aldrich) released with a programmable syringe pump (Harvard Instruments) were performed to precisely investigate the correlation between dye dispersion, CSF pulsations, and infusion parameters. Trypan blue was chosen for its intense blue color needed in optical front tracking. In a clinical setting, acute high-volume injection via a catheter is administered in a lying down (horizontal) position for safety. Thus, the horizontal patient position was chosen for our experiments.

The method of moments (MoM) was used to determine the apparent tracer dispersion velocity from video data because it minimizes the sensitivity of optically acquired concentration profiles to uneven lighting conditions, scattering effects, and uncertainty in concentration inference from intensity data (Birk et al., 2005). The MoM has been previously used to quantify dispersion in an annular tube (Mondal and Mazumder, 2005; Nagarani and sebastian, 2013; Debnath et al., 2019). Our experiment was insensitive to variations in lighting conditions, uncertainty of intensity quantification, and cardiac cycle synchronization; thus, we did not need gated acquisition. This is due to the advantage that MoM only tracks the slope of the evolution of moments over relative prolonged periods of time (=mean and variance of the entire intensity profile), so it is insensitive to perturbations of absolute concentrations or experimental variability. We first determined the speed of caudocranial motion of the infusate by tracking shifts in the center of gravity of the tracer profile. The first moment m 1 t of the observed concentration profiles at time t is calculated as follows: m 1 t = ∫ x o x m C x , t   x   d x ∫ x o x m C x , t   d x = x ¯ t , (1) where x is the position along the neuroaxis and x 0 marks the infusion point. The extreme limit of the ROI equal to full length from caudal to cranial aspects measured x m = 48 cm. In each time frame, t, the first moment, m 1 t , gives the location of the center of gravity, x ¯ t , of the tracer concentration profile at a time point, C x , t . The speed of caudocranial motion is equal to the change of the first moment with time.

We further computed the second moment or area moment of inertia for each time point as as follows: m 2 t = ∫ x o x m x t − x ¯ t 2   C x , t   d x ∫ x o x m C x , t   d x = σ t 2 . (2)

The second moment can be interpreted as the mean spread of the visible concentration profile around its center, x ¯ t . The increase in its variance with time, σ t 2 , is directly proportional to the apparent dispersion or diffusion rate, D e x p , of the tracer molecule as follows: D e x p = 1 2   ∆ σ t 2 ∆ t . (3)

The coefficient of the apparent dispersion can then be readily determined as the rate of change in the second moment of the concentration curves by plotting their variance as a function of time. Supplementary Appendix SD shows details.

In MR imaging, it is convenient to characterize natural CSF pulsations via the instantaneous total CSF volume, V(t), cervical stroke volume, v c , and CSF angular pulse frequency, ω . We estimated the average CSF flow velocity, U_rms, as follows: U r m s = ∫ 0 t V t / A 2 T d t   V t = V 0 + v c 2 − v c 2 cos ω   t , A = V L , ω = 2 π   N b p m . (4)

The average CSF flow velocity, U_rms, correlates a clinical MR quantity to parameters used in our formal flow analysis. Thus, two stroke volume settings of 0.5 mL/beat to 1 mL/beat in the frequency range of 0–120 beat/min enabled us to cover a wide range of CSF flow velocity averages covering the entire physiological range (0.09–0.93 cm/s).

The mean CSF flow velocity, U r m s , in the cervical region was computed by integrating and averaging the squared ratio of pulsatile volumetric flow rate, d V t d t , divided by the hydraulic cross-sectional area of the spinal CSF SAS, A, under cosine profile assumption as in Eq. 4, where V₀ is the initial volume of CSF in cm³ in the system and f is the frequency in beats per minute, with period T in seconds. U r m s = ∫ 0 t V t / A 2 T d t   V t = V 0 + v c 2 − v c 2 cos ω   t , A = V L , ω = 2 π f .

Initially, the experiments were repeated three times with the same settings. The repeats essentially gave the same moments, thus confirming the reproducibility and robustness of data acquisition. The results of the Shapiro–Wilk test indicated that all datasets had a normal distribution. The ANOVA test was used in the regression analyses using IBM SPSS v26 to confirm that regression models described the dependent variables. Confidence intervals were calculated for all regression results. The Durbin–Watson statistic for all regression results was 0.716; therefore, the datasets have autocorrelation. Furthermore, calculated tolerances indicated the multicollinearity effect in the datasets. A p-value of 0.05 was considered statistically significant.

Figure 2A shows tracer concentration profiles inferred from calibrated RGB pixel intensity as described in Materials and methods. Figure 2B shows the first and second moment for the selected experiment. The first and second moments for all the experiments are given in Supplementary Appendices SC, SD. The close agreement between three repetitions (=variance of the second moment ± 0 .29 cm²/min) confirms the reproducibility and acceptably low experimental variability of the experimental setting. The first moment gives information about the caudocranial velocity, and the second moment determines the dispersion coefficient of the tracer as shown in Eq. 3. The binary mask was used to crop images to the width of the spinal section covered by detectable tracer intensity. The size of this mask along the neuraxis (x-direction) is referred to as the dispersion width of the tracer at time t, DD(t). The next section shows the systematic analysis of infusion and physiological parameters on the speed of tracer dispersion.

We studied the effect of catheter diameter on the speed and size of the dispersion front. Infusion experiments are performed using three different catheters with inner diameters of 0.2 mm (needle N1), 1.0 mm (N2), and 3.2 mm (the widest needle N3). Infusion lasted 1 min over the course of 54 experiments. The dispersion width after 10 min, DD₁₀, was measured from the point of the needle tip to the tip of the dye front in the caudal and cranial directions. The duration of 10 min was chosen because this initial time window is critically important for assessing acute risks associated with high-volume IT injection. High local toxicity has been implicated with granuloma formation (Kaye et al., 2004; Allen et al., 2006). We also varied infusion volumetric flow rates (IVF = 0.5, 1.0, and 2.0 mL/min).

Figure 3A shows that high-caliber catheters promote slower tracer spread during infusate injection (1 min). We observed the formation of an injection jet, which was less pronounced in high-caliber catheters leading to shorter dispersion width after 1 min, DD₁. For identical infusion flow rate (IVFs from 0.5, 1.0, and 2.0 mL/min), the large caliber needle (N3) generated the shortest dispersion width. Figure 3B shows the effect of inner catheter diameter on the extent of the tracer spread observed 10 min after the infusion as a function of infusion flow. All experimental results (N = 54) were also fitted into a linear regression model that can be used to estimate the initial neuraxial coverage (=size of the region dosed by solute) as a function of catheter lumen and infusion volumes (IVF). The regression model D D   10   min = a N D + b I V F + c (5) can serve as a guideline in the clinical practice to estimate the length of initial dispersion as a function of catheter lumen and IVF (Figure 3B), where N D is the inner needle diameter and a , b , a n d   c are −1.4150, 2.8786, and 19.8668, respectively, with an R² value of 0.7916.

The realization of a fixed infusion rate with thinner catheters requires higher infusion pressure resulting in higher exit velocities ( V E x t ) at their tips, leading to the creation of an injection jet in the close vicinity of the catheter tip. Thus, catheter N1 at IVF = 2.0 mL/min has the largest V E x t and kinetic energy ( K E ) (as shown in Supplementary Appendix SA). Also, the exit velocity and kinetic energy of catheter N1 are higher than those of needles 2 and 3 with the same IVF. In thin catheters, a higher infusion impulse is delivered during the injection phase that translates into a wider dispersion length in the observed tracer front for the same injection flow rate.

Experimental settings of the injection flow rate, injection volume, and catheter specifications were varied to explore optimal conditions for targeting the cervical section or brain area. Fifty-four experiments (N = 54) were conducted to systematically characterize tracer targeting toward the cranial compartment as a function of infusion parameters (phase-1, shown in Figure 4A). The results in Figure 4B show that higher infusion flow rates IVF accelerate the speed of the tracer front advancing in the cranial direction from the infusion catheter tip. The distance the tracer front moves in the cranial direction from the infusion catheter tip for the CNS model in the supine position as observed 10 min after the infusion is represented by S_DD. This effect at high-volume injections is due to the increased insertion kinetic energy (Supplementary Appendix SA).

The speed of caudocranial motion, CCV, is calculated as the change of the first moment with respect to time as shown in Figure 2B. All experiments showed a shift of the first moment toward the cranium. We determined the velocity of apparent caudocranial advancement of the tracer front, C C V , by recoding its positions for different time points. The results of Figure 4C show that the caudocranial advancement of the tracer front from the injection site toward the cranium slightly increased with oscillation frequency, f , as indicated by the regression C C V f   c m / min = a   f + b (6) and U r m s , as root-mean-square velocity, as indicated by the regression C C V U r m s   c m / min = a   U r m s + b (7) with constant terms a  and  b as 0.3014 and 0.0012 for Eq. 6 and 0.3752 and 0.0012 for Eq. 7, respectively. At 1/mL stroke volume, this trend was similar but with a wider variability between runs (Figure 4D).

The rapid expansion of the drug front follows the amount of infused tracer qualitatively, as expected, because the main driver of the initial tracer spread is the injection impulse of fresh infusate concentrated in the relative narrow spatial confinement in the lumbar injection zone.

Once infusion stops (t > 1 min), further tracer spread is no longer propelled by injection impulse, but by natural oscillations of the CSF. In our closed, deformable model, induced CSF flow in the spinal SAS is also oscillatory with zero net flux. This is also approximately the situation in vivo since net flows due to bulk CSF production are much smaller than those due to oscillations. The apparent dispersion coefficient of tracer spread in the oscillatory CSF flow was determined experimentally as a function of amplitudes (cervical CSF stroke volume and mean CSF flow velocities) and frequencies in a series of dynamic tracer infusion experiments. Two stroke volume settings of 0.5–1 mL/beat in the frequency range of 0–120 beat/min enabled us to induce a wide range of CSF flow velocity averages covering a wide physiological range (U_rms = 0.09–0.93 cm/s). The root-mean-square velocity of the CSF, U r m s , in the cervical region was computed as described in Materials and methods with plots as the function of frequency for a stroke volume of 0.5 mL/stroke and 1.0 mL/stroke in Figures 5A, B. Tracer dispersion reached the cervical region in less than 10 min, and it spread quickly throughout the spinal SAS. Figures 5C, D show dispersion coefficients   D exp determined by MoM (see Materials and methods) as a function of frequency for stroke volume of 0.5 mL/stroke and 1.0 mL/stroke experiments, respectively. Experimentally obtained (=MoM) dispersion coefficients, D_exp, summarized in Figure 5 confirm our previous finding (Hsu et al., 2012; Hsu and Linninger, 2013) that CSF amplitude and frequency are critical factors of IT drug dispersion with the range of D_exp = 1.50–2.94 cm²/min for a 0.5 mL/beat stroke volume, and D_exp = 2.57–4.70 cm²/min for the 1 mL/beat stroke volume as shown in Table E2 and E3, respectively, of Supplementary Appendix SE. As the frequency increases, trypan blue disperses faster and ascends in cranial direction.

To derive functional correlations, Figure 5E lists experimental dispersion coefficient, D_exp, as a function of natural CSF oscillations in terms of amplitude, U r m s , and frequency, f , on a double logarithmic scale according to Eq. 8. The trend in double logarithmic scale, log ⁡ ⁡ ( D exp ) ∼ κ log (U_rms), had a constant slope of κ = 0.7419. A quadratic frequency dependence was also incorporated as an f-dependent offset, with λ f = a 0 + a 1   f + a 2 f 2 with a 0 = − 1.0386 , a 1 = 0.0055 , and  a 2 = − 5.8858   x   10 − 5 , and drawn parametrically in Figure 5E. Taking advantage of the logarithmic form of Eq. 8, all parameters could be estimated simultaneously by linear regression with the coefficient of determination for the best fit of R² = 0.9273. log ⁡ ⁡ ( D exp ) = λ f + κ   log ⁡ ⁡ U r m s . (8)

We also cast all experimental results into a dimensionless form for comparison to prior transport theories. We depart from the commonly known Taylor theory (Aris, 1960) and Aris (Salerno et al., 2020), which describes the increase of the effective diffusivity over molecular diffusivity, D₀, as a function of the dimensionless Peclet number, P e = U r . m . s * h D 0 , as in Eq. 9. Taylor dispersion in stationary flow has λ = 1 / 48 and κ = 2 . Hoagland and Prud’Homme (1985) modified Taylor’s approach for laminar oscillatory flow within rigid cylindrical walls citing the coefficients, λ = 0.0104 and κ = 1.88 . D e x p D 0 = 1 + λ   P e κ . (9)

Figure 5E shows a slope in our experiments ( κ = 0.7419), which does not agree with the dispersion models by Taylor (1953) and Hoagland and Prud’Homme (1985). This discrepancy can be illustrated more clearly with the help of Eq. 10, which was obtained by rearranging Eq. 9. Accordingly, dimensionless effective dispersion increases, ∆ D = D exp − D 0 D 0 , were plotted in Figure 5F based on the raw dimensional data of Figure 5E. Again, the logarithmic correlations permit linear fitting for quantifying the dependence of effective dispersion increase, log ∆ D , on the dimensionless oscillatory flow amplitude (Peclet number) and dimensionless frequency (Womersley number, α ). Here, the Womersley number is α = D H 2 ω υ , D H = 0.5 c m is the hydraulic diameter, ω = 2 π f T , f is the frequency, T = 1   sec is the period of oscillation, and υ = 7 x 10 − 7 m²/sec is the kinematic viscosity. All values are listed in Supplementary Appendix SE. log ⁡ ⁡ ∆ D = λ α + K   log ⁡ ⁡ P e . (10)

The best-fit parameter set ( Κ = 0.74, φ 0 = 2.4786 , φ 1 = 0.3130 , and  φ 2 = − 0.0229 ) in Eq. 11 has a coefficient of determination of R² = 0.9278. log ⁡ ⁡ ∆ D = φ 0 + φ 1   α + φ 2 α 2 + K   log ⁡ ⁡ Pe . (11)

We compared our experimental dispersion coefficient trends to prior work, including the work of Hoagland and Prud’Homme (1985) and Dutta and Leighton (2002). Figure 6A shows that several versions of Taylor–Aris dispersion (TAD) grossly overestimate the speed of biodispersion in the human spinal SAS. Figure 6B shows the zoomed view of the experimental portion of Figure 6A for better visibility. The Taylor (1953), Hoagland and Prud’Homme (1985), and Dutta and Leighton (2002) versions do not match our experiments even qualitatively because the slopes of the Pe dependence ( κ = 1.88 − 2.0 ) are quite different ( κ = 0.7419, Κ = 0.74) in addition to the estimated value ranges being largely off.

We also compared our experimental values to theoretical dispersion coefficient estimates by Stockman (2007), who used Lattice Boltzmann simulations. Stockman used the maximum average speed in a period of oscillation of 220.2 cm/min (in lattice unit of 1 × 10⁻² Lu/ts), and the molecular diffusion of the tracer used was D₀ = 0.0024 cm²/min. The range of the determined dispersion coefficients was D = 0.6356–1.7735 cm²/min, which is much smaller than that found in our experiment by a factor of over 2.5. The Stockman data are also plotted in Supplementary Appendix SF.

In summary, Stockman’s (2007) results underestimate the speed of tracer dispersion, and the results of Hoagland and Prud’Homme (1985), Dutta and Leighton (2002), and Taylor (1953) overestimate the dispersion coefficient by three orders of magnitude or more compared to our in vitro data.

We previously observed that annular phantoms without microanatomical features underestimate the actual dispersion after IT (Barton, 1983; Banyte et al., 2013). To test the significance of microanatomical features on tracer dispersion, we also fabricated a spinal model without nerve roots and compared tracer propagation during infusion to the more anatomically realistic model with nerve roots. Tracer dispersion in the model with nerve roots was found to be always much more rapid than in a system lacking nerve roots under the same condition of oscillatory flow (Figure 3B). Complex flow and observable mixing patterns are absent in idealized annular models lacking spinal microanatomy (Aris, 1960; Tsangaris and Athanassiadis, 1985; Salerno et al., 2020), thus failing to boost the effective dispersion of IT-injected tracers as this is the case when microanatomical features are present. The results of this study provide further evidence for the significant impact of microanatomical features on the spatial and temporal dispersion patterns of IT-administered solutes shown previously (Tangen et al., 2015; Haga et al., 2017).

Experimental data for tracer infusion experiments served as the input for a reduced-order pharmacokinetic model of IT administration. The full description of this mechanistic drug administration model is beyond the scope of this manuscript but can be found elsewhere (Tangen et al., 2017a; Linninger et al., 2023). In brief, tracer biodispersion after lumbar intrathecal injection was simulated via a convection–diffusion process distributed along the neuraxis in Eq. 12. The effect of geometry-induced mixing due to natural CSF pulsations was captured via the effective diffusivity, D e f f , determined in experiments described in Materials and methods and Results. Sensitivity of predictions to infusion settings was incorporated via a single source term, V inj x , t . The effect of injection flow rate, volume, and duration on the forces between infusate, spinal CSF interacting with deformable subarachnoid spaces (=dura) was accounted for by biomechanical fluid structure interaction (FSI). Accordingly, high-volume injections generate a non-zero caudocranial convection term, U f s i , during infusion (phase-1). After the infusion stops (phase-2), the tracer continues to spread away from the lumbar injection site in accordance with effective dispersion as shown in Eq. 12.

Tracer concentration profiles along the neuraxis as a function of time were predicted with mechanistic pharmacokinetic simulations. Pharmacokinetic tracer profile simulations with CNS dimensions and infusion settings used in the experiments took less than one CPU minute to converge, generating asymmetrical profiles with peak concentration decreases with time (Figure 7A). The comparison with the experimental run shows a qualitative match both in the spatial and the temporal dimension as shown in Figure 7A. Figure 7B shows the maximum concentration is attained in a short interval of mean path length. The preliminary results of the proposed order model show that the use of experimentally obtained effective dispersion coefficients can effectively predict drug dispersion after IT administration using a reduced-order pharmacokinetic simulation at low computational cost. It is worth noting that the simulation of biodistribution of active drugs into the CNS and the systemic circulation required additional information on biochemical parameters, denoted by the sink term R c , x , in the following equation: ∂ C x , t ∂ t = ∇ →   D e f f ∇ → C x , t − U f s i x , t , V inj   ∇ → C x , t + n inj x , t − R c , x , (12) where C is the concentration, D e f f is the effective diffusion coefficient, t is the time, x is the position along the neuroaxis, V inj is the volume of injection, U f s i is the caudocranial convection, and n inj is the number of moles. The reaction rate expression might incorporate drug half-life, tissue uptake, and clearance (Tangen et al., 2015; Linninger et al., 2023).

IT infusion studies under physiological conditions require deformability of the CSF spaces to accommodate realistic pulse pressure propagation and fluid motion inside a closed spinal SAS. Deformable fluid-filled spaces in the present bench test analog of the human CNS is fully enclosed between the soft dura/parenchyma surfaces and a distensible vascular interface (cranium vault with a distensible vascular balloon) without fluid exchange across the open system boundaries. The proposed configuration approximates the anatomy and fluid–structure interaction dynamics of the spinal SAS, so we have confidence that it reproduces the complex geometry-induced CSF mixing patterns that were reported in the previous work using the direct numerical simulation of spinal CSF flow (Tangen et al., 2015). The CNS model also incorporates microanatomical features of the spinal cord, epidural space, and peripheral nerve root bundles, especially, which are critical geometric aspects implicated with enhanced mixing (Yiallourou et al., 2012; Pahlavian et al., 2014) of fluid layers in the spinal SAS, leading to accelerated drug dispersion.

Moreover, transparent borders of the see-through human CNS replica enabled dynamic optical tracking of tracer concentration profiles during and after IT solute administration. Our experiments characterized in detail two stages occurring during high-volume IT infusion. The initial phase during which the drug is injected lasts only a few minutes in clinical settings. The duration of 10 min was chosen because it is suitable for assessing acute risks associated with high-volume IT administration (local toxicity and granuloma formation).

The axial dispersion during the injection phase (phase-1) correlated with infusion volume, IVF, and catheter diameters. Thinner catheters (inner diameter d = 0.2 mm) generate wider and fast initial solute dispersion as is attainable with high-caliber catheters (di = 3.2 mm inner diameter). Based on our series of experiments, a simple formula in Eq. 5 predicts the initial dispersion width, DD₁₀, the distribution length after 10 min of injection, as a function of injection volume/flow rate and catheter diameter. The formula in Eq. 5 may serve to estimate the initial volume of distribution, peak concentrations, and initial neuronal tissue exposure during the acute infusion phase of high-volume drug administration as a function of catheter lumen and infusion flow rate. The peak local toxicity risk may be elevated in infusion protocols generating narrower initial spread (i.e., high-volume infusion with large caliber catheters). It can also be used to estimate the expected local volume of drug action and local drug concentration to assess the risk of granuloma formations (Allen et al., 2006).

Drug pumps are often deployed for the local treatment of chronic pain or spasticity but are less suitable for targeting the brain, which was one of the main objectives of the current study. Drug pumps operate under much lower flow rate settings than possible with acute administration. Accordingly, the effect of the injection impulse may be negligible. Chronic IT administration with drug pumps are, thus, predominately governed by conditions of natural oscillation (phase-2) for its entire time course as discussed as follows. For more recent developments for effective infusion protocols from drug pumps, we refer to the work by Yaksh et al. (2012) and Hildebrand et al. (2019).

After the injection ceases, the tracer spreads due to natural CSF pulsations. Oscillatory fluid flow around microanatomical features create a geometry-induced mixing pattern, which breaks the laminar flow field by introducing eddies and vortices around nerve roots. Localized and interspersed eddies were observed in the flow directly upstream and downstream of the cylindrical peripheral nerve root bundles suspended in the flow. Moreover, trabeculae can substantially enhance this effect as reported in Eq. 12.

In all infusion experiments, slow but steady caudocranial advancement of the tracer front from the lumbar injection site toward the cranium was observed. We offer two explanations for the experimentally observed caudocranial transport.

First, lumbar injection divides the space available for the solute to disperse into a smaller distal volume containing the sacral compartment and a larger cranially facing domain stretching from the catheter tip to the thoracic, the cervical, and the cranial SAS. We observed that tracers initially spread equally in both sections, which is consistent with a diffusive process. Subsequently, the advancing dye fronts fill the closed sacral domain faster due to its smaller size. Once the closed sacral region is occupied, the tracer profile center of gravity begins shifting toward the head. The asymmetric caudocranial tracer profile develops due to a boundary effect that was confirmed with a mechanistic diffusive transport model. Systematic experiments with varying CSF conditions (amplitude and frequency) enabled the determination of a formula for the velocity of the caudocranial shift as a function of CSF pulsations and frequency as shown in Eq. 6.

The graded biomechanical deformation profile of the spinal SAS responsible for CSF pulse attenuation is a secondary factor. The CSF pulse amplitude was shown to diminish gradually from the cervical toward the sacral region (Tangen et al., 2015). The CSF pulsations and average oscillatory CSF velocities, U_rms, become larger toward the cervical region compared to the lumbar region where the pulse amplitude is almost zero. Thus, effective dispersion tends to become faster in the caudocranial direction. These two effects cannot be studied in rigid, open systems.

The infusion experiments were conducted in a closed deformable CNS model with no net CSF generation or removal. This supports the notion that bulk CSF flow or absorption is not necessary for caudocranial drug dispersion to occur. Rather, experiments suggest that caudocranial transport of infused solutes can simply result from the asymmetry of the CSF spaces and graded CSF pulse amplitudes.

Recently, an interesting analytical solution for drug dispersion in the SAS was developed by Lawrence et al. (2019). Dispersion predicted inside an idealized annular geometry also found solute transport controlled by convection but not Taylor dispersion.

Tracer dispersion after high-volume injection was very rapid, reaching the cervical region in less than 10 minutes, and it spread quickly throughout the spinal SAS. The apparent dispersion coefficient was robustly determined experimentally as a function of CSF amplitudes and frequencies. An empirical correlation (Eq. 8) between apparent diffusion coefficient and CSF pulsations, a function of CSF amplitude and oscillation, was established. This formula can be of clinical interest to predict tracer dispersion for IT drug administration.

We provide a simple guideline for estimating the volume of distribution of the drug during the injection phase as a function of catheter caliber and injection volume based on the experimental data and model in Eq. 5. Table 1 summarizes the expected size of the injection front (dispersion length after 10 min) from the lumbar injection site. It can also be used to get an idea about the advancement from the injection site toward the cranium, since the moving front of the tracer profile advances at least half of DD₁₀.

For tracer dispersion in the oscillatory CSF flow (phase-2), Eq. 8 quantifies effective dispersion as a function of stroke volume and pulse frequency. Table 2 derived from these data enables the estimation of the effective dispersion coefficient based on physiological parameters. For drug molecules with different molecular diffusion coefficients (i.e., drugs with substantially different properties of our tracer), it may be used as a first approximation when no data are available or when its Pe number is in the same range as for trypan blue with a diffusion coefficient of D₀ = 1.938 × 10⁻⁶ cm²/min.

The analytical Taylor–Aris approach was shown to grossly overestimate the bio dispersion process. All versions of the TAD models showed that the spread of the tracer in human spine does not follow the trend proposed by Taylor dispersion. Effective dispersion by geometry-induced mixing seems to follow a different trend than that described by TAD. Moreover, TAD does not consider the effect of injection volume, catheter geometry, and placement as discussed during phase-1 of the infusion experiments. Also, TAD assumes a constant flow velocity throughout the channel. Finally, poor predictions by TAD theories for IT tracer dispersion and the very large Peclet number ranges seem to suggest that the molecular diffusivity (D₀ in the denominator of Pe) may not be an ideal dimensionless scaling characteristic.

Several authors have implicated the eccentricity of idealized cross-sectional areas of CSF-filled spaces in the spinal SAS as a key factor for accelerated drug dispersion (Loth et al., 2001; Sánchez et al., 2018; Chu et al., 2019; Moral-Pulido et al., 2023). The insensitivity to centric or eccentric alignment in our experiments does not seem to support the notion of eccentricity as a significant factor for the speed of IT dispersion. Moreover, we could not observe stagnation or recirculation zones in our experiments.

There are safety limits to high-volume injections adding the CSF amount during drug administration. We have experienced in rat that no more than 10% of the CSF volume can be safely injected over a period of a couple of minutes (Venugopal et al., 2017). Higher injection impulse may also contribute to the possibility of high shear rates that nerve roots experience near the catheter tip, which may again pose an additional risk that requires clinical investigation.

We did not investigate the dispersion of active drug molecules because we believe that tracers were ideal to quantify the physical transport phenomena with optical methods. We also did not consider biological functions such as the uptake in spinal or cerebral tissue because the focus was on determining dispersion coefficients as a function of physical flow phenomenon, amplitude and frequency of the CSF. Uptake studies would require a different setting. The in vivo study by Tangen et al. (2017b) could be a starting point, but parametric studies of CSF amplitude and frequency in humans are not practical, if not unethical. Therefore, the replica of the human CSF conditions seems a good choice for the parametric dispersion studies pursued here.

The biomechanical stress–strain relation of the epidural spaces is a function of several poorly understood factors. These include viscous resistance exerted by venous blood volume and fatty tissues, possible elastic deformation resistance of nerve roots bulging into peripheral distal spaces, and the biomechanical properties of the dura membrane including the stiffening effect of trabeculae and ligaments. Our system does neither consider CSF generation in the choroid plexus nor reabsorption (elimination) in the subarachnoid villi. These additional features could be incorporated in a more refined head model, which was not needed in our study focusing on spinal dynamics. The current model was able to induce cervical CSF displacements (stroke volume 0–1 mL/beat) within the physiological range but was not designed to faithfully reproduce the biomechanical compliance of the spinal compartment. Accordingly, we also did not attempt to measure the absolute pressure changes and gradients that occur during infusion. We plan to study the rise in the line pressure or the increase in pressure in the CSF of the spinal SAS as a function of injection parameters in future work.

Another related limitation pertains to the practice of fluid removal before high-volume injection (initial CSF tapping). There was no attempt made to interrogate the CNS model regarding the biomechanical response of the CSF spaces subjected to tapping.