Elastography Validity Criteria Definition Using Numerical Simulations and MR Acquisitions on a Low-Cost Structured Phantom

MR Elastography is a novel technique enabling the quantification of mechanical properties in tissue with MRI. It relies on a three-step process that includes the generation of a mechanical vibration, motion capture using dedicated MR sequences, and data processing involving inversion algorithms. If not properly tuned to the targeted application, each of those steps may impact the final outcome, potentially causing diagnostic errors and thus eventually treatment mismanagement. Different approaches exist that account for acquisition or reconstruction errors, but simple tools and metrics for quality control shared by both developers and end-users are still missing. In this context, our goal is to provide an easily deployable workflow that uses generic validity criteria to assess the performance of a given MRE protocol, leveraging numerical simulations with an accessible experimental setup. Numerical simulations are used to help both determining sets of relevant acquisition parameters and assessing the data processing's robustness. Simple validity criteria were defined, and the overall pipeline was tested in a custom-built, structured phantom made of silicone-based material. The latter have the advantage of being inexpensive, easy to handle, facilitate the fabrication of complex structures which geometry resembles the anatomical structures of interest, and are longitudinally stable. In this work, we successfully tested and evaluated the overall performances of our entire MR Elastography pipeline using easy-to-implement and accessible tools that could ultimately translate in MRE standardized and cost-effective procedures.


INTRODUCTION
Magnetic Resonance Elastography (MRE) non-invasively assesses and quantifies mechanical properties of tissue in vivo [1,2]. Various MRE techniques have been developed, involving the combination of interdependent processes (mechanical vibration generation, motion-sensitized MR sequences, image processing algorithms) all of which may impact the qualitative and quantitative outcomes. As a result, MRE, despite being available for decades, has not entered clinical routine application yet; a consensus on precision and accuracy of this technique still needs to be defined. Indeed, one of the major drawbacks of MRE is the lack of a gold-standard that could enable to compare elastograms with each other. Rheometry may provide a reference measurement of bulk mechanical properties of various materials [3][4][5][6][7], but the precision of common rheometers or dynamic mechanical analysis (DMA) systems depend on multiple factors, like the type of instrument used or the sample's size and experimental environment [6,8,9]. Furthermore, performing reliable rheometry over a range of frequencies comparable with those employed in MRE may not be possible with common rheometers and specific highfrequency measurement systems may be required [8,9]. In addition, biomechanical properties of ex vivo tissue specimens that lack perfusion are very different than those of perfused in vivo body parts, making a direct comparison between classical rheometers and MRE difficult. In this context, biases and errors intrinsic to the technique should be assessed and controlled to avoid the generation (and further the interpretation) of unreliable elastograms. Since MRE relies on the encoding of the tissue displacement into the phase of the complex MR signal, metrics that take into account phase errors are good candidates to measure raw data quality. Phase-to-noise ratio (PNR) can be used to evaluate the impact of noise on the reconstructed displacement. Given that the phase standard deviation in radians is equal to the inverse of the magnitude SNR [10], PNR is proportional to both the SNR and the encoded phase shift. PNR can be computed in numerous ways, either by including the specific details of a given experiment [11] or independently of them, but at the expense of extra acquisitions [12]. Displacementto-noise ratio (DNR) has been similarly used as a metric, either by exploiting the knowledge of the encoding sequence details and SNR [13], or by estimating a motion SNR from the displacement noise [14,15]. The octahedral shear strain (OSS) and the derived OSS-SNR have been introduced as indicators of reconstruction accuracy showing better performance than the motion SNR [15]. However, OSS is obtained from spatial derivatives of the motion field, which makes it dependent on the processing [16]. The ratio of local shear wavelength λ to pixel size a has also been introduced as an indicator of reconstruction accuracy [17,18], and values were reported for different methods [19,20]. From these studies, the range of optimal λ/a appears to be broad (∼5- 20), although best performances are generally achieved in the range 6.7-10, especially at low SNR.
Numerical simulations can be used as a validation tool by generating displacement fields (forward problem) for given boundary conditions, visco-elasticity distributions, and vibration sources. The displacement fields are then fed into the MRE reconstruction pipeline, which solves the inverse problem and retrieves the mechanical properties for comparison with groundtruth values. Simulations are generally based on analytical formulations [5,14,17,21] or finite element methods (FEM), either with custom-developed tools [22][23][24] or dedicated commercial softwares [20,[24][25][26][27][28][29][30]. Simulations are very useful to assess the error linked to reconstruction algorithms, however they do not suffice to reflect all potential limitations of an MRE experiment (transducer, B 0 and B 1 inhomogeneities, SNR, motion sensitivity, susceptibility issues, and heterogeneity of the material at very small scales). Validation is thus typically obtained from self-made or commercial phantoms mimicking biological tissues, which shape and properties depend on a targeted application. A thorough overview of materials used for ultrasound elastography is given in [31], and here we will focus on those commonly used in MRE. Water-based phantoms are largely spread and include natural or synthetic polymers. Natural polymers comprise gelatin-based phantoms, easy to make, that are characterized by a linearly elastic behavior, low viscosity and frequency dependency [32]. However, they are sensitive to temperature changes [33] and also to the fabrication process [27,34,35]. Agar can be used alone [1,2,5,20,36,37] or added to gelatin [38] with a multi-step heating and cooling procedure at precise temperature steps. Over a period of about 1 year post fabrication, such phantoms have shown a general decrease of storage modulus following an initial increase phase (1-5 months). Synthetic polymers like polyvinyl alcohol cryogels (PVA) are obtained by performing repeated cycles of freezing and thawing of boiled and stirred PVA solutions. Their mechanical properties depend both on the number of cycles and the PVA concentration [14,[39][40][41][42]. This fabrication process enables phantoms with anisotropic mechanical properties [21], yet inhomogeneous temperature changes within the entire volume cause inhomogeneities leading to uneven inclusion-background interfaces [39,43]. Commercially available phantoms (CIRS Inc., USA), made from patented Zerdine hydrogel, are more expensive than water-based phantoms [44]. They contain well-defined inclusions but are essentially non-viscous and the measured values often differ from those declared by the manufacturer [19,20,[45][46][47][48]. Polyvinyl chloride (PVC, also known as plastisol) phantoms are fabricated by heating-up a mixture of PVC and a softening agent to high temperatures (160-177 • C) while constantly stirring before pouring into a heat-resistant mold [8,49]. Arunachalam and colleagues have assumed that a 90-day period was sufficient to reach stable mechanical properties [9]. They also reported weak MR signal and very low loss moduli of the material, causing increased reflections. Finally, silicone has been used for Digital Image Elasto-Tomography (DIET) [50,51], and for mechanical imaging using pressure sensors [52], but little so far for MRE [53,54]. Compared to the other materials described above, silicone simplifies considerably the phantom fabrication process. Typically, silicones are created by mixing two components, and the final stiffness can be controlled by the addition of a softener or hardener. Besides, many silicones are room temperature vulcanizing (RTV), which eases the fabrication process and enables the realization of homogenous structures with smooth interfaces.
Our objective is to establish simple and easily transferable validity criteria, setting the ground for basic procedures in MRE quality control and efficient deployment of the technique in the clinics. In the presented work, we assess the limitations of MRE for a given protocol, comparing reconstruction performance from both fully simulated and acquired datasets in a low-cost silicone-based phantom.

Sample Preparation and Phantom Fabrication
Different samples were obtained by mixing the two components of a silicone rubber (Eurosil 4 A+B, Schouten SynTec, the Netherlands) together with an oil softener (Eurosil Softener, Schouten SynTec, the Netherlands). An empirical density of 1,000 kg/m 3 yielded accurate volume estimation of the final material. The silicone components were stirred manually (1-5 min depending on total volume) in a plastic container until a homogeneous mixture was obtained. Different stiffness properties were obtained by fixing the weight ratio between the silicone bi-components (A and B) and varying the softener composition, expressed as 1A : 1B : Sx. For rheological measurements, silicone mixtures were poured into small cylindrical plastic containers (diameter 25 mm) to create seven samples of about 6 mm height, with softener compositions of S1.0, S1.5, S1.7, S1.9, S2.1, S2.3, and S2.5. All samples were made on the same day.
An MRE phantom was made that consisted of four cylindrical inclusions of various diameter in a homogeneous, rectangular background (Figure 1). The silicone mixture for the background (1A : 1B : S1.9) was first poured into a 3D-printed mold (Fortus 450 m, Stratasys, USA) with 4 removable 3D printed cylinders ( Figure 1A). As the 3D print may produce rough surfaces, which promotes the formation of micro air bubbles while pouring uncured silicone, a plastic sheet was placed to cover the inner surface of the mold ( Figure 1B). Furthermore, an aerosol release agent (Pol-Ease 2500, Schouten SynTec, Netherlands) was applied on the surface of the cylinders to facilitate the extraction of the phantom from the mold. We allowed the silicone to fully cure for 1 day before extracting the cylindrical inserts, and filled the remaining volume with a softer mixture 1A : 1B : S2.1 to create the inclusions. The final size of the phantom was 182 × 108 × 90 mm 3 with inclusions of diameter 10, 20, 30, and 40 mm referred to as I1, I2, I3, I4 according to their size. All samples were stored away from direct light at room temperature, in a closed plastic container together with a two-way humidity control bag ("Humidipak, " D'Addario, USA) to ensure 45-50% humidity.

Silicone Stability
Stability over time of the prepared seven cylindrical silicone samples was investigated using a rotational rheometer (MCR series, 302 and 702, Anton Paar, Austria) to extract their complex shear modulus G * = G' + iG" (G': storage modulus, G": loss modulus). For each measurement, the upper rotating plate (flat model PP25) which diameter is identical to the samples, was positioned with an automatic set point of 5-N normal force to enable a sufficient compression and rotation. A sweep of the shear strain amplitude in the range 0.01-10% was performed at 10 Hz to verify the linear behavior of each material. Subsequently, a shear strain amplitude of 1% was used for G * measurement with a linear frequency sweep from 100 to 1 Hz. For each sample, two consecutive measurements of the frequency sweep were acquired. The same measurements were repeated at several time points, respectively at 1, 2, 4, 7, 15, and 22 weeks post fabrication.
An example of the frequency sweep measurement obtained with the silicone samples is given in Supplementary Figure 1. An increased standard deviation for the two consecutive measurements can be observed for frequencies >20 Hz (Supplementary Figure 1A). Similarly, larger G" than G' are observed above 35-40 Hz, even though the material kept its solid structure. Therefore, we assumed that the measurements corresponding to 10 Hz provided reliable information and we used those values to evaluate the sample stability over time.

Forward Problem
Numerical simulations were performed using COMSOL Multiphysics 5.4 (COMSOL AB, Sweden). The geometry of the simulated phantom was defined as a rectangular solid with dimensions identical to the MRE silicone phantom described above, including the same four cylindrical inclusions ( Figure 1C). The vibration area was defined as a circular region of 16 mm diameter on the front face of the phantom, at half distance between I3 and I4, similarly to the physical vibration setup used later.
For the simulation's physics definition, the Solid Mechanics interface from the Structural Mechanics module was used. The front face of the phantom was set as a "Free" boundary, while all the other faces were assigned to a "Fixed Constraint" node to reproduce the condition of the containing mold. A "Prescribed Displacement" node was assigned to the vibration area with an imposed sinusoidal displacement perpendicular to the surface with a 200-µm amplitude and frequency f vib of 57 Hz. For both types of regions (background and inclusions) the material properties were defined using a "Linear Elastic Material" domain node with a "Viscoelasticity" domain, without any particular contact definition between the two materials.
In order to define G * = G' + iG" of the different regions, we employed an isotropic Standard Linear Solid model (see Supplementary Material). The simulation consisted in a "Time Dependent" study over a time interval of 12 vibration periods T = 1/f vib , sufficient for the waves to propagate in the whole volume, with a time step equal to T/20, for a total simulation time of 41 min (Intel Xeon CPU E5-2680 v4 @ 2.4 GHz, 14 cores, RAM 128 GB). The exported simulated data consisted in the displacement field components u, v, w (in X, Y, and Z direction, respectively) at time points tau = i * T/N (with N the number of wave phases sampled during MRE acquisitions, five in this case) of the last vibration period. A spatial grid of 9 contiguous axial slices was defined for three different isotropic voxel sizes of 2.5, 3.0, and 3.5 mm 3 , respectively. The central slice (i.e. slice 5) was set at 22 mm from the surface ( Figure 1C).

Synthetic MR Data
Displacement field components were converted into synthetic complex MRE data via a custom MATLAB (R2020a, The MathWorks Inc., USA) script. Phase-encoding was artificially achieved by multiplying the displacement data by the fractional encoding factor (Equation 1 from [12]) of the MRE encoding scheme, with a motion-encoding-gradient (MEG) of 54 mT/m: with ϕ the vibration-encoded phase, γ the gyromagnetic ratio, u and f vib the amplitude and frequency of the vibration at a given coordinate in space, A MEG and f MEG the amplitude and frequency of the sinusoidal motion encoding gradient MEG (t). Magnitude of the synthetic data was arbitrarily set to one everywhere. Since a reference dataset without MEG is required for our reconstruction, a separate dataset was generated with identical geometry, magnitude set to 1 and phase to 0 everywhere (no additional noise was considered).

MR Imaging
MRE was performed on a clinical 3 T Magnetom Prisma scanner (Siemens Healthineers, Germany) using a fractional encoding gradient-echo sequence [12] and a 20-channel head coil. To constrain the bulk motion, the phantom was kept inside its plastic

MRE Reconstruction
Both the synthetic and experimental raw MRE data were processed using a custom software performing a 3D curl-based direct inversion [12]. All the elastograms were then exported for further processing using MATLAB scripts. In this paper we focus only on G', the MRE output information most commonly used by numerous research groups.
A mask was drawn manually around the phantom boundaries before unwrapping the phase data. Since the reconstruction software requires to calculate multiple spatial derivatives, the generated elastograms are by default cropped by two voxels in 3D with respect to the initially defined volume of interest (VOI). In addition, we chose to discard two more pixels for each slice in order to take into account possible calculation error propagation that often translates in G' underestimation at boundaries. Elastograms of the five central slices were used for the analysis of synthetic data. Alternatively, elastograms of the three central slices were used with experimental data as more discrepancy was observed on surrounding slices.

Regions of Interest (ROIs)
For synthetic data, maps of the G parameter were exported at the same coordinates as the displacement field data, in order to obtain a map of the inclusion positions on the same spatial interpolation grid. This way, ROIs were automatically available for the location of the inclusions inside the background. For the experimental data, ROIs were defined as circles centered in the inclusions, segmented from T 2 maps. The respective masks were co-registered with the MRE datasets. To avoid partial volume effects due to the larger voxel sizes employed in MRE, the mask for the background was eroded by two voxels surrounding the inclusions. Mean, median, and standard deviation values of G' were calculated within the different ROIs. All the segmentation steps were based on data quality considerations detailed below.

MRE Validity Criteria
In order to evaluate the quality and reliability of raw MRE data, simple metrics were chosen, easily transposable to any MRE method. In particular, SNR was calculated voxel-wise for each set of raw data and converted into phase and displacement amplitude uncertainty maps using Equation 1. The mean SNR and standard deviation of the reference dataset (acquired without motion encoding gradients) were used to discard voxels with low SNR and therefore large phase error on final MRE maps. This threshold was empirically set to 2 standard deviations below the mean SNR: voxels on the elastograms were discarded when at least one dataset of the three encoded directions did not have a sufficient SNR. Additionally, due to derivative calculations on the reconstructed data, we discarded two more voxels on each side of the voxels excluded by the SNR threshold by dilating the SNR mask in all three directions.
To calculate the ratio of shear wavelength λ to voxel size a for the three resolutions, λ was estimated via Equation 2 [56]: where f vib is the vibration frequency of 57 Hz, ρ the density of the material (assumed to be 1,000 kg/m 3 ), and G' is the estimated storage modulus. For G' values, we considered the median values obtained with the 3-mm resolution within the background (silicone softener composition S1.9) and within the largest inclusion I4 (silicone softener composition S2.1).

MRE Validation
Reconstruction algorithms can have an impact both on the calculated stiffness and the apparent objects geometry. To validate the MRE pipeline, the contrast-to-noise ratio CNR i between the i-th inclusion and the background [19,20,57,58] was calculated using Equation 3: where m i , m B are the mean values of G' within the i-th inclusion and the background, and σ i , σ B are the corresponding standard deviations of G'. Here, the median was used instead of the mean to allow for non-Gaussian distributions. Root-mean-square-error RMSE was calculated for the voxels in the different regions according to Equation 4: where G ′ 0,j is the storage modulus defined in the simulation and G ′ j is the reconstructed value in a voxel j. For MRE data, G ′ 0,j in the background was taken as the median value at 3 mm resolution, whereas G ′ 0,j values in the inclusions were estimated as G ′ 0,j in the background divided by the ratio G ′ S1.9 /G ′ S2.1 measured via rheometry.

Statistical Analyses
For each considered resolution, the probability distribution was tested using the Kolmogorov-Smirnov method. Differences between G' values in the inclusions and those in the background were assessed with the Kruskal-Wallis test (α = 0.05), followed by a Mann-Whitney test for the two-by-two comparisons between the four inclusions and the background. The significance level was adjusted for multiple comparisons using the Bonferroni correction, and a p-value of 0.0125 was considered as significant.

Simulations
Analysis of the reconstructed elastograms for the simulated phantom is shown in Figure 3. As described above, the reconstruction does not provide any value within a 2-voxel rim around the object. The inclusions are circled in red and the sections (A-A and B-B) show G' values across the five central slices (Figure 3A).
Despite having identical compositions, large inclusions, compared to small ones, are better discriminated from the background even though they all showed significant differences (p < 0.0125) at all resolutions (see boxplots in Figure 3B). This trend is confirmed by the CNR values, greater in the inclusions I4 and I3 than in I2 and I1. The reconstruction software yields more accurate values in the background and in I4 and I3 than in I2 and I1, close to the ground-truth defined in the simulation (Figure 4). In fact, the error is similar in the background and in the inclusions I4 and I3, while it increases more than 2-fold for the smaller inclusions I2 and I1. In general, Figures 3, 4 indicate that a smaller voxel size corresponds to a greater standard deviation of G' within each region and to a larger error in the background.
The wavelength-to-voxel size ratios λ/a are calculated for each spatial resolution and summarized in Table 2. This ratio decreases with increasing voxel size, with all values included between 10 and 15.7.

Phantom MRI and MRE
The employed phantom fabrication technique produced cylindrical shapes with sharp contours and overall proper adherence between the different silicone regions at the inclusions boundaries (cf. Figure 1). The values of both T 1 and T 2 of the two silicone compositions differ by <10 ms, with minor variations over time post fabrication (Supplementary Figure 2). The inclusions shape and position were identified on T 2 maps, which presented less B 1 inhomogeneity at the center of the phantom than T 1 maps (Supplementary Figure 3), and were thus used for ROI segmentation.
The SNR and the error values both on the MR phase and the corresponding displacement amplitude are reported in Table 3. These errors were rather small and did not contribute in the overall data thresholding. An example of the reconstructed G' maps is represented in Figure 5 for the central slice at week 4 post fabrication. The number of voxels near the edges of the phantom with underestimated values of G' is between 2 and 5 voxels greater than what is observed in synthetic data. In general, the voxels discarded by our SNR threshold are located near the boundaries of the phantom and at its center, in correspondence with partial volume effects and the location of the vibration source.   Table 2). Stars: contrast-to-noise ratio CNR between each inclusion and the background.
As observed with synthetic data, the standard deviation within a given region decreases with increasing voxel size, and the values within the inclusions stand out of the background more clearly while inclusion size increases. This effect led to a higher CNR for the largest inclusions although CNR values in the phantom were smaller than those obtained from simulated data. As in synthetic data, G' in all the inclusions and for all the resolutions was significantly different from the background's (p < 0.0125). As expected, the background silicone with a lower softener content of S1.9 has greater storage modulus values than the inclusions made with a softener composition of S2.1. The mean values in the background are however lower than their median, because of the skewed G' distribution and the contribution of underestimated voxels near the phantom's edges.
The RMSE values for G' measured via MRE in the different regions for each considered resolution are illustrated in Figure 7. The error is smallest for the inclusion I4, and largest for I1, with all values greater than observed for synthetic data. As for simulations, we observed an increased error in the background All values are expressed as mean ± standard deviation.   for the finest resolution, as opposed to the inclusions where the errors increased with increasing voxel size. The evolution of G' values over time (Figure 8) shows a higher initial variation than in the case of the rheometry samples, with stability reached at around 4-weeks post fabrication. Furthermore, ratios of G' between the two silicone compositions S1.9 and S2.1 measured via MRE are similar to those measured via rotational rheometry. The wavelength-to-voxel size ratios λ/a estimated from MRE data are summarized in Table 4. This ratio increases over time for all the regions following the initial stiffening of the silicone material, with λ/a values ranging from 9.7 and 17.1.

Silicone Based Phantoms for MRE
The proposed silicone-based preparation shows great promise as a low-cost, easy-to-reproduce, and durable option for quality control and calibration in MRE experiments.
Rheometry monitoring of the complex shear modulus showed good long-term stability after a 2-to 4-week delay, across all prepared samples. During this initial phase, an overall initial FIGURE 9 | Diagram summarizing our pipeline for MRE validation using simulations (yellow), a low-cost silicone phantom (pink), MRE experiments (blue) and the relationships between different steps (arrows). Steps are represented as dashed frames, processes as large arrows, and outputs as filled boxes.
stiffening of maximum 1.6 kPa was observed. Remarkably, despite the exceeding of the maximum recommended softener composition from the manufacturer, our results showed that G' could be changed and remained stable over time. From this, we obtained a range of realistic stiffnesses covering G' values from 2 to 10 kPa, similar to that observed in biological tissues.
The fabricated MRE silicone phantom exhibited similar evolution of G' over time, while maintaining its G' ratios (background/inclusions) (Figure 8). Yet, the overall variations in G' measured with MRE (∼1.3 kPa) were generally larger than G' obtained with rheometry (∼0.7 kPa). This may be caused by variations in the aging process due to their size differences. Furthermore, such differences may also be expected since rheometry is not a gold standard and operates in a different frequency regime compared to MRE. Nevertheless, our results indicate that rheometry is a good candidate to provide an estimation of the expected values measured via MRE and consists in a good basis for fine-tuning future phantom preparations.
Here, we used two silicone compositions and a simple geometry for the purpose of the study. With appropriate modifications to our fabrication procedure, it should be possible to create more complex, multilayer geometries or anatomicallyshaped features in an easy way by exploiting the good shape conservation and interface contact between silicone mixtures once cured. A potential difficulty in the fabrication process may arise from the formation of micro air bubbles. Gradually pouring TABLE 4 | Measured G' in the largest inclusion (I4, softener composition S2.1) and background (BG, softener composition S1.9) at 3-mm resolution for different time elapsed since fabrication and corresponding estimated ratio of wavelength to voxel size λ/a for the considered voxel sizes a. the silicone over a plastic sheet covering the internal faces of the mold helped mitigating this problem, yet some susceptibility artifacts due to air bubbles at the edges of the inclusions (see Figure 5 in the proximity of I3) remained. To solve this issue, one can envision either using materials with smooth and solid surfaces for all the molding parts, such as PMMA or silicone molds with anti-adherent coating [59], scanners with lower magnetic field strengths, or a vacuum chamber for degassing during the fabrication process. Viscosity is another important aspect of the silicone material, and in our case it was greater than what is reported for other commonly used MRE phantoms. As biological tissue, our samples exhibit a dispersive behavior (Supplementary Figures 1A,B). The subsequent wave attenuation has the advantage to prevent successive reflections at boundaries and hence the presence of standing waves especially in small objects but may prevent from obtaining sufficient wave propagation inside large objects. Finally, our results showed T 1 and T 2 values in the order of 1,030 ms and 240 ms, respectively, for both compartments of the MRE phantom (differences never exceeded 7.2 ms for both T 1 and T 2 ), and with little variation over time (overall <0.9 and 2.9% in relative terms). Further studies may consider tuning the silicone relaxation times by the addition of paramagnetic substances compatible with the silicone mixture during fabrication. This would help further the sequence optimization for a faster transfer to in vivo application, in particular when other markers than stiffness are required for segmentation for example [59].

Simulations
We showed that synthetic MRE data produced via numerical simulations can be used with our reconstruction pipeline. The obtained elastograms displayed distinct (and expected) differences in G' values between the defined inclusions and the background. By varying voxel size, we could investigate the resolvability of the latter inclusions as a function of their size, as well as the accuracy of the reconstruction algorithm. It was possible to accurately retrieve G' values within large homogeneous regions while observing some limitations for inclusions smaller than 3-cm diameter. For larger voxel size and/or for a smaller inclusion diameter, our results showed a decreased standard deviation of G' as well as an overestimation of G' within the inclusions. Both findings may be partly explained by the local homogeneity assumption made in the inversion algorithm, and by the calculation of the derivatives within a local neighborhood. With large voxel sizes or small inclusions, more voxels are affected by the discontinuity in the distribution of the mechanical properties at the inclusion boundaries. The higher dispersion, especially visible in the background region when using small voxel sizes, can be explained by a sub-optimal λ/a ratio. The theoretical λ/a for the 2.5-mm resolution belongs to the higher end of the suggested optimal interval of 5-20 found in the literature. Furthermore, for this same resolution (data not shown), G' maps exhibit inhomogeneities that match the wave propagation pattern and ultimately lead to larger RMSE in the background. Considering that our simulations do not add any noise and ensure that no standing waves occur, one can conclude that such effects are exclusively due to the reconstruction pipeline and that large λ/a ratios should be avoided. All together, these results suggest that numerical simulations can provide a reliable prediction of reconstructed results for different experimental conditions in a simple, yet heterogeneous object. Possible challenges in designing accurate simulations of more complex geometries lie in the difficulties of realistic modeling, in particular of the boundary conditions and wave propagation at the various interfaces. Although simulations of a body organ may require a detailed knowledge of mechanical properties and their distribution at a fine spatial scale, their feasibility may be increased from the availability of structural datasets acquired via different imaging techniques (MRI, CT) in combination with appropriate, simplified physical assumptions. Finally, we believe that simulations could help improving the overall MRE performances by tuning either the MRE parameters, or by addressing specific aspects of the reconstruction pipeline. In this study we chose one single frequency and three spatial resolutions, but simulations can easily be performed using various frequencies and a large span of spatial resolutions.

MRE Experiments
Overall, experimental results appear very similar to the reconstructed simulation data (Figures 3, 5, 6). Some differences can be seen on the edges of the phantom exhibiting uneven contours and susceptibility artifacts at the air-silicone interface. Like the simulations, MRE demonstrated good accuracy for retrieving G' across large, homogeneous regions (background, I4, and I3), and limitations for inclusions smaller than 3 cm diameter. Those limitations are tied to the specific MRE protocol chosen in this study. We chose a vibration frequency of 57 Hz in the range commonly used in human studies for various diseases, both focal and diffuse (e.g. in the brain, breast, liver, muscle, prostate, heart), along with acquisition times that ensure patient comfort. Yet, if required, mapping smaller objects with confidence could be achieved by improving the inversion algorithm, in particular by assuming anisotropy and medium heterogeneity. Further development in the experimental protocol could also help in that direction and provide more accurate MRE results, typically by adapting the mechanical vibration settings along with the acquisition parameters or using a retrospective data resampling [60]. We also noticed biases introduced by the reconstruction process, such as the underestimation of G' values near the boundaries of the phantom (Figure 5) or the greater error in the background at the finest resolution (Figure 7). That said, confidence in the estimation of errors from the MRE side must be balanced with the fact that RMSE calculation requires an actual ground truth, here estimated from the simulation data.

Validity Criteria for MRE
The λ/a ratio was always in the range suggested by the literature, gradually increasing with time from the initial stiffening phase in our phantom. Interestingly, we note ( Table 4) that a relatively small change in voxel size (e.g. 3 mm vs. 2.5 mm) affects the ratio λ/a more than an increase in G' of almost 1.2 kPa, as could be measured from week 1 to week 4. Therefore, if decimation or interpolation of voxels are always possible after the acquisition, an appropriate native resolution remains essential for relevant reconstruction outcomes. More concretely, it shows that the quest for finer resolutions does reach a limit where it starts to impede MRE performance rather than improving it.
In our study, we used a fixed vibration frequency along with a variable resolution conversely to most studies that investigated the λ/a ratio by employing a wide range of vibration frequencies and a fixed pixel size. The range of chosen resolutions may not be sufficient to define the optimal λ/a ratio but is sufficient to show how much it impacts reconstruction, without being influenced by wave penetration or motion efficiency. Although this approach may lead to smoothing effects from filtering kernels in the processing of the raw MRE data, it is more relevant for mono-frequency MRE approaches and in organs where wave attenuation may vary significantly in a small range of frequencies.
Considering SNR, the computed values over time were high ( Table 3) and translated into rather small phase and amplitude errors, defeating their relevance to discard unreliable voxels. In a more general case, however, phase uncertainty may be used to eliminate regions where the encoding of the vibration produces a net phase lower than the corresponding error. This may occur because of low SNR, either due to short T * 2 in general, a strong vibration, highly attenuated wave amplitude or insufficient encoding efficiency. It is particularly important to be able to distinguish the two latter conditions since a lack of waves may be due to physical properties or transducer settings, whereas the encoding efficiency may be improved acting on the sequence parameters. The motion encoding efficiency, however, is not always available to the end-user, and providing this information in terms of encoded radians per displacement at the scanner user interface would allow useful data quality considerations and an estimation of PNR. Using SNR only as a thresholding metric, on the other hand, was successfully applied. In our study, it appeared particularly useful in regions where large motion amplitudes led to strong losses of MR signal, i.e. in voxels near the transducer location. Such a metric is simple and has the advantage of not requiring knowledge of the sequence details. We strongly recommend to systematically calculate the SNR and to use it as a basic quality control check to prevent any potential diagnostic errors that will in turn lead to inadequate patient treatment.

CONCLUSIONS
The proposed work presents an effective and adaptable solution for testing, assessing, or comparing MRE methods performances, which is summarized in Figure 9. Synthetic MRE was successfully obtained from simulated displacement fields. We demonstrated the fabrication of a simple, low-cost, structured MRE phantom with controlled mechanical properties close to biological tissue that are stable over time, replicating the groundtruth phantom modeled in our simulations. Furthermore, the material is quite dispersive and hence suitable for studying multifrequency MRE. Such a phantom can be easily adapted to various applications and deployed worldwide to study the performance of a given MRE experiment and reconstruction pipeline. Our simulations successfully matched experimental MRE results on the phantom. Similar results, trends, and limitations regarding the accuracy of both G' values and the inclusion shape were observed. As a consequence, simulations can be used to predict or estimate performance in MRE experiments in silico and showcase the possible intrinsic method limitations. This step could be cost and time effective for parameters fine tuning in anticipation of the optimization phases at the scanner. Finally, we showed how the interpretation of the reconstructed maps can be improved by discarding unreliable voxels. In particular, besides phase uncertainty based criteria which may be useful in certain cases and the importance of λ/a ratio, an SNR-based discrimination should be applied in general for direct thresholding to discard unreliable elastogram regions.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.