An analytical solution to the dispersion‐by‐inversion problem in magnetic resonance elastography

Magnetic resonance elastography (MRE) measures stiffness of soft tissues by analyzing their spatial harmonic response to externally induced shear vibrations. Many MRE methods use inversion‐based reconstruction approaches, which invoke first‐ or second‐order derivatives by finite difference operators (first‐ and second‐FDOs) and thus give rise to a biased frequency dispersion of stiffness estimates.


| INTRODUCTION
Magnetic resonance elastography (MRE) measures viscoelastic shear properties of biological soft tissues in vivo using externally induced vibrations at different frequencies to capture the tissue's local harmonic response by motion-sensitive magnetic resonance imaging (MRI). MRE has been proven to be sensitive to many diseases including liver fibrosis, 1 liver tumors, 2-4 renal dysfunction, 5 and neuroinflammation. 6 MRE generates parameter maps from time-harmonic shear wavefields using inverse problem solutions. 7 Inverse methods in MRE published in the literature differ by their treatment of noise, boundary conditions, and tissue properties including heterogeneity, compressibility, anisotropy, and viscosity. [8][9][10][11][12][13] Nonetheless, most MRE inversion methods are similar in that they yield parameter maps with a spatial resolution above the diffraction limit, thus providing super-resolution. To achieve this, fine features of the curvature of shear waves are analyzed, often by applying finite-difference operators (FDOs) to the wave images. 14 FDOs can retrieve the spatial gradient of the wave image with pixel-wise resolution. 15 Inversion techniques that solve the wave equation in MRE (Helmholtz equation) rely on second-order FDOs (second-FDOs), or the Laplacian, whereas other techniques rely on the phase gradient of propagating plane waves and use only first-order FDOs (first-FDOs). First-FDO methods provide wave numbers, which are inversely proportional to shear wave speed c. 16 By contrast, second-FDO inversion (Helmholtz inversion or direct inversion) yields shear modulus, which is related to c 2 . 7 Despite their different physical units, shear modulus and shear wave speed are used in the literature as parameters of tissue stiffness at vibration frequency. 17 Henceforth, the term stiffness is used when discussing the qualitative changes in c and c 2 that occur with frequency (dispersion).
Tissue-intrinsic dispersion of stiffness is determined by the material's viscoelastic behavior. A purely elastic material has no dispersion, yielding constant stiffness values over frequency. However, it has been recognized that intrinsic dispersion recovered by FDO methods is subject to bias. [18][19][20] Specifically, we and others have shown that stiffness is likely to be overestimated at low spatial support due to the discretization of shear waves while it is often underestimated at high spatial support due to noise. 18,19 Because spatial support in MRE is the number of pixels per wavelength (wave number times pixel size or wavelength divided by pixel size), multifrequency MRE examinations automatically result in a variation of spatial support due to the occurrence of different wavelengths. Consequently, FDO-based inversion approaches can result in an artificial dispersion of stiffness over frequency, which is superposed to the intrinsic dispersion caused by tissue viscosity.
This well-known dispersion-by-inversion problem in MRE has never been analytically solved in a closed form for first-and second-order FDO approaches with consideration of both discretization and noise. Instead, earlier investigators proposed spatial support of 6 to 9 pixels per wavelength to avoid discretization artifacts 18,19 without taking into account SNR, kernel width (or type of the stencil 21 ), and order of FDOs. We here provide this information, validate our analytical solution of the dispersion-by-inversion problem by numerical simulations as well as phantom data, and derive a simple rule of thumb for minimizing FDO-related dispersion artifacts for future MRE examinations.

| THEORY
In this section, we first outline the analytical responses of FDOs to data containing either discretized harmonic waves or pure Gaussian noise. We then continue by combining both signals into master equations of FDO-related dispersion functions for defining the range of unbiased MRE inversions.

| Discretization
Finite differences can be approximated by Taylor series expansions of O (h n ) with n > 1 and h = N ⋅ d, where d is the pixel size and N ≥ 1 is the width of the kernel measured in number of pixels (eg, with N = 1, derivatives are computed across adjacent pixels without gap). O(h 2 )-expansions of central first-and second-order derivatives in one-dimensional image space (x) can be represented on the basis of Dirac distributions (x) 15 : 1 FDO and 2 FDO denote central first-and second-order finite difference kernels, respectively, which are typically used in MRE inversion techniques. 13,15 Central difference operators are often applied in MRE because other nonsymmetric stencils neglect information from one direction. Furthermore, O (h n ) operators with n = 2 have been proven relatively stable whereas higher order approximations of n > 2 result in unwanted oscillations in the presence of noise due to Runge's phenomenon. 22 Henceforth, we restrict our attention to 1 FDO and 2 FDO approximations given by the first terms on the right hand sides in Equations (1) and (2) (neglecting the remainders). These central difference FDOs are used in many MRE inversion techniques including multifrequency inversions, MDEV ( 2 FDO), 23 and k-MDEV ( 1 FDO). 16 (1)
Given infinite support, Equation (1) or (2) applied to harmonic waves of wavenumber k 0 provide their first or second spatial derivatives and thus result in phase-shifted waves, which are scaled by k 0 or k 2 0 , respectively. At finite support, the harmonic responses (k harmonic ) of 1 FDO and 2 FDO can be obtained by Fourier transformation of Equations (1) and (2) with respect to x, yielding: with A is the key parameter of the resolution of wave images in MRE. Similar to the Nyquist limit, it determines when aliasing of wave numbers occurs due to low spatial support and FDO kernel sizes. As seen in Figure 1, A < 2 (ie, dk 0 N > 0.5), this is the spatial support below 2 pixels per wavelength signifies the limit of ambiguous derivatives because Equations (3) and (4) are symmetric around A = 2.
As such, A determines the artifacts caused by discretization in FDO-based MRE inversion. The second important variable is SNR. Larger pixels and wider derivative kernels minimize the noise susceptibility of the apparent wave number k obtained by FDO. Henceforth, we account for the expectation values of FDO when applied to noise in order to find the optimal A. We assume that Z is the noisy part of a complex harmonic function (signal) with normalized amplitude leading to an SNR of 1/E (|Z|). 24 Thus, according to Equation (7), SNR reads The corresponding wavenumbers obtained by FDO that were applied to pure noise Z are given as

| Discretization and noise
Apparent wave numbers n k with n = 1,2 for first-and secondorder FDOs, respectively, are the Euclidian norm of harmonic signal and noise: Figure 1 illustrates the effect of overestimation of n k in noisy signals, which in turn leads to the underestimation of stiffness (c) at high support (high A). In Figure 2 The position of these intercepts on the d ⋅ k 0 -axis depends on N such that higher spatial support is required for larger FDO kernel sizes. Figure 2C shows master curves of normalized shear wave speed for first-and second-FDOs versus A, into which all curves from Figure 2A,B collapse. Inflection points where both master curves equal 1, ie, ically. Therefore, we insert the harmonic signals of Equations (3) or (4) and noise given by Equations (11) or (12) into Equation (13) that is then solved for the inflection point A 0 as follows: These equations cannot be solved directly for A 0 but are solved for SNR instead: Both functions, plotted in Figure 3 in the range of A 0 = 1, 2, … , 100, represent rather simple powerlaws in that range of A-values. Consequently, good approximations can be obtained by series expansion. Series expansion up to second and third orders lead to numerical approximations of SNR ≈  (16) and (17), respectively, as shown in Figure 3. The corresponding functions A 0 (SNR) are readily obtained by At a given SNR, these A-values indicate kernel sizes and pixel widths and provide values near ground truth in FDObased MRE inversion.

| Experimental setup
Multifrequency MRE was investigated using an elastic, homogenous phantom made available to the MRE study group  (18) and (19), respectively. These points correspond to the analytical solutions shown in Figure 3

| Experimental data evaluation
The experimentally acquired concentric wave patterns were unwrapped and Fourier-transformed in time to extract complex-valued wave images at driving frequency f as shown in Figure 4. These wave images were mapped onto cylindrical coordinates and averaged over the azimuthal angle, resulting in a single wave profile u * (r, f ) along the radial coordinate r (see Figure 4). The obtained one-dimensional profiles were used to demonstrate FDO inversion in experimental data. These kernels were applied to u * (r, f ) using Matlab's conv function. SNR was determined by taking the amplitude of the fundamental spatial frequency over the averaged amplitudes of all higher harmonic frequencies, which we treated as noise.
Ground truth shear wave speed c 0 was reconstructed without differential operators in order to avoid stiffness biases due to discretization and noise artifacts as theoretically analyzed in the previous section. Therefore, c 0 values were determined by fitting Bessel functions to u * (r, f ) assuming that the shear waves are polarized along the infinitely large cylinder axis [26][27][28] : J 0 denotes the Bessel function of first kind, u 0 is the wave amplitude, is the phase of the shear wave at sample radius r = R, and O * is a complex-valued offset term that accounts for compression waves. k * 0 , the ground truth complex-valued wave number, was the desired result of the fit and was converted to shear wave speed by c 0 = f ∕Re k * 0 .

| RESULTS
Numerical results of FDOs applied to simulations of complex oscillations u * = exp ik 0 r with added complex noise, F I G U R E 3 Approximated and full functions of signal-to-noise ratio (SNR) versus A 0 according to Equations (18)  which was generated using Matlab's randn function, are shown in Figure 2C. The analytical solutions given in Equation (13) agree with numerical simulations. For this master curve, second-FDO results are given as wave speed c∕c 0 instead of c∕c 0 2 in order to demonstrate the alignment of first-and second-FDO results at low support, that is, in the range of discretization bias (c∕c 0 > 1). Figure 2C also demonstrates the excellent prediction of inflection point A 0 by Equations (18) and (19), allowing the further use of these approximates for minimization of inversion biases in FDObased MRE reconstruction. Figure 5 shows ground truth shear wave speed values of the phantom over frequencies. The values obtained by analytical fits display no significant dispersion over frequency, as illustrated by the linear regression function of 0.0004 m/s/Hz slope and c 0 = 1.72 m/s (R 2 = 0.15, P value = .08). Figure 6 demonstrates the relative overestimation as well as undershot of c∕c 0 (first-FDO) and c∕c 0 2 (second-FDO) due to discretization and noise in experimental data. Colors encode different stencil widths from N = 1 to 5. Analytical curves were calculated by Equation (13) based on Equations (3) and (4) for k harmonic and Equations (11) and (12) for k noise using the experimental pixel size d = 1 mm and SNR values given in Table 1 Table 1 versus A. Analytical curves (black for second-FDOs, blue for first-FDOs) were calculated based on c 0 and signal-to-noise ratio (SNR), both given in Table 1. Variations in analytical curves are mainly due to variations in SNR over frequency. Vertical arrows indicate the solutions of A 0 according to Equations (18) and (19) taking the mean SNR of all frequencies from Table 1 variations in wave profiles u * (r, f ) between different vibration frequencies f. SNR values ranged from 2.6 at 150 Hz to 23.2 at 95 Hz with a mean of 8.9 ± 5.7. Experimental results are summarized in Table 1.

| DISCUSSION
This study adds to the ongoing quest for quantitative and reproducible MRE biomarkers that are consistent with ground truth values. 17 Ground truth in MRE is normally provided by shear oscillatory rheometry, which, however, has well-known limitations regarding its sensitivity to surface texture, plate pressure, and sample geometry. Moreover, MRE biomarkers involve inverse problem solutions, which are often ill posed due to noise and unknown boundary conditions, making MRE susceptible to inversion bias. 15 Many inversion methods in MRE methods use FDOs and are known to be biased by discretization and noise. 29 Nevertheless, FDOs are frequently used in MRE due to their simple implementation, computational power, and their robustness against local stiffness changes (heterogeneity) and variation in tissue geometry and morphology. As a result, competing over-and underestimates of FDO signals occur at the same time, giving rise to an inflection point (A 0 ) at which the stiffness dispersion curve coincides with ground truth values. This synchrony of overand undershooting effects is the reason why second-FDOs, which are highly sensitive to noise, require smaller A 0 values than first-FDOs. However, considering SNR to be typically in the range of 1 to 100, the relative difference in A 0 between first-and second-FDOs is relatively small (eg, A 0 = 9.2 versus 12.4 for SNR = 10 and N = 1). In practice, approximated values of A 0 rather than precise values should be taken into account. For instance, when A 0 is more than an integer factor larger than recommended values, N should be increased in order to decrease the noise sensitivity of the FDO used.
It should be noted that many MRE inversion methods treat noise using k-space filters or wavelets prior to FDOs of constant kernel widths 7,14 or apply Savitsky-Golay filters in order to combine second derivatives with denoising by polynomial regression. 7 Effectively, these strategies can provide the same smoothing as an adaptation of N but naturally vary in technical details and are thus difficult to analyze with an explicit estimation of the error as we do here. Therefore, we focused, in theory, on unsmoothed raw wave data u * (r, f ) assuming white Gaussian noise and analyzed the effect of noise suppression by FDO kernel size. Clearly, noise, as defined as unwanted signal in MRE, 15 can have different sources including compression waves 30 or slice jittering 31 making noise anisotropic and non-Gaussian.  Furthermore, SNR given by in this study is normally addressed by other methods such as octahedral shear strain, 24 shear strain invariant noise, 32 or wavelet-based methods, 29 all of which provide values different from those obtained by Equation (10). Our method of normalizing the fundamental vibration signal with higher harmonic noise agrees remarkably well with Gaussian SNR and, thus, also with the predicted dispersion functions. However, it should be mentioned that this strategy is only feasible in homogenous materials of single spatial frequency while limited in vivo. Our phantom data are not meant to reproduce in vivo soft tissue properties but to present ground truth wave speed c 0 values without frequency dependence. A drawback of the absence of viscosity is that undamped waves resulted in multidirectional k-vectors, which we treated by directional filtering as proposed in, 33 although classical Helmholtz inversion (second-order FDO inversion) does not require such directional filtering. Finally, our study was restricted to one-dimensional analysis. Again, the aim of this study was to develop the basic layout of an inversion-bydispersion analysis in a mathematically consistent manner.
Higher-dimensional FDO methods can be used to efficiently suppress noise 34 or to address stiffness heterogeneities by employing finite volumes. 35 In summary, we have introduced an analytical approach to the well-known dispersion-by-inversion bias in MRE. Using closed-form solutions, numerical simulations, and phantom experiments, we demonstrate that FDO-based stiffness estimates can be affected by simultaneous under-and overestimation of values due to the presence of noise and wave signal discretization, respectively. Our results further show that second-order FDOs are more sensitive to noise and are more markedly biased by discretization than firstorder FDOs. A lumped parameter A was introduced, which defines spatial support of the shear waves including stencil width of FDOs. Analytical FDO dispersion functions were analytically derived, and simple rules of thumb are provided, which might help to minimize FDO-related inversion artifacts in future examinations. This study adds to research on inversion artifacts in MRE and may contribute to more reliable viscoelasticity reconstruction in future multifrequency MRE examinations.