Mixed-dimensional fluid–structure interaction simulations reveal key mechanisms of cerebrospinal fluid dynamics in the spinal canal

Abstract

Cerebrospinal flow dynamics (CSF) plays a critical role in structural disorders of the central nervous system (CNS) and in the design of effective procedures for intrathecal drug delivery. Medical imaging techniques have only partially characterized CSF dynamics. Computational models have the potential to offer a high-resolution description of CSF flow and advance our mechanistic understanding. However, anatomically-accurate computational models of CSF dynamics in the spinal canal have largely ignored the compliance of the spinal tissues, which is critical to understand the pulse wave velocity and the craniocaudal decay of CSF pulsations. Here, we propose a mixed-dimensional fluid-structure interaction method that enables high-fidelity simulations of CSF dynamics on anatomically-accurate models of the spinal canal, considering the tissue compliance effects emerging from the dura mater and epidural fat. Our mixed-dimensional approach bypasses a critical computational bottleneck that emerges from the multiscale geometry of spinal tissues. Our results show that accurate modeling of tissue compliance is critical to capture key elements of CSF dynamics. This work opens new possibilities to control and optimize intrathecal drug delivery and to understand structural abnormalities of the CNS.

Introduction

Central nervous system (CNS) disorders are a serious concern to global healthcare systems. Most of these conditions are incurable, and often require frequent therapeutic interventions to alleviate symptoms or slow down the disease progression [1, 2]. A major challenge in treating CNS disorders is drug delivery. The blood-brain barrier (BBB) blocks foreign agents from entering the CNS, including 98\(\%\) of drug molecules [3]. Hence, conventional routes of drug administration such as intravenous, oral and subcutaneous are ineffective in CNS therapeutics. Intrathecal administration is a promising approach to deliver drugs to the CNS. In this delivery modality, the therapeutic agent is injected in the spinal subarachnoid space (SAS), the region of the CNS that is located between the dura mater and the spinal cord. The SAS is filled with cerebrospinal fluid (CSF), which flows around the brain and the spinal cord and acts as a carrier of the drug. Intrathecal injection circumvents the BBB because the drug is delivered directly in the CNS [4,5,6]. Intrathecal injection has been used to deliver anesthetic agents for over a century [5, 7]. These injections are usually given at the lower lumbar levels to avoid spinal cord damage. However, the therapeutic targets are located throughout the CNS, including the spinal cord nerve roots and the supraspinal regions of the brain. The desired distribution of the drug depends on the therapeutic goals, which may range from localized delivery to broader dispersion along the spinal cord or even to supraspinal regions. A thorough understanding of the CSF dynamics and drug transport in the intrathecal space is essential to address the diverse targeting requirements in intrathecal delivery. Beyond drug delivery, investigating CSF motion is also important to mechanistically understand various CNS disorders that have been correlated with abnormal CSF motion [8,9,10]. A deeper mechanistic understanding of CSF dynamics within disease-specific environments, could enable better management and treatment strategies for CSF-related neurological disorders.

While the study of CSF dynamics has been dominated by invasive and non-invasive imaging techniques, computational modeling has emerged as a complementary tool that enables controlled manipulation of various parameters, such as, material properties and anatomical features. Earlier computational studies on CSF dynamics have considered CSF flow in multiple pathological situations [11,12,13,14,15,16,17] and have studied the impact of anatomical structures such as nerve roots [18, 19], trabeculae [20, 21] and denticulate ligaments [22] on the CSF flow. Recent computational models of the cranium have incorporated fluid–structure interaction (FSI) to capture the interactions between the parenchyma and SAS [23,24,25,26]. However, computational models of spinal CSF flow have assumed that the SAS is rigid, ignoring the key role of tissue compliance on CSF dynamics. Ignoring tissue compliance poses fundamental limitations in the predictive capabilities of the model, including the inability to accurately capture the timings of flow rate waveforms and the attenuation of CSF pulsations down the spine. This limitation is illustrated in [27], where fluid flow simulations were compared to 4D phase-contrast magnetic resonance images (PCMRI). Another fundamental limitation of fluid flow simulations of the spinal canal that ignore tissue compliance is that modeling pulsatile flow requires adding an artificial outlet boundary at the sacral end. This is a consequence of the incompressibility of CSF. Because CSF cannot be compressed or expanded, when we add or remove fluid from the spinal canal through a time-dependent boundary condition at the cranial end of the spine, the volume of the SAS needs to change accordingly. A common way to circumvent this problem is opening an outlet boundary in the spinal canal, but this is inaccurate and fundamentally changes the flow dynamics. More recent literature has attempted to model the effects of tissue compliance on CSF flow by imposing a pre-defined boundary motion on the walls of the SAS. Reference [28] has used this approach to simulate the effects of intrathoracic pressure changes caused by respiration. Similarly, in [29], the boundary motion was inferred from the phase lag and attenuation of the velocity peak amplitude of CSF velocity fields acquired from cine MRI. In reference [30], the authors used MRI-based nonuniform volumetric flow rates to prescribe radial displacements of the dura in 1 mm increments along the spinal canal. At each slice of the spinal canal, the difference between volumetric flow rate at the slice above and below determines the time-dependent radial displacement of the boundary. However, the approach based on pre-defined boundary conditions requires 4D MRI data and usually makes strong assumptions such as radially uniform displacement at each cross section of the spine. Reference [31] describes a computational model, where the dura mater is rigid, but the spinal cord is mobile and deformable. They found significant alterations in the flow with respect to flow simulations in a SAS that is time independent, but they also included an artificial outlet boundary in the caudal end of the spine, presumably because the potential compression and expansion of the spinal cord is not large enough to accommodate the fluid volume changes that occur during pulsatile flow in physiological conditions. Reference [32] used an anatomically idealized spinal SAS model with elastic dura mater showing that CSF steady streaming can be a major driving force of CSF bulk motion. However, in addition to the idealized geometry, the model is based on several approximations and uses a simplified linear elastic model for the outer wall of the SAS. To our knowledge fully coupled FSI simulations of the CSF and surrounding soft tissues have not been reported in the literature. One reason for the absence of such simulations is that accurate models of the mechanical response of the tissues surrounding the SAS requires specialized computational tools. The dura mater is a thin (\(\sim \!200\,\upmu\)m) membrane with relatively large Young’s modulus (\(\sim \!4\) MPa) [33], while epidural tissue is thicker (\(\sim \!4\) mm) and has a much smaller Young’s modulus (\(\sim \!10^{-2}\) MPa). The combination of the two supporting tissues offers a complex mechanical response. Modeling dura mater as a volumetric solid requires very fine meshes, with at least 4-5 linear elements through the thickness, to avoid membrane locking. Even if the refinement is highly localized, it will inevitably propagate into the SAS and the epidural tissue, leading to a very large number of degrees of freedom.

Here, we propose a mixed-dimensional FSI model that addresses the aforementioned computational challenges and allows us to perform the first simulations of the entire spinal canal accounting for tissue compliance through a high-fidelity model of the soft tissues surrounding the SAS. In our mixed-dimensional approach, the dura mater is modeled as a curved membrane using surface kinematics, and epidural tissue is modeled as a volumetric solid. We perform simulations on a three-dimensional, anatomically accurate geometry of full-length human spine with closed caudal end. On the foramen magnum (FM), we impose a pulsatile velocity boundary condition obtained from MRI measurements. The flow from the pulsatile inlet boundary condition is accommodated by the deformations of the dura mater and epidural tissues. Using our FSI simulations, we describe the flow patterns, CSF velocities, pressure variations and the tissue displacements comparing them to experimental measurements. We show that the FSI simulation captures flow patterns that are more realistic than those from rigid-wall simulations that ignore tissue compliance. Finally, we examine the influence of tissue displacements on the cyclic mean velocities. Our simulations indicate that modeling of FSI is essential to capture the CSF flow patterns that are observed in vivo.

Biophysical model

Anatomical model

Figure 1a illustrates a schematic cross-sectional view of the spinal canal. We model CSF flow within the SAS (blue) and its interactions with the surrounding solid tissues, including the dura mater (orange) and epidural fat (gray). In this study we disregard the presence of nerve roots (green), which may alter the flow patterns and provide structural support to the spinal cord (white). We also neglect CSF absorption within the SAS, which may occur around the nerve root sleeves because its influence is negligible relative to the dominant pulsatile component. We treat the spinal cord as rigid because the combination of its material properties and geometry make it much less deformable than dura mater or epidural fat. We assume that the epidural space solely contains fat, ignoring the presence of the venous plexus. The distribution of epidural fat along the spinal canal is uneven and discontinuous [34], but its precise geometry is not well characterized in the literature. Thus, in the absence of better data, we assume a uniform epidural tissue thickness of 4 mm throughout the spine.

[IMAGE OMITTED: SEE PDF]

We use triangulated surface-mesh files of the dura mater and spinal cord from [35] to construct the geometries of the SAS and the dura mater; see Fig. 1b. The geometries represent a full length anatomically accurate spine from the foramen magnum at the base of the skull to the coccyx. In Fig. 1b, we label the spine according to its vertebral levels, from C1 to S5, where C, T, L, S represent, respectively, the cervical, thoracic, lumbar, and sacral regions. In the sagittal plane, the spine displays significant curvature in three regions, whereas in the coronal plane, it appears nearly flat. The spinal cord geometry (gold colored in Fig. 1b) is included from foramen magnum down to the conus medullaris, where it tapers rapidly. We disregard the spinal cord geometry below the conus medullaris because its diameter is negligible and would not affect the CSF dynamics. We generate the geometry of the epidural tissue by uniformly extruding the dura mater surface outward by 4 mm along its surface normal; see Fig. 1b and the gray-shaded region in Fig. 1c.

Computational domains and boundaries

The computational domain \(\Omega\) consists of three subdomains such that \(\Omega = \Omega ^f \cup \Omega ^m \cup \Omega ^s\). To avoid going into distracting mathematical details, we treat these subdomains as open or closed sets as convenient. Here, \(\Omega ^f\) denotes the SAS, which is filled with CSF, \(\Omega ^m\) represents the dura mater and \(\Omega ^s\) is the epidural tissue. Unless otherwise specified all four domains evolve with time. The three subdomains \(\Omega ^f\), \(\Omega ^m\) and \(\Omega ^s\) are spatially arranged as shown in Fig. 1b. Because the spinal cord is assumed to be rigid, it is not part of our computational domain. The SAS and the dura mater meet at a surface that we denote by \(\Gamma ^{fm}\). The surface separating the dura mater and the epidural tissue is denoted by \(\Gamma ^{ms}\). We call \(\Gamma ^k\) the boundary of \(\Omega ^k\) for \(k=f,\, m,\, s\). The unit outward normal to \(\Gamma ^k\) is \(\varvec{n}^k\). Due to the topological arrangements of the subdomains it follows that \(\varvec{n}^f = -\varvec{n}^m \text { on } \Gamma ^{fm}\) and \(\varvec{n}^m = -\varvec{n}^s \text { on } \Gamma ^{ms}\).

Cerebrospinal fluid flow

The strong form of the mass conservation and momentum balance equations in Eulerian description is

$$\begin{aligned} \rho ^f \dot{\varvec{v}}^f&- \nabla \cdot \varvec{\sigma }^f = \varvec{f}^f&\quad&\text {in } \Omega ^f\times [0,T], \end{aligned}$$

(1)

$$\begin{aligned} \nabla \cdot \varvec{v}^f&= 0&\quad&\text {in } \Omega ^f\times [0,T]. \end{aligned}$$

(2)

Here, \(\rho ^f\) is the CSF density, \(\varvec{v}^f(\varvec{x},t)\) is the velocity field, \(\varvec{x}\) is a spatial point, t denotes time, [0, T] is our time interval of interest, \({\dot{\varvec{v}}}^f= \partial _t \varvec{v}^f + \varvec{v}^f \cdot \nabla \varvec{v}^f\) denotes the material derivative with respect to the fluid motion, and \(\varvec{f}^f\) is a force per unit volume. We assume CSF can be modeled as a Newtonian fluid, thus

$$\begin{aligned} \varvec{\sigma }^f = -p^f\varvec{I} + \mu ^f (\nabla \varvec{v}^f + \nabla ^T \varvec{v}^f), \end{aligned}$$

where \(\varvec{\sigma ^f}\) is the Cauchy stress tensor, \(p^f(\varvec{x},t)\) is the pressure field, \(\varvec{I}\) is the identity tensor \(\mu ^f\) is the dynamic viscosity, and \(\nabla ^T\) denotes the transpose operator of \(\nabla\).

Soft tissues

We model dura mater and epidural tissue as nonlinear hyperelastic materials with finite deformation kinematics. The elastodynamics equations in Lagrangian description read

$$\begin{aligned} \rho _0 \frac{\partial ^2 \varvec{u}}{\partial t^2} \bigg |_{\varvec{X}} - \nabla _{\varvec{X}} \cdot (\varvec{FS}) = \varvec{f}_0 \qquad \text {in } \Omega ^s_0 \cup \Omega ^m_0\times [0,T]. \end{aligned}$$

(3)

Here, \(\rho _0\) is the solid mass density in the reference configuration, \(\varvec{u}(\varvec{X},t)\) is the displacement field, \(\varvec{X}\) denotes a material particle, \(\varvec{F}= \nabla _{\varvec{X}} \varvec{u} + \varvec{I}\) is the deformation gradient, \(\varvec{S}\) is the second Piola-Kirchhoff stress tensor, and \(\varvec{f}_0\) is a force per unit volume of the reference configuration. The domains \(\Omega ^s_0\) and \(\Omega ^m_0\) represent, respectively, the reference configurations of \(\Omega ^s\) and \(\Omega ^m\). In Eq. (3), \(\vert _{\varvec{X}}\) indicates that the time derivative is taken holding \(\varvec{X}\) constant, and \(\nabla _{\varvec{X}}\) indicates that the gradient operator is with respect to material particles. Because we will use different kinematics and different material models for dura mater and epidural tissue, we will use the displacement representation \(\varvec{u}=\varvec{u}^s\) in \(\Omega ^s_0\), and \(\varvec{u}=\varvec{u}^m\) in \(\Omega ^m_0\), where \(\varvec{u}^m=\varvec{u}^s\) on \(\Gamma ^{ms}\).

Kinematics and material model of the epidural tissue

We treat epidural tissue as a volumetric solid modeled as a compressible Neo-Hookean material. We adopt the strain energy density \(\psi ^s\) with dilatational penalty from [36], given by,

$$\begin{aligned} \psi ^s = \frac{\mu ^s}{2} \left( J^{-2/3} \text {tr}(\varvec{C}) - 3 \right) + \frac{\kappa ^s}{2} \left( \frac{1}{2} (J^2 -1) -\ln J \right) . \end{aligned}$$

(4)

In Eq. (4), \(\varvec{C}=\varvec{F^TF}\) is the right Cauchy Green strain tensor, \(J = \text {det}(\varvec{F})\), \(\operatorname {tr}(\cdot )\) denotes the trace operator, while \(\mu ^s\) and \(\kappa ^s\) are, respectively, the shear and bulk modulus of the tissue. The shear and bulk moduli can be obtained from the Young’s modulus (\(E^s\)) and the Poisson ratio (\(\nu ^s\)) as \(\mu ^s = \frac{E^s}{2(1+\nu ^s)},\, \kappa ^s = \frac{E^s}{3(1-2\nu ^s)}\). For hyperelastic constitutive models, the second Piola-Kirchhoff stress tensor is obtained from the strain energy function as

$$\begin{aligned} & \varvec{S}^s = 2\frac{\partial \psi ^s}{\partial \varvec{C}} = \mu ^s J^{-2/3} \left( \varvec{I} - \frac{1}{3} \varvec{C}^{-1} \text {tr}(\varvec{C})\right) \nonumber \\ & \quad + \frac{1}{2} \kappa ^s \left( J^2 - 1 \right) \varvec{C}^{-1}. \end{aligned}$$

(5)

Kinematics and material model of dura mater

Unlike epidural tissue, the dura mater is very thin, and modeling it as a volumetric solid poses two challenges. First, we need sufficient mesh resolution through the dura thickness to avoid locking issues, which leads to an excessive proliferation of degrees of freedom. Second, even if meshing challenges are overcome the resulting mesh is highly non-uniform leading to linear systems with exceedingly large condition numbers. We avoid these challenges by using surface kinematics to describe the dura mater. We first express \(\Omega ^m\) as

$$\begin{aligned} \Omega ^m = \Gamma ^m \times [0,\zeta ], \end{aligned}$$

(6)

where \(\zeta\) is the thickness of the dura mater, and \(\Gamma ^m\) is the mid-plane surface of dura. Because we assume that \(\zeta\) is small, in our computations, we will identify \(\Gamma ^m\) with \(\Gamma ^{ms}\) or \(\Gamma ^{fm}\) as needed to obtain solid displacements or impose the fluid-solid interface conditions.

Because the bending stiffness of dura mater scales with the cube of its thickness, we will neglect its resistance to bending. In our model, dura mater will behave as a membrane with stress restricted to in-surface components. Although such mechanical formulation would be unstable in its own, the restriction \(\varvec{u}^m=\varvec{u}^s\) renders a stable formulation. We derive the in-surface stress tensor for dura mater following the approach for thin biological membranes described in [37].

The reference configuration of dura mater is represented by the surface \(\Gamma ^m_0\). We parametrize \(\Gamma ^m_0\) with the mapping \(\varvec{X} = \varvec{X}(\varvec{\xi })\), where \(\varvec{\xi }=\{\xi ^\alpha \}_{\alpha =1,\,2}\) are curvilinear coordinates on the surface; see Fig. 2. The tangent vectors in the reference configuration are \(\varvec{G}^\alpha = \partial \varvec{X}/\partial \xi ^\alpha\). The unit normal vector to \(\Gamma ^m_0\) is obtained as \(\varvec{N}=\varvec{G}^1\times \varvec{G}^2/||\varvec{G}^1\times \varvec{G}^2||\).

[IMAGE OMITTED: SEE PDF]

The deformed configuration of the dura mater is represented by the surface \(\Gamma ^m\), which we parametrize with the mapping \(\varvec{x} = \varvec{x}(\varvec{\xi }, t)\). Using the displacement field \(\varvec{u}^m\), the tangent vectors for the deformed surface are \(\varvec{g}^\alpha = \partial \varvec{x} / \partial \xi ^\alpha\); see Fig. 2.

To derive our constitutive model, we first decompose the total deformation gradient as

$$\begin{aligned} \varvec{F} = \varvec{F}_S + \varvec{F}_N \;\text { where }\; \varvec{F}_S = \varvec{g}_\alpha \otimes \varvec{G}^\alpha , \quad \varvec{F}_N = \lambda _N \varvec{n} \otimes \varvec{N}. \end{aligned}$$

(7)

In Eq. (7), \(\lambda _N\) is the stretch in the direction normal to the undeformed surface and the repetition of Greek index \(\alpha\) indicates summation. The surface deformation gradient \(\varvec{F}_S\) is the surface projection of the deformation gradient \(\varvec{F}\) and can be obtained using the projection tensor \(\varvec{I}_S\) as

$$\begin{aligned} \varvec{F}_S = \varvec{F}\varvec{I}_S \text { where } \varvec{I}_S = \varvec{I} - \varvec{N} \otimes \varvec{N} = \varvec{G}_{\alpha } \otimes \varvec{G}^{\alpha }. \end{aligned}$$

(8)

Using Eq. (7), the right Cauchy Green strain tensor, \(\varvec{C}\) can be split as

$$\begin{aligned} \varvec{C} = \varvec{C}_S + \varvec{C}_N \; \text { where } \; \varvec{C}_S = \varvec{F}_S^T \varvec{F}_S, \quad \varvec{C}_N = \lambda _N^2 \varvec{N} \otimes \varvec{N}. \end{aligned}$$

(9)

We model dura mater as an incompressible Neo-Hookean material with strain energy function

$$\begin{aligned} \psi ^m = \frac{\mu _m}{2}(\text {tr}(\varvec{C})-3) - \varsigma (J-1), \end{aligned}$$

(10)

where \(\varsigma\) is the Lagrange multiplier for volumetric incompressibility. The second Piola-Kirchhoff stress tensor is

$$\begin{aligned} \varvec{S}^m = 2\frac{\partial {\psi ^m}}{\partial \varvec{C}} = \mu ^m \varvec{I} - \varsigma J \varvec{C}^{-1}. \end{aligned}$$

(11)

To obtain the value of the Lagrange multiplier, we impose the restriction that the normal component of \(\varvec{S}^m\) be zero. The normal component of \(\varvec{S}^m\) can be computed as \(\varvec{S}^m_N = \varvec{S}^m \varvec{N}\otimes \varvec{N}\). From Eq. (9), we have \(\varvec{N}\otimes \varvec{N}=\lambda _N^{-2}\varvec{C}_N\), thus, \(\varvec{C}^{-1}\varvec{N}\otimes \varvec{N}= \lambda _N^{-2}\varvec{N}\otimes \varvec{N}\). Using this relation, we obtain

$$\begin{aligned} \varvec{S}^m_N = (\mu ^m - \varsigma J\lambda _N^{-2})\varvec{N}\otimes \varvec{N} \end{aligned}$$

(12)

Imposing \(\varvec{S}^m_N=0\), we obtain \(\varsigma =\lambda _N^{2}J^{-1}\mu ^m\). If we use \(J=1=J_A\lambda _N\), where \(J_A = ||\varvec{g}^1 \times \varvec{g}^2|| / ||\varvec{G}^1 \times \varvec{G}^2||\) is the surface stretch, we obtain the value of the Lagrange multiplier as \(\varsigma =\mu ^m/J_A^2\).

After substituting the in-surface stress condition, we obtain the second Piola stress tensor with only surface contributions

$$\begin{aligned} \varvec{S}^m = \mu ^m (\varvec{I} - J_A^{-2} \varvec{C}^{-1}) =\varvec{S}^m_S. \end{aligned}$$

(13)

Results

In this section, we present the flow predictions from our FSI simulation and compare them with those from a rigid-wall simulation and in vivo data from the literature.

Problem definition and parameters

Prestress in tissue

We assume that the tissues that we model are in equilibrium with the traction from the fluid and other surrounding tissues at all times [38]. Thus, the geometry obtained from magnetic resonance imaging (MRI) is not stress-free. In the absence of better information, here, we assume a prestress \(\varvec{\sigma }_0\) resulting from a mean pressure of \(p_0=600\) Pa exerted by the CSF [39] and the outer surrounding tissues

$$\begin{aligned} \varvec{\sigma }_0 = -p_0\varvec{I} \text { at } t=0. \end{aligned}$$

(14)

Boundary conditions

We obtain the CSF flow rate waveform at the C2 level (\(Q_{MRI}\)) from in vivo measurements reported in [35]. We assume a cardiac cycle period of \(T_C=0.8\) s and adjust the waveform to ensure zero net volumetric flow over the cycle; see Fig. 3a. The stroke volume associated with the flow rate waveform (defined as the net CSF volume entering the spinal canal during systole) is 0.385 mL. Since the foramen magnum is located close to the C2 level, we assume that the flow rate waveform at the foramen magnum equals \(Q_{MRI}\). Thus, we apply a time-dependent velocity Dirichlet boundary condition corresponding to the flow rate waveform \(Q_{MRI}\) at the foramen magnum. We apply a flat, blunt-shaped velocity profile at the inlet boundary (Fig. 3b). We impose no-slip boundary conditions at the spinal cord wall.

We restrict the motion of tissues at the foramen magnum boundary in all directions (Fig. 3c). Additionally, a zero-displacement Dirichlet condition is applied at the caudal end to represent rigid support from the coccyx.

On the outer wall of the epidural tissue, we employ a traction boundary condition corresponding to the uniform pressure \(p_0\). An alternative approach is to use the Robin boundary condition \(\varvec{\sigma }^s \varvec{n}^s = -k\varvec{u} -c\varvec{v}-p_{0}\varvec{n}^s\) where k and c are parameters that represent the material behavior of the the tissues surrounding the epidural fat on the exterior side; see [40]. Because reliable experimental data for k and c are lacking, we varied k over a plausible range (with \(c=0\)) and found that the overall agreement with in vivo is best for \(k=0\) (see Supplementary material). Hence, we adopt \(k=0\) for the simulations reported in this work.

Initial conditions and material parameters

Our simulations were initiated with the conditions

$$\begin{aligned} \varvec{v}(\varvec{x},0)&= \varvec{0}; \quad \varvec{x} \in \Omega ^f, \end{aligned}$$

(15)

$$\begin{aligned} p(\varvec{x},0)&= 600 \text { Pa}; \quad \varvec{x} \in \Omega ^f, \end{aligned}$$

(16)

$$\begin{aligned} \varvec{u}(\varvec{x},0)&= \varvec{0}; \quad \varvec{x} \in \Omega ^s. \end{aligned}$$

(17)

Because the above initial conditions are unrealistic, we ran the simulations until the flow field became periodic in time. We found that the flow became periodic with 92.4% accuracy after 5 cardiac cycles. We report results for the \(6^{\text {th}}\) cycle. To identify time in our results, we will use \(t_\%\in [0,100]\), which represents the percentage of the cardiac cycle in the 6th cycle. Table 1 lists the material properties used in this study.

[IMAGE OMITTED: SEE PDF]

Rigid wall simulations

To better understand the effect of tissue compliance, we also performed a rigid-wall simulation for comparison purposes. In our rigid-wall simulations, all solid tissues are assumed to be undeformable and do not need to be included in the computational domain. The rigid-wall simulation was performed on the domain \(\Omega ^f_{\text {rw}}\), which is slightly different from the initial configuration of \(\Omega ^f\). The motivation to modify the fluid domain is that rigid fluid domain is incompatible with the pulsatile inlet velocity that we impose at the foramen magnum because the fluid is incompressible. The most common approach in the literature to avoid this incompatibility is to open an outlet in the sacral region of the domain; see \(\Gamma _O^f\) in the inset in Fig. 3. We impose a zero pressure boundary condition on \(\Gamma _O^f\), which allows the fluid to exit the domain to satisfy mass conservation.

[IMAGE OMITTED: SEE PDF]

Flow patterns

To analyze the flow patterns along the spine, we define \(d=-z\) as the distance from the foramen magnum. At a distance d from the foramen magnum, we define the cross sections of the SAS (\(\Gamma ^f_d\)) as the intersection of \(\Omega ^f\) with the plane \(z=-d\). The flow rate across the cross section defined by d is given by \(Q_d (t) = \int _{\Gamma ^f_d}v_z(t) \text {d} a\).

In the rigid-wall simulation, the flow rate waveforms across all sections are identical to the inflow waveform. However, in the FSI simulations and in the in vivo measurements they vary as we move along the spine because the fluid domain is deforming with time. Figure 4a shows the flow rate waveforms at various sections. The red shaded region represents positive flow rates, indicating cranial flow, while the blue shaded region corresponds to negative flow rates, indicating caudal flow. The waveforms obtained from our FSI simulations vary noticeably across the sections, consistent with in vivo measurements from MRI taken at levels C2-C3, C7-C8, and T10-T11 [35]; see Fig. 4b.

The flow rate waveforms predicted in the FSI simulation capture two key features that are consistent with in vivo data: phase differences and craniocaudal decay in pulsation amplitude. These two aspects are analyzed separately in what follows.

Flow rate waveform phase differences

The phase differences in the FSI simulation can be quantified by tracking two key characteristics of the flow rate waveforms along the spine. The first one is the systolic peak: the time when the flow rate at a section reaches its maximum during a cycle. The second one is the flow reversal time: The time at which the flow direction reverses at the section. By recording the systolic peak along the spine, we can calculate the pulse wave speed (PWS), which is defined as the ratio of the distance traveled by the systolic peak to the time it takes to travel that distance. The PWS measures the compliance of the tissue and has been reported in multiple experimental papers based on MRI data. The PWS in our FSI simulation is spatially non-uniform, with an average value of 3.4 m/s, which is within the range of measured values in the literature; see Table 2. At the section defined by \(z=-d\), the phase difference, \(\Delta \phi (d)\), between the flow rate waveform \(Q_d\) at that section and the inflow waveform (\(Q_{MRI}\)) can be estimated as \(\Delta \phi (d) = 2\pi \frac{d}{\text {PWS} \times T_C}\), where \(T_C\) is the time period of the cardiac cycle.

[IMAGE OMITTED: SEE PDF]

[IMAGE OMITTED: SEE PDF]

The phase difference in the flow rate waveform at an intermediate location of the spine is also related to the time delay in flow reversal (\(Q_d=0\)) relative to the foramen magnum. Figure 4b depicts in vivo flow rate data from [35] showing that the flow reversal time changes from one section to another. Using the data from our FSI simulation, we produced a spatio-temporal distribution of the flow rates from the FSI simulation; see Fig. 4c. The star markers in Fig. 4c represent in vivo flow reversal times, while the black solid line defines the space-time curve of flow reversal predicted by our simulation. Our FSI simulation captures the spatial variation in flow reversal timing more accurately than the rigid-wall simulation. In the rigid-wall simulation, flow reversal occurs simultaneously at all locations, and the data would define a vertical straight line in Fig. 4c (not shown). Note that the flow rate waveforms in the lumbo-sacral regions are distorted with respect to the inflow waveform due to the pulse wave reflections near the caudal end. Hence, the descriptions of systolic peak and flow reversal in these regions are not meaningful.

Craniocaudal decay of the flow rate

Due to the tissue compliance, the pulse wave energy associated with the inflow waveform attenuates as it travels along the spine. To quantify the craniocaudal decay in CSF pulsations, we calculate the peak flow rate, which represents the amplitude of the flow rate waveform at a cross section.

The peak flow rate in the cranial and caudal direction is defined as \(Q_{d,\text {cra}}^{\text {peak}}=\max _{t_\%\in [0,100]} Q_d(t)\) and \(Q_{d,\text {cau}}^{\text {peak}}=\min _{t_\%\in [0,100]} Q_d(t)\) respectively. Figure 5a compares the peak flow rates obtained from our FSI simulation (solid lines), with those obtained from our rigid-wall simulation (dashed lines) and with the MRI observations (stars). In the rigid-wall simulation, \(Q_{d}^{\text {peak}}\) remains constant as d varies and always matches the flow rate imposed by the inlet boundary condition. However, in the results from the FSI simulation and the in vivo measurements, the peak flow rates varies along the spine. In the FSI results, the asymmetry between the peak cranial and peak caudal flow rate associated with the inflow waveform is small at the intermediate sections. For example, at \(d=10\) cm, the peak caudal and cranial volume flow rates are nearly equal (\(\approx 1.8\) mL/s), and match the values reported in vivo (1.95 ml/s); see Fig. 5a. Overall, this plot illustrates the inability of the rigid-wall model to capture the craniocaudal decay of the flow rate.

To further study the CSF flow velocities, we compare the peak velocities from our FSI simulations with those from rigid wall simulations and with in vivo experiments. We define the caudal and cranial peak velocity at distance d from the foramen magnum as

$$\begin{aligned} & v^{\text {peak}}_{d,\text {cau}} =\max _{t_\%\in [0,100]} \frac{ \int _{\Gamma _d^f} \varvec{v}\cdot \varvec{n}_d \text {d}a }{ \int _{\Gamma _d^f} \text {d}a }, \nonumber \\ & \quad v^{\text {peak}}_{d,\text {cra}} =\min _{t_\%\in [0,100]} \frac{ \int _{\Gamma _d^f} \varvec{v}\cdot \varvec{n}_d \text {d}a }{ \int _{\Gamma _d^f} \text {d}a }, \end{aligned}$$

(18)

where \(\varvec{n}_d\) is the normal to \(\Gamma _d^f\) in the cranial direction. Figure 5b shows \(v^{\text {peak}}_d\) from our simulations. Because the spine is closed at the caudal end in the FSI simulation, the peak velocities in the lower lumbar and sacral regions predicted by the FSI simulation are significantly smaller than those from the rigid-wall simulation. In the lumbar region, the value of \({v}_{d,\text {cau}}^{\text {peak}}\) obtained from the FSI simulation ranges from 2 cm/s to 0.25 cm/s, while in the sacral region is less than 0.25 cm/s. These results agree well with the experimental literature [47] (see stars in Fig 5b), where peak velocities of \(-1.07\pm\)0.49 cm/s and \(-0.32\pm\)0.33 cm/s are reported at the L2 and S1 levels, respectively.

[IMAGE OMITTED: SEE PDF]

Intraspinal pressure

We analyze the time variation of the average pressure at cross sections of the spine. We call this quantity the intraspinal pressure (ISP), which at the cross section defined by \(z=-d\), is computed as

$$\begin{aligned} \text {ISP}_d(t) = \frac{\int _{\Gamma ^f_d}p(t) \text {d} a}{\int _{\Gamma ^f_d} \text {d} a}. \end{aligned}$$

(19)

In the FSI simulation, the absolute pressure depends on the pressure initial condition and the prestress in the tissue. Hence, we only analyze the temporal and spatial variations of the pressure field. Figure 6 shows how the \(\hbox {ISP}_d\) varies with d and time in our FSI simulation. During the systolic phase, the pressure is higher in the cervical regions, and decreases toward the caudal end. At \(t_\%=50\), the pressure field becomes nearly uniform in space. In the latter half of the cycle, the pressure in the cervical regions begins to drop, resulting in pressure gradient in the cranial direction. By the early systole of the next cycle (\(t_\%\approx 5\)) the pressure gradient reverses, with the pressure higher in the cervical regions relative to the lower spinal levels.

[IMAGE OMITTED: SEE PDF]

In the FSI simulation, at any time in the cardiac cycle, the pressure difference between the foramen magnum and the caudal end does not exceed 60 Pa. By comparison, the maximum pressure difference is 140 Pa in the rigid-wall simulation. For validation, we refer to a study where the pressure was measured invasively using a catheter at the lumbar region. The amplitude of pressure at L2 level in the FSI simulation is 35 Pa (Fig. 6), which agrees with the clinical range of 4-10 mm \(\hbox {H}_2\)O (\(\approx\) 39-98 Pa) [48].

The pressure waveform at different levels exhibits a phase lag relative to the pressure waveform at the foramen magnum. This phase lag increases caudally, reaching approximately \(\pi\) at the caudal end. However, in the rigid-wall simulations, the pressure waveforms at all locations are in-phase.

In contrast to the rigid-wall simulations, where the pressure waves cannot be captured due to infinite speed, in the FSI simulations, we observe several wave reflections near the caudal end (Fig. 6). Successive reflections occur at \(t_\% \approx 40\), and \(t_\% \approx 60\). However, such wave reflections cannot be captured through MRI because they produce only small flow deviations, which are below the detection threshold.

Tissue displacements

The pulsatile motion of CSF in the spinal canal causes deformations in dura and epidural fat. Here, we are interested in two quantities: (i) percentage change in cross sectional area of SAS, and (ii) maximum tissue displacements. We calculate the maximum percentage change in area of cross-section of SAS due to the tissue deformations over the cardiac cycle as

$$\begin{aligned} {\% \Delta A^{SAS}_d = \frac{\max \limits _{t_\%\in [0,100]} \bigg (\int _{\Gamma ^f_d(t_\%)} \text {d} a \bigg ) - \int _{\Gamma ^f_d(t=0)} \text {d} a}{\int _{\Gamma ^f_d(t=0)} \text {d} a} \times 100.} \end{aligned}$$

(20)

The maximum tissue displacement is calculated as \(u^{\text {peak}}_d\! =\max _{\varvec{x}\in \Gamma _d^f} \max _{t_\%\in [0,100]} {\Vert \varvec{u}\Vert }\). Figure 7a plots \(\%\Delta A_d^{SAS}\) and \(u_d^{\text {peak}}\) with respect to distance from the foramen magnum. The maximum \(\%\Delta A^{SAS}_d\) occurs at C5 level and takes the value \(10\%\), which corresponds to an absolute change of \(\sim\)22 \(\text {mm}^2\) in the cross section. Secondary maxima occur near the T9 and L3 levels. We observe large tissue displacements, with a maximum of 1.5 mm occurring at the T4 level, i.e., \(d\approx\) 19 cm; see Fig. 7a. At the T4 level, however, \(\%\Delta A^{SAS}_d\) is very small. In general, the large tissue displacements occur predominantly with rigid body motions of the cross sections, and do not contribute significantly to \(\%\Delta A^{SAS}_d\).

[IMAGE OMITTED: SEE PDF]

Figure 7b, highlights the cross section where the rigid body displacements in the tissue are the largest (\(d\approx 19\) cm). The dashed lines indicate the boundary of dura mater in the undeformed configuration. We display the z-component of the velocity using a colormap to indicate axial flow direction. In the anterior side of spine, the SAS narrows during the systole (\(t_\%=12.5\)), and widens during diastole (\(t_\%=62.5\)). Conversely, in the posterior spine, the SAS widens during systole and narrows by mid-diastole.

The black arrows in the Fig. 7b illustrate the in-plane velocity vectors. The in-plane velocities in the FSI simulation are larger than those in the rigid-wall simulation due to the motion of the dura mater. For instance, at \(t_\%=12.5\), the dura displaces in the anterior direction, therefore pushing the CSF in the anterior direction. Whereas, at \(t_\%=75\), the dura displaces towards the posterior direction, pushing the CSF from anterior to posterior side.

To our knowledge, there is no experimental literature on dura deformations, hence we could not compare our model results in terms of tissue deformations. However, since we did not account for the mechanical support for the epidural fat, the large displacements observed (up to 1.5 mm) may be an overestimate. In actual physiology, the distribution of epidural fat is highly discontinuous and is backed by the rigid bone. Modeling tissue support around the epidural fat would require precise knowledge of the support locations.

Steady-streaming

Pulsatile flow in a deformable channel results in a non-zero mean flow over a cycle. This phenomenon is referred to as steady-streaming [49, 50]. In the radionuclide scanning studies [51, 52], the authors observed that a tracer injected in the lumbar region moves upward and reaches the basal cisterns in 15-20 min. This observation indicates net upward motion of the tracer over one cardiac cycle with velocities in the order of 1 cm/min, which is significantly faster than the velocities expected from diffusion. The precise reason for this rapid net upward motion of tracers is multifaceted and poorly understood. Recent studies [32, 53] on steady-streaming in simplified geometrical models of the spine have shown that, under certain assumptions, the cumulative effects of nonlinear convective accelerations result in a CSF bulk motion consistent with the observations from radionuclide scanning studies [51, 52]. Here, we explore the steady-streaming in our simulations due to the nonlinear convective acceleration and tissue deformation in an accurate three-dimensional patient specific geometry.

Steady-streaming has been quantified in the literature using the Eulerian and Lagrangian descriptions. The Eulerian steady-streaming velocity is the cyclic time averaged velocity at fixed points in space. When the fluid domain moves with time, as in our FSI simulation, the exact definition of Eulerian steady-streaming velocity becomes inapplicable because some regions outside of the initial configuration can have active mean velocities. The concept of Eulerian steady-streaming velocity can be extended to moving domains in multiple ways, none of them completely satisfactory. Here, we define the Eulerian steady-streaming velocity as the cyclic time averaged velocity as observed from points fixed with respect to the moving domain,

$$\begin{aligned} \langle \varvec{v} ({\widehat{\varvec{x}}}) \rangle _{E} = \frac{1}{T_c}\int _0^{T_c} \varvec{v}(\widehat{\varvec{x}},t) \,\text {d}t. \end{aligned}$$

(21)

The Lagrangian steady-streaming velocity is the time averaged velocity of a fluid particle along its trajectory over a cardiac cycle,

$$\begin{aligned} \langle \varvec{v}(\varvec{X}) \rangle _{L} = \frac{1}{T_c}\left[ \varvec{X}(t+T_c) - \varvec{X}(t)\right] . \end{aligned}$$

(22)

We carried out the calculations in Eqs. (21) and (22) using PyVista [54]. Both the Eulerian and Lagrangian steady-streaming velocity fields in our FSI simulation are nearly identical (not shown here), therefore we only report the steady-streaming fields in Eulerian description.

Figures 8a and 8b show the z-component of steady-streaming velocities from our rigid-wall simulation and our FSI simulation, respectively. Negative and positive velocity represent, respectively, caudal and cranial flow. The steady-streaming velocities in the rigid-wall simulation are smaller than those in the FSI simulation and do not exhibit any sensible pattern. In the FSI simulation, steady-streaming velocities show a bi-directional flow, with upward (cranial) steady-streaming on one side (anterior or posterior) and downward steady-streaming on the other side. These flow directions alternate along the length of the spine; see Fig. 8b. We name the regions with alternating patterns as I, II and III; see Fig. 8b. The inset figures (top and bottom) in Fig. 8b show the 3D glyphs of the steady-streaming velocity inside II, and near the intersection between II, III. The glyph size is scaled to the z-component of the steady-streaming velocity. The flow in II is bi-directional with upward steady-streaming in the posterior side, and downward in the anterior side. In the lateral side of spine in II, we observe that the flow is predominantly azimuthal, with velocities of nearly the same order of magnitude as the axial steady-streaming velocities.

The flow patterns in II agree with those reported in [32, 53], where the alternating bi-directional flow regions are identified as vortices divided by separation lines. In our simulation, we do not observe clear vortex type behavior; see bottom inset in Fig. 8b. The black line with arrows illustrates the flow direction near the intersection between II and III. The caudal flow in anterior side of II deviates azimuthally as it approaches III and only a part of the flow returns back into the II, while rest of the flow passes into III. Moreover in [32, 53], the authors show that the eccentricity of the spinal cord made with dura mater is a key geometrical feature that determines the location of separation lines. However in our simulation the eccentricity of spinal cord changes significantly with time, and we observe the separation lines near C5-C6 and T8-T9 levels, which are located slightly above the levels where the eccentricity of spinal cord in initial configuration reverses.

[IMAGE OMITTED: SEE PDF]

To put the magnitudes of steady-streaming velocity in context, in our FSI simulation the spatial maximum of \(\langle \varvec{v} \rangle _{E}\) is 0.54 cm/s in the caudal direction, which occurs near the T5 level. Because the pulsatile velocities in the FSI simulation are smaller at the lumbar and sacral levels, we observe smaller steady-streaming velocities (order of 0.001 cm/s) relative to those in the cervical and thoracic levels. To our knowledge, there is no experimental study that has characterized steady-streaming velocities due to their small value, hence we could not be compare them with experiments. Comparing our steady-streaming velocities with the radionuclide scanning observations at this stage is premature because tracer transport needs to account for several factors such as binding and degradation.

To further quantify the effects of tissue deformation on steady-streaming, we evaluate the steady-streaming flow rate across a section in the upward or downward direction. The steady-streaming flow rates in the upward direction at the section \(z=-d\) are defined, respectively, as

$$\begin{aligned} Q^{\text {SS}}_d = \int _{\Gamma _d^f} \max (\langle v_z \rangle _{E},0) \text {d}a. \end{aligned}$$

(23)

[IMAGE OMITTED: SEE PDF]

Figure 9 compares \(Q^{\text {SS}}_d\) in rigid-wall simulations with those in FSI simulations. The \(Q^{\text {SS}}_d\) in FSI simulation is \(\sim\)10 times larger than in the rigid-wall simulation. In the FSI data, \(Q^{\text {SS}}_d\) exhibits local maxima at \(d\approx 14\) cm and \(d\approx 40\) cm, where the tissue displacements are largest; and local minima at \(d\approx 31\) cm, where the steady-streaming flow direction is non-axial.

Understanding the Eulerian and Lagrangian steady streaming velocities of CSF can reveal key factors influencing drug transport. For instance, identifying regions with minimal steady streaming such as separation lines can aid in optimizing drug delivery procedures. Additionally, the Eulerian cyclic mean velocities are used as frozen velocity field in [28] when simulating drug transport that requires computations of oscillatory flow fields across thousands of cycles.

Conclusion

We proposed a mixed-dimensional computational FSI model to simulate the CSF flow dynamics and tissue mechanics in a 3D patient specific model of the spinal canal. Due to the different length scales of dura mater and epidural fat, we discretized dura mater as a 2D membrane, and epidural fat as a 3D solid. To emphasize the critical role of tissue compliance in the CSF dynamics, we compared our FSI simulations with rigid-wall fluid flow simulations. Our FSI simulation revealed several key flow mechanisms that the rigid-wall simulations failed to capture: (1) Phase lag between the flow rates at different sections of the spine, corresponding to a finite pulse wave speed of 3.4 m/s that agrees with in vivo measurements, (2) Craniocaudal decay in flow pulsations with small velocities in the caudal region, (3) Strong and well-defined bi-directional flow pattern of steady-streaming velocities qualitatively compatible with previous experimental [51, 52] and theoretical [32, 53] studies.

Our model could be further developed in several ways. First, incorporating the nerve roots, which alter the details of the CSF flow dynamics. Second, accounting for the compliance of the spinal cord. Considering the deformation of the spinal cord and its mechanical connection to the dura mater through the nerve roots would likely reduce pulse wave speed, and increase the phase lags in flow rate and pressure waveforms. Third, using a more detailed anatomical and physiological description of the epidural fat, including its interactions with surrounding tissues and the presence of the venous and arterial plexus. The arterial and venous plexus affect the CSF motion by transmitting the intra-thoracic or intra-abdominal pressure fluctuations. The effects of arterial/venous plexus could be modeled by using a time-dependent traction boundary condition on the outer wall of the SAS, with experimental pressure data.

We anticipate that our computational method will open new possibilities to better understand structural abnormalities of the CNS. Perhaps more importantly, it can enable control and optimization of intrathecal drug delivery techniques which will lead to more effective therapies for CNS disorders.

Data availability

The data supporting the findings of this study are available within the paper.

References

Scheltens P, Blennow K, Breteler MM, De Strooper B, Frisoni GB, Salloway S, Flier WM. Alzheimer’s disease. Lancet. 2016;388(10043):505–17. https://doi.org/10.1016/S0140-6736(15)01124-1.

Hampel H, Mesulam M-M, Cuello AC, Farlow MR, Giacobini E, Grossberg GT, Khachaturian AS, Vergallo A, Cavedo E, Snyder PJ, et al. The cholinergic system in the pathophysiology and treatment of Alzheimer’s disease. Brain. 2018;141(7):1917–33. https://doi.org/10.1093/brain/awy132.

Pardridge WM. The blood-brain barrier: bottleneck in brain drug development. NeuroRx. 2005;2:3–14. https://doi.org/10.1602/neurorx.2.1.3.

De Andres J, Hayek S, Perruchoud C, Lawrence MM, Reina MA, De Andres-Serrano C, Rubio-Haro R, Hunt M, Yaksh TL. Intrathecal drug delivery: advances and applications in the management of chronic pain patient. Front Pain Res. 2022;3:900566. https://doi.org/10.3389/fpain.2022.900566.

Kakkis E, McEntee M, Vogler C, Le S, Levy B, Belichenko P, Mobley W, Dickson P, Hanson S, Passage M. Intrathecal enzyme replacement therapy reduces lysosomal storage in the brain and meninges of the canine model of MPS I. Mol Genet Metab. 2004;83:163–74. https://doi.org/10.1016/j.ymgme.2004.07.003.

Papisov MI, Belov VV, Gannon KS. Physiology of the intrathecal bolus: the leptomeningeal route for macromolecule and particle delivery to CNS. Mol Pharm. 2013;10:1522–32. https://doi.org/10.1021/mp300474m.

Pilcher C, Meacham WF. The chemotherapy of intracranial infections: iii. The treatment of experimental staphylococcic meningitis with intrathecal administration of penicillin. J Am Med Assoc. 1943;123:330–2. https://doi.org/10.1001/jama.1943.02840410012004.

Wagshul ME, Eide PK, Madsen JR. The pulsating brain: a review of experimental and clinical studies of intracranial pulsatility. Fluids Barriers CNS. 2011;8:1–23. https://doi.org/10.1186/2045-8118-8-5.

Lindstrøm EK, Ringstad G, Mardal K-A, Eide PK. Cerebrospinal fluid volumetric net flow rate and direction in idiopathic normal pressure hydrocephalus. NeuroImage Clin. 2018;20:731–41. https://doi.org/10.1016/j.nicl.2018.09.006.

Schubert JJ, Veronese M, Marchitelli L, Bodini B, Tonietto M, Stankoff B, Brooks DJ, Bertoldo A, Edison P, Turkheimer FE. Dynamic 11C-PiB PET shows cerebrospinal fluid flow alterations in Alzheimer disease and multiple sclerosis. J Nucl Med. 2019;60:1452–60. https://doi.org/10.2967/jnumed.118.223834.

Yiallourou TI, Kröger JR, Stergiopulos N, Maintz D, Martin BA, Bunck AC. Comparison of 4D phase-contrast MRI flow measurements to computational fluid dynamics simulations of cerebrospinal fluid motion in the cervical spine. PLoS ONE. 2012;7:52284. https://doi.org/10.1371/journal.pone.0052284.

Rutkowska G, Haughton V, Linge S, Mardal K-A. Patient-specific 3D simulation of cyclic CSF flow at the craniocervical region. Am J Neuroradiol. 2012;33:1756–62. https://doi.org/10.3174/ajnr.a3047.

Roldan A, Wieben O, Haughton V, Osswald T, Chesler N. Characterization of CSF Hydrodynamics in the presence and absence of tonsillar ectopia by means of computational flow analysis. Am J Neuroradiol. 2009;30:941–6. https://doi.org/10.3174/ajnr.a1489.

Linge S, Haughton V, Løvgren AE, Mardal K-A, Helgeland A, Langtangen HP. Effect of tonsillar herniation on cyclic CSF flow studied with computational flow analysis. Am J Neuroradiol. 2011;32:1474–81. https://doi.org/10.3174/ajnr.a2496.

Linge SO, Mardal K-A, Helgeland A, Heiss JD, Haughton V. Effect of craniovertebral decompression on CSF dynamics in Chiari malformation type I studied with computational fluid dynamics. J Neurosurg Spine. 2014;21:559–64. https://doi.org/10.3171/2014.6.spine13950.

Linge S, Mardal K-A, Haughton V, Helgeland A. Simulating CSF flow dynamics in the normal and the Chiari I subarachnoid space during rest and exertion. Am J Neuroradiol. 2013;34:41–5. https://doi.org/10.3174/ajnr.a3282.

Helgeland A, Mardal K-A, Haughton V, Reif B. Numerical simulations of the pulsating flow of cerebrospinal fluid flow in the cervical spinal canal of a Chiari patient. J Biomech. 2014;47:1082–90. https://doi.org/10.1016/j.jbiomech.2013.12.023.

Stockman HW. Effect of anatomical fine structure on the flow of cerebrospinal fluid in the spinal subarachnoid space. J Biomech Eng. 2006;128(1):106–14. https://doi.org/10.1115/1.2132372.

Stockman HW. Effect of anatomical fine structure on the dispersion of solutes in the spinal subarachnoid space. J Biomech Eng. 2007;129(5):666–75. https://doi.org/10.1115/1.2768112.

Gupta S, Soellinger M, Boesiger P, Poulikakos D, Kurtcuoglu V. Three-dimensional computational modeling of subject-specific cerebrospinal fluid flow in the subarachnoid space. J Biomech Eng. 2009;131:021010. https://doi.org/10.1115/1.3005171.

Gupta S, Soellinger M, Grzybowski DM, Boesiger P, Biddiscombe J, Poulikakos D, Kurtcuoglu V. Cerebrospinal fluid dynamics in the human cranial subarachnoid space: an overlooked mediator of cerebral disease. I Comput Model J Royal Soc Interface. 2010;7:1195–204. https://doi.org/10.1098/rsif.2010.0033.

Pahlavian SH, Yiallourou T, Tubbs RS, Bunck AC, Loth F, Goodin M, Raisee M, Martin BA. The impact of spinal cord nerve roots and denticulate ligaments on cerebrospinal fluid dynamics in the cervical spine. PLoS ONE. 2014;9:91888. https://doi.org/10.1371/journal.pone.0091888.

Gholampour S. FSI simulation of CSF hydrodynamic changes in a large population of non-communicating hydrocephalus patients during treatment process with regard to their clinical symptoms. PLoS ONE. 2018;13(4):0196216. https://doi.org/10.1371/journal.pone.0196216.

Gholampour S, Fatouraee N. Boundary conditions investigation to improve computer simulation of cerebrospinal fluid dynamics in hydrocephalus patients. Commun Biol. 2021;4(1):394. https://doi.org/10.1038/s42003-021-01920-w.

Gholampour S, Yamini B, Droessler J, Frim D. A new definition for intracranial compliance to evaluate adult hydrocephalus after shunting. Front Bioeng Biotechnol. 2022;10:900644. https://doi.org/10.3389/fbioe.2022.900644.

Causemann M, Vinje V, Rognes ME. Human intracranial pulsatility during the cardiac cycle: a computational modelling framework. Fluids Barriers CNS. 2022;19(1):84. https://doi.org/10.1186/s12987-022-00376-2.

Heidari Pahlavian S, Bunck AC, Loth F, Shane Tubbs R, Yiallourou T, Robert Kroeger J, Heindel W, Martin BA. Characterization of the discrepancies between four-dimensional phase-contrast magnetic resonance imaging and in-silico simulations of cerebrospinal fluid dynamics. J Biomech Eng. 2015;137(5):051002. https://doi.org/10.1115/1.4029699.

Kuttler A, Dimke T, Kern S, Helmlinger G, Stanski D, Finelli LA. Understanding pharmacokinetics using realistic computational models of fluid dynamics: biosimulation of drug distribution within the CSF space for intrathecal drugs. J Pharmacokinet Pharmacodyn. 2010;37:629–44. https://doi.org/10.1007/s10928-010-9184-y.

Tangen KM, Hsu Y, Zhu DC, Linninger AA. CNS wide simulation of flow resistance and drug transport due to spinal microanatomy. J Biomech. 2015;48(10):2144–54. https://doi.org/10.1016/j.jbiomech.2015.02.018.

Khani M, Xing T, Gibbs C, Oshinski JN, Stewart GR, Zeller JR, Martin BA. Nonuniform moving boundary method for computational fluid dynamics simulation of intrathecal cerebrospinal flow distribution in a cynomolgus monkey. J Biomech Eng Trans Asme. 2017;139:081005. https://doi.org/10.1115/1.4036608.

Cheng S, Fletcher D, Hemley S, Stoodley M, Bilston L. Effects of fluid structure interaction in a three dimensional model of the spinal subarachnoid space. J Biomech. 2014;47:2826–30. https://doi.org/10.1016/j.jbiomech.2014.04.027.

Sanchez AL, Martinez-Bazan C, Gutierrez-Montes C, Criado-Hidalgo E, Pawlak G, Bradley W, Haughton V, Lasheras JC. On the bulk motion of the cerebrospinal fluid in the spinal canal. J Fluid Mech. 2018;841:203–27. https://doi.org/10.1017/jfm.2018.67.

Cavelier S, Quarrington RD, Jones CF. Tensile properties of human spinal dura mater and pericranium. J Mater Sci Mater Med. 2022;34:4. https://doi.org/10.1007/s10856-022-06704-0.

De Andrés J, Reina MA, Machés F, De Sola RG, Oliva A, Prats-Galino A. Epidural fat: considerations for minimally invasive spinal injection and surgical therapies. J Neurosurg Rev. 2011;1:45–53.

Sass LR, Khani M, Natividad GC, Tubbs RS, Baledent O, Martin BA. A 3D subject-specific model of the spinal subarachnoid space with anatomically realistic ventral and dorsal spinal cord nerve rootlets. Fluids Barriers CNS. 2017;14:36. https://doi.org/10.1186/s12987-017-0085-y.

Simo JC, Hughes TJ. Computational inelasticity vol. 7. New York: Springer; 2006. https://doi.org/10.1007/b98904.

Tepole AB, Kabaria H, Bletzinger K-U, Kuhl E. Isogeometric Kirchhoff-love shell formulations for biological membranes. Comput Methods Appl Mech Eng. 2015;293:328–47. https://doi.org/10.1016/j.cma.2015.05.006.

Hsu M-C, Bazilevs Y. Blood vessel tissue prestress modeling for vascular fluid-structure interaction simulation. Finite Elem Anal Des. 2011;47:593–9. https://doi.org/10.1016/j.finel.2010.12.015.

Rangel-Castillo L, Gopinath S, Robertson CS. Management of intracranial hypertension. Neurol Clin. 2008;26(2):521–41. https://doi.org/10.1016/j.ncl.2008.02.003.

Moireau P, Xiao N, Astorino M, Figueroa CA, Chapelle D, Taylor CA, Gerbeau J-F. External tissue support and fluid-structure simulation in blood flows. Biomech Model Mechanobiol. 2012;11:1–18. https://doi.org/10.1007/s10237-011-0289-z.

Lui AC, Polis TZ, Cicutti NJ. Densities of cerebrospinal fluid and spinal anaesthetic solutions in surgical patients at body temperature. Can J Anaesth. 1998;45:297–303. https://doi.org/10.1007/bf03012018.

Bloomfield I, Johnston I, Bilston L. Effects of proteins, blood cells and glucose on the viscosity of cerebrospinal fluid. Pediatr Neurosurg. 1998;28:246–51. https://doi.org/10.1159/000028659.

Runza M, Pietrabissa R, Mantero S, Albani A, Quaglini V, Contro R, Camann WR. Lumbar dura mater biomechanics: experimental characterization and scanning electron microscopy observations. Obstet Anesth Dig. 1999;19:141–56. https://doi.org/10.1097/00000539-199906000-00022.

Shoham N, Gefen A. The biomechanics of fat: from tissue to a cell scale. Boston: Springer; 2016. p. 79–92.

Sonnabend K, Brinker G, Maintz D, Bunck AC, Weiss K. Cerebrospinal fluid pulse wave velocity measurements: in vitro and in vivo evaluation of a novel multiband cine phase-contrast MRI sequence. Magn Reson Med. 2021;85:197–208. https://doi.org/10.1002/mrm.28430.

Kalata W, Martin BA, Oshinski JN, Jerosch-Herold M, Royston TJ, Loth F. MR measurement of cerebrospinal fluid velocity wave speed in the spinal canal. IEEE Trans Biomed Eng. 2009;56:1765–8. https://doi.org/10.1109/tbme.2008.2011647.

Chun S-W, Lee H-J, Nam K-H, Sohn C-H, Kim KD, Jeong E-J, Chung SG, Kim K, Kim D-J. Cerebrospinal fluid dynamics at the lumbosacral level in patients with spinal stenosis: a pilot study. J Orthop Res. 2017;35:104–12. https://doi.org/10.1002/jor.23448.

Tumani H, Petereit H, Gerritzen A, Gross C, Huss A, Isenmann S, Jesse S, Khalil M, Lewczuk P, Lewerenz J, et al. S1 guidelines “lumbar puncture and cerebrospinal fluid analysis’’(abridged and Translated Version). Neurol Res Pract. 2020;2:1–28. https://doi.org/10.1186/s42466-020-0051-z.

Wang D, Tarbell J. Nonlinear analysis of flow in an elastic tube (artery): steady streaming effects. J Fluid Mech. 1992;239:341–58. https://doi.org/10.1017/S0022112092004439.

Riley N. Steady streaming. Annu Rev Fluid Mech. 2001;33:43–65. https://doi.org/10.1146/annurev.fluid.33.1.43.

Di Chiro G. Observations on the circulation of the cerebrospinal fluid. Acta Radiol Diagn. 1966;5:988–1002. https://doi.org/10.1177/02841851660050P242.

Greitz D, Hannerz J. A proposed model of cerebrospinal fluid circulation: observations with radionuclide cisternography. Am J Neuroradiol. 1996;17(3):431–8.

Coenen W, Gutierrez-Montes C, Sincomb S, Criado-Hidalgo E, Wei K, King K, Haughton V, Martinez-Bazan C, Sanchez AL, Lasheras JC. Subject-specific studies of CSF bulk flow patterns in the spinal canal: implications for the dispersion of solute particles in intrathecal drug delivery. Am J Neuroradiol. 2019;40:1242–9. https://doi.org/10.3174/ajnr.a6097.

Sullivan CB, Kaszynski A. PyVista: 3d plotting and mesh analysis through a streamlined interface for the visualization toolkit (VTK). J Open Sour Softw. 2019;4:1450. https://doi.org/10.21105/joss.01450.

Anderson JD. Governing equations of fluid and structural mechanics. John Wiley & Sons, Ltd.; 2013. p. 1–35.

Donea J, Huerta A, Ponthot J-Ph, Rodríguez-Ferran A. Arbitrary Lagrangian-Eulerian methods. Encyclopedia of computational mechanics. Wiley; 2004.

Bazilevs Y, Calo VM, Hughes T, Zhang Y. Isogeometric fluid-structure interaction: theory, algorithms, and computations. Comput Mech. 2008;43:3–37. https://doi.org/10.1007/s00466-008-0315-x.

Johnson AA, Tezduyar TE. Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces. Comput Methods Appl Mech Eng. 1994;119:73–94. https://doi.org/10.1016/0045-7825(94)00077-8.

Stein K, Tezduyar T, Benney R. Mesh moving techniques for fluid-structure interactions with large displacements. J Appl Mech. 2003;70:58–63. https://doi.org/10.1115/1.1530635.

Geuzaine C, Remacle J-F. Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities. Int J Numer Meth Eng. 2009;79:1309–31. https://doi.org/10.1002/nme.2579.

Alnæs M, Blechta J, Hake J, Johansson A, Kehlet B, Logg A, Richardson C, Ring J, Rognes ME, Wells GN. The fenics project version 1.5. Arch Numer Softw. 2015. https://doi.org/10.11588/ans.2015.100.20553.

Logg A, Mardal K-A, Wells G. Automated solution of differential equations by the finite element method: the FEniCS book vol. 84. Berlin: Springer; 2012.