Patient-Specific 3D Simulation of Cyclic CSF Flow at the Craniocervical Region

BACKGROUND AND PURPOSE: Flow simulations in patient-specific models of the subarachnoid space characterize CSF flow in more detail than MR flow imaging. We extended previous simulation studies by including cyclic CSF flow and patient-specific models in multiple patients with Chiari I. We compared simulation results with MR flow measurements. MATERIALS AND METHODS: Volumetric high resolution image sets acquired in 7 patients with Chiari I, 3 patients who had previous craniovertebral decompression, and 3 controls were segmented and converted to mathematical models of the subarachnoid space. CSF flow velocities and pressures were calculated with high spatial and temporal resolution during simulated oscillatory flow in each model with the Navier-Stokes equations. Pressures, velocities, and bidirectional flow were compared in the groups (with Student t test). Peak velocities in the simulations were compared with peak velocities measured in vivo with PCMR. RESULTS: Flow visualization for patients and volunteers demonstrated nonuniform reversing patterns resembling those observed with PCMR. Velocities in the 13 subjects were greater between C2 and C5 than in the foramen magnum. Chiari patients had significantly greater peak systolic and diastolic velocities, synchronous bidirectional flow, and pressure gradients than controls. Peak velocities measured in PCMR correlated significantly (P = .003; regression analysis) despite differences between them. CONCLUSIONS: In simulations of CSF, patients with Chiari I had significantly greater peak systolic and diastolic velocities, synchronous bidirectional flow, and pressure gradients than controls.

C omputer simulations created with CFD demonstrate properties of CSF flow that may not be visualized with MR imaging. 1,2 The simulations use 3D mathematical models of subarachnoid space geometries. CSF flow velocities and pressures are then computed on discretized representations (meshes) of these models by solving Navier-Stokes equations. CSF flow-velocity vectors and pressure scalars are computed at each designated computational point and time-step. CSF flow is displayed in 3D with high spatial and temporal resolution for the selected time period. Spatial and temporal resolutions may, in principle, be as fine as desired, the upper limit being dictated by computer resources. Compared with MR imaging, all of the selected subarachnoid space is covered at every chosen point in time rather than just selected sections or times.
Patient-specific computer simulations allow detailed 3D flow and pressure analysis to be undertaken in actual subarachnoid spaces. In a previous study, 3 models of the subarachnoid spaces in 1 healthy adult and in 1 patient with Chiari I were created from segmented MR images. Displayed flow patterns matched those observed with PCMR in the same patients. The study showed variations in pressure and velocity through the upper cervical spinal canal and greater velocities and pressure gradients in the patient with Chiari I than in the control subject. This pilot study, however, showed flow simulations only for constant caudal or cranial flow, not for oscillating flow. Comparison of computed velocities with MR measurements was performed in 6 patients.
The goal of the present study was to create a group of patient-specific models for simulation of cyclic craniocervical CSF flow.

Patients and Subjects
The institutional review board at the University of Wisconsin School of Medicine and Public Health approved this retrospective study of anonymized anatomic and phase-contrast MR studies in a registry of image data at the university. Cases were selected without reference to clinical or demographic data.

MR Imaging Data Collection
Imaging was acquired on 1 of 3 1.5T scanners similarly equipped. Imaging included localizer images, sagittal T1-weighted images (TR/TE ϭ 520/10, section thickness ϭ 3 mm, matrix ϭ 256 ϫ 256, FOV ϭ 16 cm, NEX ϭ 3), and sagittal T2-weighted images (TR/TE ϭ 3600/116, section thickness ϭ 3 mm, matrix ϭ 512 ϫ 512, FOV ϭ 16 cm) of the cervical spine, and other head and axial spine images as clinically indicated. Sagittal phase-contrast MR images and axial images at the foramen magnum and at other levels in the upper cervical spinal canal were acquired. The images were gated to the cardiac cycle by electrocardiogram. The plane of acquisition for the PCMR flow images was selected according to a standard protocol. For sagittal flow images, the midline sagittal plane was chosen. For the axial flow images, a plane transverse to the axis of the spinal canal immediately below the tonsilar tips and other axial planes below were selected. For the phase-contrast flow measurement with the axial PCMR, the acquisition parameters were flip angle (␣) ϭ 20°, TR/TE ϭ 20/5 ms, section thickness ϭ 5 mm, FOV ϭ 180 mm, matrix ϭ 256 ϫ 256, and velocity encoding ϭ 10 cm/s. Flow was calculated from axial PCMR images with software resident in the scanner. 4 T2-weighted images were inspected for evidence of syrinx and tonsilar herniation and other abnormal findings. Tonsilar herniation was measured with resident software as the distance from the tip of the tonsils to a line defining the foramen magnum. 4 A syrinx was diagnosed on the basis of a well-defined region with the signal intensity of clear fluid within the spinal cord.
In each patient, a 3D high-resolution image set was acquired with a volumetric multiecho 3D radial acquisition with fat and water separation balanced steady-state free precession. 5 The sequence imaged a 20 ϫ 20 ϫ 20 cm 3 volume, with 256 ϫ 256 ϫ 256 voxels, acquired by 10,000 projections in a total scan time of 5 minutes, resulting in a spatial resolution of 0.8 ϫ 0.8 ϫ 0.8 mm 3 and high contrast between CSF and the surrounding tissues. Other scan parameters included flip angle ϭ 15°, receiver bandwidth ϭ Ϯ125 kHz, and TR ϭ 2.9 ms.

Patient-Specific Models
Volumetric images were processed and segmented, and meshes were created using VMTK (http//www.vmtk.org). 6 During this process the surface was smoothed with a passband of 0.1. The investigator performed some editing to facilitate the subsequent simulations, including adjustment of the resolution of the model to achieve the desired anatomic detail 7 required and to eliminate undesired structures such as blood vessels in the subarachnoid space. 8 The models were cut rostral to the tip of the occipital bone above foramen magnum and at the C5-C6 intervertebral disk level (Fig 1).

CSF Flow Simulations
Flow simulations were performed on the VMTK models by solving the Navier-Stokes equations for a Newtonian fluid over 5 flow cycles in FEniCS (http://www.fenicsproject.org) using the Finite Element Method in combination with the Incremental Pressure Correction Scheme. [8][9][10] A temporal resolution of 0.005 seconds was chosen. Flow was cycled at 1 Hz with a sine wave pattern.
Boundary conditions at the inlet and outlet were specified by equal top-bottom time-dependent pressure difference in all models. Local differences in pressure gradients resulting from geometric differences between the models were disregarded. No-slip boundary conditions were set at the dura and at the cord boundaries. Other conditions assumed for the simulation included initial stationary fluid, kinematic fluid viscosity of 0.7 ϫ 10 6 m 2 /s, fluid attenuation of 1 g/cc, and rigid walls. Flow in the caudal direction was assigned a negative sign.
The computed CSF flow velocities and pressures were visualized and analyzed in Paraview (http://www.paraview.org). Sagittal and axial sections were reviewed in cine mode. For comparison of different subjects and groups, snapshot images with a sagittal section and with axial sections at each cervical level were displayed for selected time points in the cardiac cycle (Fig 2). CSF pressures or velocities were displayed on the selected sections for each patient. Peak velocities for craniad and caudad flows were calculated with the "rescale to data-range" function and were tabulated for each level and each subject. One investigator assessed size and duration of synchronous bidirectional flow by parsing the images through the cardiac cycle. Magnitude of synchronous bidirectional flow was calculated as the maximal difference between flow velocities in the superior and inferior directions within 1 axial level; duration was quantified as the time span for which the magnitude of both caudad and craniad flow exceeded 0.1 cm/s.

Statistics
Difference in means between groups for pressure, velocity, and bidirectional flow were tested with the 2-tailed Student t test. The peak velocities in the simulations were compared with peak velocities in PCMR images at the same level, where available. Significance was calculated by the regression analysis program in Excel (Microsoft, Botthell, Washington).

Patients and Subjects
From the registry of 34 cases, 13 subjects who had highresolution volumetric images, sagittal T2-weighted images, PCMR images, and limited demographic data were selected (Tables 1-3).

MR Imaging Data Collection
On the basis of the tonsil position and clinical information tabulated in the registry, the subjects were classified as patients with Chiari I (n ϭ 7), controls (n ϭ 3), or post-craniovertebral decompression for a Chiari I (n ϭ 3). The patients with Chiari I included 3 females and 4 males, ranging in age from 2-40 years. The patients with Chiari I had 5-to 13-mm tonsilar herniation. Three controls subjects included 2 presumedhealthy volunteers and 1 female patient with abdominal symptoms who had an incidental giant cisterna magna, which prompted a CSF flow study to exclude an arachnoid cyst. One volunteer was female and the other 2 did not have sex tabulated. The postoperative patients ranged from 3-54 years of age. Two of these patients had syrinxes. Axial velocity measurements were obtained at 5 levels in 1 subject and at 1 level in 5 other patients. In 1 patient with Chiari I, a small lipoma was present in the foramen magnum, assumed to be incidental.

Patient-Specific Models
The meshes created for each subject contained 90,000 to 600,000 computational points. Resolution varied from 0.1-1.0 mm, depending on the detail required in the region. The lipoma in the subarachnoid space in 1 patient was removed in the editing process. Blood vessels and spinal nerve roots in the subarachnoid space were eliminated in editing. The segmentation processes produced technically satisfactory results, except in the regions with narrow subarachnoid spaces, where manual methods were needed. Syrinxes in 2 cases were edited to replace the fluid signal intensity within the cord with signal intensity similar to cord tissue, to facilitate segmentation. Minimal smoothing was required to achieve the final model that was used in FEniCS.

CSF Flow Simulations
Simulations obtained in each case showed a similar oscillating flow and pressure pattern in the upper cervical spinal canal. The time resolved images showed reversal of CSF flow during the cardiac cycle in all subjects. Velocities peaked at the middle (t ϭ 0.5) and at the end of the cycle (t ϭ 1.0) when the pressure gradient nulled and flow reversed when the pressure gradients reached maximum (t ϭ 0.25 and 0.75). In Paraview cine and snapshot sagittal and images, velocities had a nonuniform distribution through the subarachnoid space in all subjects in the study.

Velocity Characteristics for the Controls and Patients
Paraview snapshot images of peak velocities in each patient with Chiari I (Fig 3), control (Fig 4), and postoperative patient (Fig 5) showed greater velocities in the C5 end of the model than in the foramen magnum. Velocities peaked at C5 (in 5 patients with Chiari I) or between C2 and C4 (2 patients with Chiari I). At the level of C1, velocities anterior to the spinal cord were higher than posterior to the cord. Systolic velocities, on average, had slightly higher magnitudes, with peak veloci-    In healthy models, peak systolic velocities ranged from 6.8 cm/s to 7.9 cm/s ( Table 4). Peak systolic velocities increased below the foramen magnum in the 3 volunteers, reaching a maximum at C3, C4, or C5. Peak diastolic veloci-ties, slightly slower than systolic, ranged from 6.5 cm/s to 7.5 cm/s. Patients with Chiari I had bands of higher velocity near the tonsils. In patients with Chiari I, peak systolic velocities ranged from 7.5 cm/s to 12.1 cm/s. The 27% difference with respect to controls was significant (P ϭ .008, 2-tailed Student t test). Diastolic velocities ranged in the patients with Chiari I from 6.9 to 11.8 cm/s, significantly higher than the control group (P ϭ .027, 2-tailed Student t test). Temporal velocity variations were much more evident in patients with Chiari I than in volunteers. Five of the patients and 2 of the postoperative patients had large flow jets anterolateral to the spinal cord, stretching from C1 to C5 (Fig 6). The magnitude and intensity of the jets changed with time. These were most prominent at times of flow-direction reversal and less prominent at times of peak systolic flow. They appeared to have less uniformly increasing velocities along the spinal canal.
Postoperative patients had peak systolic velocities from 8.3 to 15.6 cm/s and peak diastolic velocities from 7.8 to 15.7 cm/s. Postoperative patients had no significant differences in velocities compared with patients with Chiari I (P ϭ .8).
Synchronous bidirectional flow appeared in all subjects at times of flow-direction reversal. It was most prominent along the lateral boundaries of subarachnoid space and typically reached maximum at the same levels at which maximum systolic/diastolic velocities were found ( Table 4). The magnitude of bidirectional flow ranged from 3.8 to 5.3 cm/s in controls and from 4.5 to 7.5 cm/s in patients with Chiari I. The 34% difference between patients with Chiari I and volunteers was significant (P ϭ .03, 2-tailed Student t test). The duration of synchronous bidirectional flow did not differ significantly between groups.

Pressure Characteristics
Pressures varied through the cardiac cycle in the 3 groups (Fig 7). Pressure gradients varied nonlinearly with level. At C1, differences between volunteers and patients with Chiari I were significant (P Ͻ .01, 2-tailed Student t test). No significant difference was found between volunteers or Chiari I and postoperative subjects.

Comparison of Computed Velocities to In Vivo Velocity Measurements
In PCMR images, velocities were higher in the C4 and C5 region than at C1, as in simulations. CSF velocities were greater in the patients with Chiari than in the controls, as in the simulations. Synchronous bidirectional flow was visualized in the patients with Chiari I. It was not visualized in the control subjects, who had small-magnitude bidirectional flow. Peak velocities measured in PCMR correlated significantly (P ϭ .003; regression analysis), despite differences between them (Fig 8A and B).

Discussion
Computations showed that CSF velocities peak between C2 and C5 in patients with Chiari I and controls. CSF velocities are greater in the patients with Chiari I than in healthy subjects. Pressure differences at the level of the tonsils in patients with Chiari I are greater than what is observed at the same level in controls. Synchronous bidirectional flow appears in all subjects at times of flow-direction reversal. The magnitude of synchronous bidirectional flow is greater in patients with Chiari I than in controls. It shows agreement between CSF velocities computed with CFD methods and CSF velocities measured with PCMR.
The PCMR and simulation results agree with previous studies. 4,[11][12][13][14] The greater spatial and temporal variations observed in Chiari I models have been reported in clinical 12 and CFD studies. [1][2][3] Jets observed in our study in the anterolateral locations in patients with Chiari I correspond to published observations. 3,4,12 Patients with Chiari I and, to a lesser extent, healthy subjects had synchronous bidirectional flow in agreement with previous reports. 1,12 The magnitude of systolic and diastolic flow in patients with Chiari I and in controls agrees well with previous reports. 15,16 The accuracy of CFD studies depends on the appropriateness of the conditions assumed. The suitability of the modeling techniques used to simulate CSF flow, 8 as well as for more general flow problems, 10,17 was extensively tested. To determine whether flow velocities had reached a steady state for the time period used, the peak velocities in the fourth and the fifth cycles at specific times in the cycle were compared. Between cycles, velocities varied typically by 1%-2% and up to 6%, indicating that flow was nearly stabilized at the fourth cycle.
To determine the adequacy of the 0.005-second time resolu-  tion, simulations at 0.005 seconds and at 0.001 seconds were compared, with differences usually of 1%-2% and not exceeding 6%. To test the adequacy of the spatial resolution selected, the resolution of coarsest case was increased with 3 times as many computational points. On inspection the maximal velocities differed by less than 10%. To test the validity of applying a uniform pressure gradient, simulations were performed with the gradient increased by 20% posterior to the cord, with minimal differences resulting. To determine the validity of assuming a symmetric pulse wave for CSF flow, simulations with this wave were compared with a simulation in which the pressure wave had a short systole and longer diastole. The asymmetric wave resulted in higher maximal flux and greater bidirectional flow but little change in peak velocities. Assuming rigid walls of the subarachnoid space, though unrealistic, facilitates the simulations. Assessing the effect of wall movement requires additional simulations. We selected oscillating intracranial pressure as the condition producing CSF flow. Whether alternative boundary conditions result in similar conclusions will require additional study. The agreement of the CSF simulation results with PCMR measurements was imperfect, either due to inaccuracy in PCMR measurement of CSF flow 18 or to errors resulting from the selected boundary conditions.