Reverberant magnetic resonance elastographic imaging using a single mechanical driver

Abstract Reverberant elastography provides fast and robust estimates of shear modulus; however, its reliance on multiple mechanical drivers hampers clinical utility. In this work, we hypothesize that for constrained organs such as the brain, reverberant elastography can produce accurate magnetic resonance elastograms with a single mechanical driver. To corroborate this hypothesis, we performed studies on healthy volunteers (n = 3); and a constrained calibrated brain phantom containing spherical inclusions with diameters ranging from 4–18 mm. In both studies (i.e. phantom and clinical), imaging was performed at frequencies of 50 and 70 Hz. We used the accuracy and contrast-to-noise ratio performance metrics to evaluate reverberant elastograms relative to those computed using the established subzone inversion method. Errors incurred in reverberant elastograms varied from 1.3% to 16.6% when imaging at 50 Hz and 3.1% and 16.8% when imaging at 70 Hz. In contrast, errors incurred in subzone elastograms ranged from 1.9% to 13% at 50 Hz and 3.6% to 14.9% at 70 Hz. The contrast-to-noise ratio of reverberant elastograms ranged from 63.1 to 73 dB compared to 65 to 66.2 dB for subzone elastograms. The average global brain shear modulus estimated from reverberant and subzone elastograms was 2.36 ± 0.07 kPa and 2.38 ± 0.11 kPa, respectively, when imaging at 50 Hz and 2.70 ± 0.20 kPa and 2.89 ± 0.60 kPa respectively, when imaging at 70 Hz. The results of this investigation demonstrate that reverberant elastography can produce accurate, high-quality elastograms of the brain with a single mechanical driver.


Introduction
Magnetic resonance elastography (MRE) is an imaging technique for determining the mechanical properties of tissues noninvasively and in vivo (Muthupillai et al 1995), and the current gold standard imaging method for diagnosing liver fibrosis (Yin et al 2007, Venkatesh et al 2013. Given the success of MRE in this application, several studies are now focused on investigating MRE's utility in other clinical applications. These include improving the differential diagnosis of breast cancer (McKnight et al 2002, Sinkus et al 2007, Patel et al 2021, identifying tears in skeletal muscles, (Dresner et al 2001), detecting pulmonary disease in lungs (Mariappan et al 2014), and diagnosing prostate cancer (Brock et al 2015). MRE could also prove helpful in assessing the progression of Alzheimer's disease (Murphy et al 2011), Parkinson's disease (Lipp et al 2013), multiple sclerosis (Wuerfel et al 2010), brain integrity and microstructural changes in health and disease (Sack et al 2013), and evaluating normal pressure hydrocephalus (Streitberger et al 2011, Freimann et al 2012. Li and colleagues reviewed the crucial elements common to all successful magnetic resonance elastographic imaging systems (Li et al 2014): appropriate mechanical stimulation of the organ under investigation, acquiring wave images with a good signal-to-noise ratio (SNR), and computing robust estimates of shear modulus from the measured wave fields. In this work, we focus on the third element, i.e. efficiently computing shear modulus.
Researchers have proposed different approaches for computing shear modulus that vary in accuracy and computational efficiency (Doyley 2012). One method computes shear modulus directly from local estimates of wavelengths (Manduca et al 1996). Although this is computationally efficient, estimating the wavelength in complex organs, such as the brain, is difficult because waves reflecting from the skull and internal structures superimpose to create complex shear wave fields (Muthupillai et al 1995). The local frequency estimator (LFE) is an alternative method for estimating shear modulus (Knutsson et al 1994, Kruse et al 2000, Hu 2020). Like the local wavelength estimation (LWE) approach, the LFE method is computationally efficient but produces erroneous shear modulus when applied to complex wave fields (Hiscox et al 2016). Hu et al (2020) proposed an enhanced local frequency estimator (ELFE) that used directional filters to eliminate undesirable reflections. They demonstrated that ELFE produced more accurate shear modulus estimates than the conventional LFE approach and reduced far-field artifacts (i.e. artifacts in regions far from the wave source). Researchers have proposed several direct inversion algorithms to overcome challenges incurred when estimating shear modulus with either LWE or the LFE approach by algebraically solving for the complex shear modulus from the Helmholtz equation (Manduca et al 2001, Oliphant et al 2001, Papazoglou et al 2008, Barnhill et al 2018. Although these direct inversion schemes are fast and accurate, they are more susceptible to noise than the LFE, LWE, or ELFE methods (Hu 2020) investigated the impact of measurement range on two shear modulus estimation approaches, the algebraic-inversion-of-differential-equation (AIDE) and the local frequency estimator. Using the wavelength-to-pixel size ratio performance metric, they revealed that AIDE incurred significant errors when the wavelength-to-pixel ratio was less than 10. In contrast, the LFE method incurred errors only when the wavelength-to-pixel ratio was less than 2, showing its superiority over the AIDE method. To improve the robustness of the direct inversion method (Barnhill et al 2018), developed a multi-frequency inversion approach that incorporates first-order gradients and combines shear modulus estimates from a narrow range of frequencies. Researchers have also used filtering schemes to enhance performance (Scott et al 2020). However, excessive filtering degrades spatial resolution. An artificial neural network has recently been used to reconstruct the shear modulus distribution (Scott et al 2020). Neural networks should provide fast and reliable shear modulus estimates once the neural network is sufficiently trained; however, their performance in different clinical scenarios has yet to be revealed. Iterative inversion methods offer the opportunity to model heterogeneous, viscoelastic tissues appropriately (Van Houten et al 1999, 2001a, Doyley et al 2000. This inversion approach is robust but computationally expensive, requiring several hours to compute highresolution elastograms. Our long-term objective is to integrate MRE into our clinical workflow. More specifically, to develop methods to provide accurate MR shear modulus elastograms at the MR console when imaging the brain. This paper revisits the local wavelength estimation approach by considering the shear modulus estimation problem as a reverberant problem. More specifically, we seek to estimate the local wavelength of complex wave fields using a technique known as reverberant elastography. Reverberant elastography uses multiple point sources to generate complex wave fields (Parker et al 2017); the resulting wave fields' autocorrelation function provides reliable local wavelength estimates. However, utilizing many mechanical drivers to produce complex wave fields can hamper clinical utility. More specifically, performing MRE with multiple drivers is impractical for many clinical applications and could prove uncomfortable for patients. In this work, we hypothesize that we can produce reliable elastograms with a single mechanical driver in constrained organs such as the brain, where complex wave fields are generated naturally. This hypothesis was based on the observation that the skull has many surfaces that act as point mechanical sources (Clayton et al 2012, Smith et al 2020. To corroborate this hypothesis, we performed studies on a constrained gelatin phantom and healthy volunteers (n = 3). We used a similarity metric to quantify the degree of reverberance induced in different displacement components. We evaluated the performance (accuracy) of reverberant elastograms relative to those computed using the subzone inversion method (Van Houten et al 2001a, Doyley et al 2010). Also, we assessed the feasibility of recovering shear modulus from a single component of the reverberant wave field, because doing so would overcome a limitation of the subzone inversion method that requires the entire 3D displacement field to produce reliable elastograms. Currently, the 3D displacement field is acquired by applying a phase contrast pulse sequence (Weaver et al 2001) three times, one for each component. Reverberant elastography could reduce MRE acquisition time (by a factor of three) if it produced accurate elastograms from a single displacement component.

Reverberant shear wave elastography
The general principles of reverberant elastography have been previously described (Ormachea et al 2018, Zvietcovich et al 2019, Ormachea and Parker 2021; therefore, in this section, we provide a summary of the approach. Complex wave fields are produced when plane waves that originate from multiple point sources or reflected from various angles superimpose. The particle velocity V (ε, t) of the complex wave field is given by Parker et al (2017): where t and ε represent the time and position in the complex (reverberant) wave field, respectively; k represents the wavenumber, and ω 0 the angular frequency. The subscript q denotes the random unit vector, and the subscript l denotes a unit vector, n , ql parallel to a disk formed by orthogonal basis vectorsAE andf (Parker et al 2017). The random variable v ql describes the particle velocity magnitude. For isotropic mediums, the autocorrelation of the wave field in a plane transverse to the detected motion vector is given by (Parker et al 2017): where B vv represents the 2D autocorrelation of V (ε, t), b represents the expected value of the squared particle velocity magnitude, k is the wavenumber, and j 0 and j 1 are spherical Bessel functions of the first kind of order 0 and 1, respectively. Since n ql and v ql are independent of position, an ersatz form of equation (1) can be written as: where  ( ) AE represents the spatially varying phase, and V o is related to the root mean squared amplitude of the field. The wavenumber, k, is computed from the ensemble average of the reverberant field as follows: where A is a scaling constant, which was determined empirically to be one in this study using known measurement of shear wave speed. The bracket denotes the average value over a homogeneous kernel.
In this study, we acquired MR motion data over four-time points or equally spaced phase offsets. The Fourier transform of these temporal data provided complex motion at the vibration frequency. We computed the phase angle of the motion data at each pixel in the imaging field of view (FOV). Assuming phase varies in three dimensions: x y z k x k z k z c , , , 5 where k x , k y and k z are components of the wave vectors in the x, y and z coordinate directions, respectively, and c o is a constant. Substituting equation (5) into equation (4) gives an approximate estimate of the wavenumber: We applied the two-dimensional unwrapping algorithm described in Zhao et al (2018) to phase maps acquired from each coordinate direction. We computed local estimates of wave number by using the singular value decomposition method (Strang 2016) to fit a plane to the unwrapped phase maps within threedimensional overlapping kernels. Shear wave speed (c) was computed from local estimates of wavenumber as follows: where f represents the shear wave frequency. Local estimates of shear modulus, m, were computed from the estimated shear wave speed as follows (Parker et al 2010(Parker et al , 2011: where, ρ represents tissue density (1 g cm −3 ). In this study, shear wave speed was assumed to be isotropic. The shear modulus of the tissue or phantom under investigation was estimated by applying the reverberant method to each phase direction and then computing the average of the resulting images (i.e. the composite shear modulus elastograms).

Subzone elastography
The overlapping subzone inversion method computes shear modulus elastograms from MR-measured internal tissue displacements by combining the finite element method and the Newton-Raphson iterative scheme, previously described (Van Houten et al 1999, 2001b, Doyley et al 2003. This inversion approach seeks the distribution of mechanical parameters (in our case, lambda (λ) and shear (μ) modulus) that minimize the difference between internal tissue displacements calculated with the finite element model and those measured with MR. To reduce the memory required to solve the three-dimensional inverse elasticity problem on high- where U z (μ z , λ z ) and U z m represent vectors of the calculated and measured displacements at the nodal coordinates of each subzone, respectively. μ z and λ z represent shear and lambda modulus at the nodal where J z (μ z , λ z ) is the n × n Jacobian matrix, and a is a positive number that was used to improve the condition of the Hessian matrix,

Phantom experiments
This study aimed to show that reverberant elastography produces reliable shear modulus estimates from complex wave fields induced in constrained objects, with multiple reflecting surfaces using a single mechanical driver. To demonstrate this, we compared the performance (accuracy and contrast-to-noise ratio (CNR)) of reverberant elastograms to subzone elastograms. We also explore the feasibility of recovering shear modulus from a single component of the reverberant wave field. Independent mechanical testing provided absolute shear modulus values.

Phantom fabrication
We fabricated a brain-shaped phantom ( Louis, MO, USA) in a highly controlled and repeatable manner as described in Doyley et al (2003). The phantom contained three spherical gelatin inclusions with diameters of 18 mm, 12 mm and 4 mm. Table 1 gives the percentage by weight of the gelatin, water, and copper sulfate used to fabricate the surrounding background and inclusions.

Elastographic imaging
We performed all elastographic imaging in a whole body 3T MRI scanner (Prisma, Siemens, Erlangen, Germany) with a 20-channel head coil. A pneumatic actuator with a passive driver (Resoundant, Inc., Rochester, MN, USA) was used to induce shear waves in each phantom, as illustrated in figure 1(a). We used two vibration frequencies (50 Hz and 70 Hz) separately during elastographic imaging, with motionencoding gradients matched in the period to vibration and with a variable number of gradients depending on frequency. It took approximately six minutes to acquire each 3D data set. Table 2 summarizes the actuator amplitude, gradient amplitude, and the number of gradients used at each frequency. The singleshot echo-planar imaging (EPI) sequence (Johnson et al 2014, Chaze et al 2019) measured the resulting time-varying harmonic tissue displacements. We configured the MR scanner with echo and repetition times of 76.0 ms and 8640 ms, respectively. Forty axial slices were acquired for the phantom with a 153 mm × 153 mm × 60 mm field-of-view (1.6 mm isotropic voxel size), with four images with relative phase offsets. For brain, eighty axial slices with four phase offsets were acquired with a 240 mm × 240 mm × 120 mm field-ofview (2.5 mm isotropic voxel size). We used the Fourier transform method described in Sinkus et al (2000) to estimate the complex, three-dimensional displacement field. To remove low-frequency longitudinal waves and high-frequency noise, we applied a two-dimensional bandpass filter in all directions. The cutoff spatial frequency related to the wavenumber k of the filter was determined from pre-selected low (c l ) and high (c h ) shear wave speed values. The corresponding filter cutoffs were k l = (2 × k × π × f )/c h , and k h = (2 × π × f )/c l . For this study c h was 3.5 m s −1 , and c l was 0.3 m s −1 . These values were selected using phantoms with known shear modulus and an average shear modulus of the whole brain reported in the literature. Two groups of modulus elastograms were computed from each data set; one was computed using the reverberant shear modulus estimation method and the other with the subzone reconstruction method. The reverberant method was performed using 6.4 mm × 6.4 mm × 6.4 mm overlapping kernels. All subzone reconstructions were conducted on a finite element mesh consisting of 70 000 nodes and 415 000 elements (created using MATLAB version R2022b). A spatial filtering weight of 20% and subzone radius size of 9 mm and a regularization value of 1 e −7 were also employed during subzone reconstructions. A homogeneous trial solution (shear and lambda moduli of 0.33 kPa and 33 kPa, respectively) was assumed at the start of all subzone reconstructions. Reconstructions were terminated either after 100 global iterations or when the relative error of the global objective function did not decrease significantly (2%) for 10 consecutive iterations, whichever condition occurred first. In general, it took 7 h to compute subzone elastograms and 3 min to compute reverberant elastograms on an Intel Xeon Gold 6330 CPU computer system (20 cores) running at 2 GHz (Dell Technologies, Round Rock, Texas USA).   (Fung 1981): where v is the poison's ratio, assumed to be 0.495 in this work. The ground truth for the quantitative evaluations was computed from the average of five statistically independent measurements. The actual shear modulus was 8.15 ± 0.05 kPa for the inclusions and 3.34 ± 0.04 kPa for the surrounding gel, which was consistent with previously reported values for brain tissues (

Performance metrics
Elastograms from each reconstruction method were visually inspected for quality, and four quantitative metrics were used to evaluate performance. The octahedral shear strain signal-to-noise ratio (OSS-SNR) (McGarry et al 2011) was used to evaluate the quality of the measured displacement fields. The quality of the modulus elastograms was evaluated quantitatively by computing the contrast-to-noise ratio (CNR) performance metric. CNR was defined on a logarithmic scale as follows (Techavipoo and Varghese 2005): where μ b and μ i represent the mean shear modulus chosen from regions-of-interest (ROIs) in the background and inclusion, respectively, while σ b and σ i represent the standard deviation of the shear modulus in the corresponding ROIs. The mean error (ME) performance metric was used to evaluate the accuracy of the recovered modulus compared to mechanical testing as follows: here μ e and μ ref represent the estimated and from mechanical testing shear modulus, respectively.
The similarity (χ) between the measured and theoretically derived autocorrelation (see equation (2)) was used to quantify each displacement component's reverberance level. For a given window, the similarity metric was computed as follows where x i and y i represent vectors of the measured and theoretically computed autocorrelation function, and N is the number of elements in each vector. We used a threshold to classify wave fields as complex (reverberant) or directed. This threshold was determined empirically by simulating complex wave fields created with increasing numbers of incident shear waves (i.e. from 10 to 100 with an increment of 10), as described in (Zvietcovich et al 2019). For each wave field, we computed the autocorrelation function of the simulated wave field and those computed theoretically using equation (2). We observed that the similarity metric increased rapidly as the number of incident shear waves increased from 10 to 50, then plateaued to 80% with a further increase in incident waves (not shown). 80% represents the threshold for differentiating fully reverberant (n = 100) from non-reverberant wave fields (N < 50). Therefore, in this study, pixels in the complex wave field whose similarity matrix (χ) exceeded 80% were deemed reverberant.

Clinical study
To evaluate the performance of reverberant elastograms under standard clinical conditions, we performed brain MRE imaging on three healthy volunteers. We acquired two datasets from each volunteer using an imaging protocol approved by the University of Delaware Institutional Review Board. In one acquisition, we induced 50 Hz shear waves within the brain. Without moving the volunteers, we increased the frequency of the induced shear waves to 70 Hz. In both acquisitions, an inflatable pillow (see figure 1(b)) positioned at the base of the skull generated shear waves within the brain. Table 3 provides the scanning parameters used in clinical imaging; all other imaging parameters are consistent with the phantom experiment. We used the FMRIB Software Library (FSL), i.e. FLIRT registration and FAST segmentation toolboxes (Smith et al 2004), to register and segment brain elastograms.

Results
The proceeding subsection reports the results of experiments conducted on heterogeneous phantoms and volunteers to evaluate the performance of reverberant elastograms relative to subzone elastograms.  when imaging at 50 Hz ((a), (c), (e)) and 70 Hz ((b), (d), (f)). The percentage of reverberant pixels (i.e. pixels whose similarity matrix exceeded 80%) in the X, Y, and Z displacement components was 89%, 88% and 91%, respectively, when imaging at 50 Hz. Similarly, the percentage of reverberant pixels in the X, Y, and Z displacement components was 92%, 91% and 88%, respectively, when imaging 70 Hz. Although the measured wave field appears to be more reverberant when imaging at 70 Hz, statistical evaluation of the wave fields computed with the Kolmogorov-Smirnov test showed no statistically significant difference (p > 0.05) in wave fields obtained at 50 and 70 Hz.  computed from individual displacements) and subzone elastograms was 6.8% higher and 1.9% lower respectively than its actual value when imaging at 50 Hz ( figure 5(a)). For reverberant elastograms, the mean shear modulus of the 18 mm inclusion computed from the X, Y and Z components of displacements was 1.3%, 3.3% and 14.1% higher than the actual modulus. The mean shear modulus of the 12 mm diameter inclusion ( figure 5(b)) estimated from subzone elastograms was 2.6% lower than the actual value when imaging at 50 Hz. Composite reverberant elastograms and elastograms computed from the X, Y and Z displacement components underestimated the shear modulus of the 12 mm diameter inclusion by 14.1%, 16.6%, 15.3 and 4.3%, respectively. We observed a similar trend for the estimated shear modulus of the background gel figure 5(c). In this case, the composite reverberant inversion was more accurate than the subzone inversion, and the most accurate reverberant elastograms were estimated from the Y displacement component. Figure 5(d) shows the CNR computed from the elastograms shown in figures 4(a)-(j). The CNR of subzone elastograms was comparable to those produced with the reverberant inversion method. The CNR of composite reverberant elastograms and those computed from the X displacement yielded the highest and lowest CNR, respectively, when imaging at 50 or 70 Hz.

Clinical studies
3.2.1. Quantitative assessment displacement fields induced within the brain when imaging at 50 Hz and 70 Hz Figures 6(a)-(r) shows representative examples of MR magnitude images and MR elastograms obtained from a healthy brain. Like the phantom studies, there was no statistically significant difference (p > 0.05) in wave fields obtained when imaging at 50 and 70 Hz. The levels of reverberance incurred in the brain were comparable to those incurred in the phantoms. More specifically, the percentage of reverberant pixels in each displacement  component (X, Y, and Z) was 90%, 88% and 92% when imaging at 50 Hz and 91%, 87% and 93% when imaging at 70 Hz. The values were similar (88%-94%) for the other two brain data.

Assessment of shear modulus values of different brain tissues from reverberant and subzone elastograms
Representative reverberant (figures 6(g)-(l)) and subzone figures 6(m)-(r) nt coronal positions within the brain when elastographic imaging was performed at 50 Hz (figures 6(g), (h), (i), (m), (n), (o)) and 70 Hz (figures 6(j), (k), (l), (p), (q), (r)). The global brain stiffnesses (GBT) estimated from subzone and reverberant elastograms shown in figure 6 were 2.5 ± 0.96 kPa and 2.3 ± 0.89 kPa, respectively, when imaging at 50 Hz. When imaging at 70 Hz, the GBT estimated from subzone and reverberant elastograms were 2.5 ± 0.82 kPa and 2.89 ±. 83 kPa, respectively. Figure 7 shows representative elastograms obtained from three different slices when the reverberant method was applied to individual displacement components when imaging at 50 Hz (figures 7(a)-(i)) and 70 Hz (figures 7(j)-(r)). The GBT estimated from X, Y, and Z motion component elastograms were 2.1 ± 0 .88 kPa, 2.4 ± 0.98 kPa, and 2.3 ± 0.99 kPa, respectively, when imaging at 50 Hz. When imaging at 70 Hz, the GBT estimated from the X, Y, and Z displacement components were 2.4 ± 0.82 kPa, 2.7 ± 0.75 kPa, and 2.89 ± 0.92 kPa, respectively. Figures 8(a), (b) shows box plots of the mean shear modulus of the whole brain, white matter, and gray matter estimated from composite reverberant and subzone elastograms for the volunteers employed in this study for 50 Hz and 70 Hz. At both frequencies, the global shear modulus of the brain was consistent with values reported in Murphy et al (Murphy et al 2013), Ingolf et al (Sack et al 2011), Matthew et al (Murphy et al 2011). and Arani et al (Arani et al 2015). The shear modulus of the white and gray matter was consistent with results previously reported (Zhang et al 2011, Johnson et al 2013a, 2013b. Figures 8(c), (d) shows box plots of the mean shear modulus of the whole brain, white matter, and gray matter estimated from reverberant elastograms computed from individual displacement components for 50 Hz and 70 Hz. Table 4 summarizes the shear modulus of the entire brain, white matter, and gray matter computed from reverberant and subzone elastograms for the volunteers employed in this study.

Discussion
This study introduced an analog of the reverberant shear wave elastographic imaging technique previously developed for ultrasound elastography (Parker et al 2017), now adapted for MRE. Reverberant elastography assumes a superposition of waves generated by multiple point sources, and reflections will generate a complex 3D shear wave field. Under these conditions, the two-dimensional autocorrelation function of the complex wave field and the phase gradient provides fast and accurate estimates of the wavenumber from which shear wave speed and modulus are derived. Reverberant elastography has been demonstrated with ultrasound (Parker et al 2017) and optical coherence tomography (Zvietcovich et al 2019, Ge et al 2022) using multiple point drivers. However, to our knowledge, this is the first reported study demonstrating its utility with MRE and its feasibility in the brain using a single mechanical driver. We used the goodness of fit performance metric to quantify the degree of reverberance in different displacement components. We also studied three healthy volunteers to assess how reverberant elastography performs within the clinical setting. The primary findings of this investigation were as follows: 1. A single mechanical source generates a complex wave field (reverberant) in constrained phantoms, quantified by assessing how well the measured autocorrelation function matched the theoretically derived one. Applying the reverberant shear modulus estimation method to the resulting shear wave fields produced elastography with variable accuracy (83.2%-98.7%) and CNR ranging from 63.1 to 73 dB (figure 5).
2. The performance (accuracy and CNR) of subzone and composite reverberant elastograms were comparable. However, the performance of reverberant elastograms degraded marginally when reverberant elastograms were computed from a single displacement component ( figure 5).
3. The reverberant reconstruction method produced meaningful elastograms when applied to one displacement component, but performance varied based on the displacement component used. We observed the best performance when reconstructions were performed using all displacement components. This raises concerns about the usefulness of reverberant elastograms produced from a single displacement component obtained from a partially reverberant field.
4. The mean shear modulus of the whole brain, white, and gray matter estimated from composite reverberant and subzone elastograms obtained from healthy volunteers were consistent with previously reported shear modulus estimates of the healthy brain.
In ultrasound-reverberant elastography, elastograms acquired at higher frequencies were superior to those obtained at lower frequencies (Parker et al 2017, Ormachea et al 2019; in this study, the quality of elastograms produced at 50 and 70 Hz was similar. Figure 2 demonstrates the impact of attenuation on the displacements when images are acquired at 50 Hz and 70 Hz. Despite the reduction in shear wave amplitude, the displacement  amplitudes at 70 Hz were high enough to produce reasonable elastograms. Although we increased the amplitude of shear waves when imaging at 70 Hz, however, since the mechanical driver is pneumatic, the amplitude of the shear wave at the passive driver is likely to be much less than at the active driver. Therefore, future studies will measure the amplitude shear waves at both the active and passive drivers during imaging. The composite elastograms and those computed from the individual displacement components overestimated the modulus of the 18 mm diameter inclusion (see figure 5(a)). The 12 mm diameter inclusion shear modulus was underestimated in all cases (see figure 5(b)). The cause of this behavior is unclear, but it may be due to using a suboptimal kernel or a partially reverberant field. We expect that using larger kernels will increase the CNR but cause the modulus of smaller inclusion to be underestimated. Employing a larger kernel would also reduce the background noise observed in the elastograms (see figure 5(c)). In contrast, smaller kernels should provide more accurate results, although with a lower CNR. Further work is needed to fully understand the kernel size's impact on the experiment's performance and to determine the underlying cause of the observed behavior.
Model-based inversion approaches need all three displacement components to compute shear modulus precisely (Doyley 2012). Researchers have demonstrated theoretically that computing shear modulus with fewer components reduces the accuracy of the resulting elastograms (Skovoroda et al 1994). Sampling error and poor ultrasound penetration, especially in patients with high body mass index (Zhao et al 2014) cause the diagnostic performance of ultrasound elastography for detecting and staging liver fibrosis to be lower than that achieved with MRE (Li et al 2021). Figure 5 demonstrates that reverberant shear wave elastography can create reasonably accurate elastograms (on the order of 16.8% error) from a single displacement component. Although shear elastograms computed from individual displacement components may differ in appearance (see figure 7), the average modulus computed from each brain region was similar (see figure 8). Using only a single displacement component when performing reverberant MR elastography is not recommended due to inconsistent results obtained with different displacement components. We recommend that composite reverberant elastograms be computed from all three displacement components to obtain the best results.
Researchers have demonstrated that longitudinal waves make it difficult to reconstruct shear modulus from a single displacement field (Honarvar et al 2013). The presence of longitudinal waves could be one reason the performance of reverberant elastograms varied with different components of displacements. The curl operator is typically used to minimize the impact of longitudinal waves. The bandpass filter employed in this study was equivalent to computing the vector curl of the complex wave field (Sinkus et al 2005). Studies performed on a constrained gelatin phantom and the brains of healthy volunteers demonstrated that reverberant shear wave fields could be generated using standard elastographic imaging equipment (i.e. a single mechanical driver). Besides mechanical sources, other factors, such as the position of the external reflectors, size of the mechanical sources, shear wave attenuation, vibration amplitude, etc, will dictate the degree of reverberance induced within tissues. The impact of these variables on reverberant fields induced in soft tissue is beyond the scope of this work.
Nevertheless, intuitively, we expect to use fewer mechanical sources to generate a fully reverberant shear wave field in constrained organs such as the brain compared to unconstrained organs. Using multiple shakers can direct waves in different directions. Furthermore, the vibration from these sources can be out of phase to reduce the formation of standing waves. Therefore, we plan to investigate the clinical performance of reverberant brain elastograms obtained using a mechanical actuation system like that described in (Anderson et al 2016, Li et al 2021, which employs multiple mechanical sources. A limitation of this study is that we did not perform any imaging on an unconstrained phantom or organ. To address this, we plan to compare the performance of the shear modulus estimation method on both constrained and unconstrained phantoms. Although figures 6 and 7 demonstrated that the reverberant elastography method produces clinically useful elastograms, the spatial resolution of acquisition may be another factor that may affect the shear modulus. In this study, we used a 2.5 mm isotropic voxel size. Curtis et al (Johnson et al 2014). used 2 mm, and (Zhang et al 2011). used 3 mm in their studies. Reconstructions from only one component of motion displayed consistent shear modulus values (see table 4). We are currently conducting studies on a larger cohort of patients to evaluate the clinical performance of the reverberant elastograms; the results we will report in a future publication. Although the reconstructions from reverberant elastography have a lower spatial resolution, the reconstruction times are much faster. The shorter reconstruction times and accurate results make the technique well-suited for clinical studies in the brain.

Conclusion
This study corroborated the hypothesis that reverberant elastography produces reliable shear modulus elastograms of constrained organs, such as the brain, with a single mechanical driver. Studies performed on a brain-shape phantom demonstrated that the accuracy of reverberant elastograms computed using all three displacement components was comparable to those calculated using the subzone inversion method. The clinical study results were consistent with those of the phantom study (i.e. the performance of reverberant and subzone elastograms was comparable) and sufficiently encouraging enough to warrant further evaluation with a larger cohort of subjects.