Determination of spinal tracer dispersion after intrathecal injection in a deformable CNS model

Background: Traditionally, there is a widely held belief that drug dispersion after intrathecal (IT) delivery is confined locally near the injection site. We posit that high-volume infusions can overcome this perceived limitation of IT administration. Methods: To test our hypothesis, subject-specific deformable phantom models of the human central nervous system were manufactured so that tracer infusion could be realistically replicated in vitro over the entire physiological range of pulsating cerebrospinal fluid (CSF) amplitudes and frequencies. The distribution of IT injected tracers was studied systematically with high-speed optical methods to determine its dependence on injection parameters (infusion volume, flow rate, and catheter configurations) and natural CSF oscillations in a deformable model of the central nervous system (CNS). Results: Optical imaging analysis of high-volume infusion experiments showed that tracers spread quickly throughout the spinal subarachnoid space, reaching the cervical region in less than 10 min. The experimentally observed biodispersion is much slower than suggested by the Taylor–Aris dispersion theory. Our experiments indicate that micro-mixing patterns induced by oscillatory CSF flow around microanatomical features such as nerve roots significantly accelerate solute transport. Strong micro-mixing effects due to anatomical features in the spinal subarachnoid space were found to be active in intrathecal drug administration but were not considered in prior dispersion theories. Their omission explains why prior models developed in the engineering community are poor predictors for IT delivery. Conclusion: Our experiments support the feasibility of targeting large sections of the neuroaxis or brain utilizing high-volume IT injection protocols. The experimental tracer dispersion profiles acquired with an anatomically accurate, deformable, and closed in vitro human CNS analog informed a new predictive model of tracer dispersion as a function of physiological CSF pulsations and adjustable infusion parameters. The ability to predict spatiotemporal dispersion patterns is an essential prerequisite for exploring new indications of IT drug delivery that targets specific regions in the CNS or the brain.


ELECTRONIC SUPPLEMENTARY MATERIALS Appendix A: Realization of cerebrospinal fluid pulsation due to vascular expansion
Pulsatile fluid motion was induced via an inflatable cerebrovascular compartment generating cerebrospinal fluid (CSF) flow patterns in the cervical, thoracic, and lumbar regions in a deformable spinal subarachnoid space (SAS) model that covered a physiological range of infusion volumetric flow rate (IVF) obtained previously in subject-specific magnetic resonance imaging CSF flow measurements.
Induction of CSF motion.Induction of pulsatile cerebrovascular expansion is controlled by a distensible balloon connected to a peristaltic piston pump, which puts in motion artificial CSF inside the spinal SAS.
The peristaltic piston pump has two settings for the displacement volume moved with each stroke (0 to 1.0 mL) and allows for change in pulsation frequency (0 to 300 bpm).The pump has a set duration of stroke time.The vascular volumetric strain mainly serves to displace fluid from the cranial to the spinal CSF.Dilation of the vascular balloon can be adjusted using a peristaltic piston pump (Walchem E-Class Metering Pump) to control the range of expansion and contraction to generate desired levels of CSF stroke volumes (i.e., the additional volume of CSF introduced into the extramedullary space of the central nervous system model with each beat) of up to 1.0 ml/beat at oscillations from 0-180 beats per minute in the cervical region.
Injection system.Precise settings for bolus or continuous infusion volumes and timing were controlled with a syringe pump (Infusion Pump: NE-4000 Syringe Pump).Bolus volumes in the range of 1.0 mL and 2.0 mL of the trypan blue dye (0.1% and 0.4%-vol aq., Sigma-Aldrich) were selected.The exit velocity and kinetic energies for different needles and infusion volumetric flow rate is shown in table A1.
During and after the infusion, CSF oscillations were maintained with the expanding vascular component (=piston pump) at 40, 72, and 120 beats per minute.These conditions were set to match the natural dynamics of pulsatile CSF found with cine phase-contrast magnetic resonance imaging of the human cranium and spinal canal.Another approach, as shown in Fig. C1.B, was also used which is the averaging over a selected pixels in the vertical columns (radial direction).The same results were observed in both cases.The RGB data obtained using the two approaches are plotted in Fig. R2 below.The curve overlap which shows that there is no noticeable difference in both methods.Grayscale formula • The RGB extraction from the frame.The RGB ranges from 0 to 255.
• With the use of the gray scale formular, intensities were obtained.
• The concentration profiles were plotted for each frame and scaled to have the maximum concentration as 0.4 (%M).
• Area under the curves is constant to ensure that mass is conserved.
• Second moments were computed, and the half of the gradient of the second moment plot gives the dispersion coefficient

Conversion formula from intensity to tracer concentration needed in Method of Moments (MoM).
Note that knowledge of the absolute concentration is not needed in the MoM, although this can be done by calibration [A1].
Where  is the intensity curve   = maximum value of the intensity curve   = minimum value of the intensity curve   = the intensity level  at each position i   = the inferred intensity (concentration) at a point along the neuraxis Fig. C3 shows the evolution of the tracer using the grayscale formula for experiment 40bpm and 72bpm at 1.0mL/Stroke.Figure C5 shows the plots of the first moment for all the experimental data.To have clean image data, a white florescence light was used with a white canvas for background control.
All experiments were run three times to ensure reproducibility.The low variation in the results obtained for each three experiments shows that the study is reproducible.
Fig. C3.The evolution of the tracer in the domain using the grayscale formula for experiment 40bpm and 72bpm both at 1.0mL/Stroke.It must be noted that the white offset has been done to ensure that the concentration profile pattern is correct.Also, the area under the curves is constant to ensure accuracy of the analysis.Binary (0-1) Image analysis.Alternative method not requiring quantitative intensity information at all.
This method could be used on qualitative (binary) immunobiological essays.
• The RGB extraction from the frame • With the use of the gray scale formular as shown in Eq. (C1), intensities were obtained.
• Binary method was applied to assign the concentration of 0.4 (%M) form region where there is tracer (i.e., where the value of the gray data is 1).0 (%M) was assigned to the region where gray data is 0.
• The concentration profiles were plotted for each frame • Area under the curves is constant to ensure that mass is conserved.
• Second moments were computed, and the half of the gradient of the second moment plot gives the dispersion coefficient.
Fig. C5 shows the evolution of the tracer using the binary method for experiment 40bpm and 72bpm at 1.0mL/Stroke.Womersley number ().We determined the dispersion coefficient experimentally with the method of moment and listed each dispersion coefficient in the context of the dimensionless experimental settings.
We plotted the experimental dispersion coefficient as a function of frequency and stroke volume.We also plotted the experimental dispersion coefficient as a function of Peclet number and Womersley number.The description of how to compute the Peclet number is given in Eq. (E1).The major term in Eq. ( E1) is the root-mean-square velocity ( .. ) of the solute.Considering that the domain of study is deformable and there is no flow out of the domain, the velocity of the solute is attenuated from the cervical region to the lumbar region.In our study, the mean velocity used for computation is that at the cervical region which is the maximum velocity in the domain.Table E1 shows the parameter used for the computation of the dimensionless quantities.Where V is the total volume change in cm 3 , V0 is the initial volume of CSF in cm 3 in the system,  is the angular frequency, Nb.p.m is the frequency in beat per minute, and t is the time in seconds.
The volume change was modeled as a function of time using cosine function which later converted to the instantaneous volume change through the differentiation of the volume change.Using cross-section area of the spine at the cervical region, the instantaneous volume change was converted into the velocity of the CSF at the cervical region.Finally, the root-mean-square of the obtained velocity was obtained which was used in the computation of the Peclet number.
It must be noted that the root-mean-square velocity varies along the length of the spine with its value as zero at the sacral which is as result of the deformation of the spine.In our analysis of the TAD, we used the highest  .. in the spine which is obtained at the cervical region.This is to estimate the highest dispersion that could be obtained in the system at different stroke volumes and different frequencies.
The Womersley number f is the frequency, T = 1 sec, is the period of oscillation, and υ = 7x10 −7 m 2 /sec, is the kinematic viscosity.
Table E2 and table E3 show the dispersion coefficients data for stroke volume of 0.5mL/Stroke and 1.0mL/Stroke respectively.The tables also show the data for dispersion coefficient using Taylor approach as shown in Eq. ( 9) in the main text.
From table E2 and table E3, some of the experiments were repeated for some frequencies.Even though variance exists, what we can see is that the trend is clear.An empirical correlation Eq. (E3) between apparent diffusion coefficient and CSF pulsations, a function of CSF amplitude and oscillation was established as a function of the root-mean-square velocity of the CSF, measured in the cervical region, as well as the frequency of the CSF pulsations.Where () =  0 +  1  +  2  2 . 0 = −1.0386, 1 = 0.0055,   2 = −5.8858x 10 −5 are constant terms, and  = 0.7419.  is the root-mean-square velocity of the CSF,  is the frequency and   is the experimental dispersion coefficient for the dimensional model.
made of the set and converted to a greyscale matrix.The difference of the greyscale images with tracer present to the baseline reference is saved as a third matrix and saved as a 1 or 0 to create a mask of the region of interest (ROI) in the original images where tracer is present.These frames (matrix pages) are now the masks respective of each frame time point.With the mask frame matrix, the software matches each greyscale difference frame to the original, unaltered image to select only the pixel red green blue (RGB) values of tracer in the original, unaltered frames using the mask ROI.A waterfall diagram of the final sequence of frames from the video recordings from the experiment is inferred from the software, which represents the mass flux of the tracer for each time point of phase-2 (Fig.B1).A calibration of known concentrations of trypan blue matched to the RGB values of images of the known concentrations provides a sure fit of a concentration of trypan blue to a given RGB value within an image of tracer in the experiments.Extraction of the concentration value saved to a fourth matrix occurs only for the pixels of the original frames located at the mask ROI.Then, the software saves each frame of the pixel ROI determined from the mask matrix as an RGB to a fifth matrix.Each frame of the fifth matrix is computed for the mean RGB value for each row of the matrix.The semi-automatic program in MATLAB calculates the likely concentration with the polynomial fit.Normalization is performed to ensure the mass flux represented by the integral area of the curve remains constant for each curve, and the results are plotted and saved to a '.mat' file for further analysis.

Fig. B1 .
Fig. B1.The waterfall represents the mass flux of the tracer for each time point of phase-2.

Fig
Fig. C1.(A) When the RGB data in the center pixels are analyzed.(B) When RGB data in the pixels between the selected points a and b are analyzed and averaged in radial direction (C) Sample of the raw image data obtained on which pixels and RGB data are obtained.

Fig. C2 .
Fig. C2. Figure showing the plot of the RGB intensities using both methods (Fig. C1 A and Fig. C1 B) described above

Fig. C4 .
Fig. C4.First moment plots for all the experimental data.The first moment, M1, indicates that the tracer moves in cranial direction from the injection position (50 cm) towards the thoracic region position (44 cm).

Fig. C5 .
Fig.C5.The evolution of the tracer in the domain using the grayscale formula (Binary) for experiment 40bpm and 72bpm both at 1.0mL/Stroke.The area under the curves is constant to ensure accuracy of the analysis.

Fig. D2 .
Fig. D2.Second moment plots for all the experimental data which shows the gradients of the second moment with which the dispersion coefficients are obtained.

Table E1 :
The table of the parameter used for the computation.Where  .. is the root-mean-square velocity.The values for the Peclet number are shown in tableE2and table E3 for different velocities.The  .. of the CSF was obtained using the different stroke volumes at different frequencies.The stroke volume model used is shown in Eq. (E2).For the stroke volume of c mL/stroke, the stroke volume conversion to Total volume is shown in Eq. (E2)