Evaluation of MR Elastography as a Noninvasive Diagnostic Test for Spontaneous Intracranial Hypotension

Patient Selection

After we obtained institutional review board approval and written informed consent, brain MRE examinations were performed on 15 patients with SIH who were diagnosed with a CVF from September 2022 to April 2023. Dynamic myelogram images diagnosing the CVF were obtained and interpreted by one of our board-certified neuroradiologists who specializes in CSF leaks. The CVF was subsequently verified by an additional neuroradiologist (I.T.M.). All patients met The International Classification of Headache Disorders-3 criteria for SIH. The number of days between diagnostic myelography and MRE was recorded. MRE was performed before catheter embolization treatment. Patient age, sex, and symptom duration were recorded. Brain MRIs were reviewed for the Bern score.⁹

Image Acquisition

MRE data and T1-weighted anatomic images were collected on 3T MR imaging scanners (GE Healthcare). All MRE examinations were performed at 60 Hz with 3-mm resolution.

Data from 44 control participants (age range, 56–89 years) were included from a previously published study.¹⁰ This cohort was scanned on a 3T GE Healthcare scanner and confirmed to be negative for amyloid-PET and cognitive impairment. In this cohort, an inversion recovery echo-spoiled gradient echo (the pulse sequence was used for anatomic imaging with parameters including an imaging matrix of 256 by 256 pixels, TR/TE of 6.3/2.8 ms, flip angle of 11°, inversion time of 400 ms, FOV of 27 cm, and bandwidth of 31.25 kHz). The MRE data were acquired using a spin-echo echo-planar imaging technique at 60-Hz vibration. The acquisition parameters were the following: TR/TE of 3600/62 ms, FOV of 24 cm, imaging matrix of 72 × 72 reconstructed to 80 × 80, 48 contiguous 3-mm-thick axial slices, one 18.2 -ms motion-encoding gradient on each side of the refocusing radiofrequency pulse, motion-encoding in x, y, and z directions, and 8 phase offsets spaced evenly during a period corresponding to 60-Hz motion.

The remaining 21 control participants (range, 19–53 years of age) and 15 patients with SIH (range, 35–70 years of age) were scanned on a compact GE Healthcare 3T system.^11,12 The controls in this cohort did not have any neurologic conditions or increased intracranial pressure. Anatomic images were obtained using a T1-weighted MPRAGE acquisition with parameters of TR/TE of 6.1/2.5 ms, inversion time of 600 ms, flip angle of 8°, FOV of 160 × 160 mm, image matrix of 256 × 256, section thickness of 1.2 mm, and acquisition time of 4 minutes 18 seconds. MRE acquisition was performed using a flow-compensated spin-echo echo-planar imaging method, and vibrations were applied by a pneumatic actuator with a frequency of 60 Hz. The acquisition parameters were TR/TE of 4000/59 ms, FOV of 24 × 24 cm, image matrix of 80 × 80, 48 contiguous 3-mm-thick axial slices, one 16.7 -ms 8-G/cm motion-encoding gradient on each side of the refocusing pulse, motion-encoding in x, y, and z directions, and 8 phase offsets spaced evenly over 1 period corresponding to a motion of 60 Hz.¹³

Viscoelastic Property Map Estimation

Stiffness and damping ratio maps were computed by neural network inversion.¹⁴ The shear modulus of a viscoelastic material at a given frequency is given by the complex valued quantity G = G’ + iG″. From this modulus, stiffness is computed as 2|G|²/(G’+|G|), and the damping ratio is the defined as .^15,16 Training data were generated with a finite difference model of the wave equations, assuming harmonic motion in linear viscoelastic materials using a random distribution of shear moduli (stiffness range, 1–5 kPa; damping ratio range, 0–0.5), which was described in detail in previous studies.^14,17 The true mechanical properties were spatially varied to relax the tissue homogeneity assumption. From the training set, random patches were selected and randomized by the application of phase cycling, noise, and masking. Separate neural networks were trained to estimate the stiffness and damping ratio using Keras and TensorFlow backend.^18,19 Training was performed with an Adam optimizer using a batch size of 100 and 1000 batches per epoch.²⁰ Two learning rates were used (0.001 and 0.0001), and learning was stopped at each rate when the mean squared error did not improve for 3 consecutive epochs.

To apply the Neural Network Intelligence (NNIs) in vivo, first, we computed tissue probability maps for each participant using the unified segmentation algorithm in SPM12 software (http://www.fil.ion.ucl.ac.uk/spm/software/spm12) and tissue priors from the Mayo Clinic Adult Life span Template (MCALT).^21,22 T1-weighted images were registered and resliced to the MRE magnitude image (T2-weighted), and segmentation was performed according to both channels to obtain tissue probability maps in MRE space. A brain mask was computed to indicate voxels with a combined probability of white and gray matter tissues greater than the probability of CSF. The inverse deformation field was applied to the MCALT lobar atlas to obtain these region assignments in MRE space.

Displacement data were filtered to reduce interslice phase discontinuities,²³ unwrapped by using a graph cuts algorithm,²⁴ and an adaptive curl operator was applied to reduce the effects of longitudinal waves on shear property estimation. Each examination was processed in 3 subregions to avoid processing across the major dural folds, which act as shear wave sources. These subregions, defined as the intersection of relevant lobar atlas assignments and the brain mask, were the left cerebrum plus the corpus callosum, right cerebrum plus the corpus callosum, and the cerebellum plus the brainstem. Because the corpus callosum is included in 2 subregions, final estimates were computed as the average of the 2 subregion estimates. This procedure was repeated for the stiffness and damping ratio separately.

Voxelwise Mapping

Using the previously computed deformation field, we warped mechanical property maps into the MCALT space using nearest-neighbor interpolation. Then at each voxel, a linear model was fit to stiffness estimates at each voxel with predictors, including a constant, age, age², sex, and a categoric variable for the diagnosis of SIH. These predictors were included because the brain stiffness has been shown to be associated with age and sex,¹⁰ and these effects could obscure the SIH effect. A quadratic term for age was included because of the wide age range of study participants, resulting in a significant nonlinear age effect on both mechanical properties, as shown in the Online Supplemental Data. A t test was performed on the SIH predictor, and false discovery rate (FDR) correction was applied according to the method of Storey, with Q < 0.05 considered significant. An additional categoric predictor for the scanner was evaluated, but the scanner predictor was excluded because none of the voxels were found to have a significant effect after the FDR correction.

Pattern Analysis

We used a previously described pattern analysis to summarize each participant’s MRE findings.²⁵ For computing an individual’s pattern scores with leave-one-out cross-validation, first their viscoelastic map was held aside and voxelwise modeling was performed in the remaining participants (with the same predictors as above). The map of the held-out individual was corrected for age, sex, and the mean viscoelastic property on the basis of the modeling result. Then a correlation coefficient was calculated between the corrected viscoelastic property map and a map of the expected SIH effect and converted to a z score via the Fisher transformation. For each participant, the z scores were calculated separately for the stiffness and damping ratio maps and used as features in a support vector machine (SVM) classifier with a linear kernel and box constraint of 1. Using leave-one-out cross-validation, we evaluated the SVM model for its ability to distinguish the SIH group from the control group as measured by accuracy and the area under the receiver operating characteristic curve. Including the scanner effect slightly improved the estimated area under the receiver operating characteristic curve and accuracy, but the scanner effect was not statistically significant. Therefore, we excluded it and opted to report the more conservative result.