- Email: [email protected]
Computational fluid dynamics modelling of cerebrospinal fluid pressure in Chiari malformation and syringomyelia
Journal of Biomechanics 46 (2013) 1801–1809
Contents lists available at SciVerse ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Computational fluid dynamics modelling of cerebrospinal fluid pressure in Chiari malformation and syringomyelia Elizabeth C. Clarke a,n, David F. Fletcher b, Marcus A. Stoodley c, Lynne E. Bilston d a Murray Maxwell Biomechanics Laboratory, Kolling Institute of Medical Research, Sydney Medical School, University of Sydney, Level 10, Kolling Building 6, RNS Hospital, St Leonards, NSW 2065, Australia b School of Chemical and Biomolecular Engineering, University of Sydney, Australia c The Australian School of Advanced Medicine, Macquarie University, Australia d Neuroscience Research Australia and Prince of Wales Clinical School, University of New South Wales, Australia
art ic l e i nf o
a b s t r a c t
Article history: Accepted 20 May 2013
The pathogenesis of syringomyelia in association with Chiari malformation (CM) is unclear. Studies of patients with CM have shown alterations in the CSF velocity profile and these could contribute to syrinx development or enlargement. Few studies have considered the fluid mechanics of CM patients with and without syringomyelia separately. Three subject-specific CFD models were developed for a normal participant, a CM patient with syringomyelia and a CM patient without syringomyelia. Model geometries, CSF flow rate data and CSF velocity validation data were collected from MRI scans of the 3 subjects. The predicted peak CSF pressure was compared for the 3 models. An extension of the study performed geometry and flow substitution to investigate the relative effects of anatomy and CSF flow profile on resulting spinal CSF pressure. Based on 50 monitoring locations for each of the models, the CM models had significantly higher magnitude (p o0.01) peak CSF pressure compared with normal. When using the same CSF input flow waveform, changing the upper spinal geometry changed the magnitude of the CSF pressure gradient, and when using the same upper spinal geometry, changing the input flow waveform changed the timing of the peak pressure. This study may assist in understanding syringomyelia mechanisms and relative effects of CSF velocity profile and spinal geometry on CSF pressure. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Syringomyelia Chiari malformation Cerebrospinal fluid (CSF) Computational fluid dynamics (CFD) Magnetic resonance imaging (MRI)
1. Introduction Chiari malformation Type I (CM) is an abnormality characterised by protrusion of the cerebellar tonsils through the foramen magnum. Approximately two thirds (Pinna et al., 2000) of patients with CM also develop syringomyelia, where fluid filled cysts (syrinxes) develop in the spinal cord. Many patients report headaches, pain, and motor and sensory disturbances. Surgery is the only available treatment but syrinxes may persist or enlarge after surgery, and it is not clear why some patients with CM develop syringomyelia and others do not. Understanding the mechanisms of syringomyelia has the potential to improve treatment for these patients. Early theories of syrinx development in CM proposed that herniation of the cerebellar tonsils causes increased resistance across the foramen magnum thereby preventing spinal cord extracellular fluid outflow or driving fluid into the central canal (Heiss et al., 1999; Oldfield et al., 1994; Olivero and Dinh, 1992). However a syrinx may be under high pressure and it is not clear how the CSF
n
Corresponding author. Tel.: +61 299 264821; fax: +61 299 265266. E-mail address: [email protected] (E.C. Clarke).
0021-9290/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jbiomech.2013.05.013
would accumulate in the central canal against a pressure gradient. Another theory proposes that syringomyelia develops passively as a result of the Venturi effect, where narrowed regions of the CSF pathway result in higher velocity and lower pressure (Greitz, 2006), however clinical studies of CM patients have reported significant reductions in syrinx size following surgery without significant reduction in CSF velocity near the syrinx (Mauer et al., 2011), and the velocity of the fluid flow is low, so such an effect is likely to be small. Experimental studies have demonstrated that fluid flows into the cord and syrinx via the perivascular spaces (Rennels et al., 1990; Rennels et al., 1985; Stoodley et al., 1996). It is possible that the pressure gradient and/or timing of arterial pulsation-dependent CSF flow could be a driving force of syrinx fluid from the perivascular spaces into the central canal. A recent computational simulation demonstrated that even a moderate shift in the timing of the CSF pressure pulse relative to arterial pulse could influence the mass flow rate of CSF into the perivascular spaces (Bilston et al., 2010). More recently, a follow up study by Elliott et al. demonstrated that this mechanism requires that the spinal cord have non-zero compliance with the vascular bed (Elliott et al., 2011). While MRI scanning can provide structural anatomical, velocity and flow rate data, and has been used to demonstrate the contribution
1802
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
of arterial and venous pulsations to CSF pulsation (Bhadelia et al., 1997), it is not possible to measure non-invasively spinal CSF pressures in humans. Invasive CSF pressure measurements in patients with CM have demonstrated dissociation between cranial and spinal CSF pressures (Williams, 1981a; Williams, 1981b), and increased CSF pressure in CM patients that reduced following posterior fossa decompression surgery (Heiss et al., 1999). However these studies did not compare patients with and without syringomyelia and did not investigate any temporal changes in the CSF waveform (e.g. delays in peak CSF pressure or delays in onset of systolic flush). A very recent study by Heiss et al. has compared CSF pressures in CM patients with and without syringomyelia (Heiss et al., 2012). Subject-specific computational simulations provide an opportunity to non-invasively study CSF pressure in patients with CM, with and without syringomyelia. Previous computational simulations of the spinal subarachnoid space (SSAS) have investigated the mechanisms of syrinx development in conditions such as arachnoiditis, finding that a porous obstruction in the SSAS changes timing and magnitude of peak CSF pressure (Cheng et al., 2012). Another study used computational modelling to investigate the effect of tonsillar herniation on CSF pressure by extending normal tonsil geometry inferiorly, finding that modified tonsil geometry increased the CSF pressure gradient (Linge et al., 2011). To our knowledge this is the first study to investigate the cyclic CSF pressure profile using subject-specific models of Chiari patients with and without syringomyelia.
2. Methods 2.1. MRI scanning and flow measurements The University of NSW Human Research and Ethics Committee approved all experimental protocols, and all participants gave written informed consent to enter
the study. Three participants underwent MRI scanning of the head and neck; one patient with CM and syringomyelia (female, 58 years, 80 kg), one patient with CM but no syringomyelia (female, 60 years, 65 kg), and one healthy control subject (female, 60 years, 80 kg). A 3D sagittal isotropic (0.94 mm) T1-weighted turbo field echo MRI of the head and neck, and cardiac-gated cine phase-contrast (PC) MRI scans were performed using a 3 T MRI scanner (Philips Achieva TX, Best, Netherlands). The parameters for the 3D isotropic MRI scan were flip angle ¼ 81, matrix ¼288 288, FOV ¼270 270 169, TR/TE¼ 5.5/2.4 m s, slice thickness¼0.94 mm. For the PC-MRI scans, retrospective cardiac gating was achieved using VCG leads, and 30 cardiac phase images were collected, measured from the onset of the R wave. PC-MRI scans were performed at the following cranial/spinal levels and encoding velocities (Venc) using an (axial) scan plane aligned perpendicular to the CSF flow: 5 mm cranial to the tip of the cerebellar tonsils (Venc 10 cm/s); mid-vertebral C2 (Venc 9 cm/s); and midvertebral C5 (Venc 13 cm/s) (Fig. 1A). The spatial resolution of the PC-MRI scans was 0.694–0.977 mm2. Encoding velocities were determined during pilot studies in order to prevent aliasing. Other typical PC-MRI parameters were flip angle ¼101, matrix ¼240 176, FOV ¼ 250 250, TR/TE¼21/6.8 m s, slice thickness¼ 5 mm. CSF velocity and flow rate were measured from PC-MRI scans using the freely available analysis software Segment (Heiberg et al., 2010). For the tonsil-level scan the entire CSF cross-section was segmented and analysed for flow rate data (Fig. 1D and E), and for the C2- and C5-level scans velocity data were analysed using 10 regions of interest (2 mm diameter) spaced evenly around the SSAS (Fig. 1B and C). The tonsil-level flow rate was represented by a 6th order Fourier series for use as a flow input function for the models (Fig. 5B), and the C2 and C5 velocity data were used for model validation (Fig. 3).
2.2. Computational modelling The spinal cord, spinal canal and cerebellum were segmented manually from the 3D MRI scans using SURFdriver software. These point cloud segmentations were imported into Rhinoceros (v3.0, Robert McNeal and associates, WA) and fitted with closed non-rational B-splines of degree 3 passing through the points for smoothing. These were then fitted with a surface (loft) and capped to create volumes. The CSF volume was created by Boolean subtraction of the cerebellum and spinal cord volumes from the spinal canal volume, then trimmed 5 mm cranial to the base of the cerebellar tonsils and 95 mm caudal to the base of the tonsils (Fig. 2A). The cranial end of the CSF volume was extended upwards by 5 mm to
Fig. 1. Example MRI scans from a healthy control volunteer showing (A) the approximate alignment of the cine phase-contrast MRI scan planes (cerebellar tonsils, mid-C2 and mid-C5) aligned perpendicular to the CSF flow, (B) an example anatomical scan through the mid-C2, (C) a magnified view of a mid-C2 level phase image (as indicated in B) depicting approximate locations of the 10 regions used for validation of the models against MRI-measured CSF velocity data, (D) an example anatomical scan through the cerebellar tonsil level, and (E) a magnified view of a tonsil level phase image (as indicated in D) depicting manual segmentation of the entire CSF region.
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
Normal
With Syrinx
No Syrinx
1803
cardiac cycle time for an indication of the relative peak pressure timing between the models. Since there are known differences in both hindbrain anatomy and CSF flow profile for CM patients, we performed an additional substitution step (changing either the anatomy or flow profile while holding the other constant) in order to investigate the relative effects of the anatomy and CSF flow on resulting CSF pressure profile. For example, it is possible that the CM patient with syringomyelia developed the syrinx due to the hindbrain anatomy, so by normalising the input CSF flow rate function we aimed to explore the effects of hindbrain and upper spinal geometry on CSF pressure profiles between the CM patients with and without syringomyelia. The geometries from each of the CM participants were used with the Fourier series flow-rate function calculated from the PC-MRI scan (as described above) of the healthy participant (termed geometry substitution), and the healthy participant geometry was used with the Fourier series flow-rate functions calculated from the PC-MRI scans (as described above) for each of the CM participants (termed flow substitution). The pressure monitoring locations from the original geometries were used for analysis of these cases, and again a uniform velocity across the inlet was assumed.
3. Results The simulation velocities and MRI measured velocities were well correlated (examples shown in Fig. 3). There were no significant differences between MRI and model peak velocities for any of the 3 cases (p ¼0.1 with syringomyelia, p ¼0.7 without syringomyelia and p¼ 0.07 normal participant). There were no significant differences between MRI and model peak velocity timing for CM with syringomyelia, (p ¼0.1) or normal participant (p ¼0.5). Fig. 4 shows examples of the velocity (Fig. 4 A and B) and pressure gradients (Fig. 4C and D) at 2 time steps for a representative model (geometry substitution of no-syrinx patient with normal participant input flow rate profile). 3.1. CSF hydrodynamics in Chiari malformation
Fig. 2. (A) The CSF volumes for the 3 cases from the inlet side (cerebellar tonsil level) and (B) an example of the model mesh at the outlet showing the inflation layers.
provide an entry zone ahead of rapid changes in cross-sectional shape. The volumes were then imported into ANSYS CFX (v13, ANSYS Inc., PA). Tetrahedral meshing of the CSF volumes was performed using ANSYS with 5 layers of prismatic inflation at the walls (Fig. 2B). Proximity meshing was used to ensure that small flow paths were resolved. CSF was modelled as a Newtonian fluid with viscosity 0.8 m Pa s and density 1000 kg/m3 (Bloomfield et al., 1998). For a typical velocity of 0.1 m/s and hydraulic diameter of 5 mm the Reynolds number was calculated as 625, so laminar flow was assumed. The tonsil-level Fourier series flow rate function (using original timescale information from PC-MRI scans) was assigned as the input at the cranial end of the model, with a uniform velocity across the inlet (see Section 4). The average gauge pressure at the caudal end of the models was set to 0 Pa. Velocity monitoring points were manually placed in the models at coordinates approximately corresponding to the monitoring points from the MRI scans (Fig. 1). Pressure monitors were placed at 5 evenly spaced cranialcaudal levels along the model, with 10 points placed evenly around the SSAS. In all cases an initial steady state run was performed to establish the flow and provide a consistent starting flow-field for the transient run. A transient solution was performed and data were collected on the second cycle, with the solution converged to a normalised residual of o10−5 at each time step (5 m s) for all variables. The CFD code used is a 3D Navier–Stokes solver that uses a finite volume method applied to a vertex based discretisation. The equations were solved using a coupled solver, with the pressure–velocity coupling achieved by a modified Rhie– Chow scheme. The spatial derivatives were calculated using a bounded second order scheme and the 2nd order backwards Euler scheme was used for time derivatives. Double precision arithmetic was used to avoid any problems from rounding errors. Validation of the models was performed using the original timescales by comparing the simulated velocity profiles (magnitude and time of peak velocity) against the MRI measured velocity profiles (magnitude and time of peak velocity) at each of the 20 corresponding locations, for each participant, using paired t-tests. The time axes of output pressures were also normalised to proportion of the
Fig. 5 shows simulated CSF pressures, relative to the outlet at the caudal end of the models, in the anterior SSAS at the C2 level for the 3 models, both on the normalised timescale Fig. 5(C) and the original timescale Fig. 5(D) and the corresponding flow rate data Fig. 5(A–B). All 3 cases showed an initial peak in CSF pressure during systole. In this location, the peak pressure was 26.1 Pa in the model for the healthy participant, 33.6 Pa for the patient with syringomyelia, and 56.6 Pa in the CM patient without syringomyelia. The time of peak pressure was nearly 40%, or approximately 50 m s, earlier than normal in the patient without syringomyelia. These patterns were consistent across other pressure monitoring locations. 3.2. Flow substitution Fig. 6(A–B) shows simulated CSF pressures at one corresponding location (anterior, C2-level) for the normal model geometry with input flow functions from the 3 participants, both on the normalised timescale (A) and the original timescale (B), and the corresponding flow rate data are shown in Fig. 5(A–B). The times of peak pressures for the flow substitution models aligned with the times of peak pressures for the corresponding participant original models (i.e. in Fig. 5A with Fig. 6A and Fig. 5B with Fig. 6B). As with the original models shown in Fig. 5C the peak pressure arrived earlier in the cardiac cycle for the flow substitution model of the CM patient without syringomyelia (shown in Fig. 6A). 3.3. Geometry substitution Fig. 6(C–D) shows simulated CSF pressures at one corresponding location (anterior, C2-level) for the 3 model geometries with input flow functions from the healthy participant, both on the normalised timescale (C) and the original timescale (D). Corresponding flow rate data for the normal participant is shown in
1804
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
Model MRI data
CSF velocity (m/s)
CSF velocity (m/s)
Model MRI data
Time (% cardiac cycle)
Time (% cardiac cycle)
CSF velocity (m/s)
Model MRI data
Time (% cardiac cycle) Fig. 3. MRI-measured and simulation velocity–time data used for model validation. Each figure shows data from the C2 level (left panels) and C5 level (right panels) of the 5 regions around the subarachnoid space (left or right only shown) for the Chiari malformation patient with syringomyelia (A), the Chiari malformation patient without syringomyelia (B), and the normal participant (C). Data from the model simulations are shown as black solid lines and the data from the MRI scans are shown as grey dashed lines. The peak velocities used for validation statistics are indicated by the triangle symbols (simulation data are shown in black and MRI scan data are shown in grey).
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
1805
Fig. 4. Representative model outputs showing (A–B) CSF velocity and (C–D) CSF pressures at 2 time steps from one simulation.
Fig. 5(A–B). As with the original models shown in Fig. 5(C–D) the peak pressure was higher for the model geometry of the CM patient without syringomyelia than the model geometry of the CM patient with syringomyelia, and both CM geometries produced higher pressure than healthy geometry.
4. Discussion 4.1. Implications of results for syringomyelia development in CM The results of this study show that the model of the CM patient with syringomyelia had slightly elevated peak pressure compared with normal, but closer to the timing of peak pressure in the normal participant model. Interestingly, the CM patient without syringomyelia had earlier and higher magnitude peak pressure than both the normal participant and the patient with syringomyelia. The current study included 3 subjects, and while flow profiles were typical of their respective groups in a larger MRI study (Clarke et al., 2013), individual patient anatomy may not be, so the pressure results from these subjects may not be representative of their groups. Previous simulations (Bilston et al., 2010) have suggested
that a modest shift in the timing of the CSF pressure profile with respect to the arterial pressure profile may enhance CSF flow through the perivascular spaces into the central canal of the spinal cord, and this may contribute to formation or enlargement of a syrinx. More specifically, if the peak SSAS pressure occurs at a time in the arterial cycle when spinal arteriole pressure is low (periarterial spaces at their widest), flow into the spinal cord may be enhanced. This would encourage fluid accumulation in the spinal cord, unless fluid outflow is able to accommodate this additional influx. Similarly, if the lowest SSAS pressures occur when spinal arteriole pressure is high (peri-arterial spaces at their narrowest), this may impede backflow through the peri-arterial spaces encouraging fluid accumulation in the cord. This could lead to the “presyrinx state” observed clinically, and eventually to a syrinx. The results from the flow and geometry substitution studies suggest that the pressure magnitude is strongly influenced by the anatomy at the craniocervical junction and in the spinal SSAS. It is difficult to compare our predicted pressure magnitudes (peak 20–60 Pa) with in vivo pressure measurements due to the location/nature of reference pressure points, penetration of the dura for pressure sensors, the rigid in vitro nature of our simulations with absence of structures such as spinal nerve. An in vivo
1806
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
Normal participant Chiari with syrinx Chiari without syrinx
Normal participant Chiari with syrinx Chiari without syrinx
Normal participant Chiari with syrinx Chiari without syrinx
Normal participant Chiari with syrinx Chiari without syrinx
Fig. 5. (A) and (B) shows CSF flow rate data from the tonsil-level PC-MRI scans (segmented as shown in Fig. 1E) alongside the Fourier functions that were fit to these data points and used as input flow rates for the simulations. MRI data are shown as symbols for the normal (triangle), Chiari with syrinx (circle), and Chiari without syrinx (square) participants, and Fourier fits are shown as lines for the normal (sold black), Chiari with syrinx (solid grey), and Chiari without syrinx (dashed black) participants. Both the normalised timescale (A) and original timescale (B) are shown. (C) and (D) shows simulated CSF pressures, relative to the outlet at the caudal end of the models, at one representative location (anterior, C2-level) for the normal healthy participant (solid black line), the Chiari malformation participant with syringomyelia (solid grey line) and the Chiari malformation patient without syringomyelia (dashed black line). Both the normalised timescale (C) and original timescale (D) are shown.
study (Williams, 1981b) measured the mean amplitude of CSF under cardiac pulsation as 3.86 mm Hg (514 Pa) at the lumbar relative to ventricular level, and another (Heiss et al., 1999) measured mean cervical subarachnoid pressure as 2 mm Hg (267 Pa) at the cervical relative to the ventricular level. Our study considers only a small length of the spinal canal (pressures relative the lower cervical level). The monitoring and reference pressure points are in different locations in our study and these in vivo examples but it should be noted that our predicted pressures are significantly lower than those reported in the in vivo studies and care should be taken if using the absolute magnitudes of pressures reported in our study. The current study is however useful for making comparisons between the 3 patient groups. The findings of the current study agree with previous computational fluid dynamics studies. One study demonstrated that simulated tonsillar herniation increased spinal CSF pressure gradient (Linge et al., 2011). Other studies have used lumped parameter and coupled hydrodynamics to investigate CSF pressure and flow environments (Kim and Cirovic, 2011). A study by Elliott et al. (2011) investigated mechanisms of syringomyelia and effects of simulated shunt and bypass surgery. Interestingly they found that, while the simulation of arachnoiditis increased the steady state pressure mean and amplitude, including a syrinx lowered the mean and amplitude of the steady state pressure back towards normal. Similar stenosis modelling studies by Martin et al. found that the inclusion of a syrinx reduced the mean and amplitude of SSAS pressure (Martin et al., 2010; Martin and Loth, 2009). This observation is consistent with our findings that CSF peak pressure is closer to normal in patients with syringomyelia than without.
Several previous studies have combined MRI and modelling to investigate CSF dynamics for clinical purposes. For example, Hsu et al. (2012) developed an intricate computational model to investigate CSF pulsation effects on intrathecal drug distribution; Alperin et al. (1996) investigated perturbations in arterial and venous flow on CSF and spinal tissues, and Sweetman and Linninger (2011) developed a complex 3D fluid–structure interaction model of the entire CNS to investigate interactions between pulsating vasculature, compliant tissues, and CSF flow. To our knowledge the models presented in the current study are the first subject-specific CFD models for investigating CSF pressure in CM patients with and without syringomyelia. The findings of this study are important because the only treatment for symptomatic CM is posteria fossa decompression surgery and this is sometimes unsuccessful. The surgery certainly changes geometry and therefore likely changes the CSF pressure gradient. However, it is still not known whether this change in geometry also changes the flow-rate profile, which may also affect the CSF pressure profile. This may help to explain why decompression surgery is sometime unsuccessful and further investigation may improve treatment for these patients. 4.2. Limitations of the current study The models presented here are CFD models with solid boundaries. We hypothesised that deformation of the spinal cord in response to these pressures is likely to be very small, because the magnitudes of the pressures involved (o 100 Pa) are orders of magnitude lower than the elastic modulus of the spinal cord
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
1807
Normal geometry, normal input flow Normal geometry, Chiari with syrinx flow Normal geometry, Chiari without syrinx flow
Normal geometry, normal input flow Normal geometry, Chiari with syrinx flow Normal geometry, Chiari without syrinx flow
Normal geometry, normal input flow Chiari with syrinx geometry, normal flow Chiari without syrinx geometry, normal flow
Normal geometry, normal input flow Chiari with syrinx geometry, normal flow Chiari without syrinx geometry, normal flow
Fig. 6. (A) shows simulated CSF pressures on a normalised timescale at one representative location (anterior, C2-level) for the flow substitution models, where the same model geometry (normal participant) was used but the input flow rate function was different (the 3 flow rate functions are shown in Fig. 5B). The solid black line is derived from the normal participant geometry with normal participant input flow rate function, the solid grey line is derived from the normal geometry with flow rate function from the patient with syringomyelia, and the dashed black line is derived from the normal geometry with flow rate function from the patient without syringomyelia. Arrows indicate the change in timing of the peak pressures for the 3 models. (B) Shows the same data as in (A), but using the original timescale, rather than the normalised timescale. (C) Shows simulated CSF pressures on a normalised timescale at one representative location (anterior, C2-level) for the geometry substitution models, where the same input flow rate function (normal participant, shown in Fig. 5B) was used but model geometry was different (the 3 geometries are shown in Fig. 2 A). The solid black line is derived from the normal participant geometry with normal participant input flow rate function, the solid grey line is derived from the geometry of the patient with syringomyelia and normal participant input flow rate function, and the dashed black line is derived from the geometry of the patient without syringomyelia and normal participant input flow rate function. Arrows indicate the change in peak pressure magnitude for the 3 models. (D) Shows the same data as in (C), but using the original timescale, rather than the normalised timescale.
(approx. 1 MPa) (Bilston and Thibault, 1996). We confirmed this using fluid–structure interaction modelling for the normal participant and found no discernable difference in the pressure wave timing or magnitude, and micron-order displacement of the spinal cord. We also measured tonsil motion from the tonsil-level PC-MRI (by segmenting the tonsils and calculating their displacement from the velocity–time curves) in these subjects and did not detect any cranial-caudal tonsil motion throughout the cardiac cycle. We therefore believe these assumptions will have minimal impact on the estimated pressures. While the locations of the C2 and C5 PC-MRI scan planes were known with respect to the model geometries, the 10 velocity probes on each plane (used for validation) were placed manually. Despite this the agreement between MRI and simulations velocities is good (Fig. 3). Also probe volumes for velocity validation different sizes (point probes were used for simulations whereas 2 mm diameter probes were used for MRI scans to include some spatial averaging over multiple pixels). This difference in size is unlikely to affect the overall patterns of results in this study, but closer matching of probe volumes may have improved the validation at narrow regions of the SSAS. Overall, however the validation results match well (Fig. 3) and there were no significant differences in peak systolic velocity magnitude and no significant differences in timing of peak systolic velocity for 2 of the 3 cases.
Furthermore, there are different types of error in CSF velocity– time profiles from PC-MRI scanning depending on the type of gating used. Retrospective cardiac gating was used in order to ensure that the entire cycle was captured and to minimise sensitivity to irregular heart rates (Lenz et al., 1989). However, we did not simultaneously measure arterial and CSF waveforms but rather measured the CSF waveforms relative to the R–R interval. It is important to keep this in mind when comparing the cycle timing for the 3 models presented here. Previous studies using cine PC-MRI investigated relative timing of pulses in the epidural vessels and CSF simultaneously, finding no significant differences in the arterial and CSF peak flow timing at the C2-C3 levels (Baledent et al., 2001; Wagshul et al., 2006). The inlet flow-rate condition to these models assumed a uniform velocity across the inlet because the inlet-level MRI low-rate data was sampled by segmenting the entire SSAS (i.e. the mean flow-rate from all SSAS pixels). In reality there is some variation in flow-rate across the axial plane, but data from our larger cohort PC-MRI flow study (Clarke et al., 2013) found no significant differences in peak velocities between probes placed around the SSAS, so we believe the assumption of uniform inlet velocity would have negligible effect on the outcomes of the models. The 3 participants were chosen because they had flow-rate and velocity profiles that were representative of their respective
1808
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
3 groups in a larger cohort quantitative flow study (Clarke et al., 2013). From that study we found patients without syringomyelia had significantly earlier and faster peak systolic CSF flow than both normal volunteers and patients with syringomyelia. It follows that in the current study, the peak CSF pressure occurred earlier in the model for the patient without syringomyelia compared with the other 2 models. However a recently published similar study reported that Chiari patients had significantly later peak systolic CSF flow than normal (Bunck et al., 2012). Therefore it is important to recognise that the models presented here may not be representative of all patients in their respective groups, however the outcomes of the models may be useful for making inferences about the relative effects of upper spinal geometry and CSF flow rate on the timing and magnitude of peak CSF pressure. Also, the geometry and flow substitution models are not true subjectspecific models with matched flow and geometry, and therefore cannot be validated experimentally and are presented only as an exploration of the role of the anatomy and flow profiles on resulting SSAS pressures. It should also be noted that MRI segmentations were smoothed to assist formation of smooth model geometries, the models do not include epidural blood vessels or nerve roots, and do not simulate the pulsatile blood flow in the epidural vessels. Surface smoothing and exclusion of small structures leads to a better quality computational mesh, which was deemed preferable over working with raw data that necessarily contains imperfections in geometry. While these features may affect localised pressure magnitudes we would not expect these features to change the overall patterns observed in these models and we would expect the exclusion of these features to have similar effects across the 3 models. Also, spine and spinal cord compliance were assumed equal for the 3 participants, but it is possible that these may be different (this is not known) and this may affect the magnitude and timing of the CSF pressure profiles. For example, Martin and colleagues have demonstrated that changing spinal canal compliance and including the presence of a modelled syrinx influences the magnitude of CSF pressure and timing (Martin et al., 2010; Martin and Loth, 2009). The models presented here are simplified in terms of tissue compliance and inclusion of spinal structures, but do include the subject-specific spinal geometries and inlet flow conditions. The situation appears more complex than either study considered alone but it was interesting to note that Martin et al. found the presence of a syrinx to increase the magnitude and delay the peak CSF pressure compared with an equivalent model without a syrinx (Martin et al., 2010). In the current study we also found the model for the patient with a syrinx to have later peak CSF pressure than the patient without a syrinx but reduced peak CSF pressure. It is acknowledged that in the current study using a short section of upper spinal geometry that timing of the pressure waves may be affected mildly by fluid–structure interaction, but the complete study of CSF dynamics in the entire (closed) CNS space needs to account for the interaction of CSF filled spaces with deformable structures such as soft brain parenchyma boundary, the dura and the blood vessels.
Conflicts of Interest No conflicts of interest to declare.
Acknowledgements A research project grant from the Column of Hope provided funding for the study. ECC and LEB are supported by NHMRC research fellowships.
References Alperin, N., Vikingstad, E.M., Gomez-Anson, B., Levin, D.N., 1996. Hemodynamically independent analysis of cerebrospinal fluid and brain motion observed with dynamic phase contrast MRI. Magnetic Resonance in Medicine 35, 741–754. Baledent, O., Henry-Feugeas, M.C., Idy-Peretti, I., 2001. Cerebrospinal fluid dynamics and relation with blood flow: a magnetic resonance study with semiautomated cerebrospinal fluid segmentation. Investigative Radiology 36, 368–377. 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, 258–264. Bilston, L.E., Stoodley, M.A., Fletcher, D.F., 2010. The influence of the relative timing of arterial and subarachnoid space pulse waves on spinal perivascular cerebrospinal fluid flow as a possible factor in syrinx development. Journal of Neurosurgery 112, 808–813. Bilston, L.E., Thibault, L.E., 1996. The mechanical properties of the human cervical spinal cord in vitro. Annals of Biomedical Engineering 24, 67–74. Bloomfield, I.G., Johnston, I.H., Bilston, L.E., 1998. Effects of proteins, blood cells and glucose on the viscosity of cerebrospinal fluid. Pediatric Neurosurgery 28, 246–251. Bunck, A.C., Kroeger, J.R., Juettner, A., Brentrup, A., Fiedler, B., Crelier, G.R., Martin, B. A., Heindel, W., Maintz, D., Schwindt, W., Niederstadt, T., 2012. Magnetic resonance 4D flow analysis of cerebrospinal fluid dynamics in Chiari I malformation with and without syringomyelia. European Radiology. Cheng, S., Stoodley, M.A., Wong, J., Hemley, S., Fletcher, D.F., Bilston, L.E., 2012. The presence of arachnoiditis affects the characteristics of CSF flow in the spinal subarachnoid space: A modelling study. Journal of Biomechanics 45, 1186–1191. Clarke, E., Stoodley, M., Bilston, L., 2013. Changes in temporal flow characteristics of CSF in Chiari I Malformation with and without syringomyelia: Implications for theory of syrinx development. Journal of Neurosurgery, 118(5):1135-1140 http://dx.doi.org/10.3171/2013.2.JNS12759. Epub 2013 Mar 15. Elliott, N.S., Lockerby, D.A., Brodbelt, A.R., 2011. A lumped-parameter model of the cerebrospinal system for investigating arterial-driven flow in posttraumatic syringomyelia. Medical Engineering and Physics 33, 874–882. Greitz, D., 2006. Unraveling the riddle of syringomyelia. Neurosurgical Review 29, 251–263, discussion 264. Heiberg, E., Sjogren, J., Ugander, M., Carlsson, M., Engblom, H., Arheden, H., 2010. Design and validation of segment—freely available software for cardiovascular image analysis. BMC Medical Imaging 10, 1. Heiss, J.D., Patronas, N., DeVroom, H.L., Shawker, T., Ennis, R., Kammerer, W., Eidsath, A., Talbot, T., Morris, J., Eskioglu, E., Oldfield, E.H., 1999. Elucidating the pathophysiology of syringomyelia. Journal of Neurosurgery 91, 553–562. Heiss, J.D., Snyder, K., Peterson, M.M., Patronas, N.J., Butman, J.A., Smith, R.K., Devroom, H.L., Sansur, C.A., Eskioglu, E., Kammerer, W.A., Oldfield, E.H., 2012. Pathophysiology of primary spinal syringomyelia. Journal of Neurosurgery: Spine 17, 367–380. Hsu, Y., Hettiarachchi, H.D., Zhu, D.C., Linninger, A.A., 2012. The frequency and magnitude of cerebrospinal fluid pulsations influence intrathecal drug distribution: key factors for interpatient variability. Anesthesia and Analgesia 115, 386–394. Kim, M., Cirovic, S., 2011. A computational model of the cerebrospinal fluid system incorporating lumped-parameter cranial compartment and one-dimensional distributed spinal compartment. Journal of Biorheology 25, 78–87. Lenz, G.W., Haacke, E.M., White, R.D., 1989. Retrospective cardiac gating: a review of technical aspects and future directions. Magnetic Resonance Imaging 7, 445–455. Linge, S.O., Haughton, V., Lovgren, A.E., Mardal, K.A., Helgeland, A., Langtangen, H.P., 2011. Effect of tonsillar herniation on cyclic CSF flow studied with computational flow analysis. AJNR American Journal of Neuroradiology 32, 1474–1481. Martin, B.A., Labuda, R., Royston, T.J., Oshinski, J.N., Iskandar, B., Loth, F., 2010. Spinal subarachnoid space pressure measurements in an in vitro spinal stenosis model: implications on syringomyelia theories. Journal of Biomechanical Engineering 132, 111007. Martin, B.A., Loth, F., 2009. The influence of coughing on cerebrospinal fluid pressure in an in vitro syringomyelia model with spinal subarachnoid space stenosis. Cerebrospinal Fluid Research 6, 17. Mauer, U.M., Gottschalk, A., Mueller, C., Weselek, L., Kunz, U., Schulz, C., 2011. Standard and cardiac-gated phase-contrast magnetic resonance imaging in the clinical course of patients with Chiari malformation Type I. Neurosurgical Focus 31, E5. Oldfield, E.H., Muraszko, K., Shawker, T.H., Patronas, N.J., 1994. Pathophysiology of syringomyelia associated with Chiari I malformation of the cerebellar tonsils. Implications for diagnosis and treatment. Journal of Neurosurgery 80, 3–15. Olivero, W., Dinh, D., 1992. Chiari I malformation with traumatic syringomyelia and spontaneous resolution: case report and literature review. Neurosurgery 30, 758–760. Pinna, G., Alessandrini, F., Alfieri, A., Rossi, M., Bricolo, A., 2000. Cerebrospinal fluid flow dynamics study in Chiari I malformation: implications for syrinx formation. Neurosurgical Focus 8, E3. Rennels, M.L., Blaumanis, O.R., Grady, P.A., 1990. Rapid solute transport throughout the brain via paravascular fluid pathways. Advances in Neurology 52 431–439
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
Rennels, M.L., Gregory, T.F., Blaumanis, O.R., Fujimoto, K., Grady, P.A., 1985. Evidence for a ‘paravascular’ fluid circulation in the mammalian central nervous system, provided by the rapid distribution of tracer protein throughout the brain from the subarachnoid space. Brain Research 326, 47–63. Stoodley, M.A., Jones, N.R., Brown, C.J., 1996. Evidence for rapid fluid flow from the subarachnoid space into the spinal cord central canal in the rat. Brain Research 707, 155–164. Sweetman, B., Linninger, A.A., 2011. Cerebrospinal fluid flow dynamics in the central nervous system. Annals of Biomedical Engineering 39, 484–496. 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. Journal of Neurosurgery 104, 810–819.
1809
Williams, B., 1981a. Simultaneous cerebral and spinal-fluid pressure recordings 1. Technique, physiology, and normal results. Acta Neurochirurgica 58, 167–185. Williams, B., 1981b. Simultaneous cerebral and spinal-fluid pressure recordings 2. Cerebrospinal dissociation with lesions at the foramen magnum. Acta Neurochirurgica 59, 123–142.
Contents lists available at SciVerse ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Computational fluid dynamics modelling of cerebrospinal fluid pressure in Chiari malformation and syringomyelia Elizabeth C. Clarke a,n, David F. Fletcher b, Marcus A. Stoodley c, Lynne E. Bilston d a Murray Maxwell Biomechanics Laboratory, Kolling Institute of Medical Research, Sydney Medical School, University of Sydney, Level 10, Kolling Building 6, RNS Hospital, St Leonards, NSW 2065, Australia b School of Chemical and Biomolecular Engineering, University of Sydney, Australia c The Australian School of Advanced Medicine, Macquarie University, Australia d Neuroscience Research Australia and Prince of Wales Clinical School, University of New South Wales, Australia
art ic l e i nf o
a b s t r a c t
Article history: Accepted 20 May 2013
The pathogenesis of syringomyelia in association with Chiari malformation (CM) is unclear. Studies of patients with CM have shown alterations in the CSF velocity profile and these could contribute to syrinx development or enlargement. Few studies have considered the fluid mechanics of CM patients with and without syringomyelia separately. Three subject-specific CFD models were developed for a normal participant, a CM patient with syringomyelia and a CM patient without syringomyelia. Model geometries, CSF flow rate data and CSF velocity validation data were collected from MRI scans of the 3 subjects. The predicted peak CSF pressure was compared for the 3 models. An extension of the study performed geometry and flow substitution to investigate the relative effects of anatomy and CSF flow profile on resulting spinal CSF pressure. Based on 50 monitoring locations for each of the models, the CM models had significantly higher magnitude (p o0.01) peak CSF pressure compared with normal. When using the same CSF input flow waveform, changing the upper spinal geometry changed the magnitude of the CSF pressure gradient, and when using the same upper spinal geometry, changing the input flow waveform changed the timing of the peak pressure. This study may assist in understanding syringomyelia mechanisms and relative effects of CSF velocity profile and spinal geometry on CSF pressure. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Syringomyelia Chiari malformation Cerebrospinal fluid (CSF) Computational fluid dynamics (CFD) Magnetic resonance imaging (MRI)
1. Introduction Chiari malformation Type I (CM) is an abnormality characterised by protrusion of the cerebellar tonsils through the foramen magnum. Approximately two thirds (Pinna et al., 2000) of patients with CM also develop syringomyelia, where fluid filled cysts (syrinxes) develop in the spinal cord. Many patients report headaches, pain, and motor and sensory disturbances. Surgery is the only available treatment but syrinxes may persist or enlarge after surgery, and it is not clear why some patients with CM develop syringomyelia and others do not. Understanding the mechanisms of syringomyelia has the potential to improve treatment for these patients. Early theories of syrinx development in CM proposed that herniation of the cerebellar tonsils causes increased resistance across the foramen magnum thereby preventing spinal cord extracellular fluid outflow or driving fluid into the central canal (Heiss et al., 1999; Oldfield et al., 1994; Olivero and Dinh, 1992). However a syrinx may be under high pressure and it is not clear how the CSF
n
Corresponding author. Tel.: +61 299 264821; fax: +61 299 265266. E-mail address: [email protected] (E.C. Clarke).
0021-9290/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jbiomech.2013.05.013
would accumulate in the central canal against a pressure gradient. Another theory proposes that syringomyelia develops passively as a result of the Venturi effect, where narrowed regions of the CSF pathway result in higher velocity and lower pressure (Greitz, 2006), however clinical studies of CM patients have reported significant reductions in syrinx size following surgery without significant reduction in CSF velocity near the syrinx (Mauer et al., 2011), and the velocity of the fluid flow is low, so such an effect is likely to be small. Experimental studies have demonstrated that fluid flows into the cord and syrinx via the perivascular spaces (Rennels et al., 1990; Rennels et al., 1985; Stoodley et al., 1996). It is possible that the pressure gradient and/or timing of arterial pulsation-dependent CSF flow could be a driving force of syrinx fluid from the perivascular spaces into the central canal. A recent computational simulation demonstrated that even a moderate shift in the timing of the CSF pressure pulse relative to arterial pulse could influence the mass flow rate of CSF into the perivascular spaces (Bilston et al., 2010). More recently, a follow up study by Elliott et al. demonstrated that this mechanism requires that the spinal cord have non-zero compliance with the vascular bed (Elliott et al., 2011). While MRI scanning can provide structural anatomical, velocity and flow rate data, and has been used to demonstrate the contribution
1802
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
of arterial and venous pulsations to CSF pulsation (Bhadelia et al., 1997), it is not possible to measure non-invasively spinal CSF pressures in humans. Invasive CSF pressure measurements in patients with CM have demonstrated dissociation between cranial and spinal CSF pressures (Williams, 1981a; Williams, 1981b), and increased CSF pressure in CM patients that reduced following posterior fossa decompression surgery (Heiss et al., 1999). However these studies did not compare patients with and without syringomyelia and did not investigate any temporal changes in the CSF waveform (e.g. delays in peak CSF pressure or delays in onset of systolic flush). A very recent study by Heiss et al. has compared CSF pressures in CM patients with and without syringomyelia (Heiss et al., 2012). Subject-specific computational simulations provide an opportunity to non-invasively study CSF pressure in patients with CM, with and without syringomyelia. Previous computational simulations of the spinal subarachnoid space (SSAS) have investigated the mechanisms of syrinx development in conditions such as arachnoiditis, finding that a porous obstruction in the SSAS changes timing and magnitude of peak CSF pressure (Cheng et al., 2012). Another study used computational modelling to investigate the effect of tonsillar herniation on CSF pressure by extending normal tonsil geometry inferiorly, finding that modified tonsil geometry increased the CSF pressure gradient (Linge et al., 2011). To our knowledge this is the first study to investigate the cyclic CSF pressure profile using subject-specific models of Chiari patients with and without syringomyelia.
2. Methods 2.1. MRI scanning and flow measurements The University of NSW Human Research and Ethics Committee approved all experimental protocols, and all participants gave written informed consent to enter
the study. Three participants underwent MRI scanning of the head and neck; one patient with CM and syringomyelia (female, 58 years, 80 kg), one patient with CM but no syringomyelia (female, 60 years, 65 kg), and one healthy control subject (female, 60 years, 80 kg). A 3D sagittal isotropic (0.94 mm) T1-weighted turbo field echo MRI of the head and neck, and cardiac-gated cine phase-contrast (PC) MRI scans were performed using a 3 T MRI scanner (Philips Achieva TX, Best, Netherlands). The parameters for the 3D isotropic MRI scan were flip angle ¼ 81, matrix ¼288 288, FOV ¼270 270 169, TR/TE¼ 5.5/2.4 m s, slice thickness¼0.94 mm. For the PC-MRI scans, retrospective cardiac gating was achieved using VCG leads, and 30 cardiac phase images were collected, measured from the onset of the R wave. PC-MRI scans were performed at the following cranial/spinal levels and encoding velocities (Venc) using an (axial) scan plane aligned perpendicular to the CSF flow: 5 mm cranial to the tip of the cerebellar tonsils (Venc 10 cm/s); mid-vertebral C2 (Venc 9 cm/s); and midvertebral C5 (Venc 13 cm/s) (Fig. 1A). The spatial resolution of the PC-MRI scans was 0.694–0.977 mm2. Encoding velocities were determined during pilot studies in order to prevent aliasing. Other typical PC-MRI parameters were flip angle ¼101, matrix ¼240 176, FOV ¼ 250 250, TR/TE¼21/6.8 m s, slice thickness¼ 5 mm. CSF velocity and flow rate were measured from PC-MRI scans using the freely available analysis software Segment (Heiberg et al., 2010). For the tonsil-level scan the entire CSF cross-section was segmented and analysed for flow rate data (Fig. 1D and E), and for the C2- and C5-level scans velocity data were analysed using 10 regions of interest (2 mm diameter) spaced evenly around the SSAS (Fig. 1B and C). The tonsil-level flow rate was represented by a 6th order Fourier series for use as a flow input function for the models (Fig. 5B), and the C2 and C5 velocity data were used for model validation (Fig. 3).
2.2. Computational modelling The spinal cord, spinal canal and cerebellum were segmented manually from the 3D MRI scans using SURFdriver software. These point cloud segmentations were imported into Rhinoceros (v3.0, Robert McNeal and associates, WA) and fitted with closed non-rational B-splines of degree 3 passing through the points for smoothing. These were then fitted with a surface (loft) and capped to create volumes. The CSF volume was created by Boolean subtraction of the cerebellum and spinal cord volumes from the spinal canal volume, then trimmed 5 mm cranial to the base of the cerebellar tonsils and 95 mm caudal to the base of the tonsils (Fig. 2A). The cranial end of the CSF volume was extended upwards by 5 mm to
Fig. 1. Example MRI scans from a healthy control volunteer showing (A) the approximate alignment of the cine phase-contrast MRI scan planes (cerebellar tonsils, mid-C2 and mid-C5) aligned perpendicular to the CSF flow, (B) an example anatomical scan through the mid-C2, (C) a magnified view of a mid-C2 level phase image (as indicated in B) depicting approximate locations of the 10 regions used for validation of the models against MRI-measured CSF velocity data, (D) an example anatomical scan through the cerebellar tonsil level, and (E) a magnified view of a tonsil level phase image (as indicated in D) depicting manual segmentation of the entire CSF region.
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
Normal
With Syrinx
No Syrinx
1803
cardiac cycle time for an indication of the relative peak pressure timing between the models. Since there are known differences in both hindbrain anatomy and CSF flow profile for CM patients, we performed an additional substitution step (changing either the anatomy or flow profile while holding the other constant) in order to investigate the relative effects of the anatomy and CSF flow on resulting CSF pressure profile. For example, it is possible that the CM patient with syringomyelia developed the syrinx due to the hindbrain anatomy, so by normalising the input CSF flow rate function we aimed to explore the effects of hindbrain and upper spinal geometry on CSF pressure profiles between the CM patients with and without syringomyelia. The geometries from each of the CM participants were used with the Fourier series flow-rate function calculated from the PC-MRI scan (as described above) of the healthy participant (termed geometry substitution), and the healthy participant geometry was used with the Fourier series flow-rate functions calculated from the PC-MRI scans (as described above) for each of the CM participants (termed flow substitution). The pressure monitoring locations from the original geometries were used for analysis of these cases, and again a uniform velocity across the inlet was assumed.
3. Results The simulation velocities and MRI measured velocities were well correlated (examples shown in Fig. 3). There were no significant differences between MRI and model peak velocities for any of the 3 cases (p ¼0.1 with syringomyelia, p ¼0.7 without syringomyelia and p¼ 0.07 normal participant). There were no significant differences between MRI and model peak velocity timing for CM with syringomyelia, (p ¼0.1) or normal participant (p ¼0.5). Fig. 4 shows examples of the velocity (Fig. 4 A and B) and pressure gradients (Fig. 4C and D) at 2 time steps for a representative model (geometry substitution of no-syrinx patient with normal participant input flow rate profile). 3.1. CSF hydrodynamics in Chiari malformation
Fig. 2. (A) The CSF volumes for the 3 cases from the inlet side (cerebellar tonsil level) and (B) an example of the model mesh at the outlet showing the inflation layers.
provide an entry zone ahead of rapid changes in cross-sectional shape. The volumes were then imported into ANSYS CFX (v13, ANSYS Inc., PA). Tetrahedral meshing of the CSF volumes was performed using ANSYS with 5 layers of prismatic inflation at the walls (Fig. 2B). Proximity meshing was used to ensure that small flow paths were resolved. CSF was modelled as a Newtonian fluid with viscosity 0.8 m Pa s and density 1000 kg/m3 (Bloomfield et al., 1998). For a typical velocity of 0.1 m/s and hydraulic diameter of 5 mm the Reynolds number was calculated as 625, so laminar flow was assumed. The tonsil-level Fourier series flow rate function (using original timescale information from PC-MRI scans) was assigned as the input at the cranial end of the model, with a uniform velocity across the inlet (see Section 4). The average gauge pressure at the caudal end of the models was set to 0 Pa. Velocity monitoring points were manually placed in the models at coordinates approximately corresponding to the monitoring points from the MRI scans (Fig. 1). Pressure monitors were placed at 5 evenly spaced cranialcaudal levels along the model, with 10 points placed evenly around the SSAS. In all cases an initial steady state run was performed to establish the flow and provide a consistent starting flow-field for the transient run. A transient solution was performed and data were collected on the second cycle, with the solution converged to a normalised residual of o10−5 at each time step (5 m s) for all variables. The CFD code used is a 3D Navier–Stokes solver that uses a finite volume method applied to a vertex based discretisation. The equations were solved using a coupled solver, with the pressure–velocity coupling achieved by a modified Rhie– Chow scheme. The spatial derivatives were calculated using a bounded second order scheme and the 2nd order backwards Euler scheme was used for time derivatives. Double precision arithmetic was used to avoid any problems from rounding errors. Validation of the models was performed using the original timescales by comparing the simulated velocity profiles (magnitude and time of peak velocity) against the MRI measured velocity profiles (magnitude and time of peak velocity) at each of the 20 corresponding locations, for each participant, using paired t-tests. The time axes of output pressures were also normalised to proportion of the
Fig. 5 shows simulated CSF pressures, relative to the outlet at the caudal end of the models, in the anterior SSAS at the C2 level for the 3 models, both on the normalised timescale Fig. 5(C) and the original timescale Fig. 5(D) and the corresponding flow rate data Fig. 5(A–B). All 3 cases showed an initial peak in CSF pressure during systole. In this location, the peak pressure was 26.1 Pa in the model for the healthy participant, 33.6 Pa for the patient with syringomyelia, and 56.6 Pa in the CM patient without syringomyelia. The time of peak pressure was nearly 40%, or approximately 50 m s, earlier than normal in the patient without syringomyelia. These patterns were consistent across other pressure monitoring locations. 3.2. Flow substitution Fig. 6(A–B) shows simulated CSF pressures at one corresponding location (anterior, C2-level) for the normal model geometry with input flow functions from the 3 participants, both on the normalised timescale (A) and the original timescale (B), and the corresponding flow rate data are shown in Fig. 5(A–B). The times of peak pressures for the flow substitution models aligned with the times of peak pressures for the corresponding participant original models (i.e. in Fig. 5A with Fig. 6A and Fig. 5B with Fig. 6B). As with the original models shown in Fig. 5C the peak pressure arrived earlier in the cardiac cycle for the flow substitution model of the CM patient without syringomyelia (shown in Fig. 6A). 3.3. Geometry substitution Fig. 6(C–D) shows simulated CSF pressures at one corresponding location (anterior, C2-level) for the 3 model geometries with input flow functions from the healthy participant, both on the normalised timescale (C) and the original timescale (D). Corresponding flow rate data for the normal participant is shown in
1804
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
Model MRI data
CSF velocity (m/s)
CSF velocity (m/s)
Model MRI data
Time (% cardiac cycle)
Time (% cardiac cycle)
CSF velocity (m/s)
Model MRI data
Time (% cardiac cycle) Fig. 3. MRI-measured and simulation velocity–time data used for model validation. Each figure shows data from the C2 level (left panels) and C5 level (right panels) of the 5 regions around the subarachnoid space (left or right only shown) for the Chiari malformation patient with syringomyelia (A), the Chiari malformation patient without syringomyelia (B), and the normal participant (C). Data from the model simulations are shown as black solid lines and the data from the MRI scans are shown as grey dashed lines. The peak velocities used for validation statistics are indicated by the triangle symbols (simulation data are shown in black and MRI scan data are shown in grey).
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
1805
Fig. 4. Representative model outputs showing (A–B) CSF velocity and (C–D) CSF pressures at 2 time steps from one simulation.
Fig. 5(A–B). As with the original models shown in Fig. 5(C–D) the peak pressure was higher for the model geometry of the CM patient without syringomyelia than the model geometry of the CM patient with syringomyelia, and both CM geometries produced higher pressure than healthy geometry.
4. Discussion 4.1. Implications of results for syringomyelia development in CM The results of this study show that the model of the CM patient with syringomyelia had slightly elevated peak pressure compared with normal, but closer to the timing of peak pressure in the normal participant model. Interestingly, the CM patient without syringomyelia had earlier and higher magnitude peak pressure than both the normal participant and the patient with syringomyelia. The current study included 3 subjects, and while flow profiles were typical of their respective groups in a larger MRI study (Clarke et al., 2013), individual patient anatomy may not be, so the pressure results from these subjects may not be representative of their groups. Previous simulations (Bilston et al., 2010) have suggested
that a modest shift in the timing of the CSF pressure profile with respect to the arterial pressure profile may enhance CSF flow through the perivascular spaces into the central canal of the spinal cord, and this may contribute to formation or enlargement of a syrinx. More specifically, if the peak SSAS pressure occurs at a time in the arterial cycle when spinal arteriole pressure is low (periarterial spaces at their widest), flow into the spinal cord may be enhanced. This would encourage fluid accumulation in the spinal cord, unless fluid outflow is able to accommodate this additional influx. Similarly, if the lowest SSAS pressures occur when spinal arteriole pressure is high (peri-arterial spaces at their narrowest), this may impede backflow through the peri-arterial spaces encouraging fluid accumulation in the cord. This could lead to the “presyrinx state” observed clinically, and eventually to a syrinx. The results from the flow and geometry substitution studies suggest that the pressure magnitude is strongly influenced by the anatomy at the craniocervical junction and in the spinal SSAS. It is difficult to compare our predicted pressure magnitudes (peak 20–60 Pa) with in vivo pressure measurements due to the location/nature of reference pressure points, penetration of the dura for pressure sensors, the rigid in vitro nature of our simulations with absence of structures such as spinal nerve. An in vivo
1806
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
Normal participant Chiari with syrinx Chiari without syrinx
Normal participant Chiari with syrinx Chiari without syrinx
Normal participant Chiari with syrinx Chiari without syrinx
Normal participant Chiari with syrinx Chiari without syrinx
Fig. 5. (A) and (B) shows CSF flow rate data from the tonsil-level PC-MRI scans (segmented as shown in Fig. 1E) alongside the Fourier functions that were fit to these data points and used as input flow rates for the simulations. MRI data are shown as symbols for the normal (triangle), Chiari with syrinx (circle), and Chiari without syrinx (square) participants, and Fourier fits are shown as lines for the normal (sold black), Chiari with syrinx (solid grey), and Chiari without syrinx (dashed black) participants. Both the normalised timescale (A) and original timescale (B) are shown. (C) and (D) shows simulated CSF pressures, relative to the outlet at the caudal end of the models, at one representative location (anterior, C2-level) for the normal healthy participant (solid black line), the Chiari malformation participant with syringomyelia (solid grey line) and the Chiari malformation patient without syringomyelia (dashed black line). Both the normalised timescale (C) and original timescale (D) are shown.
study (Williams, 1981b) measured the mean amplitude of CSF under cardiac pulsation as 3.86 mm Hg (514 Pa) at the lumbar relative to ventricular level, and another (Heiss et al., 1999) measured mean cervical subarachnoid pressure as 2 mm Hg (267 Pa) at the cervical relative to the ventricular level. Our study considers only a small length of the spinal canal (pressures relative the lower cervical level). The monitoring and reference pressure points are in different locations in our study and these in vivo examples but it should be noted that our predicted pressures are significantly lower than those reported in the in vivo studies and care should be taken if using the absolute magnitudes of pressures reported in our study. The current study is however useful for making comparisons between the 3 patient groups. The findings of the current study agree with previous computational fluid dynamics studies. One study demonstrated that simulated tonsillar herniation increased spinal CSF pressure gradient (Linge et al., 2011). Other studies have used lumped parameter and coupled hydrodynamics to investigate CSF pressure and flow environments (Kim and Cirovic, 2011). A study by Elliott et al. (2011) investigated mechanisms of syringomyelia and effects of simulated shunt and bypass surgery. Interestingly they found that, while the simulation of arachnoiditis increased the steady state pressure mean and amplitude, including a syrinx lowered the mean and amplitude of the steady state pressure back towards normal. Similar stenosis modelling studies by Martin et al. found that the inclusion of a syrinx reduced the mean and amplitude of SSAS pressure (Martin et al., 2010; Martin and Loth, 2009). This observation is consistent with our findings that CSF peak pressure is closer to normal in patients with syringomyelia than without.
Several previous studies have combined MRI and modelling to investigate CSF dynamics for clinical purposes. For example, Hsu et al. (2012) developed an intricate computational model to investigate CSF pulsation effects on intrathecal drug distribution; Alperin et al. (1996) investigated perturbations in arterial and venous flow on CSF and spinal tissues, and Sweetman and Linninger (2011) developed a complex 3D fluid–structure interaction model of the entire CNS to investigate interactions between pulsating vasculature, compliant tissues, and CSF flow. To our knowledge the models presented in the current study are the first subject-specific CFD models for investigating CSF pressure in CM patients with and without syringomyelia. The findings of this study are important because the only treatment for symptomatic CM is posteria fossa decompression surgery and this is sometimes unsuccessful. The surgery certainly changes geometry and therefore likely changes the CSF pressure gradient. However, it is still not known whether this change in geometry also changes the flow-rate profile, which may also affect the CSF pressure profile. This may help to explain why decompression surgery is sometime unsuccessful and further investigation may improve treatment for these patients. 4.2. Limitations of the current study The models presented here are CFD models with solid boundaries. We hypothesised that deformation of the spinal cord in response to these pressures is likely to be very small, because the magnitudes of the pressures involved (o 100 Pa) are orders of magnitude lower than the elastic modulus of the spinal cord
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
1807
Normal geometry, normal input flow Normal geometry, Chiari with syrinx flow Normal geometry, Chiari without syrinx flow
Normal geometry, normal input flow Normal geometry, Chiari with syrinx flow Normal geometry, Chiari without syrinx flow
Normal geometry, normal input flow Chiari with syrinx geometry, normal flow Chiari without syrinx geometry, normal flow
Normal geometry, normal input flow Chiari with syrinx geometry, normal flow Chiari without syrinx geometry, normal flow
Fig. 6. (A) shows simulated CSF pressures on a normalised timescale at one representative location (anterior, C2-level) for the flow substitution models, where the same model geometry (normal participant) was used but the input flow rate function was different (the 3 flow rate functions are shown in Fig. 5B). The solid black line is derived from the normal participant geometry with normal participant input flow rate function, the solid grey line is derived from the normal geometry with flow rate function from the patient with syringomyelia, and the dashed black line is derived from the normal geometry with flow rate function from the patient without syringomyelia. Arrows indicate the change in timing of the peak pressures for the 3 models. (B) Shows the same data as in (A), but using the original timescale, rather than the normalised timescale. (C) Shows simulated CSF pressures on a normalised timescale at one representative location (anterior, C2-level) for the geometry substitution models, where the same input flow rate function (normal participant, shown in Fig. 5B) was used but model geometry was different (the 3 geometries are shown in Fig. 2 A). The solid black line is derived from the normal participant geometry with normal participant input flow rate function, the solid grey line is derived from the geometry of the patient with syringomyelia and normal participant input flow rate function, and the dashed black line is derived from the geometry of the patient without syringomyelia and normal participant input flow rate function. Arrows indicate the change in peak pressure magnitude for the 3 models. (D) Shows the same data as in (C), but using the original timescale, rather than the normalised timescale.
(approx. 1 MPa) (Bilston and Thibault, 1996). We confirmed this using fluid–structure interaction modelling for the normal participant and found no discernable difference in the pressure wave timing or magnitude, and micron-order displacement of the spinal cord. We also measured tonsil motion from the tonsil-level PC-MRI (by segmenting the tonsils and calculating their displacement from the velocity–time curves) in these subjects and did not detect any cranial-caudal tonsil motion throughout the cardiac cycle. We therefore believe these assumptions will have minimal impact on the estimated pressures. While the locations of the C2 and C5 PC-MRI scan planes were known with respect to the model geometries, the 10 velocity probes on each plane (used for validation) were placed manually. Despite this the agreement between MRI and simulations velocities is good (Fig. 3). Also probe volumes for velocity validation different sizes (point probes were used for simulations whereas 2 mm diameter probes were used for MRI scans to include some spatial averaging over multiple pixels). This difference in size is unlikely to affect the overall patterns of results in this study, but closer matching of probe volumes may have improved the validation at narrow regions of the SSAS. Overall, however the validation results match well (Fig. 3) and there were no significant differences in peak systolic velocity magnitude and no significant differences in timing of peak systolic velocity for 2 of the 3 cases.
Furthermore, there are different types of error in CSF velocity– time profiles from PC-MRI scanning depending on the type of gating used. Retrospective cardiac gating was used in order to ensure that the entire cycle was captured and to minimise sensitivity to irregular heart rates (Lenz et al., 1989). However, we did not simultaneously measure arterial and CSF waveforms but rather measured the CSF waveforms relative to the R–R interval. It is important to keep this in mind when comparing the cycle timing for the 3 models presented here. Previous studies using cine PC-MRI investigated relative timing of pulses in the epidural vessels and CSF simultaneously, finding no significant differences in the arterial and CSF peak flow timing at the C2-C3 levels (Baledent et al., 2001; Wagshul et al., 2006). The inlet flow-rate condition to these models assumed a uniform velocity across the inlet because the inlet-level MRI low-rate data was sampled by segmenting the entire SSAS (i.e. the mean flow-rate from all SSAS pixels). In reality there is some variation in flow-rate across the axial plane, but data from our larger cohort PC-MRI flow study (Clarke et al., 2013) found no significant differences in peak velocities between probes placed around the SSAS, so we believe the assumption of uniform inlet velocity would have negligible effect on the outcomes of the models. The 3 participants were chosen because they had flow-rate and velocity profiles that were representative of their respective
1808
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
3 groups in a larger cohort quantitative flow study (Clarke et al., 2013). From that study we found patients without syringomyelia had significantly earlier and faster peak systolic CSF flow than both normal volunteers and patients with syringomyelia. It follows that in the current study, the peak CSF pressure occurred earlier in the model for the patient without syringomyelia compared with the other 2 models. However a recently published similar study reported that Chiari patients had significantly later peak systolic CSF flow than normal (Bunck et al., 2012). Therefore it is important to recognise that the models presented here may not be representative of all patients in their respective groups, however the outcomes of the models may be useful for making inferences about the relative effects of upper spinal geometry and CSF flow rate on the timing and magnitude of peak CSF pressure. Also, the geometry and flow substitution models are not true subjectspecific models with matched flow and geometry, and therefore cannot be validated experimentally and are presented only as an exploration of the role of the anatomy and flow profiles on resulting SSAS pressures. It should also be noted that MRI segmentations were smoothed to assist formation of smooth model geometries, the models do not include epidural blood vessels or nerve roots, and do not simulate the pulsatile blood flow in the epidural vessels. Surface smoothing and exclusion of small structures leads to a better quality computational mesh, which was deemed preferable over working with raw data that necessarily contains imperfections in geometry. While these features may affect localised pressure magnitudes we would not expect these features to change the overall patterns observed in these models and we would expect the exclusion of these features to have similar effects across the 3 models. Also, spine and spinal cord compliance were assumed equal for the 3 participants, but it is possible that these may be different (this is not known) and this may affect the magnitude and timing of the CSF pressure profiles. For example, Martin and colleagues have demonstrated that changing spinal canal compliance and including the presence of a modelled syrinx influences the magnitude of CSF pressure and timing (Martin et al., 2010; Martin and Loth, 2009). The models presented here are simplified in terms of tissue compliance and inclusion of spinal structures, but do include the subject-specific spinal geometries and inlet flow conditions. The situation appears more complex than either study considered alone but it was interesting to note that Martin et al. found the presence of a syrinx to increase the magnitude and delay the peak CSF pressure compared with an equivalent model without a syrinx (Martin et al., 2010). In the current study we also found the model for the patient with a syrinx to have later peak CSF pressure than the patient without a syrinx but reduced peak CSF pressure. It is acknowledged that in the current study using a short section of upper spinal geometry that timing of the pressure waves may be affected mildly by fluid–structure interaction, but the complete study of CSF dynamics in the entire (closed) CNS space needs to account for the interaction of CSF filled spaces with deformable structures such as soft brain parenchyma boundary, the dura and the blood vessels.
Conflicts of Interest No conflicts of interest to declare.
Acknowledgements A research project grant from the Column of Hope provided funding for the study. ECC and LEB are supported by NHMRC research fellowships.
References Alperin, N., Vikingstad, E.M., Gomez-Anson, B., Levin, D.N., 1996. Hemodynamically independent analysis of cerebrospinal fluid and brain motion observed with dynamic phase contrast MRI. Magnetic Resonance in Medicine 35, 741–754. Baledent, O., Henry-Feugeas, M.C., Idy-Peretti, I., 2001. Cerebrospinal fluid dynamics and relation with blood flow: a magnetic resonance study with semiautomated cerebrospinal fluid segmentation. Investigative Radiology 36, 368–377. 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, 258–264. Bilston, L.E., Stoodley, M.A., Fletcher, D.F., 2010. The influence of the relative timing of arterial and subarachnoid space pulse waves on spinal perivascular cerebrospinal fluid flow as a possible factor in syrinx development. Journal of Neurosurgery 112, 808–813. Bilston, L.E., Thibault, L.E., 1996. The mechanical properties of the human cervical spinal cord in vitro. Annals of Biomedical Engineering 24, 67–74. Bloomfield, I.G., Johnston, I.H., Bilston, L.E., 1998. Effects of proteins, blood cells and glucose on the viscosity of cerebrospinal fluid. Pediatric Neurosurgery 28, 246–251. Bunck, A.C., Kroeger, J.R., Juettner, A., Brentrup, A., Fiedler, B., Crelier, G.R., Martin, B. A., Heindel, W., Maintz, D., Schwindt, W., Niederstadt, T., 2012. Magnetic resonance 4D flow analysis of cerebrospinal fluid dynamics in Chiari I malformation with and without syringomyelia. European Radiology. Cheng, S., Stoodley, M.A., Wong, J., Hemley, S., Fletcher, D.F., Bilston, L.E., 2012. The presence of arachnoiditis affects the characteristics of CSF flow in the spinal subarachnoid space: A modelling study. Journal of Biomechanics 45, 1186–1191. Clarke, E., Stoodley, M., Bilston, L., 2013. Changes in temporal flow characteristics of CSF in Chiari I Malformation with and without syringomyelia: Implications for theory of syrinx development. Journal of Neurosurgery, 118(5):1135-1140 http://dx.doi.org/10.3171/2013.2.JNS12759. Epub 2013 Mar 15. Elliott, N.S., Lockerby, D.A., Brodbelt, A.R., 2011. A lumped-parameter model of the cerebrospinal system for investigating arterial-driven flow in posttraumatic syringomyelia. Medical Engineering and Physics 33, 874–882. Greitz, D., 2006. Unraveling the riddle of syringomyelia. Neurosurgical Review 29, 251–263, discussion 264. Heiberg, E., Sjogren, J., Ugander, M., Carlsson, M., Engblom, H., Arheden, H., 2010. Design and validation of segment—freely available software for cardiovascular image analysis. BMC Medical Imaging 10, 1. Heiss, J.D., Patronas, N., DeVroom, H.L., Shawker, T., Ennis, R., Kammerer, W., Eidsath, A., Talbot, T., Morris, J., Eskioglu, E., Oldfield, E.H., 1999. Elucidating the pathophysiology of syringomyelia. Journal of Neurosurgery 91, 553–562. Heiss, J.D., Snyder, K., Peterson, M.M., Patronas, N.J., Butman, J.A., Smith, R.K., Devroom, H.L., Sansur, C.A., Eskioglu, E., Kammerer, W.A., Oldfield, E.H., 2012. Pathophysiology of primary spinal syringomyelia. Journal of Neurosurgery: Spine 17, 367–380. Hsu, Y., Hettiarachchi, H.D., Zhu, D.C., Linninger, A.A., 2012. The frequency and magnitude of cerebrospinal fluid pulsations influence intrathecal drug distribution: key factors for interpatient variability. Anesthesia and Analgesia 115, 386–394. Kim, M., Cirovic, S., 2011. A computational model of the cerebrospinal fluid system incorporating lumped-parameter cranial compartment and one-dimensional distributed spinal compartment. Journal of Biorheology 25, 78–87. Lenz, G.W., Haacke, E.M., White, R.D., 1989. Retrospective cardiac gating: a review of technical aspects and future directions. Magnetic Resonance Imaging 7, 445–455. Linge, S.O., Haughton, V., Lovgren, A.E., Mardal, K.A., Helgeland, A., Langtangen, H.P., 2011. Effect of tonsillar herniation on cyclic CSF flow studied with computational flow analysis. AJNR American Journal of Neuroradiology 32, 1474–1481. Martin, B.A., Labuda, R., Royston, T.J., Oshinski, J.N., Iskandar, B., Loth, F., 2010. Spinal subarachnoid space pressure measurements in an in vitro spinal stenosis model: implications on syringomyelia theories. Journal of Biomechanical Engineering 132, 111007. Martin, B.A., Loth, F., 2009. The influence of coughing on cerebrospinal fluid pressure in an in vitro syringomyelia model with spinal subarachnoid space stenosis. Cerebrospinal Fluid Research 6, 17. Mauer, U.M., Gottschalk, A., Mueller, C., Weselek, L., Kunz, U., Schulz, C., 2011. Standard and cardiac-gated phase-contrast magnetic resonance imaging in the clinical course of patients with Chiari malformation Type I. Neurosurgical Focus 31, E5. Oldfield, E.H., Muraszko, K., Shawker, T.H., Patronas, N.J., 1994. Pathophysiology of syringomyelia associated with Chiari I malformation of the cerebellar tonsils. Implications for diagnosis and treatment. Journal of Neurosurgery 80, 3–15. Olivero, W., Dinh, D., 1992. Chiari I malformation with traumatic syringomyelia and spontaneous resolution: case report and literature review. Neurosurgery 30, 758–760. Pinna, G., Alessandrini, F., Alfieri, A., Rossi, M., Bricolo, A., 2000. Cerebrospinal fluid flow dynamics study in Chiari I malformation: implications for syrinx formation. Neurosurgical Focus 8, E3. Rennels, M.L., Blaumanis, O.R., Grady, P.A., 1990. Rapid solute transport throughout the brain via paravascular fluid pathways. Advances in Neurology 52 431–439
E.C. Clarke et al. / Journal of Biomechanics 46 (2013) 1801–1809
Rennels, M.L., Gregory, T.F., Blaumanis, O.R., Fujimoto, K., Grady, P.A., 1985. Evidence for a ‘paravascular’ fluid circulation in the mammalian central nervous system, provided by the rapid distribution of tracer protein throughout the brain from the subarachnoid space. Brain Research 326, 47–63. Stoodley, M.A., Jones, N.R., Brown, C.J., 1996. Evidence for rapid fluid flow from the subarachnoid space into the spinal cord central canal in the rat. Brain Research 707, 155–164. Sweetman, B., Linninger, A.A., 2011. Cerebrospinal fluid flow dynamics in the central nervous system. Annals of Biomedical Engineering 39, 484–496. 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. Journal of Neurosurgery 104, 810–819.
1809
Williams, B., 1981a. Simultaneous cerebral and spinal-fluid pressure recordings 1. Technique, physiology, and normal results. Acta Neurochirurgica 58, 167–185. Williams, B., 1981b. Simultaneous cerebral and spinal-fluid pressure recordings 2. Cerebrospinal dissociation with lesions at the foramen magnum. Acta Neurochirurgica 59, 123–142.