## Abstract

Abnormal cerebrospinal fluid (CSF) flow is suspected to be a contributor to the pathogenesis of neurodegenerative diseases such as Alzheimer's through the accumulation of toxic metabolites, and to the malfunction of intracranial pressure regulation, possibly through disruption of neuroendocrine communication. For the understanding of transport processes involved in either, knowledge of *in vivo* CSF dynamics is important. We present a three-dimensional, transient, subject-specific computational analysis of CSF flow in the human cranial subarachnoid space (SAS) based on *in vivo* magnetic resonance imaging. We observed large variations in the spatial distribution of flow velocities with a temporal peak of 5 cm s^{−1} in the anterior SAS and less than 4 mm s^{−1} in the superior part. This could reflect dissimilar flushing requirements of brain areas that may show differences in susceptibility to pathological CSF flow. Our methods can be used to compare the transport of metabolites and neuroendocrine substances in healthy and diseased brains.

## 1. Introduction

The cerebrospinal fluid (CSF) not only provides mechanical protection to the central nervous system, but also acts as a transport medium for nutrients, neuroendocrine substances and for the removal of toxic metabolites, preserving the chemical environment of the brain (Davson *et al*. 1987). As outlined in the second part of this work (Holman *et al*. 2010), it is hypothesized that the disruption of CSF flow may be linked to neurodegenerative diseases such as Alzheimer's through disturbed regulation of intracranial pressure, through accumulation of toxic metabolites or through a combination of both (Segal 2000; Stopa *et al*. 2001; Kivisakk *et al*. 2003; Silverberg *et al*. 2003; Abbott 2005; Johanson *et al*. 2005).

The accumulation of any substance transported through the CSF is governed to a large extent by the pulsatile dynamics of this fluid. Surprisingly, little detail is known about CSF flow. This can be attributed to the fact that potentially accurate invasive flow measurements necessarily alter its motion, and that non-invasive studies are limited in resolution and accuracy, except in areas of well-defined shape with CSF motion along a predominant axis such as in the aqueduct of Sylvius, in the superior part of the spinal subarachnoid space (SAS) and, to a certain extent, in some planes within the various cranial cisterns. Magnetic resonance imaging (MRI) velocimetry is currently the most widely used technique to study the CSF flow non-invasively (Feinberg & Mark 1987; Enzmann & Pelc 1992; Poncelet *et al*. 1992).

Using three-dimensional computational fluid dynamics (CFD), it is possible to reconstruct the velocity field in the CSF space. This has been done for parts of the ventricular system, the spinal SAS and the inferior cranial SAS, but not for the remainder of the cranial SAS that poses serious modelling challenges owing to its complex anatomy spanning multiple length scales (Jacobson *et al*. 1996; Aroussi *et al*. 2000; Loth *et al*. 2001; Stockman 2006; Kurtcuoglu *et al*. 2007*a*; Gupta *et al*. 2009). This includes the spongy framework of trabeculae, which are delicate tissue filaments that protrude through the meninges into the SAS (Killer *et al*. 2003*a*). CSF is produced in the choroid plexus and is absorbed through the arachnoid granulations (AGs) (Haroun *et al*. 2007; Weller *et al*. 2009), extra-cranial lymphatic pathways (Johnston 2003) and possibly through additional intraparenchymal routes (Greitz *et al*. 1997). The drainage through the lymphatic pathways is well validated in mammals other than man, but there is no consensus on the relative importance of lymphatic CSF drainage compared with absorption through the AGs in humans, which is the traditionally accepted route (Koh *et al*. 2005; Weller *et al*. 2009). Despite the paucity of evidence supporting the AGs as the sole CSF absorption route, we only take into account this pathway by way of convenience, since more quantitative data on absorption through the AGs are available than on the lymphatic pathway. We show in §4 that this simplification will not influence the reported results substantially. With an average diameter of approximately 300 µm (Upton & Weller 1983), the AGs are well below the resolution limit of current clinical MRI units, and thus have to be modelled based on *ex vivo* (Grzybowski *et al*. 2007) and *in vitro* (Grzybowski *et al*. 2006; Holman *et al*. 2010) data.

In the work at hand, we use CFD to compute and analyse the flow in the cranial SAS based on *in vivo* anatomic and velocimetric MRI data of a healthy male volunteer (figure 1). We take into account the distribution of AGs based on *ex vivo* data gathered from 35 human brains (Grzybowski *et al*. 2007) and estimate the CSF outflow rate through the AGs again using *in vivo* MRI flow measurements. We model the influence of the trabecular framework using a previously described anisotropic porous model (Gupta *et al*. 2009). In that work, we investigated CSF flow in the superior spinal and inferior cranial SAS and the fourth ventricle. These areas are smaller than the here covered superior cranial SAS, they do not contain any AGs and they have less intricate anatomic features.

We show in the work at hand that there are large spatial variations in the velocity distribution in the cranial SAS that will influence the transport behaviour of toxic metabolites and neuroendocrine and other substances released into the CSF. Assuming there are physiological reasons for these variations, we hypothesize that not all regions of the brain will be affected with the same severity by disrupted CSF flow.

## 2. Material and methods

### 2.1. SAS anatomy and CSF flow acquisition

T2-weighted MRI scans were performed on a 25-year-old healthy male volunteer using a whole-body 3T MRI scanner (Achieva, Philips Medical Systems, Best, The Netherlands). The in-plane resolution of the scans was 0.45 × 0.45 mm and the slice spacing was 0.6 mm. The details of the MRI acquisition sequence are presented in Gupta *et al*. (2009). The acquired images were manually segmented using Amira 4.1 (Mercury Computer Systems, San Diego, CA, USA) to obtain a three-dimensional representation of the SAS (figures 1 and 2). The voxel-based segmented structures were then converted to non-uniform rational B-spline surfaces in order to allow for efficient generation of a high-quality computational grid.

MRI velocimetry was used to acquire CSF velocity data at the inferior boundaries of the cranial SAS, i.e. in the pontine and cerebellomedullary cisterns. A phase contrast velocity mapping sequence combined with an ECG-triggered T1-weighted transient field echo read-out scheme was applied (Gupta *et al*. 2009). In order to maximize the velocity signal-to-noise ratio, the phase-encoding velocity was adjusted so as to conform to the maximum velocity at the boundary. The phase-encoding velocity was set to 2 cm s^{−1} for the measurements in the cerebellomedullary cistern and to 4 cm s^{−1} for measurements in the pontine cistern. The measured profiles were manually phase-unwrapped, filtered with a filter mask of 5 × 5 pixels (Lim 1990) and then smoothed using cubic spline interpolation. Figure 3 shows the boundary velocity profiles at the measurement locations within the pontine and cerebellomedullary cisterns during one cardiac cycle.

### 2.2. Distribution and hydraulic conductivity of arachnoid granulations

The detailed characterization of the AGs' topography and morphology is currently not possible with any *in vivo* measurement technique. Grzybowski *et al.* (2007) characterized the AG distribution in the superior sagittal sinus based on *en face* digital images of 35 human brain samples. For the present analysis, areas of AGs that were present in more than 10 per cent of the analysed brains in Grzybowski *et al*. (2007) were mapped to the current anatomic geometry. The surface area occupied by the AGs, *A*_{AG}, amounted to 114 mm^{2}, which is approximately 0.29 per cent of the total brain surface area. It is to be noted that while CSF drainage may take place through several pathways, in our model all CSF absorption occurs exclusively through the AGs clearing into the superior sagittal sinus.

The most important parameter that governs CSF outflow through the granulations is their hydraulic conductivity (Albeck *et al*. 1991; Eklund *et al*. 2007; Grzybowski *et al*. 2007). Hydraulic conductivity, *L*_{P,AG}, is a measure of the amount of CSF that flows through a unit cross-sectional area of the granulation surface as a result of a unit pressure drop across that surface,
2.1
where Δ*P*_{AG} is the total pressure drop across the AG surface and is the net outflow through the AGs. Alternatively, the resistance to cerebrospinal outflow can be defined as
2.2

We modelled the AGs as differential pressure valves that respond to the pressure difference Δ*P*_{AG} between the SAS and the superior sagittal sinus, allowing only unidirectional flow out of the SAS. Estimates of Δ*P*_{AG} across the granulation membranes are published in Davson *et al*. (1987) and Grzybowski *et al*. (2006) as Δ*P*_{AG} = 3.15 mmHg (420 Pa).

### 2.3. Computational model

We have reconstructed CSF motion in the cranial SAS using a finite-volume (FV) CFD approach (Versteeg & Malalasekra 1996). The fundamental basis of this approach is the discretization of the governing fluid flow partial differential equations known as Navier–Stokes (NS) equations into a system of coupled algebraic equations that can be solved iteratively. While FV CFD is a computationally expensive technique, particularly when dealing with complex geometries, it is well validated and its ability to handle complex flows has been demonstrated in various research and industrial applications. In order to incorporate the effect of the trabecular morphology of the SAS into the calculations, the SAS was modelled as a porous medium by extension of the NS equations to the NS/Brinkman equations (Brinkman 1948). The permeability *k* of the SAS varies anisotropically throughout the domain. For its estimation, the porous SAS is assumed to feature an idealized morphology with straight cylindrical pillars representing the trabeculae extending perpendicularly from the pia to the arachnoid mater. This simplification allows for the application of a closed-form mathematical solution for the permeability in an arbitrary direction within the SAS for a given trabecular radius and porosity (Van der Westhuizen & Du Plessis 1994). For the current computations, we used an SAS porosity value of *ε* = 0.99 (Tada & Nagashima 1994) and a trabecular radius of *r* = 15 µm (Killer *et al*. 2003*b*). A detailed description of the derivation of the porous SAS model is given in Gupta *et al*. (2009).

CFD analysis requires the discretization of the SAS domain into a large number of small cells, within which pressure and velocity are approximately uniform. These cells define the computational grid. The presence of sulci and the wide range of length scales within the cranial SAS necessitate a very fine computational grid for accurate CFD computations. We modelled CSF as an incompressible, Newtonian fluid with the same density and viscosity as that of water at 37°C (Bloomfield *et al.* 1998; Kurtcuoglu *et al*. 2007*b*). The computations were performed using a commercial FV CFD solver Fluent 6.3 (Fluent Inc., Lebanon, NH, USA) on a high-quality computational grid consisting of 14 million cells. A discussion of the influence of the choice of computational grid on the reported results is presented in the electronic supplementary material, figure S7.

Velocity profiles obtained using MRI as described above were imposed at boundaries within the pontine and cerebellomedullary cisterns (figure 3). The measured profiles were first interpolated onto the computational grid and then interpolated in time in accordance with the temporal step size used for the transient computations. In order to take into account brain motion, an artificial transient volumetric flow rate, , was imposed uniformly at the entire brain surface, with flow direction normal to the surface. This approach lets us avoid the exceedingly complex and computationally costly transient deformation of the computational grid, while still accounting for spatially averaged CSF motion induced by brain deformation. This artificial volumetric flow rate was determined as 2.3 where and are the transient volumetric flow rates into the SAS domain at the boundaries in the pontine and cerebellomedullary cisterns, respectively, and is the net outflow through the AGs given by equation (2.4). The volumetric flow rates and were obtained by integrating the spatial velocity profiles (acquired using velocimetric MRI as described above) over the respective boundary. The calculation of is shown in §3. The obtained spatially averaged transient brain surface deformation rate, calculated by dividing by the brain surface area, is shown in figure 4. The boundary conditions used in the current computational model are further summarized in table 1.

If we neglect absorption of CSF at locations other than the AGs, the net CSF produced within a cardiac cycle will equal the total outflow through the AGs, , which can be expressed as
2.4
where *T* is the time period of one cardiac cycle. Using equation (2.4), evaluates to 0.47 ml min^{−1}. With a pressure drop, Δ*P*_{AG}, of 3.15 mmHg across the granulation membrane and an AG surface area of 114 mm^{2}, the hydraulic conductivity of the granulations *L*_{P,AG} is 130.9 µl min^{−1} mmHg^{−1} cm^{−2} (equation (2.1)). Correspondingly, the resistance, *R*, of the AGs to CSF outflow evaluates to 6.7 mmHg ml^{−1} min^{−1} (equation (2.2)). This high value is in good agreement with the results reported in Albeck *et al*. (1991) and Eklund *et al*. (2007) and is indicative of the fact that the mesothelial cell layer lining the AG surface offers high resistance to CSF outflow, thus limiting CSF drainage into the venous blood.

## 3. Results

We will refer to the space surrounding the parietal and superior frontal lobes of the cerebral cortex as the superior cranial SAS (figure 2). The space surrounding the frontal lobe and the temporal lobe including chiasmatic and interpeduncular cisterns will be referred to as the anterior cranial SAS. Inferior to the interpeduncular cistern is the pontine cistern, where CSF enters the cranial SAS through the left and right foramina of Luschka and through the spinal SAS. The posterior cranial SAS surrounds the occipital lobe, including the quadrigeminal cistern and the region below it. It also surrounds the cerebellum and extends up to the cerebellomedullary cistern.

Figures 5 and 6 and the movies in the electronic supplementary material visualize the CSF velocity and pressure field. With respect to velocity distribution, the cranial SAS can be divided into three distinct regions: (i) posterior cranial and superior cranial SASs, exhibiting low velocities with peak values less than 4 mm s^{−1} throughout the cardiac cycle, except in the vicinity of the cerebellomedullary cistern, where the peak velocity is 8.3 mm s^{−1}; (ii) anterior cranial SAS, excluding chiasmatic and interpeduncular cisterns, exhibiting higher velocities than the posterior cranial and superior cranial spaces with absolute maximum velocity up to 1.3 cm s^{−1}; (iii) chiasmatic, interpeduncular and pontine cisterns with maximum velocity up to 5 cm s^{−1}. Table 2 presents various flow parameters at selected cross sections within the SAS. The stroke volume, SV, is obtained in two steps, first by determining the area under the transient volumetric flow-rate curve of the respective cross section for caudio-cranial (positive) and cranio-caudal (negative) flow rates, and then by taking the mean of the absolute value of these two areas. For the AGs, the stroke volume is calculated as the product of net volumetric flow rate through the granulations, , and the time period, *T*, of the cardiac cycle.

The peak Reynolds number, *Re*_{peak}, is calculated based on the peak velocity, *ϑ*_{peak}, and the hydraulic diameter, *D*_{h} = 4 ⋅ area/perimeter, of the respective cross section and is defined as
where *ρ* and *μ* are the CSF density and dynamic viscosity, respectively. The Reynolds number, *Re*, in general, is defined as the ratio of the inertial to the viscous forces present at a particular location within a given system. It is used as a measure of flow characteristics, where higher values indicate turbulent behaviour, lower values imply laminar flow and very low values indicate diffusive behaviour. The limits between these regimes are system dependent, but can be approximated as *Re* < 1 for the diffusive range, 1 < *Re* < 1000 for the laminar range and potentially turbulent dynamics for larger values. The maximum Reynolds number in the cranial SAS, *Re*_{peak} = 386, was observed at the pontine boundary and lies well within the laminar range. Low Reynolds numbers with a peak value of 72 were observed in the posterior cranial SAS at the cerebellomedullary cistern boundary, and with a peak value of around 20 at plane B (see insert in table 2) in the posterior cranial SAS. The flow is primarily laminar in this region, with pockets of diffusive zones with *Re* < 1 such as in the quadrigeminal cistern.

The ratio of transient inertial forces to viscous forces can be studied by considering the Womersley parameter, *α*, which is defined as
where *T* is the time period of the cardiac cycle. For small Womersley parameters (*α* < 1), the transient inertial forces are low, and hence the flow remains in phase with the driving pressure gradient. For larger Womersley parameters (*α* > 1), the transient inertial forces are high enough and the velocity profiles do not immediately follow the driving pressure gradient. As can be seen in table 2, the Womersley parameter is quite high in the entire treated domain, indicating strong transient inertial effects on the CSF flow field. Figure 5 shows the pressure distribution, *P*, at selected points within the cardiac cycle relative to the CSF outlet pressure *P*_{CSF,outlet} = Δ*P*_{AG} + *P*_{SSS}, where *P*_{SSS} is the blood pressure in the superior sagittal sinus and Δ*P*_{AG} is the pressure drop across the AGs' surface. The absolute CSF pressure, *P*_{abs}, at any point within the domain can be determined as *P*_{abs} = *P*+Δ*P*_{AG} + *P*_{SSS}. The flow dynamics within a cardiac cycle can be described as follows.

At the beginning of the cardiac cycle, which is chosen to coincide with the R-peak of the electrocardiogram, CSF flows in the cranio-caudal direction in the anterior cranial SAS and in the pontine and cerebellomedullary cisterns. CSF flows in the cranio-caudal direction in the posterior cranial space surrounding the cerebellum and in the caudio-cranial direction in the SAS surrounding the occipital lobe. Thus, at the beginning of the cardiac cycle, CSF flows from the posterior cranial to the anterior cranial space with a net flow from the cranial SAS to the ventricular space and to the spinal SAS. This is caused, presumably, by the expansion of the lateral and third ventricles and the expansion of blood vessels in the SAS and beneath the pia mater, or by the compression of the SAS. The pressure, *P*, varies from −0.15 mmHg in the anterior cranial space to 0.04 mmHg in the posterior cranial space.

At around 20 per cent of the cardiac cycle, the flow in the entire space turns caudio-cranial. This is probably caused by the compression of blood vessels in or beneath the SAS or by the expansion of the SAS itself. In our computational model, both these effects are merged into a combined distensibility boundary condition at the brain surface and can, thus, not be distinguished (§2.3). The flow in this phase of the cardiac cycle features the highest velocities of the entire cycle. The pressure varies from 0 to 0.20 mmHg, with maximum value within the anterior cranial space.

At around 40 per cent of the cycle, CSF flows from the anterior cranial to the posterior cranial space with a net flow from the ventricular system to the cranial SAS. The reason for this is ventricular contraction, forcing CSF into the pontine cistern via the two foramina of Luschka. CSF in the pontine cistern is further carried upwards towards the superior cranial space, with the inertia of the CSF arriving from the spinal SAS. This phase of the cycle features very low velocities (less than 2 mm s^{−1}) in the superior cranial and the posterior cranial SAS and a high inertia flow in the pontine cistern and in the anterior cranial SAS. The overall pressure varies from 0 to 0.11 mmHg, with maximum pressure still in the anterior cranial space.

The entire flow almost stagnates at around 60 per cent of the cardiac cycle. Peak velocities in the entire treated domain are less than 5 mm s^{−1} at this time and the pressure variation is small.

At approximately 80 per cent of the cardiac cycle, the flow closely resembles that at the beginning of the cycle. The absolute maximum velocity during this phase is 3 cm s^{−1}, occurring within the pontine cistern. The velocities in the rest of the domain are very low. There is a decrease in the pressure in the anterior cranial space, while the pressures within the superior cranial and the posterior cranial SAS remain similar to those at 60 per cent of the cycle.

## 4. Discussion and limitations

The pressure in the superior cranial SAS remains positive throughout the entire cardiac cycle, facilitating the outflow of CSF through the AGs. Our computations indicate the hydraulic conductivity of the granulations, *L*_{P}_{,}_{AG}, to be approximately 130.9 µl min^{−1} mmHg^{−1} cm^{−2} for a physiological pressure drop of 3.15 mmHg across the granulation membrane. This compares relatively well to a value of 92.49 ± 11.79 µl min^{−1} mmHg^{−1} cm^{−2} obtained using an *in vitro* model of cells cultured from human AGs that we present in the second part of this work (Part II. *In vitro* arachnoid outflow model) (Holman *et al*. 2010).The difference between the computationally derived and the *in vitro* hydraulic conductivity can be explained by several factors, including the fact that the computations do not consider CSF outlets other than the AGs, necessitating a higher hydraulic conductivity to obtain the same absorption rate. A detailed discussion is given in Part II (Holman *et al*. 2010).

Limiting CSF outflow to the AGs clearing into the superior sagittal sinus may also induce errors in the calculation of the CSF flow field. There is, in particular, strong evidence for extracranial lymphatic drainage (Johnston *et al*. 2000; Johnston 2003), but also indications for drainage of CSF through granulations present in transverse and sigmoid sinuses (Haroun *et al*. 2007) as well as parenchymal (Greitz *et al*. 1997; Levine 1999) and spinal (Edsbagge *et al*. 2004) absorption of CSF. However, the impact of CSF absorption on the flow field is likely to be insignificant, because the outflow stroke volume of 6.5 µl (table 2) is much smaller than the stroke volumes in the rest of the domain. The accurate choice of CSF drainage locations and absorption rate distribution between them will become important when substance transport within the SAS is to be calculated.

Large variations in the CSF velocity field indicate substantial spatial differences in transport characteristics. For flow in the third cerebral ventricle, there are indications that the ventricle shape is optimized with regards to substance transport between the pituitary gland and the hypothalamus (Kurtcuoglu *et al*. 2007*b*). It is conceivable that, similarly, the flow variations in the SAS reflect an optimization with regards to efficient transport of toxic metabolites and other substances used for communication through the CSF space. This would suggest that different areas of the brain would be affected dissimilarly by a disruption of CSF flow. As an example, a brain region requiring very efficient flushing of toxic metabolites may experience accumulation of waste products that could lead to or accelerate neurodegenerative processes, whereas, in other areas, the disruption may have little to no effect.

For the current computations, we used an SAS porosity value of *ε* = 0.99 (Tada & Nagashima 1994) and a trabecular radius of *r* = 15 µm (Killer *et al*. 2003*b*). The physiological variations in trabecular morphology and topography are likely to influence the CSF flow field within the SAS. In order to determine the extent of this influence, we performed a sensitivity analysis with a twofold increase in trabecular density and radius, respectively. The results are summarized in table 3. Doubling the trabecular density decreases the longitudinal permeability by 75 per cent, transverse permeability by 68 per cent and the porosity by 1 per cent. It further changes the total pressure variation to 0.42–0.69 mmHg within one cardiac cycle, from the value of −0.15 to 0.20 mmHg in the simulation case with the original trabecular density and radius. Similarly, doubling the trabecular radius decreases longitudinal permeability by 76 per cent, transverse permeability by 62 per cent and porosity by 3 per cent. This results in a total pressure variation of −0.36 to 0.60 mmHg_{.} In conclusion, trabecular radius and density both significantly influence the pressure drop across the SAS domain.

The total nominal volume of the reconstructed cranial SAS is approximately 127 ml. This volume increases and decreases cyclically owing to the compression and expansion of the blood vessels in or beneath the SAS and owing to distensibility of the parenchyma itself. This volume change is taken into account in our computational model by means of homogeneously distributed in and outflow of CSF over the entire brain surface. While there is general agreement that parenchyma deformation is not uniform owing to finite brain asymmetries and owing to the brain's anisotropic material properties (Toga & Thompson 2003), accurate *in vivo* measurement of these deformations is still very challenging (Soellinger *et al*. 2007, 2009). Consequently, our model cannot reproduce any local flow variations resulting from the non-uniform deformation of the SAS. Additional local flow variations may be induced by blood vessels in or beneath the SAS, which are not all taken into account in our model.

## 5. Conclusion

In the work at hand, we have described quantitatively and qualitatively the cranial CSF flow dynamics in an adult human using a computational model based on *in vivo* MRI data. We have included details of the SAS down to the resolution limit of the employed MRI unit and have taken into account the sub-resolution trabecular structures using an anisotropic porous media model. We have included the outflow of CSF through the AGs based on an estimation of their flow resistance, and have obtained a prediction of their hydraulic conductivity that agrees reasonably well with corresponding *in vitro* results. While two-dimensional models of CSF flow in the cranial SAS exist, the work at hand is the first to report transient, three-dimensional *in vivo* velocity and pressure distributions. It is also the first to report the effect of changes in trabecular morphology and topography on the cranial CSF pressure gradient.

## Acknowledgements

The financial support of the ETH Zurich Research Commission and the Swiss National Science Foundation through *SmartShunt—The Hydrocephalus Project* are kindly acknowledged. We would also like to thank the Ohio Lions Foundation and the North American Neuro-Ophthalmology Society for support of this work.

## Footnotes

- Received January 25, 2010.
- Accepted February 24, 2010.

- © 2010 The Royal Society