Three-dimensional cerebrospinal fluid flow within the human central nervous system

This work describes a three-dimensional (3D) numerical simulation of the flow field of the complete enclosed Central Nervous System (CNS) including the ventricular system, the spinal cord and spinal sub-arachnoid space (SAS). Previous works on the topic consider only parts of the complete system imposing artificial boundary conditions at internal cross sections. The computational domain was constructed from MRI data using Materialise Software Mimics. In this work pulsatile velocity inlets in the lateral ventricles, due to the cardiac cycle, were used to simulate the dynamic nature of the CSF, whilst pressure outlets were used to model the areas of CSF re-absorption. A porous medium formulation (Darcy flow) was considered in the SAS to account for the effect of the arachnoid trabeculae within these areas. The simulation was run using the commercial CFD code Fluent using the laminar solver and transient simulation. A maximum CSF velocity was found to be in the region 11.8 mm/s, with a peak pressure drop through the aqueduct of the order of 2.8 Pa corresponding to a calculated peak Reynolds number of 12. CSF pressure at the exits of Magendie and Luschke were found to vary over the cardiac cycle, with pressure at the exits of Luschke being higher than Magendie for large periods of the cycle. CSF was seen to enter the SAS as a laminar jet from the exits of Magendie and Luschke. By only considering the cardiac cycle a very slow CSF motion within the spinal SAS was observed with magnitudes significantly reduced after a depth of 50mm down the column from the exit of Magendie. This result suggests that pulsating wall motion in the region of the spinal cord, due to respiratory effects, needs to be considered in order to predict experimentally observed flow recirculation within the spinal sac. 287 words.


Introduction. Cerebrospinal fluid (CSF) is a Newtonian fluid with similar
properties to water and is primarily secreted by the choroid plexus within the cerebral ventricles of the brain. It flows within the Central Nervous System (CNS) coating the spinal cord and skull, protecting the delicate organs of the body from injury, as well as carrying nutrients. The motion of CSF has been debated within the medical literature for a number of years. Understanding its complex motion will allow greater clarity in this area of medicine. This in turn will lead to improvements on techniques that use the CSF flow motion to deliver medical drugs to specific areas of the CNS, as is the case of the current treatment via CSF of Leptomeningeal dissemination, in which tumour cells migrate around the CSF. Besides having good understanding of the CSF flow motion, it could be possible to treat some severe forms of mental illness, such as epilepsy and Parkinsons disease, via drugs delivery through the CSF. Using measurement techniques such as magnetic resonance imaging (MRI) clinicians have demonstrated the pulsatile nature of the CSF and its relationship with the cardiac cycle, however despite its usefulness there are large limitations with MRI when studying fluid motion within the human body, and especially with the CSF dynamics within the CNS. Detailed MRI studies can only be conducted over a small area of interest. At the same time, the type of data required from a scan must be decided in advance. The CNS is a large and interacting system, with CSF moving through its passage in a complex and unusual way that is not fully understood. At some points the flow is almost stagnant, at others it is rapidly accelerating. All of this information cannot be captured in one MRI investigation. There exists therefore, a role for numerical simulations to address the limitations of MRI and resolve the ambiguities surrounding CSF motion. There have been some attempts to simulate CSF motion, however early numerical simulations were limited in either the anatomical geometry of the flow domain, generally presented as a 2D section through the system, or by modelling individual parts of the CNS [1,2,3,4,5,6]. Only limited attempts have been made to examine the CSF flow dynamic within the wider context of the CNS [7,8,9]. Limiting computational analysis to small anatomical parts of the CNS, such as the aqueduct or ventricular system without considering the influence of the remaining part of flow system will not lead to a thorough understanding of CSF dynamics at that anatomical part of the CNS. This point has been highlighted by Kurtcuolgu et al (2005) [4] who argued that meaningful analysis of CSF motion of all parts of the CNS needs to be coupled together in one numerical simulation. This argument was analysed in a previous paper by the authors [10] who showed that static boundary conditions at the exits of Magendie and Luschke should not be applied in CFD models of the ventricles as these boundaries are transient and are influenced by the fluid dynamics of the whole system. This paper seeks to address these concerns by presenting a simplified three dimensional MRI generated CFD numerical simulation of the CNS, including the ventricular system and subarachnoid spaces (SAS). The investigation uses a previously reported method to simulate the expansion of the choroids plexus [10] to examine the pulsatile nature of the CSF. The work also examines the effect of porous arachnoid trabeculae within the SAS on the fluid motion of the CSF within the CNS.
2. Methodology. A three dimensional (3D) model of the CNS was produced computationally using MRI data with 1mm sided voxels from cranial vault to sacrum. The MRI manipulation software Mimics (Materialise NV) was used for the threshold separation technique used to capture regions of CNS in each MRI slice and produce a 3D surface model of the entire volume ( Figure 1). The modelling software Mimics was used to generate the flow exit apertures, the inlets and to connect the ventricles to the SAS. The inlets at the choroid plexus were approximated to be two discs of diameter 6mm embedded 2mm into the each of the lateral ventricles. The cerebellum and hindbrain were simplified to the shape of a cone due to the complexity of interpreting the boundaries on the MRI scan in this region. This was connected to the top of the spinal cord, allowing the CSF exiting the foramen of Magendie to flow into the spinal SAS. The exits of Luschke are known to be tube like structures which project at approximately 90 • to the ventricular system. These tubes (in the model 4mm diameter) extend from the fourth ventricle and finish at the point where the Cranial SAS begins. The areas of CSF re-absorption (flow exits), known as the Pacchionian granulations are composed of arachnoid villi positioned along the superior sagittal sinus, a large venous space running from front to back along the centre of the inside of the skull. These granulations were approximated as two circular exits, of diameter 2mm, at the top-centre of the model, where a pressure outlet boundary condition is specified. Figure 2 shows views of the upper section of the model, not including the spinal cord.
The surface model was converted to a volumetric model containing 640,080 tetrahedral cells using the CFD pre-processing software GAMBIT. Boundaries were applied to the model with the areas of re-absorption set as pressure outlets whilst the choroid plexus were modelled as pulsating velocity inlets.
2.1. CFD numerical simulation. The commercial CFD code FLUENT R , was used to calculate the CSF flow dynamics within the CNS. FLUENT R uses the Finite Volume Method to solve the Navier-Stokes equations; the fluid is incompressible and the SIMPLE pressure-velocity linking algorithm was used with a laminar viscous model. The choroid plexus was treated as the source of pulsing in the CSF, using a technique previously developed by Linninger et al (2005) [7]. This uses the hypothesis of Bering (1955) [11], that the pulsatile motion of the CSF is a result of the expansion and contraction of the choroid plexus, due to the pulsing of arterial blood around the site during the cardiac cycle. Velocity inlet boundary conditions were applied to the choroid plexus, with a User Defined Function (UDF) to mimic the pulsing nature of CSF. The velocity history in the model is based on the CSF pulse wave derived mathematically by Portnoy and Chopp (1981) [12] which was included in previous work by the authors [10]. Figure 3 shows the transient velocity profile at the choroid plexus, for which the average velocity corresponds to a CSF daily rate production of 500 ml per day. In order to simulate the porous medium produced by the arachnoid trabeculae the numerical simulation was set up containing a porous medium in all fluid regions of the SAS. The porosity was determined from the work of Jacobson et al (1999(b)) [13], who examined a simplified CFD model of the SAS and found that a reasonable permeability of the SAS was between 10 −7 and 10 −9 m 2 . Using the lower value of a 10 −9 m 2 , calculations based on Darcys law were performed that produced the equivalent values of porosity and volume fraction required to simulate the porous medium using CFD. The simulation was modelled with rigid walls with a non slip velocity boundary condition. The low Reynolds number within the region of greatest velocity magnitude (the aqueduct) allowed the simulation to be modelled using the laminar solver. Cell refinement, until the changes between velocities at various points in the grid between each refinement were negligible, was used to achieve grid independence. The numerical simulation was run over 5 cardiac cycles to ensure that a cyclic condition was established using a time step of 0.1ms, with readings taken at 10 ms intervals. Numerical results. The numerical simulation allows for the full 3D fluid structure of CSF motion within the CNS to be examined. In order to discuss the findings of the work, the numerical results have been divided into the following four sections: 3.1. Aqueduct of Sylvius. Velocity in the aqueduct of Sylvius is the highest in the ventricular system and fluctuates due to the pulsatility of the CSF. Figure 4 compares maximum CSF velocity in the aqueduct with velocity at the choroid plexus during one cycle. Maximum CSF velocity follows the same profile as the velocity at the choroid plexus, but there is a slight phase shift in the peak maximum velocity due to the local inertial effect. CSF velocity at the Choroids plexus peaks at 0.33s into the cycle, whilst CSF velocity magnitude peaks in the aqueduct at 0.34s into the cycle. Peak CSF velocity was calculated to be 11.8 mm/s, with a corresponding peak Reynolds number in the order of 12. The minimum peak velocity in the aqueduct is 8.05 mm/s. The peak pressure over one cycle during the simulation is shown in figure 5. The results are presented alongside the peak pressure calculations for a similar simulation in which the SAS was modelled without considering a porous structure. Peak CSF pressure within the CNS was 3.77 Pa with a minimum pressure of -0.92Pa when the SAS was modelled as a porous medium.
3.2. Ventricular system. Visualising the flow of CSF within the ventricular system is difficult due to its complex shape. Figure 6 shows contours of velocity magnitude on a cross section through the ventricles. The image was taken at time of maximum CSF velocity, 0.34s into the cycle and the velocities plot were limited to a maximum of 0.5mm/s for image clarity, thus ignoring the detail in the aqueduct.
The velocity magnitudes from figure 6 show the acceleration of the CSF through the foramen of Magendie and again at the entrance to the aqueduct. A deceleration of the fluid can also be seen at the exit of the aqueduct, as the CSF enters the larger volume of the fourth ventricle. The laminar jet from the foramen of Magendie into  the SAS can also be clearly seen. A pressure drop of 2.8Pa was calculated through the aqueduct at the time during the cycle of maximum CSF velocity magnitude.

3.3.
Exit of Magendie and entrance to the Spinal SAS. CSF dynamics in the fourth ventricle were examined over the cycle at discrete time points. These are shown in figure 7, which show CSF velocity magnitude and directional vectors of the velocity field indicating flow movement at the exit of Magendie and entrance to the SAS. In the simulation CSF velocity increases as it enters the Foramen of Magendie and exits into the Spinal SAS as a defined laminar jet. During the simulation CSF velocities at the entrance to the foramen of Magendie varied between 0.40mm/s and 0.50mm/s and between 0.20mm/s to 0.30mm/s at the entrance to the SAS. Pressure at the exits of Magendie and Luschke were also examined over the cycle and are shown in figure 8. Pressures at the exits of Luschke are greater than at Magendie for a large part of the cycle. As can be observed, a time variable pressure field at the exits of Luschke and Magendie is obtained from the simulation of the pulsating CSF flow within the CNS, contradicting the constant pressure outlet boundary condition previously used by several authors (see [1,3,7,4]) on their research analysis of the ventricular system. the spine. Figure 9 shows directional vectors of the velocity field (not indicating velocity magnitude) on cross sections taken at various heights (z direction) through the spinal column and laterally (y direction) through the CNS, to show the Foramen of Magendie and SAS. These direction vectors illustrate CSF entering the spinal SAS and flowing around the spinal cord towards the base of the spine and up, towards the outlets at the top of the skull. Further analysis of CSF within the spinal column is shown in figure 10, which shows CSF velocity magnitudes on contour lines made on cross sections through the spinal column at several distances (z) along its length at various times over one cycle. The datum (z = 0mm) for distance down the spinal column is the bottom of the exit of Magendie. Whilst the position of maximum velocity at each height down the spinal column does not change over one cycle there is a noticeable acceleration and deceleration of the CSF around the spinal cord during the cycle.
Average CSF velocities along the length of the spinal column were examined and at the time in the cycle of peak CSF velocity (0.34s) it was found to fall from a maximum of 25 µm/s at the foramen of Magendie to almost zero by 50 mm below the foramen of Magendie, inside the spinal cord. At the spinal column and cord, a rigid wall boundary condition was applied. With this condition no significant recirculation was observed within the spinal column, except near to the inlet to the spinal SAS, i.e. in the region of 50 mm down from the exit of Magendie, showing the small effect of the cardiac cycle on the flow motion within the spinal sac. This is to be expected as there is likely to be a pulsatile CSF motion within the spine due to breathing [14], which has not been considered in the present work. On the other Figure 6. Velocity magnitude (mm/s) on the plane sectioning the ventricular system showing CSF motion within the ventricular system and into the SAS during the numerical simulation. Time shown is 0.34s into the cycle.
hand, almost stagnation flow at the entire spinal cord was observed in the case of not pulsating CSF flow motion, i.e. with inlet condition given only by a constant velocity corresponding to the CSF daily rate production of 500 ml per day. The results found in this work, give a good indication of the position within the spinal cord where ventricular cardiac related pulsations cease to influence the flow of the CSF. The effect of respiration pulsating of the tissues lining the spinal column as an additional driving flow force in this region is the topic of ongoing work.
4. Discussion and conclusions. Most of the previously reported numerical simulations considering the flow motion within the CNS have been limited to 2D geometries. It is difficult to examine the full aspects of CSF flow within the CNS as large anatomical regions are ignored due to the simplification of the geometry. CSF motion has been demonstrated experimentally to be truly a three dimensional problem [9] with 3D structures dominating the flow. The 3D numerical simulations presented in this paper allow a broader examination of the flow domain than has previously been presented in the literature. The pressures reported in this paper are due to movement of the CSF and should not be compared directly to values for ICP recorded in medical procedures as this value is influenced by the base pressure which has been set equal to zero in this work. A porous medium was used to model the effect of the arachnoid trabeculae in the SAS, which affects significantly the magnitude of the pressure field (see figure 5). Peak pressure through the system was found to vary between 3.77Pa and -0.92Pa over the cycle. These values are slightly higher than in the recently presented numerical simulation data by Cheng  [6], who calculated maximum pressure within the lateral ventricles to be in the range of 2Pa. The higher pressure in our simulation is due to the flow resistance of the porous medium, which as shown in figure 5 increases the pressure in the system by approximately 80% when compared to simulations without the porous medium resistance. The pressure drop through the aqueduct calculated at the time of peak CSF velocity was found to be in the region of 2.8Pa. This value is slightly higher than previously reported, (1-2Pa) and is again due to the modelling of the porous medium within the SAS. Pressure at the exits of Luschke and Magendie were also examined, with peak pressure approximately 25% higher at the exit of Luschke, likely due to the geometry of the exit of Luschke modelled in the simulation. The effect of pulsing CSF was found to be marginal in the lateral ventricles, where the volume of CSF is large and therefore less responsive to acceleration of fluid by the displacement of the choroid plexus. Once CSF enters the third ventricle there is a noticeable acceleration of the CSF in the centre of the third ventricle. Pulsating CSF within the aqueduct fluctuates between 8.05mm/s and 11.85mm/s ( fig  4). This value is within the velocity range reported in both the medical literature ( [15,16,17,18,19]) and computational literature ( [1,7,8,6]). Peak CSF velocity within the aqueduct is 1% of cycle out of phase with the peak velocity at the choroid plexus. This can be explained by the local inertia effects in transient flow. A maximum Reynolds number in the region of 12 was calculated within the cerebral aqueduct. Simulation of pulsatile CSF motion within the spinal SAS (figures 9 and 10) shows CSF moving in three directions from its exit at the foramen of Magendie. CSF moves up towards the exits at the Pacchionian granulations, down the spinal column and finally around the spinal cord. The upward trend is to be expected as CSF flows towards the pressure outlet exits at the top of the skull and has been predicted by Greitz et al (1993) [17] who suggest that CSF motion down the length of the spinal column should resemble a meandering river, which moves upwards. The flow recirculation around the spinal cord is a consequence of the closeness of this part of the domain and the opening at the Pacchionian granulations located at the top of the skull. Part of the fluid leaving the exit of Magendie flows into spinal cord and after having an internal recirculation moves back in a spiral motion into the SAS to reach the pressure outlets at the top of the skull. Velocity within the spinal column was found to fall off significantly after approximately 50mm into the column. It is expected that stronger flow recirculation at the spinal cord will be observed if the effect of the respiratory cycle is considered and not only the cardiac cycle. This result does however suggest the likely influence limit of the cardiac cycle of spinal CSF flow. The present numerical simulation is able to predict completely the flow field within the ventricular system and cranial SAS, as well as partially describing the flow field within the spinal column. In order to improve the modelling further, accurate information regarding respiratory effects of spinal wall motion are required to improve the boundary conditions within the spinal region of the CNS