Hostname: page-component-8448b6f56d-tj2md Total loading time: 0 Render date: 2024-04-18T05:47:33.698Z Has data issue: false hasContentIssue false

A one-dimensional model for the pulsating flow of cerebrospinal fluid in the spinal canal

Published online by Cambridge University Press:  30 March 2022

S. Sincomb
Department of Mechanical and Aerospace Engineering, UC San Diego, La Jolla, CA 92093-0411, USA
W. Coenen
Grupo de Mecánica de Fluidos, Universidad Carlos III de Madrid, Leganés (Madrid) 28911, Spain
C. Gutiérrez-Montes
Grupo de Mecánica de Fluidos, Universidad de Jaén, Jaén 23071, Spain
C. Martínez Bazán
Grupo de Mecánica de Fluidos, Universidad de Granada, Granada 18071, Spain
V. Haughton
Department of Radiology, University of Wisconsin-Madison, Madison, WI 53792-3252, USA
A.L. Sánchez*
Department of Mechanical and Aerospace Engineering, UC San Diego, La Jolla, CA 92093-0411, USA
Email address for correspondence:


The monitoring of intracranial pressure (ICP) fluctuations, which is needed in the context of a number of neurological diseases, requires the insertion of pressure sensors, an invasive procedure with considerable risk factors. Intracranial pressure fluctuations drive the wave-like pulsatile motion of cerebrospinal fluid (CSF) along the compliant spinal canal. Systematically derived simplified models relating the ICP fluctuations with the resulting CSF flow rate can be useful in enabling indirect evaluations of the former from non-invasive magnetic resonance imaging (MRI) measurements of the latter. As a preliminary step in enabling these predictive efforts, a model is developed here for the pulsating viscous motion of CSF in the spinal canal, assumed to be a linearly elastic compliant tube of slowly varying section, with a Darcy pressure-loss term included to model the fluid resistance introduced by the trabeculae, which are thin collagen-reinforced columns that form a web-like structure stretching across the spinal canal. Use of Fourier-series expansions enables predictions of CSF flow rate for realistic anharmonic ICP fluctuations. The flow rate predicted using a representative ICP waveform together with a realistic canal anatomy is seen to compare favourably with in vivo phase-contrast MRI measurements at multiple sections along the spinal canal. The results indicate that the proposed model, involving a limited number of parameters, can serve as a basis for future quantitative analyses targeting predictions of ICP temporal fluctuations based on MRI measurements of spinal-canal anatomy and CSF flow rate.

JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (, which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
© The Author(s), 2022. Published by Cambridge University Press

1. Introduction

The cerebrospinal fluid (CSF) is a clear fluid that fills the ventricles of the brain as well as the subarachnoid spaces (SSASs) (Linninger et al. Reference Linninger, Tangen, Hsu and Frim2016). Normal CSF behaves as a Newtonian fluid (Ommaya Reference Ommaya1968; Bloomfield, Johnston & Bilston Reference Bloomfield, Johnston and Bilston1998) and its properties are very close to those of water (i.e. density $\rho =10^3$ kg m$^{-3}$ and kinematic viscosity $\nu = 0.7 \times 10^{-6}$ m$^2$ s$^{-1}$). The motion of CSF in the central nervous system, which has important physiological functions and plays a role in the development of different neurological diseases, has been the subject of numerous studies, as reviewed by Linninger et al. (Reference Linninger, Tangen, Hsu and Frim2016). Recent efforts include numerical simulations of the entire cranial cavity (Gholampour & Fatouraee Reference Gholampour and Fatouraee2021) and investigations of flow in the perivascular spaces of cerebral arteries (Thomas Reference Thomas2019; Carr et al. Reference Carr, Thomas, Liu and Shang2021; Coenen, Zhang & Sánchez Reference Coenen, Zhang and Sánchez2021) and along the cerebral aqueduct (Sincomb et al. Reference Sincomb, Coenen, Sánchez and Lasheras2020, Reference Sincomb, Coenen, Criado-Hidalgo, Wei, King, Borzage, Haughton, Sánchez and Lasheras2021), for example. The present paper deals with the motion along the SSAS, a slender compliant canal of length $L\simeq 60~\rm {cm}$ bounded internally by the pia mater surrounding the spinal cord and externally by the deformable dura membrane (see figure 1a). The arterial blood flow in and out of the rigid cranial vault causes the intracranial pressure (ICP) to fluctuate in time following the cardiac cycle (Du Boulay Reference Du Boulay1966; Bhadelia et al. Reference Bhadelia, Bogdan, Kaplan and Wolpert1997; Wagshul et al. Reference Wagshul, Chen, Egnor, McCormack and Roche2006), driving the pulsatile motion of CSF along the SSAS.

Figure 1. (a) Main anatomical features of the spinal canal for subject 1; (b) ICP wave form (Wagshul, Eide & Madsen Reference Wagshul, Eide and Madsen2011) CC BY 2.0 (left subpanel) and dimensionless (right subpanel); (c) dimensionless canal functions for subjects 1 and 2.

Continuous ICP monitoring is key in the assessment of surgical intervention and also for guiding therapy for cases of traumatic brain injury (TBI), normal pressure hydrocephalus (NPH) and other neurointensive states. Although the mean ICP value is often used clinically, it is of interest to also assess the pulsatile ICP variation or morphology. As shown in figure 1(b), the ICP waveform generally has three peaks associated with the cardiac cycle whose amplitudes decrease in a stepwise manner in a healthy individual (Singh & Cheng Reference Singh and Cheng2021). The waveform is altered due to the changes in intracranial volume (Unnerbäck, Ottesen & Reinstrup Reference Unnerbäck, Ottesen and Reinstrup2018), an important result in the context of disease conditions that produce an increase in mean ICP (TBI or oedema formation), which can result in the waveform becoming more rounded (Ellis, McNames & Aboy Reference Ellis, McNames and Aboy2007), or NPH, which leads to greater fluctuation amplitudes (Eide & Sorteberg Reference Eide and Sorteberg2016). The insertion of ICP sensors requires a burr hole made into the skull. Since the procedure has inherent risks, including haemorrhage and infection (Evensen & Eide Reference Evensen and Eide2020), there is interest in developing non-invasive techniques for ICP characterization. The approach postulated here exploits the close connection existing between the ICP waveform and the resulting CSF motion along the spinal canal. It is reasoned that detailed knowledge of the flow rate along the canal, obtained via phase-contrast magnetic resonance imaging (PC-MRI) techniques (Feinberg & Mark Reference Feinberg and Mark1987), can be used to infer the associated ICP waveform.

The flow in the canal fundamentally involves a fluid–structure interaction problem, which depends on detailed anatomical features of the canal determining its compliance and flow resistance (Linninger et al. Reference Linninger, Tangen, Hsu and Frim2016). As summarized by Khani et al. (Reference Khani, Sass, Xing, Keith Sharp, Balédent and Martin2018), most previous modelling efforts are based on numerical simulations with different levels of complexity. Analytic flow models involving a reduced number of parameters can be more useful in enabling inverse predictions of ICP from measurements of flow rates. One-dimensional models for pressure/flow wave propagation along the spinal canal have been developed in the past using a coaxial cylindrical tube configuration (Berkouk, Carpenter & Lucey Reference Berkouk, Carpenter and Lucey2003; Carpenter, Berkouk & Lucey Reference Carpenter, Berkouk and Lucey2003; Cirovic & Kim Reference Cirovic and Kim2012). More elaborate three-dimensional flow models assuming a thin annular canal of non-uniform width are also available (Sánchez et al. Reference Sánchez, Martínez-Bazan, Gutiérrez-Montes, Criado-Hidalgo, Pawlak, Bradley, Haughton and Lasheras2018; Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019; Gutiérrez-Montes et al. Reference Gutiérrez-Montes, Coenen, Lawrence, Martínez-Bazán, Sánchez and Lasheras2021). These previous efforts have assumed the canal section to be open, thereby neglecting the pressure loss introduced by spinal microstructures, effects of which have been quantified numerically by Tangen et al. (Reference Tangen, Hsu, Zhu and Linninger2015). Their analysis showed that most of the increase in pressure loss is associated with the arachnoid trabeculae, which are thin collagen-reinforced columns that form a web-like structure stretching across the SSAS (Mortazavi et al. Reference Mortazavi2018). Following Gupta et al. (Reference Gupta, Soellinger, Boesiger, Poulikakos and Kurtcuoglu2009), our analysis will model the complex trabeculae network as a porous medium of variable permeability. For increased generality, no specific shape will be assumed for the canal cross-section, thereby generalizing our previous analyses (Sánchez et al. Reference Sánchez, Martínez-Bazan, Gutiérrez-Montes, Criado-Hidalgo, Pawlak, Bradley, Haughton and Lasheras2018; Lawrence et al. Reference Lawrence, Coenen, Sánchez, Pawlak, Martínez-Bazán, Haughton and Lasheras2019; Gutiérrez-Montes et al. Reference Gutiérrez-Montes, Coenen, Lawrence, Martínez-Bazán, Sánchez and Lasheras2021), postulating the SSAS to be a thin annular canal surrounding the spinal cord, an assumption that necessarily fails in the sacral region, as shown in the cross-sectional views of figure 1(a).

2. Preliminary considerations

2.1. The ICP

Attention will be focused on the motion induced by the cardiac cycle, associated with the periodic temporal fluctuations of the ICP $p_c(t)$ from its mean value $\langle p_c \rangle =T^{-1} \int _t^{t+T} p_c\, {\rm d}t$, where $T \simeq 1$ s is the period of the cardiac cycle. These fluctuations can in general be represented in the form $p_c(t)-\langle p_c \rangle =\Delta p \varPi (\omega t)$, involving the mean fluctuation amplitude $\Delta p=T^{-1} \int _t^{t+T} |p_c-\langle p_c \rangle | \,{\rm d}t$, which typically takes values of the order of $\Delta p \sim 50\text {--}100$ Pa, along with a dimensionless function $\varPi (\omega t)$ describing the waveform, with $\omega = 2 {\rm \pi}/T$ denoting the relevant angular frequency. Note that the function $\varPi$ must satisfy $\int _t^{t+T} \varPi \,{\rm d}t=0$ and $T^{-1} \int _t^{t+T} |\varPi |\,{\rm d}t= 1$, for consistency with the definition of $\Delta p$.

2.2. The canal geometry

In deriving a simple one-dimensional model for the flow dynamics, the spinal canal will be modelled as a tube displaying a slowly varying shape over its length $L \simeq 50\text {--}70$ cm. The flow is to be described in terms of curvilinear coordinates $(x,y,z)$, with the streamwise distance $x$ measured from the open end, connected to the cranial cavity through the foramen magnum, with the closed sacral end corresponding to $x=L$, as indicated in figure 1(a). As seen in figure 1(a), in the stretch of canal occupied by the spinal cord the cross-sectional shape is an annulus, bounded internally by the pia mater and externally by the dura mater. In the one-dimensional model developed below the morphology of the cross-section enters only through two related quantities that vary along the spinal canal, namely, the cross-sectional area occupied by CSF at each transverse section $A(x)$ and the length of the wetted boundary $\ell (x)$, the latter including the pia mater surrounding the spinal cord. The average cross-sectional area is given by $A_o=V_{CSF}/L \simeq 1.5$ cm$^2$, where $V_{CSF}=\int _0^L A \,{\rm d} x \simeq 80$ cm$^{3}$ is the total volume of CSF in the in the SSAS (Edsbagge et al. Reference Edsbagge, Starck, Zetterberg, Ziegelitz and Wikkelso2011). Since the characteristic transverse length $A_o^{1/2}$ satisfies $A_o^{1/2} \ll L$, the flow is fundamentally slender. Fluid motion predominantly occurs in the axial direction, with streamwise velocity $u(x,y,z,t)$ driven by the streamwise pressure distribution $p(x,t)$ (Sánchez et al. Reference Sánchez, Martínez-Bazan, Gutiérrez-Montes, Criado-Hidalgo, Pawlak, Bradley, Haughton and Lasheras2018), assumed to be uniform across the canal section, as is consistent with the slender-flow approximation.

2.3. Governing equations

During each cardiac cycle, the ICP pulsation drives a small volume $\Delta V \sim 1$ cm$^{3}$ of CSF in and out of the spinal canal (Linninger et al. Reference Linninger, Tangen, Hsu and Frim2016). This oscillating flow is accommodated by the displacement of fat tissue and venous blood, which results in a periodic change $\Delta A$ of the local cross-sectional area $A$ at a given location. Since the characteristic stroke volume $\Delta V$ is much smaller than the total CSF volume $V_{CSF}$, these temporal changes are small, i.e. $\Delta A \sim (\Delta V/V_{CSF}) A_o\sim 1$ mm$^2$. To model these changes, we shall adopt the linear elastic model,

(2.1)\begin{equation} \frac{\partial A}{\partial t}=\gamma \frac{\partial p}{\partial t}, \end{equation}

involving a local compliance $\gamma (x)$ having dimensions of surface over pressure with mean value $\gamma _o =L^{-1} \int _0^L \gamma (x) \,{\rm d} x$.

To maximize the simplicity and facilitate comparisons with in vivo results, CSF motion is to be characterized with use of the local volumetric flow rate $Q(x,t)=\iint u \, {\rm d} y\,{\rm d}z$, obtained by integrating the streamwise velocity $u(x,y,z,t)$ across the canal section. Its streamwise variation is related with the temporal variation of the cross-sectional area through the integrated continuity equation ${\partial Q}/{\partial x}+{\partial A}/{\partial t}=0$, which can be rewritten in the form

(2.2)\begin{equation} \frac{\partial Q}{\partial x}+\gamma(x) \frac{\partial p}{\partial t}=0, \end{equation}

after substitution of (2.1). On the other hand, the axial component of the momentum balance equation can be simplified by neglecting convective acceleration, whose magnitude can be shown to be a factor $\Delta V/V_{CSF}$ smaller than that of the local acceleration (Sánchez et al. Reference Sánchez, Martínez-Bazan, Gutiérrez-Montes, Criado-Hidalgo, Pawlak, Bradley, Haughton and Lasheras2018), along with the contribution of streamwise derivatives to the viscous force. Furthermore, following Gupta et al. (Reference Gupta, Soellinger, Boesiger, Poulikakos and Kurtcuoglu2009), the pressure loss caused by the trabeculae network is modelled using Darcy's law, yielding (Kurtcuoglu, Jain & Martin Reference Kurtcuoglu, Jain and Martin2019)

(2.3)\begin{equation} \frac{\partial {u}}{\partial t}={-}\frac{1}{\rho} \frac{\partial {p}}{\partial x}+\nu \left(\frac{\partial^2 {u}}{\partial {y}^2}+\frac{\partial^2 {u}}{\partial {z}^2} \right)-\frac{\nu}{\kappa} {u}, \end{equation}

where $\kappa (x)$ is the SSAS permeability, whose value depends on the number and structure of the arachnoid trabeculae.

2.4. The inviscid wave model

It is illustrative to consider first the inviscid case $\nu =0$, for which integration of (2.3) across the canal yields

(2.4)\begin{equation} \frac{\partial Q}{\partial t}+\frac{A(x)}{\rho} \frac{\partial p}{\partial x}=0. \end{equation}

In writing the above equation, we have neglected the small temporal variation of the cross-sectional area $A$, as is consistent with the condition $\Delta A \ll A_o$ previously discussed. Equations (2.2) and (2.4) can be integrated with boundary conditions $p=p_c$ at $x=0$ and $Q=0$ at $x=L$ to determine the periodic variation of the pressure and flow rate along the canal. The wave nature of the flow can be emphasized by considering a canal with constant section $A(x)=A_o$ and constant compliance $\gamma (x)=\gamma _o$, for which (2.2) and (2.4) can be combined to give

(2.5)\begin{equation} \frac{\partial^2 p}{\partial t^2}=c^2 \frac{\partial^2 p}{\partial x^2} \quad {\rm and} \quad \frac{\partial^2 Q}{\partial t^2}=c^2 \frac{\partial^2 Q}{\partial x^2}, \end{equation}


(2.6)\begin{equation} c=\left(\frac{A_o}{\rho \gamma_o} \right)^{1/2} \end{equation}

is the elastic wave speed of the problem (Grotberg & Jensen Reference Grotberg and Jensen2004). For a harmonic ICP fluctuation $p_c-\langle p_c \rangle =\Delta p ({\rm \pi} /2) \cos (\omega t)$, the above wave equation can be solved to give

(2.7)\begin{equation} Q={-}\frac{\rm \pi}{2} \gamma_o L \omega \Delta p \sin(\omega t) \frac{\sin[k(1-x/L)]}{k \cos k}, \end{equation}


(2.8)\begin{equation} k=\frac{L\omega}{c}=\left(\frac{\rho \gamma_o L^2 \omega^2}{A_o}\right)^{1/2} \end{equation}

is a relevant dimensionless wavenumber. The flow rate (2.7) oscillates in phase along the entire canal, that being a fundamental limitation of the inviscid model, which is unable to reproduce the streamwise phase lag of the flow rate that has been consistently observed in MRI measurements (Yallapragada & Alperin Reference Yallapragada and Alperin2004; Wagshul et al. Reference Wagshul, Chen, Egnor, McCormack and Roche2006; Tangen et al. Reference Tangen, Hsu, Zhu and Linninger2015). As shown below, consideration of viscous pressure losses, including those associated with the trabeculae, is needed to describe both the phase lag and the rate of streamwise attenuation of the flow rate.

3. A dimensionless flow model accounting for flow resistance

To reduce the parametric dependence, it is convenient to formulate the problem in dimensionless form using $\omega ^{-1}$, $L$ and $A_o^{1/2}$ as scales for the time and for the longitudinal and transverse length scales, respectively. In looking for appropriate scales for $Q$ and $p$, one may note from (2.2) that the characteristic value $Q_c$ of the volume flux associated with an ICP fluctuation of magnitude $\Delta p$ is $Q_c=\gamma _o \omega L \Delta p$, and from (2.4) that the corresponding streamwise variations of the pressure are of order $p-p_c \sim k^2 \Delta p$, with $k$ denoting the wavenumber defined in (2.8). These scales lead to the new variables

(3.1)\begin{equation} \left.\begin{gathered} \tau=\omega t, \quad \xi=\frac{x}{L}, \quad \hat{y}=\frac{y}{A_o^{1/2}}, \quad \hat{z}=\frac{z}{A_o^{1/2}},\\ \hat{u}=\frac{u}{Q_c/A_o}, \quad \hat{Q}=\frac{Q}{Q_c}, \quad \hat{p}=\frac{p-p_c}{k^2 \Delta p}. \end{gathered}\right\} \end{equation}

Similarly, $A_o$ and $\gamma _o$ are used to define the functions $\hat {A}=A/A_o$ and $\hat {\gamma }=\gamma /\gamma _o$.

The development begins by writing the continuity equation (2.2) in the reduced form

(3.2)\begin{equation} \frac{\partial \hat{Q}}{\partial \xi}+\hat{\gamma} \left(\frac{{\rm d} \varPi}{{\rm d} \tau}+k^2 \frac{\partial \hat{p}}{\partial \tau} \right)=0. \end{equation}

With the scales selected, the momentum equation (2.3) takes the dimensionless form

(3.3)\begin{equation} \frac{\partial \hat{u}}{\partial \tau}={-}\frac{\partial \hat{p}}{\partial \xi}+\frac{1}{\alpha^2} \left(\frac{\partial^2 \hat{u}}{\partial \hat{y}^2}+\frac{\partial^2 \hat{u}}{\partial \hat{z}^2} \right)-\mathcal{R}(\xi) \hat{u}, \end{equation}

where $\alpha =(A_o \omega /\nu )^{1/2}$ is the relevant Womersley number and $\mathcal {R}(\xi )=\nu /(\kappa \omega )$ is a dimensionless resistance coefficient. The velocity must satisfy the no-slip condition $\hat {u}=0$ on the canal boundary $\varSigma$. Integrating the above equation across the canal section yields

(3.4)\begin{equation} \frac{\partial \hat{Q}}{\partial \tau}+ \hat{A}(\xi) \frac{\partial \hat{p}}{\partial \xi}={-}\frac{1}{\alpha^2} \int_{\varSigma} \hat{\tau}_f \,{\rm d}s-\mathcal{R}(\xi) \hat{Q}, \end{equation}

where $\textrm {d}s$ the element of arclength measured at a given section along the canal boundary $\varSigma$ and $\hat {\tau }_f=\partial \hat {u}/\partial n$ is the dimensionless viscous stress at $n=0$, where $n$ denotes the dimensionless distance from the wall.

3.1. Simplifications for $\alpha \gg 1$

The solution can be simplified by taking into account that the characteristic viscous time across the canal section $A_o/\nu$ is fairly large compared with the characteristic pulsation time $\omega ^{-1}$. In the associated limit $\alpha \gg 1$, the longitudinal velocity is uniform outside a thin near-wall Stokes layer of rescaled thickness $\alpha ^{-1} \ll 1$. The uniform velocity in the inviscid core varies along the canal according to $\hat {u}=\hat {Q}(\xi,\tau )/\hat {A}(\xi )$, while the accompanying pressure gradient is $\partial \hat {p}/\partial \xi =-(\partial \hat {Q}/\partial \tau + \mathcal {R} \hat {Q})/\hat {A}$. Viscous forces are important in the Stokes layer, across which the velocity $\hat {u}_{S}$ evolves from the inviscid value $\hat {Q}/\hat {A}$ to a zero value at $n=0$, as described by the reduced problem

(3.5)\begin{equation} \frac{\partial \hat{u}_{S}}{\partial \tau}-\frac{1}{\hat{A}(\xi)} \frac{\partial \hat{Q}}{\partial \tau}=\frac{\partial^2 \hat{u}_{S}}{\partial \eta^2}+\mathcal{R}(\xi) \left(\frac{\hat{Q}}{\hat{A}(\xi)}-\hat{u}_{S} \right) \left\{\begin{array}{@{}ll} \eta=0: & \hat{u}_{S}=0 \\ \eta \rightarrow \infty: & \hat{u}_{S} \rightarrow \hat{Q}/\hat{A} \end{array} \right., \end{equation}

where $\eta =\alpha n$. The solution to (3.5) determines in particular the value of $\hat {u}'_o(\xi,\tau )={\partial \hat {u}_{S}}/\partial \eta |_{\eta =0}$, which can be used to write (3.4) in the form

(3.6)\begin{equation} \frac{\partial \hat{Q}}{\partial \tau}+ \hat{A}(\xi) \frac{\partial \hat{p}}{\partial \xi}={-}\frac{\hat{\ell}(\xi)}{\alpha} \hat{u}'_o - \mathcal{R}(\xi) \hat{Q}. \end{equation}

As expected, since at leading order in the limit $\alpha \gg 1$ the structure of the Stokes layer is identical all around the canal wall, the term $-\hat {\ell } \hat {u}'_o/\alpha$ representing in (3.6) the viscous force $-\int _{\varSigma } \hat {\tau }_f\, \textrm {d}s/\alpha ^2$ is linearly proportional to the dimensionless length of the wetted boundary $\hat {\ell }(\xi )=\ell /A_o^{1/2}$.

3.2. Solution in terms of Fourier expansions

The problem can be solved for a given general periodic function $\varPi (\tau )=\sum _{n=1}^\infty \textrm {Re}(B_n \,\textrm {e}^{\mathrm {i} n \tau })$, where $\textrm {Re}$ indicates the real part, $\mathrm {i}$ is the imaginary unit, and $B_n$ are complex constants, by introducing accompanying Fourier expansions

(3.7) \begin{equation} \left.\begin{array}{l@{}} \displaystyle\hat{p}=\sum_{n=1}^\infty {\rm Re}(B_n P_n(\xi)\, {\rm e}^{\mathrm{i} n \tau}), \\ \displaystyle\hat{Q}=\sum_{n=1}^\infty {\rm Re}(B_n \mathrm{i} n Q_n(\xi) \,{\rm e}^{\mathrm{i} n \tau}), \\ \displaystyle\hat{u}_{S}=\dfrac{1}{\hat{A}(\xi)}\sum_{n=1}^\infty {\rm Re}(B_n \mathrm{i} n Q_n(\xi) f_n(\eta) \,{\rm e}^{\mathrm{i} n \tau}), \end{array} \right\} \end{equation}

for the pressure, volume flow rate and Stokes-layer velocity, with $P_n(\xi )$, $Q_n(\xi )$ and $f_n(\eta )$ representing complex functions. The function $f_n=1-\exp [-(\mathrm {i} n+\mathcal {R})^{1/2} \eta ]$ is obtained from the reduced Stokes problem

(3.8)\begin{equation} (\mathrm{i} n+\mathcal{R}) (f_n-1)=\frac{{\rm d}^2 f_n}{{\rm d} \eta^2} \left\{\begin{array}{@{}ll} \eta=0: & f_n=0 \\ \eta \rightarrow \infty: & f_n \rightarrow 1 \end{array} \right., \end{equation}

which follows from introduction of (3.7) into (3.5), thereby yielding $\textrm {d}f_n/\textrm {d}\eta (0)=(\mathrm {i} \, n+\mathcal {R})^{1/2}$ and

(3.9)\begin{equation} \hat{u}'_o=\frac{1}{\hat{A}(\xi)} \sum_{n=1}^\infty {\rm Re}[B_n \mathrm{i} n Q_n(\xi) [\mathrm{i} n+\mathcal{R}(\xi)]^{1/2} \,{\rm e}^{\mathrm{i} n \tau}]. \end{equation}

Substituting (3.7) and (3.9) into (3.2) and (3.6) leads to the first-order linear ordinary differential equations

(3.10)$$\begin{gather} \frac{{\rm d} Q_n}{{\rm d} \xi}+\hat{\gamma}(\xi) [1+k^2 P_n]=0, \end{gather}$$
(3.11)$$\begin{gather}\hat{A}(\xi) \frac{{\rm d} P_n}{{\rm d} \xi}=\left\{n^2- \mathrm{i} n \left[\frac{\hat{\ell}(\xi)}{\alpha \hat{A}(\xi)} [\mathrm{i} n+\mathcal{R}(\xi)]^{1/2} + \mathcal{R}(\xi)\right]\right\} Q_n, \end{gather}$$

which can be further combined to generate the boundary-value problem

(3.12)\begin{align} \left.\begin{gathered} \frac{{\rm d}}{{\rm d} \xi}\left[\frac{1}{\hat{\gamma}(\xi)} \frac{{\rm d} Q_n}{{\rm d} \xi}\right]+\frac{k^2}{\hat{A}(\xi)}\left\{n^2- \mathrm{i} n \left[\frac{\hat{\ell}(\xi)}{\alpha \hat{A}(\xi)} [\mathrm{i} n+\mathcal{R}(\xi)]^{1/2} + \mathcal{R}(\xi)\right]\right\} Q_n=0,\\ \frac{{\rm d} Q_n}{{\rm d} \xi}(0)+\hat{\gamma}(0)=Q_n(1)=0, \end{gathered}\right\} \end{align}

for $Q_n(\xi )$, with the value of $\textrm {d} Q_n/\textrm {d} \xi$ at $\xi =0$ following from using $P_n(0)=0$ in (3.10).

It is worth noting that the terms in the large square brackets in (3.11) represent the pressure losses associated with the Stokes layer developing on the canal boundary (the term proportional to $[\mathrm {i} n+\mathcal {R}]^{1/2}$) and with the trabeculae (the term proportional to $\mathcal {R}$). The two pressure losses have different phase, so that they have distinct effects on the resulting flow rate, as can be inferred from (3.12). Since the resistance exerted by the Stokes layer is inversely proportional to $\alpha$, it tends to have a lesser effect, especially on the first Fourier mode $n=1$, for which the pressure loss associated with the trabeculae is significantly higher. Although simplified flow descriptions neglecting the presence of the Stokes layer are worth exploring in future work, for completeness the computations presented below utilize the full equation (3.12) in evaluating the flow rate.

For a canal of uniform section, uniform compliance and uniform permeability (i.e. $\hat {A}=\hat {\gamma }=1$, $\hat {\ell }\,=\,\textrm {constant}$ and $\mathcal {R}=\textrm {constant}$) the solution reduces to $Q_n={\sin [\beta _n (1-\xi )]}/(\beta _n \cos \beta _n)$, with

(3.13)\begin{equation} \beta_n= k n \left\{1+ \left[\frac{1- \mathrm{i}}{(\alpha/\hat{\ell}) \sqrt{2 n}}\left(1-\frac{\mathcal{R}\mathrm{i} }{n}\right)^{1/2}-\frac{\mathcal{R}\mathrm{i} }{n} \right]\right\}^{1/2}. \end{equation}

It is worth noting that the inviscid limit considered earlier corresponds to $\beta _n= k n$, the limiting form of (3.13) for $\alpha \gg 1$ and $\mathcal {R}=0$. For the case $B_1={\rm \pi} /2$ and $B_n=0$ for $n>1$, associated with the harmonic ICP $p_c=\Delta p \ ({\rm \pi} /2) \cos (\omega t)$, the Fourier series for the flow rate reduces to the single term

(3.14)\begin{equation} \hat{Q}={\rm Re}\left(\frac{\rm \pi}{2} \mathrm{i} \frac{\sin[k (1-\xi)]}{k \cos k}\, {\rm e}^{\mathrm{i} \tau}\right)={-}\frac{\rm \pi}{2} \sin \tau \frac{\sin[k (1-\xi)]}{k \cos k}, \end{equation}

consistent with (2.7).

4. Illustrative sample applications

In general, numerical integration of (3.12) is needed to determine $Q_n(\xi )$. The solution depends on the anatomical characteristics of the SSAS, including its length $L$, average cross-sectional area $A_o$ and geometric functions $\hat {A}=A/A_o$ and $\hat {\ell }=\ell /A_o^{1/2}$. For a given subject, the necessary anatomical data can be determined from high-resolution MRI images, as explained for instance in Coenen et al. (Reference Coenen, Gutiérrez-Montes, Sincomb, Criado-Hidalgo, Wei, King, Haughton, Martínez-Bazán, Sánchez and Lasheras2019). The computations reported below correspond to two subjects: a healthy 27-year-old female with $L=60$ cm, $A_o=182$ mm$^2$ and $\alpha =(A_o \omega /\nu )^{1/2} = 40$ (subject 1); and a healthy 39-year-old male with $L=64$ cm, $A_o = 138$ mm$^2$ and $\alpha = 35$ (subject 2), both with cardiac period $T=2{\rm \pi} /\omega =1$ s. The functions $\hat {A}(\xi )$ and $\hat {\ell }(\xi )$ shown in figure 1(b) are obtained by the boundaries of the binary image stack resulting from segmentation of the high resolution MRI images.

The value of $\mathcal {R}=\nu /(\kappa \omega )$ depends on the permeability $\kappa$, which in turn is a function of the SSAS porosity $\epsilon$ and trabeculae transverse size. In the following, we shall assume that $\kappa$ is uniform and adopt the approximate formula $\kappa ={\rm \pi} a^2 \epsilon (1-\sqrt {1-\epsilon })^2/[24(1-\epsilon )^{3/2}]$, derived by Gupta et al. (Reference Gupta, Soellinger, Boesiger, Poulikakos and Kurtcuoglu2009) for a trabeculae network comprising cylindrical posts of radius $a$ extending normally to the arachnoid layer. For a porosity $\epsilon =0.99$, a value estimated by Tada & Nagashima (Reference Tada and Nagashima1994), and a trabeculae radius $a=15$ $\mathrm {\mu }$m (Stockman Reference Stockman2006) this approximate formula gives $\kappa =2.362 \times 10^{-8}$ m$^2$, corresponding to a dimensionless resistance factor $\mathcal {R} = 4.78$, to be used below.

Limited information is available on the spinal canal compliance $\gamma$ (Tangen et al. Reference Tangen, Hsu, Zhu and Linninger2015). Although departures of $\gamma$ from its mean value $\gamma _o$ can be expected as a result of the unequal distribution of fat tissue and epidural veins along the canal, possibly resulting in a smaller value of $\gamma$ in the cervical region (Yallapragada & Alperin Reference Yallapragada and Alperin2004), a uniform value $\gamma =\gamma _o$ is to be employed in the following sample computation, for which $\hat {\gamma }=\gamma /\gamma _o=1$. To estimate the average compliance $\gamma _o$, which determines from (2.8) the dimensionless wavenumber $k$, one may use (2.6) to relate $\gamma _o$ with the elastic wave speed $c$. Using the value $c=4.6$ m s$^{-1}$ reported by Kalata et al. (Reference Kalata, Martin, Oshinski, Jerosch-Herold, Royston and Loth2009) for healthy humans, it follows from (2.6) that $\gamma _0= 8.6\times 10^{-9}$ m$^2$ Pa$^{-1}$ and from (2.8) that $k\simeq 0.81$ for subject 1, with corresponding values $\gamma _0=6.5\times 10^{-9}$ m$^2$ Pa$^{-1}$ and $k\simeq 0.86$ for subject 2.

The modes $\hat {Q}_n$ determined from (3.12) with use made of the anatomical values reported above were used in (3.7) to compute the corresponding dimensionless flow rate $\hat {Q}$ yielding the results shown in figures 2(c) and 3(c). The computation of $\hat {Q}$ involves 10 Fourier coefficients $B_n$, corresponding to the function $\varPi (\tau )$ shown in figure 1(b), taken as representative of a healthy state ICP waveform (Di Ieva, Schmitz & Cusimano Reference Di Ieva, Schmitz and Cusimano2013). The results are compared in figures 2 and 3 with flow rate measurements acquired at 12 vertebral levels using PC-MRI with retrospective cardiac gating in in vivo experiments involving the two subjects (see Coenen et al. (Reference Coenen, Gutiérrez-Montes, Sincomb, Criado-Hidalgo, Wei, King, Haughton, Martínez-Bazán, Sánchez and Lasheras2019) for details of the data-acquisition process). To establish a quantitatively consistent comparison, the results are presented in normalized form. Thus, the stroke volume $V_{S}(\xi )=\tfrac {1}{2} \int _t^{t+T} |{Q}|\, \textrm {d} t$ measured via MRI is scaled with its entrance value (i.e. $V_{S}(0)=0.7$ cm$^3$ for subject 1 and $V_{S}(0)=0.6$ cm$^3$ for subject 2), while the model prediction $\hat {V}_{S}(\xi )=\tfrac {1}{2} \int _\tau ^{\tau +2 {\rm \pi}} |\hat {Q}| \,\textrm {d} \tau$ is correspondingly scaled with $\hat {V}_{S}(0)$, with $\hat {V}_{S}(0)=1.9$ for subject 1 and $\hat {V}_{S}(0)=1.8$ for subject 2. Similarly, the PC-MRI flow rate measurements are scaled with its characteristic value $\omega V_{S}(0)$ while the predicted flow rate $\hat {Q}$ is scaled with $\hat {V}_{S}(0)$.

Figure 2. Subject 1: selected spinal cord locations (a) with corresponding flow rate variation obtained from PC-MRI measurements (b) and from model predictions for $\mathcal {R}=4.78$ (c), $\mathcal {R}=0$ (d), $\mathcal {R}=1$ (e) and $\mathcal {R}=20$f), with associated normalized stroke volumes shown in (g).

Figure 3. Subject 2: selected spinal cord locations (a) with corresponding flow rate variation obtained from PC-MRI measurements (b) and from model predictions for $\mathcal {R}=4.78$ (c), $\mathcal {R}=0$ (d), $\mathcal {R}=1$ (e) and $\mathcal {R}=20$f), with associated normalized stroke volumes shown in (g).

With the uniform values $\hat {\gamma }=1$ and $\mathcal {R} = 4.78$ selected in the computations, the comparisons in figures 2 and 3 indicate that the model is able to describe reasonably well the main features of the pulsating flow rate, including its particular shape and the phase lag of its peak value. In addition, the streamwise decay of the stroke volume (figures 2g and 3g), a metric often used in the medical community to characterize CSF flow, is seen to agree well with the corresponding values computed from the MRI data, with root-mean-square differences between model and MRI remaining below 0.09 for both subjects. As a further consistency check, one may use the measured value of $V_{S}(0)$ along with the accompanying dimensionless prediction $\hat {V}_{S}(0)$ to obtain from the definition $V_{S}(0)=\gamma _o L \Delta p \hat {V}_{S}(0)$ an estimate for the mean fluctuation amplitude $\Delta p$. The value obtained, $\Delta p=73$ Pa for subject 1 and $\Delta p=80$ Pa for subject 2, corresponding to a peak-to-peak value of approximately 200 Pa, is well within the range of values reported in the literature (Eide & Brean Reference Eide and Brean2006), thereby providing additional confidence in the model.

To investigate the sensitivity of the model predictions to changes in trabeculae resistance, results for varying $\mathcal {R}$ are also included in figures 2 and 3. In the absence of trabeculae, i.e. for $\mathcal {R}=0$, both the rate at which the flow rate fluctuation decays along the canal and the general flow rate waveform, including the associated phase lag, are in poor agreement with the MRI measurements. Consequently, for $\mathcal {R}=0$ the modelled streamwise decay of stroke volume, shown in figures 2($g$) and 3($g$), is also seen to strongly overpredict the amount of CSF flow in the thoracic and lumbar regions compared with the MRI measurements. These quantitative findings emphasize the important effects of trabeculae, previously pointed out by Tangen et al. (Reference Tangen, Hsu, Zhu and Linninger2015). As can be seen in figures 2 and 3, the results are not very sensitive to the specific choice of $\mathcal {R}$, provided that an order-unity value is selected. Nevertheless, the computations obtained with $\mathcal {R} = 4.78$, the value determined using the permeability and other anatomical features taken from the literature (Tada & Nagashima Reference Tada and Nagashima1994; Stockman Reference Stockman2006; Gupta et al. Reference Gupta, Soellinger, Boesiger, Poulikakos and Kurtcuoglu2009), appear to give better overall agreement with the MRI measurements, providing the optimal amount of attenuation all along the canal.

5. Concluding remarks

The temporal and spatial variation of the flow rate predicted by the one-dimensional model developed here supplemented with a presumed ICP waveform has been shown to predict the main features revealed in in vivo experiments. Utilization of the model in future efforts to develop non-invasive measurements of ICP requires the solution of the inverse problem, i.e. the computation of the ICP fluctuation from the PC-MRI flow measurements. Given the uncertainty regarding the canal compliance and the trabeculae-network properties, it is unclear whether conventional parameter-fitting approaches will be successful in these future developments or whether more elaborate optimization algorithms, possibly based on machine-learning techniques, will be needed. In improving the model, these future efforts should also account for localized pressure losses associated with the presence of additional spinal microstructures, such as nerve roots and ligaments.


The work of S.C., V.H. and A.L.S. was supported by the National Institute of Neurological Disorders and Stroke through contract no. 1R01NS120343-01. The work of W.C. was supported by the Comunidad de Madrid through the contract CSFFLOW-CM-UC3M. C.M.B., C.G. and W.C. acknowledge the support of the Spanish MICINN through the coordinated project PID2020-115961RB. C.M.B. and C.G. also acknowledge the support provided by Junta de Andalucía and European Funds through grant P18-FR-4619.

Declaration of interests

The authors report no conflict of interest.



Berkouk, K., Carpenter, P.W. & Lucey, A.D. 2003 Pressure wave propagation in fluid-filled co-axial elastic tubes part 1: basic theory. Trans. ASME J. Biomech. Engng 125 (6), 852856.CrossRefGoogle ScholarPubMed
Bhadelia, R.A., Bogdan, A.R., Kaplan, R.F. & Wolpert, S.M. 1997 Cerebrospinal fluid pulsation amplitude and its quantitative relationship to cerebral blood flow pulsations: a phase–contrast MR flow imaging study. Neuroradiology 39 (4), 258264.CrossRefGoogle ScholarPubMed
Bloomfield, I.G., Johnston, I.H. & Bilston, L.E. 1998 Effects of proteins, blood cells and glucose on the viscosity of cerebrospinal fluid. Pediatr. Neurosurg. 28 (5), 246251.CrossRefGoogle ScholarPubMed
Carpenter, P.W., Berkouk, K. & Lucey, A.D. 2003 Pressure wave propagation in fluid-filled co-axial elastic tubes. Part 2: mechanisms for the pathogenesis of syringomyelia. Trans. ASME J. Biomech. Engng 125 (6), 857863.CrossRefGoogle ScholarPubMed
Carr, J.B., Thomas, J.H., Liu, J. & Shang, J.K. 2021 Peristaltic pumping in thin non-axisymmetric annular tubes. J. Fluid Mech. 917, A10.CrossRefGoogle ScholarPubMed
Cirovic, S. & Kim, M. 2012 A one-dimensional model of the spinal cerebrospinal-fluid compartment. Trans. ASME J. Biomech. Engng 134 (2), 021005.CrossRefGoogle ScholarPubMed
Coenen, W., Gutiérrez-Montes, C., Sincomb, S., Criado-Hidalgo, E., Wei, K., King, K., Haughton, V., Martínez-Bazán, C., Sánchez, A.L. & Lasheras, J.C. 2019 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. 40, 12421249.CrossRefGoogle ScholarPubMed
Coenen, W., Zhang, X. & Sánchez, A.L. 2021 Lubrication analysis of peristaltic motion in non-axisymmetric annular tubes. J. Fluid Mech. 921, R2.CrossRefGoogle Scholar
Di Ieva, A., Schmitz, E.M. & Cusimano, M.D. 2013 Analysis of intracranial pressure: past, present, and future. Neuroscientist 19 (6), 592603.CrossRefGoogle ScholarPubMed
Du Boulay, G.H. 1966 Pulsatile movements in the CSF pathways. Br. J. Radiol. 39 (460), 255262.CrossRefGoogle ScholarPubMed
Edsbagge, M., Starck, G., Zetterberg, H., Ziegelitz, D. & Wikkelso, C. 2011 Spinal cerebrospinal fluid volume in healthy elderly individuals. Clin. Anat. 24 (6), 733740.CrossRefGoogle ScholarPubMed
Eide, P.K. & Brean, A. 2006 Intracranial pulse pressure amplitude levels determined during preoperative assessment of subjects with possible idiopathic normal pressure hydrocephalus. Acta Neurochir. 148 (11), 11511156.CrossRefGoogle ScholarPubMed
Eide, P.K. & Sorteberg, W. 2016 Outcome of surgery for idiopathic normal pressure hydrocephalus: role of preoperative static and pulsatile intracranial pressure. World Neurosurg. 86, 186193.CrossRefGoogle ScholarPubMed
Ellis, T., McNames, J. & Aboy, M. 2007 Pulse morphology visualization and analysis with applications in cardiovascular pressure signals. IEEE Trans. Biomed. 54 (9), 15521559.CrossRefGoogle ScholarPubMed
Evensen, K.B. & Eide, P.K. 2020 Measuring intracranial pressure by invasive, less invasive or non-invasive means: limitations and avenues for improvement. Fluids Barriers CNS 17 (1), 34.CrossRefGoogle ScholarPubMed
Feinberg, D.A. & Mark, A.S. 1987 Human brain motion and cerebrospinal fluid circulation demonstrated with MR velocity imaging. Radiology 163 (3), 793799.CrossRefGoogle ScholarPubMed
Gholampour, S. & Fatouraee, N. 2021 Boundary conditions investigation to improve computer simulation of cerebrospinal fluid dynamics in hydrocephalus patients. Commun. Biol. 4 (1), 394.CrossRefGoogle ScholarPubMed
Grotberg, J.B. & Jensen, O.E. 2004 Biofluid mechanics in flexible tubes. Annu. Rev. Fluid Mech. 36, 121147.CrossRefGoogle Scholar
Gupta, S., Soellinger, M., Boesiger, P., Poulikakos, D. & Kurtcuoglu, V. 2009 Three-dimensional computational modeling of subject-specific cerebrospinal fluid flow in the subarachnoid space. Trans. ASME J. Biomech. Engng 131 (2), 021010.CrossRefGoogle ScholarPubMed
Gutiérrez-Montes, C., Coenen, W., Lawrence, J.J., Martínez-Bazán, C., Sánchez, A.L. & Lasheras, J.C. 2021 Modelling and direct numerical simulation of flow and solute dispersion in the spinal subarachnoid space. Appl. Math. Model. 94, 516533.CrossRefGoogle Scholar
Kalata, W., Martin, B.A., Oshinski, J.N., Jerosch-Herold, M., Royston, T.J. & Loth, F. 2009 MR measurement of cerebrospinal fluid velocity wave speed in the spinal canal. IEEE Trans. Biomed. 56 (6), 17651768.CrossRefGoogle ScholarPubMed
Khani, M., Sass, L.R., Xing, T., Keith Sharp, M., Balédent, O. & Martin, B.A. 2018 Anthropomorphic model of intrathecal cerebrospinal fluid dynamics within the spinal subarachnoid space: spinal cord nerve roots increase steady-streaming. Trans. ASME J. Biomech. Engng 140 (8), 0810121–08101215.CrossRefGoogle ScholarPubMed
Kurtcuoglu, V., Jain, K. & Martin, B.A. 2019 Modelling of cerebrospinal fluid flow by computational fluid dynamics. In Biomechanics of the Brain, pp. 215–241. Springer.CrossRefGoogle Scholar
Lawrence, J.J., Coenen, W., Sánchez, A.L., Pawlak, G., Martínez-Bazán, C., Haughton, V. & Lasheras, J.C. 2019 On the dispersion of a drug delivered intrathecally in the spinal canal. J. Fluid Mech. 861, 679720.CrossRefGoogle Scholar
Linninger, A.A., Tangen, K., Hsu, C.-Y. & Frim, D. 2016 Cerebrospinal fluid mechanics and its coupling to cerebrovascular dynamics. nnu. Rev. Fluid Mech. 48, 219257.CrossRefGoogle Scholar
Mortazavi, M.M. et al. 2018 Subarachnoid trabeculae: a comprehensive review of their embryology, histology, morphology, and surgical significance. World Neurosurg. 111, 279290.CrossRefGoogle ScholarPubMed
Ommaya, A.K. 1968 Mechanical properties of tissues of the nervous system. J. Biomech. 1 (2), 127138.CrossRefGoogle Scholar
Sánchez, A.L., Martínez-Bazan, C., Gutiérrez-Montes, C., Criado-Hidalgo, E., Pawlak, G., Bradley, W., Haughton, V. & Lasheras, J.C. 2018 On the bulk motion of the cerebrospinal fluid in the spinal canal. J. Fluid Mech. 841, 203227.CrossRefGoogle Scholar
Sincomb, S.J., Coenen, W., Criado-Hidalgo, E., Wei, K., King, K., Borzage, M., Haughton, V., Sánchez, A.L. & Lasheras, J.C. 2021 Transmantle pressure computed from mr imaging measurements of aqueduct flow and dimensions. AJNR Am. J. Neuroradiol. 42 (10), 18151821.Google ScholarPubMed
Sincomb, S., Coenen, W., Sánchez, A.L. & Lasheras, J.C. 2020 A model for the oscillatory flow in the cerebral aqueduct. J. Fluid Mech. 899, R1.CrossRefGoogle Scholar
Singh, V. & Cheng, R. 2021 Neurovascular physiology and neurocritical care. In Handbook of Clinical Neurology, vol. 176, pp. 71–80. Elsevier.CrossRefGoogle Scholar
Stockman, H.W. 2006 Effect of anatomical fine structure on the flow of cerebrospinal fluid in the spinal subarachnoid space. Trans. ASME J. Biomech. Engng 128 (1), 106114.CrossRefGoogle ScholarPubMed
Tada, Y. & Nagashima, T. 1994 Modeling and simulation of brain lesions by the finite-element method. IEEE Engng Med. Biol. Mag. 13 (4), 497503.CrossRefGoogle Scholar
Tangen, K.M., Hsu, Y., Zhu, D.C. & Linninger, A.A. 2015 CNS wide simulation of flow resistance and drug transport due to spinal microanatomy. J. Biomech. 48 (10), 21442154.CrossRefGoogle ScholarPubMed
Thomas, J.H. 2019 Fluid dynamics of cerebrospinal fluid flow in perivascular spaces. J. R. Soc. Interface 16 (159), 20190572.CrossRefGoogle ScholarPubMed
Unnerbäck, M., Ottesen, J.T. & Reinstrup, P. 2018 ICP curve morphology and intracranial flow-volume changes: a simultaneous ICP and cine phase contrast MRI study in humans. Acta Neurochir. 160 (2), 219224.CrossRefGoogle ScholarPubMed
Wagshul, M.E., Chen, J.J., Egnor, M.R., McCormack, E.J. & Roche, P.E. 2006 Amplitude and phase of cerebrospinal fluid pulsations: experimental studies and review of the literature. J. Neurosurg. 104 (5), 810819.CrossRefGoogle ScholarPubMed
Wagshul, M.E., Eide, P.K. & Madsen, J.R. 2011 The pulsating brain: a review of experimental and clinical studies of intracranial pulsatility. Fluids Barriers CNS 8 (1), 5.CrossRefGoogle ScholarPubMed
Yallapragada, N. & Alperin, N. 2004 Characterization of spinal canal hydrodynamics and compliance using bond graph technique and CSF flow measurements by MRI. Proc. Intl Soc. Magn. Reson. Med. 11, 2658.Google Scholar
Figure 0

Figure 1. (a) Main anatomical features of the spinal canal for subject 1; (b) ICP wave form (Wagshul, Eide & Madsen 2011) CC BY 2.0 (left subpanel) and dimensionless (right subpanel); (c) dimensionless canal functions for subjects 1 and 2.

Figure 1

Figure 2. Subject 1: selected spinal cord locations (a) with corresponding flow rate variation obtained from PC-MRI measurements (b) and from model predictions for $\mathcal {R}=4.78$ (c), $\mathcal {R}=0$ (d), $\mathcal {R}=1$ (e) and $\mathcal {R}=20$f), with associated normalized stroke volumes shown in (g).

Figure 2

Figure 3. Subject 2: selected spinal cord locations (a) with corresponding flow rate variation obtained from PC-MRI measurements (b) and from model predictions for $\mathcal {R}=4.78$ (c), $\mathcal {R}=0$ (d), $\mathcal {R}=1$ (e) and $\mathcal {R}=20$f), with associated normalized stroke volumes shown in (g).