Robustness of MR Elastography in the Healthy Brain: Repeatability, Reliability, and Effect of Different Reconstruction Methods

Changes in brain stiffness can be an important biomarker for neurological disease. Magnetic resonance elastography (MRE) quantifies tissue stiffness, but the results vary between acquisition and reconstruction methods.

MRE is usually performed using an external mechanical transducer that vibrates on the surface of the body. Compressional waves from the vibration are converted into shear waves when passing the tissue interface, causing tissue displacement from both compressional and shear waves. This tissue displacement is encoded into the phase of the magnetization by a modified phase-contrast MRI sequence, using motion-encoding gradients synchronized with the vibrational movement. Next, quantitative maps of tissue stiffness and viscosity are created by applying inversion algorithms to the phase images. 4 MRE is a well-established technique in the liver. However, its primary application in brain studies has been in research. 5 To use brain MRE for diagnostic purposes, the robustness and the reliability of the technique needs to be established. In order to determine true biological changes in tissue stiffness, as caused by therapy, the inherent variation in the measurements needs to be quantified.
Estimates of tissue stiffness in the brain may vary with experimental design, study timing, hardware, vibrational frequency, acquisition methods, and processing pipelines, as well as physiological variations between individuals. 1 Consequently, reproducibility of tissue stiffness estimates across sites is challenging. Moreover, a wide range of reconstruction algorithms are available, all with different underlying physical assumptions. 6 Within one experimental design, both reliability and repeatability should be high in order to track changes over time or to separate between normal and abnormal tissue stiffness. Researchers have reported brain MRE repeatability as measured by coefficients of variation (CV), with values below 10% for subregions in the healthy brain, and below 2% for the brain as a total. [7][8][9][10] In order to trust stiffness estimates, the underlying data need to be of sufficient quality, both in terms of wave propagation and data acquisition. 11 As the tissue displacement caused by the shear waves are the underlying signal prior to processing, MRE data quality is here quantified by the ratio between the shear waves and the compressional waves.
The aim of this study was to assess repeatability and test-retest reliability of MRE in the human brain, and to evaluate the effect of different reconstruction methods and varying MRE data quality on stiffness estimates.

Image Acquisition
The study was approved by the national Research Ethics Committee and the Institutional Review Board. Informed consent was obtained from healthy test subjects. The study used a test-retest design with MRI exams performed on a 3T scanner (Ingenia, Philips Medical Systems, Best, the Netherlands) using a 32-channel head coil. A T 1weighted anatomical reference scan for tissue segmentation was acquired by an 3D inversion recovery turbo field echo acquisition with flip angle = 8 , repetition time (TR) = 5.2 msec, echo time (TE) = 2.3 msec, shot interval = 3000 msec, inversion delay = 853 msec, with field of view (FOV) = 256 mm × 256 mm × 368 mm and 1 mm isotropic resolution.
The MRE was performed using a mechanical transducer placed on the side of the subject's head (Fig. 1). This device induced shear waves at 50 Hz into the brain. 12 Image acquisition was performed with a multishot gradient-echo MRE sequence 10 using a Hadamard encoding scheme 13 with bipolar 13 mT/m motionencoding gradients in three orthogonal directions at 115 Hz. Synchronization to the external wave generator was achieved using a transistor-transistor logic (TTL) trigger signal. A reference scan without motion encoding gradients was used as a phase reference. Fifteen contiguous transversal slices were scanned using an isotropic 3.1 mm FIGURE 1: The transducer setup, thoroughly described previously. 12 (a) The mechanical transducer contains an asymmetrical mass which is rotated by a timing belt on a rod connected to the rotating axis. (b) Illustration of the transducer placement in the head coil illustrated. The transducer is connected to the motor by a flexible rotating axis, shown in black, and a curved plastic head holder connects the head and the transducer. Hearing protection and extra padding to stabilize the head not shown in picture. resolution, a matrix size of 72 × 70, and FOV = 22 cm. Other scan parameters were: flip angle = 20 , TR = 384 msec, TE = 12 msec, Cartesian readout, and a sensitivity encoding factor of 2. 14 Eight equally distributed wave phases were sampled over one oscillation period at 25 Hz. The actual mechanical vibration frequency was shifted to the second index of the Fourier transform, thereby filtering out potential contributions from the frequencies 25 Hz, 75 Hz, and 100 Hz. The total MRE acquisition time was 5.5 minutes. Two MRE acquisitions were performed during the MR exam, 25 minutes apart. The same image geometry was used for both scans. No repositioning of the subject was allowed and the same operator acquired all scans.
In addition to MRE, diffusion tensor images were acquired. These images were acquired using a diffusion-weighted spin-echo, single-shot, echo-planar imaging sequence with Cartesian readout, accelerated by a sensitivity encoding factor 2. Fifty slices of 2.5 mm isotropic resolution were acquired with an acquisition matrix of 94 × 94, FOV = 240 × 240 mm, TR = 9.8 sec, and TE = 60 msec. Fifteen gradient directions and b-values 0 sec/mm 2 and 800 sec/ mm 2 were used, with a total scan duration of 6 minutes.

Image Processing
In this study the MRE phase images were unwrapped by a method based on the mathematical problem of minimum cost flow analysis. 15 Next, pixelwise Fourier transformation of the data was performed in order to obtain the tissue displacement in the frequency domain, before each component of the complex-valued displacement vector u was Gaussian filtered, using a 3D filter of width σ = 0.75 pixels with 3 × 3 × 3 pixels support. The viscoelasticity was then solved using two different state-of-the-art reconstruction methods, both of which have been used in recent scientific work. 12,[16][17][18] Both are direct inversion methods that assume incompressibility, local isotropy, and stiffness homogeneity. The following linear timeharmonic viscoelastic equations governs the wave behavior: where ρ is the tissue density, ω is the angular frequency of the transducer, u are the complex-valued wave displacements, and p is the complex hydrostatic pressure. G* is the complex-valued shear modulus, G* = G 0 + iG 00 , where G 0 is the shear storage modulus or stiffness, and G 00 is the shear loss modulus or viscosity.
The first method performs the inversion by applying the curl operator on the displacement data u. This eliminates the pressure term in Eq. (1) and separates the shear waves from the compressional waves. 19 which is then solved for G* by least squares polynomial fitting. 20 As the reconstruction is performed on the curl field, this method will be denoted curl reconstruction. The second method is a finite element method (FEM) reconstruction. Details are thoroughly described elsewhere. 21 Briefly, however, the method uses a compact FEM with divergence-free basis functions. Data leading to negative shear modulus values are removed and a weighted averaging of the shear modulus based on residual error is also performed. In this case, only first-order spatial derivatives are needed, since no curl operator is applied. The derivatives of u were computed by least squares polynomial fitting. 20 Final stiffness maps were spike-filtered for both reconstructions using a sliding 3 × 3 pixel window. If the central pixel of the sliding window carried the largest or the smallest entry of all values within the sliding window, and this entry was 3 standard deviations away from the average value of the surrounding eight values, this central value was replaced by the average value of the neighboring pixels.
Diffusion analysis was performed in nordicICE (NordicNeuroLab, Bergen, Norway) using the following preprocessing steps: automatic detection of noise threshold, noise level cutoff, and motion correction. Maps of apparent diffusion coefficient (ADC) were produced as previously described. 22

Image Registration and Analysis
Image registration was performed using MatLab (v. R2018a, MathWorks, Natick, MA) and SPM12 (v. 7487, Wellcome Trust Centre for Neuroimaging, London, UK). First, anatomical T 1weighted images were coregistered by affine transformations into the native spaces of each MRE scan. Second, the downsampled images were warped to the Montreal Neurological Institute (MNI) brain region template, 23 and the inverse deformation field was used to reorient the binary maps of the brain regions of interest (ROIs) into the space of the MRE data for each scan. The diffusion images were also warped into the same space.
The following ROIs of the brain were pulled from the MNI template: white matter, the deep gray matter regions caudate nucleus, thalamus, putamen, and hippocampus, the cortical gray matter in the frontal, occipital, parietal, and temporal regions ( Fig. 2a-c), as well as an ROI of the entire brain. To avoid artifacts from MRE reconstruction at the edge of the brain, the outermost voxels in the brain were removed. This was performed using a mask from each MRE scan's magnitude image, which was segmented in nordicICE, morphologically closed, and then eroded at a depth equaling two pixels. ADC maps from the diffusion acquisition were used to reduce any potential errors from partial volume effects. A binary mask with a cutoff ADC value of <1.2 × 10 -3 mm 2 /s was used to exclude voxels with a high content of cerebrospinal fluid. 24 The size of the masked MRE ROIs ranged from 100 to 15,000 voxels.
For the gray matter regions, tissue stiffness normalized to each subject's white matter was calculated in addition to the absolute stiffness measurements.

Statistical Analysis
Mean values of the stiffness in the ROIs were calculated for all subjects, and median values with range were calculated for the 15 subjects. Ratios between white and gray matter were calculated by the median of the ratio of values for white and gray matter based on the average value of both scans for each volunteer.
Test-retest reliability was estimated using both absolute and relative indices. Relative reliability was measured by intraclass correlation coefficients (ICCs) with 95% confidence intervals (CIs), calculated using Stata (release 16, StataCorp, College Station, TX) based on an absolute-agreement, two-way mixed-effects model. 25 Repeatability coefficients (RC) were calculated in MatLab and are defined as 1.96Á ffiffi ffi 2 p σ, where σ is equal to the standard deviation of the measurement differences between the two scans. 26,27 The comparisons of stiffness values from different reconstruction methods were assessed using a paired Wilcoxon signed-rank test. Correlations between shear-compression wave ratio and stiffness were measured using Spearman's rank correlation. A significance level of P < 0.05 was assumed after Holm-Bonferroni corrections for multiple comparisons.

Results
Fifteen healthy volunteers were examined in the study, between the ages of 21 and 33 years (median 27). Of these, six were female and nine male. Figure 3 shows stiffness maps for three subjects, reconstructed using both the curl and the FEM reconstruction method.

Repeatability and Reliability
CURL RECONSTRUCTION. Using the curl reconstruction, the median tissue stiffness in the whole-brain ROI was 1.29 kPa (range 1.01-1.39 kPa, N = 15) in the first scan and 1.28 kPa (range 1.15-1.40 kPa, N = 15) in the second scan ( Table 1). The CV for the measured tissue stiffness in the whole-brain ROI was 4.3%. The CV in the gray matter regions was 4.2-12.8% and 5.1% for the white matter ( Table 2). The ICC between the tissue stiffness of the wholebrain ROI for scans 1 and 2 was 0.74. The ICC for the subregion ROIs ranged from 0.20 in thalamus to 0.89 in the frontal cortex. The RC was 0.14 kPa (95% CI: 0.10--0.19 kPa) for the whole-brain ROI, while the lowest RC was observed for the frontal cortex (0.11 kPa) and the highest for hippocampus (0.43 kPa).
FEM RECONSTRUCTION. Using the FEM reconstruction, the median of average tissue stiffness estimates for the wholebrain ROI was 1.76 kPa (range 1.53-1.91 kPa) in the first  Table 1). The corresponding median CV was 3.8%. The CV was 4.0-9.6% for the gray matter ROIs and 4.5% for white matter ( Table 2). The ICC between the tissue stiffness of the whole-brain ROI for scans 1 and 2 was 0.60, and ranged from 0.15 in putamen to 0.75 in the occipital cortex. The RC for the whole-brain ROI was 0.17 kPa (95% CI: 0.12-0.24 kPa), while the lowest RC was observed in the occipital and the temporal cortex (0.19 kPa) and the highest RC in the putamen (0.48 kPa).
There were no significant differences in the reliability and repeatability measurements between the curl and FEM reconstruction methods: ICC (P = 0.16), CV (P = 0.82) and RC (P = 0.10). Figure 4 shows the stiffness estimates in each scan for all subjects, and Fig. 5 shows the correlation plot and the Bland-Altman plot for tissue stiffness estimates in the subregions of the brain.

FEM RECONSTRUCTION YIELDS HIGHER STIFFNESS
ESTIMATES THAN THE CURL RECONSTRUCTION. For the whole-brain ROI, the FEM reconstruction resulted in 39% higher estimated stiffness compared to the curl reconstruction (P < 0.05). For the individual subregions, no significant difference was found in stiffness estimates from the FEM and the curl reconstruction (P = 0.004) (Fig. 4b).

STIFFNESS VALUES NORMALIZED TO WHITE MATTER.
For the curl reconstruction, white matter was 8% stiffer than gray matter (cortical and deep combined) (P < 0.05). For the FEM reconstruction, white matter was 5% stiffer than gray matter (P < 0.05). Figure 4b shows the distribution of stiffness in ROIs after normalization to white matter. Normalizing to white matter did not lead to significant differences in ICC (curl: P = 0.008, FEM: P = 0.20), CV (curl: P = 0.46, FEM: P = 0.74), or RC (curl: P = 0.008, FEM: P = 0.016) regardless of reconstruction.  Stiffness estimates from scans 1 and 2, from both reconstruction methods and shear-compression wave ratios (MRE data quality). The MRE data quality is higher in the more peripheral regions of interest (ROIs) than in the deep gray matter ROIs toward the center of the brain, with the exception of hippocampus, where the data quality lie between that of the deep gray matter and cortical gray matter regions.

RELATIONSHIP BETWEEN SHEAR-COMPRESSION WAVE
RATIO AND TISSUE STIFFNESS. The median shearcompression wave ratio for the whole-brain ROI was 7.7 (range 5.9-15.5, N = 15) for the first scan and 9.5 (range 6.3-14.9) in the second scan. (Table 1). The median absolute difference between the shear-compression wave ratio of the two scans was 1.7 (range 0.03-4.7). Grouped together, the deep gray matter regions caudate nucleus, hippocampus, putamen, and thalamus had a lower shear-compression wave ratio than white matter (P < 0.05), and the cortical gray matter regions had a higher ratio than white matter (P < 0.05).
No significant correlations were found between the shear-compression wave ratio and stiffness in the whole-brain ROI for either reconstruction method (curl: P = 0.08 in scan 1 and P = 0.04 in scan 2, FEM: P = 0.07 in scan 1 and P = 0.12 in scan 2). (See Fig. 6a,b.) Further, no significant correlations were found between median shear-compression wave ratio for a subregion ROI and the median stiffness in that ROI for either reconstruction method (curl: P = 0.19, FEM: P = 0.43). (See Fig. 6c,d.) ASSOCIATION BETWEEN DATA QUALITY AND REPEATABILITY. No significant correlation was observed between difference in shear-compression wave ratio between scans and stiffness difference between scans (P = 0.09 for curl reconstruction, P = 0.12 for FEM reconstruction), nor between the lowest data quality and the difference in wholebrain stiffness (P = 0.09 for curl reconstruction, P = 0.49 for FEM reconstruction).

Discussion
This study measured tissue stiffness of the healthy brain and its subregions using two different MRE reconstruction methods. The variation in stiffness estimates due to normal biologic and technical factors was quantified, and MRE was found to be a reliable method for assessing brain stiffness. The results suggest that MRE can be used to track changes in tissue stiffness caused by disease, and to track the biomechanical effects of treatment.

Stiffness Estimates Depend on Reconstruction Method
In 15 healthy subjects, stiffness estimates depended on the reconstruction method, where the FEM reconstruction yielded a 39% higher stiffness estimate than the curl reconstruction. This is consistent with earlier results, which found 10-42% higher stiffness values in phantom regions when the data were reconstructed using the FEM compared with the curl approach. 21 As Fovargue et al suggest, this may be due to noise sensitivity. Results are further expected to vary with different acquisition strategies and processing pipelines.  Coefficients of variation (CV), intraclass correlation coefficients (ICC), and repeatability coefficients (RC) for the brain and its subregions, using curl reconstruction and the FEM reconstruction, respectively. There were no significant differences in the reliability and repeatability measurements ICC, CV, nor RC, between the curl-and FEM reconstruction methods.
Circumspection is therefore required when comparing stiffness values obtained with different techniques. Measured stiffness depends on the choice of data acquisition and image reconstruction technique. A way to address this issue it to normalize stiffness measurements to a reference tissue in each subject. 2 In studies of brain cancer patients, tissue stiffness is usually measured relative to healthy white matter. [28][29][30] This study presents both absolute and normalized measurements. Because stiffness has been shown to vary with age and sex, normalizing within each subject is even more critical with a more heterogeneous subject population. 31,32 Relationship Between MRE Data Quality and Estimated Stiffness MRE data quality was quantified by the ratio between the shear waves and the compressional waves, namely, the magnitude of the curl of the displacement field relative to the magnitude of the divergence of the displacement field. Brain tissue is nearly incompressible in vivo. As a consequence, the divergence should be close to zero, while the curl carries the shear signal used to calculate the shear modulus. 16 This study was not able to identify a relationship between MRE data quality and brain stiffness. Regardless of reconstruction method, the estimated correlation was insignificant at the 0.05 level. MRE was performed twice in each subject. For most subjects, the difference in tissue stiffness between the two scans was small, even when the MRE data quality of the two scans differed.
For a subset of subjects, data quality varied substantially between the two scans. Half of the subjects had a relative difference in shear-compression wave ratio of more than 20% between scans. Subject movement may have contributed to these differences. Moreover, despite the differences in data quality, the median difference in measured brain stiffness for these cases was only 2%. Can the difference between stiffness measurements in the first and second scan be due to differences in data quality? For example, are large differences in stiffness measurements associated with substantial differences in MRE data quality? The answer appears to be no. There was no significant correlation between the difference in data quality and differences in tissue stiffness estimates between scans.
Regions of low stiffness, such as caudate nucleus and thalamus, tended toward low shear-compression wave ratios, while regions of high stiffness, such as white matter, tended to have higher shear-compression wave ratios. In general, the shear-compression wave ratio was higher in the brain regions lying closest to the skull, and decreased toward the center of the brain. As the shear waves propagate from the skull inwards, the waves are attenuated, possibly causing a lower MRE signal for the central regions.

Factors Affecting MRE Data Quality
While all subjects were scanned using the same setup, the MRE data quality differed between subjects. The subject with the highest shear-compression wave ratio had twice that of the subject with the lowest shear-compression wave ratio. These challenges may be further exacerbated in a clinical setting. A useful MRE techniques should therefore be robust to data quality.
MRE relies on the transmission of vibrations into the brain tissue. In this study, the vibrations were transmitted by a mechanical transducer placed on the side of the subject's head. Specifically, the transducer was applied on the side of the subject's head, and firmly placed using padding on both sides of the head, leaving no room for head movement. The positioning of the transducer was controlled by visual inspection of images from the localizer scan.
In practice, a central challenge is producing consistent vibrations. Effective transmission requires a firm contact between the transducer and the subject's head. To ensure this contact, the transducer was modified with a curved piece of plastic that lay flush with the side of the test subject's head, slightly behind and above the temple. This piece of plastic was 3D-printed to match the average curvature of the human head, and lined with a thin silicon pad. This is in contrast to the only available commercial MRE hardware, which uses a passive acoustic driver underneath the subject's head. 33 Using a mechanical transducer positioned on the side of the head has potential drawbacks and advantages relative to using a passive driver beneath the head. The setup used in this study may be susceptible to small movements of the transducer relative to the head. For example, there may be less friction with the transducer for subjects with long or smooth hair. However, an advantage of positioning the transducer on the side of the head rather than on a passive driver underneath the head, avoiding that vibrations are dampened by the weight of the head. 12 The mechanical vibration inside the applied transducer is caused by a cable with a rotating central axis. Bends in the FIGURE 4: (a): Summary plot of mean stiffness in whole-brain ROI from scans 1 and 2 in all 15 subjects, with the two scans taken $25 minutes apart, shown for both the curl and the FEM reconstruction. Tissue stiffness from the FEM reconstruction is higher than for the curl reconstruction for all subjects and scans. (b) Median stiffness estimates (average of both scans for each subject), ribbon showing first and third quartile range, from both the curl reconstruction and the FEM reconstruction for all the investigated brain regions. The FEM reconstruction yields higher stiffness values than the curl reconstruction for all regions, while the relationship between regions show similar trends for the two reconstructions, for most of the regions. White matter is measured to be stiffer than gray matter using the curl reconstruction, while white and gray matter show similar stiffness using the FEM reconstruction. (c) Median stiffness estimates normalized to values in white matter, from the curl reconstruction and the FEM reconstruction for all the ROIs. The difference between the reconstructions is more apparent after normalization, as the white matter is stiffer relatively to the other regions with curl reconstruction compared to FEM reconstruction.
axis can lead to upper harmonics in the vibration. This effect was minimized by placing the cable as straight as possible from the motor to the scanner isocenter. A different room layout with the MRI scanner placed differently relative to the motor could affect the data quality. The ideal setup would be one where the axis would go straight without bends from the motor to the isocenter.

Repeatability
The RCs found in this study suggest that in order to track whole-brain stiffness changes in a patient over time, this change would need to be larger than 0.14 kPa for the curl reconstruction, and larger than 0.17 kPa for the FEM reconstruction. For the regions deep in the brain, with higher measurement errors, changes would have to be larger than 0.4 kPa. Earlier studies in patients found a 10% reduction of the whole brain jG*j in patients with Alzheimer's disease compared to healthy controls, 34 while a 10% reduction of the whole brain G 0 was found for patients with normal pressure hydrocephalus, 35 and for patients with multiple sclerosis, 36 compared to healthy controls. As this study has shown that a significant change in whole-brain stiffness would have to be 11% and 10% for the curl and FEM reconstruction, respectively, this could pose a challenge. However, for several of the brain disorders, for example, brain tumors, a potential stiffness change could be local, and a localized stiffness change could be easier to detect than a global stiffness change of the whole brain. For brain tumors, a large heterogeneity of tumor biomechanics has been shown, with G 0 ranging from -60% to +70% of the value in the patient's normal-appearing white FIGURE 6: (a) The shear-compression wave ratio vs. the stiffness in the whole-brain ROI for the curl reconstruction, showing the two scans for each subject linked with a line. (b) Similar for the FEM reconstruction. (c) Shear-compression wave ratio for each subregion vs. the stiffness in that region for the curl reconstruction, where every entry corresponds to the average of the two scans, and shown for all subjects (N = 15). (d) Similar for the FEM reconstruction. The only significant correlation between the global shearcompression wave ratio and whole-brain ROI stiffness was observed for the FEM reconstruction (Spearman's rho = 0.49, P < 0.05 for scan 1 and rho = 0.59, P < 0.05 for scan 2).
matter. When planning a study in patients, the effect size and repeatability needs to be taken into account to be sure to obtain adequate statistical power.
In addition to assessing repeatability, test-retest reliability was estimated using both relative and absolute indices. The results showed moderate reliability in terms of ICC. The ICC, unlike the CV, depends on the underlying distribution of subjects. Hence, a more homogeneous distribution will result in a lower ICC, as ICC is a measure of how much the measurements for one subject align compared to those of other subjects. For a homogeneous subject group, such as the one used in this study, the test-retest reliability is better assessed by repeatability coefficients, which provide the reliability assessments in kPa, the same units as the tissue stiffness estimates.

Limitations
A general challenge for work in this area is the lack of a gold standard for in vivo tissue stiffness measurements. Limited knowledge of the underlying true values makes comparing reconstruction methods with different stiffness estimates difficult.
The applicability of the results from this study is limited by the relatively small sample comprised of healthy volunteers. Future work should include a larger cohort, and include patients with neurological disease in addition to healthy volunteers. Furthermore, the MRE-derived shear modulus is a frequency-dependent quantity, so measures performed at 50 Hz will only be valid at 50 Hz. 37 This frequency was chosen as it balances resolving power and penetration, and is commonly used in MR elastography. 38 In this study the test and the retest scan were performed during the same session. The goal was to minimize the factors that could affect the result. A future study where subjects were repositioned between scans or scanned on different days would be clinically interesting and add to the understanding of factors affecting repeatability, as larger variability can be included both in subjects' biological variation, but also introducing uncertainty in varying the placement of equipment. For a full evaluation of brain MRE reliability, reproducibility of stiffness estimates between imaging systems and manufacturers is warranted.

Conclusion
This study found MRE of the human brain to be a robust technique in terms of repeatability and reliability. The RC values suggest that in order to track stiffness changes in a patient over time, the change needs to be on the order of 0.2 kPa, and up to 0.5 kPa for deeper-lying regions with higher measurement error. The estimated tissue stiffness was higher when using the divergence-free FEM reconstruction than using the curl-based reconstruction. Data quality was higher in the more peripheral brain regions than in its central regions. Caution is warranted when comparing stiffness values obtained with different techniques, and normalizing stiffness values to healthy matter is recommended in a patient population.