Several groups have investigated the use of magnetic resonance elastography (MRE) in neurological diseases including brain tumors (1), normal pressure hydrocephalus (2), and Alzheimer’s disease (3). In MRE, mechanically induced motion in tissue is estimated from a time-encoded series of phase contrast images (4). For steady-state MRE, relevant motion information is encoded within the first temporal harmonic of this image series. Given an estimate of this quantity, quantitative spatial maps of tissue stiffness can be constructed (5). Recently, Trzasko et al. (6) described a statistical signal processing framework that optimally estimates this quantity, and robust graph cuts-based optimization procedure for solving this problem. Although tissue stiffness map accuracy is fundamentally dependent on the quality of the temporal harmonic estimate from which it is derived, the practical advantage of the advanced harmonic estimation tool for brain MRE is still unknown. The goal of this work is to quantitatively investigate the improvements in stiffness estimation accuracy provided by the robust estimation tool in a previously-described skull phantom (7), and to test its potential benefit for brain MRE in vivo.In a standard steady-state MRE acquisition, the signal observed by receiver channel $$$c\in[0,C)$$$ during phase offset $$$t\in[0,T)$$$ at spatial position $$$x$$$ (out of $$$N$$$ pixels) for a single motion-encoding direction can be modeled as $$$G[x,t,c] = M[x,c]\text{exp}(j\text{Re}\{\eta[x]\zeta[t]\}) + Z[x,t,c]$$$, where $$$M$$$ is the complex-valued background signal, $$$\eta$$$ is the (complex) harmonic quantity, $$$\zeta[t] = \text{exp}(j2\pi{t}/T)$$$, and $$$Z[x,t,c]$$$ is complex Gaussian noise. This signal can be expressed in ($$$NT\times{C}$$$) matrix form as $$$G=H(\eta)M+Z$$$, where $$$[H(\eta)]_{(x,t),(x,t)} = \text{exp}(j\text{Re}\{\eta[x]\zeta[t]\})$$$ is a block-diagonal operator. Robust harmonic estimation comprises identifying $$$\eta$$$ which minimizes a cost functional constructed of a spatial smoothness penalty and the marginalized (negative) log-likelihood for $$$G$$$, i.e., $$$J(\eta) = \lambda P(\eta) - \text{trace}\left\{(H(\eta)[H(\eta)]^{\dagger})G\Lambda^{-1}G^{*}\right\} $$$, where $$$\Lambda = \mathbb{E}[Z^{*}Z]$$$. As shown in [6], this problem can be solved via iterative $$$\alpha$$$-expansion (i.e., graph cuts) [8], which provides robustness against phase wrapping.

To test stiffness accuracy, MRE was performed on a geometrically-accurate skull phantom with dynamical mechanical analysis (DMA) tested magnitude complex shear modulus (|G*|) of 4.0±0.3 kPa. MRE acquisition used a modified spin-echo echo planar imaging (EPI) sequence, a vibration frequency of 40Hz; TR/TE=4000/82ms; FOV=24 cm; 72x72 image matrix; 48 contiguous 3-mm-thick axial slices; two 40Hz motion-encoding gradient pairs for all 3 motion-encoding directions; and $$$T=4$$$ phase offsets spaced evenly over one vibration period. To test the efficacy of this technique in vivo, this experiment was repeated in a patient with a highly vascular well-delineated acoustic neuroma after IRB approval and written consent. RHE was implemented in C++/OpenMP using the max-flow library from (9). This procedure was run for 250 iterations using $$$\lambda=1e6$$$ (phantom), $$$\lambda=1e7$$$ (patient), and $$$P(\eta)$$$ defined with 6 cardinal neighbors. $$$\alpha$$$-expansions were performed with initial $$$|\Delta_{i}|=\sqrt{2}\pi$$$, and decimated in magnitude every 20 iterations; $$$\angle\{\Delta_{i}\}$$$ was changed every iteration according to the following schedule: $$$\{\pi/4,3\pi/4,5\pi/4,7\pi/4\}$$$. On a dual 3.0 GHz quad-core server with 32 GB memory, our code takes ~20 min to complete 3D estimations for all three motion directions. Quantitative stiffness map estimation was performed by direct inversion of the Helmholtz equation using the curl of the estimated 3D harmonic displacement fields (10).

The first 2 rows of Fig. 1 and 2 show the curled wave fields and stiffness (|G*|) maps calculated using the standard and robust harmonic estimation (RHE) methods in the geometrically-accurate brain phantom and in vivo, respectively. The bottom row in figure 1 shows a line profile (red dotted line in the curl images) for both the phase difference (red curve) and the RHE (blue curve) data sets, along with a corresponding magnitude MR image. The bottom row of figure 2 shows a magnified view of the tumor and adjacent tissue. Observe that the waves in RHE – compared to the standard result – are more sinusoidal (as would be expected), and that the stiffness estimate becomes more uniform. The line profiles also show apparent improvement in wave information in thinner anatomical regions of the phantom. Fig. 3 shows a map of relative error between the DMA and MRE –based stiffness estimates – again note that RHE enables provision of more accurate quantitative stiffness information overall.In this work, we adapted and applied a graph cut-based harmonic estimation method for brain MRE, and demonstrated the performance benefits of this tool both in a geometrically-accurate skull phantom and in vivo. The phantom data suggests that this technique improves the accuracy of brain MRE, and the in vivo data demonstrates that better tumor delineation is achieved when RHE processing is used.