Magnetic resonance elastography in nonlinear viscoelastic materials under load

Characterisation of soft tissue mechanical properties is a topic of increasing interest in translational and clinical research. Magnetic resonance elastography (MRE) has been used in this context to assess the mechanical properties of tissues in vivo noninvasively. Typically, these analyses rely on linear viscoelastic wave equations to assess material properties from measured wave dynamics. However, deformations that occur in some tissues (e.g. liver during respiration, heart during the cardiac cycle, or external compression during a breast exam) can yield loading bias, complicating the interpretation of tissue stiffness from MRE measurements. In this paper, it is shown how combined knowledge of a material’s rheology and loading state can be used to eliminate loading bias and enable interpretation of intrinsic (unloaded) stiffness properties. Equations are derived utilising perturbation theory and Cauchy’s equations of motion to demonstrate the impact of loading state on periodic steady-state wave behaviour in nonlinear viscoelastic materials. These equations demonstrate how loading bias yields apparent material stiffening, softening and anisotropy. MRE sensitivity to deformation is demonstrated in an experimental phantom, showing a loading bias of up to twofold. From an unbiased stiffness of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4910.4 \pm 635.8$$\end{document}4910.4±635.8 Pa in unloaded state, the biased stiffness increases to 9767.5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm \,$$\end{document}±1949.9 Pa under a load of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx $$\end{document}≈ 34% uniaxial compression. Integrating knowledge of phantom loading and rheology into a novel MRE reconstruction, it is shown that it is possible to characterise intrinsic material characteristics, eliminating the loading bias from MRE data. The framework introduced and demonstrated in phantoms illustrates a pathway that can be translated and applied to MRE in complex deforming tissues. This would contribute to a better assessment of material properties in soft tissues employing elastography.


Introduction
Diagnosis and therapy planning often depend on developing an accurate understanding of the state of a patient's disease. Complex alterations occur in diseased tissue, depending on a range of factors such as changes to protein isoforms, alterations in protein density, extracellular matrix reorganisation, inflammation, etc. These modifications fundamentally impact tissue mechanical properties, triggering changes in stiffness, elasticity and viscosity (Yeh et al. 2002;Yin et al. 2009;Baiocchini et al. 2016). The ability to detect changes in tissue mechanical properties in vivo has the potential to significantly benefit clinical medicine, enabling a more thorough assessment of pathology as well as monitoring of treatment progression. This potential has lead to the development of a number of approaches for the noninvasive assessment of tissue mechanics.
Manual palpation is a classical approach of inspecting tissue for changes in elasticity in accessible organs (Weiss 2003, Nguyen et al. 2014. Advancements in medical imaging have led to the development of elastography-a modern quantitative technique for assessing the mechanical properties of internal organs. Elastography imaging is usually achieved using Ultrasound or Magnetic Resonance Imaging (MRI) (see Carey and Carey 2010 for a review), enabling the determination of tissue properties in vivo. Elastography imaging focuses on observing waves as they traverse through tissues, thus allowing quantification of strain response due to an excitation. In the case of magnetic resonance elastography (MRE), low-amplitude harmonic waves (10-1000, 30-80 Hz in vivo) are produced on the body and subsequently imaged in three-dimensions (Manduca et al. 2001, Fovargue et al. 2018b. Harmonic wave motion is then related to local tissue properties, in most cases, through linear viscoelastic wave equations (Glaser et al. 2012). MRE has been successfully employed for a range of tissues-including brain (Klatt et al. 2007;Green et al. 2008;Schregel et al. 2012), liver (Huwart et al. 2006, Klatt et al. 2007) and breast (Sinkus et al. 2000, 2005a, Houten et al. 2003showing contrast between diseased and healthy tissue as well as providing resolution of apparent tissue anisotropy ). More recently, it has also been applied to heart (Elgeti et al. 2008, Kolipaka et al. 2009, Robert et al. 2009, Kolipaka et al. 2010, Couade et al. 2011, Kolipaka et al. 2012, Elgeti and Sack 2014, skeletal muscles (Green et al. 2012, and other soft tissues, reviewed in Manduca et al. (2001), Mariappan et al. (2010) and Fovargue et al. (2018b).
While interpreting stiffness measures in organs, such as the liver, has shown diagnostic value, a challenge yet to be addressed is the impact of macroscale deformations on wave propagation behaviour. Broadly, by looking at the 1D wave equation ∂ 2 u/∂t 2 = c 2 ∂ 2 u/∂ x 2 , term c = √ T /ρ depends on tension T and linear density ρ and is scaling the spatial second derivative of displacement u. The resultant force term on the right hand side balances the acceleration term on the left hand side, leading to the phase velocity increasing with higher tension. Elastic, hyperelastic and viscoelastic nonlinear materials would exhibit this effect. This is relevant in the context of tissues that undergo natural deformations (e.g: the heart muscle during the cardiac cycle (Asner et al. 2017), the liver during respiration (Kang et al. 2012), etc.) and tissues that undergo imposed deformations (e.g. breast during MRE scans Sinkus et al. 2005aSinkus et al. , 2007. As many tissues are known to exhibit significant nonlinearity (Nash and Hunter 2000;Taber 2004;Liu et al. 2006;Holzapfel and Ogden 2009;Gao et al. 2010;Nordsletten et al. 2011b;Goenezen et al. 2012) under physiological conditions, both natural and imposed deformation can result in loading bias-an apparent stiffness that is a conglomer-ate measure of the intrinsic material properties and current kinematic state. For instance, the nonlinearity arising from compressive strain was quantified in ex vivo bovine liver samples using MRE measurements (Clarke et al. 2011). Hence, mixtures of patient-specific kinematics and intrinsic tissue properties are non-trivial and presumably result in a reduction in clinical specificity.
In this paper, the harmonic perturbation of Cauchy's equations of motion for a general incompressible nonlinear viscoelastic material is examined, in order to understand the links between intrinsic properties and kinematics in MRE. From this analysis, a set of governing equations are derived, illustrating the general impact of large deformation on the motion of harmonic waves. Importantly, these equations enable the elimination of loading bias and therefore retrieval of intrinsic properties when large-scale deformation and constitutive behaviour are known. Theoretical test cases are presented, showing the ability of deformation to yield apparent increase, decrease and anisotropy in materials. In order to demonstrate the applicability of these perturbed equations, experiments were performed in polyvinyl alcohol (PVA) phantoms. Rheological tests were performed to characterise material behaviour, and loaded phantoms were imaged using MRE. Results show the capacity of these equations to correct for the complex dynamics of wave motion in deformed materials. To the authors' knowledge, this is the first cross-validation study linking mechanical characterisation of a material in rheology with wave motion through the deformed material in MRE. The developments presented provide a key step towards improving the estimation of material properties in soft tissues, using elastography.
As a starting point, Sect. 2.1 presents perturbation theory employed within Cauchy's equations of motion. The perturbed equations are analysed in the context of a general viscoelastic material under load, further leading to the derivation of apparent stiffness moduli influenced by large deformations and material constitutive law. For exemplification purposes, an idealised case of harmonic wave motion traversing a uniformly compressed Neo-Hookean material is presented. Experiments are employed to test the theoretical ground developed. In Sect. 2.2, PVA phantoms are tested in a rheological setup, to determine the material governing law. In Sect. 2.3, MRE data are acquired in PVA phantoms at different uniaxial compression levels. The knowledge on rheological behaviour and deformation is integrated into MRE data analysis, in order to undo the loading bias. The experimental results are presented in Sect. 3. Particularly, the estimation of the intrinsic stiffness of PVA is presented in Sect. 3.3. The theoretical framework presented in this study could be extended to in vivo MRE data analysis.

Harmonic wave motion in deformed nonlinear viscoelastic materials
MRE relies on the analysis of small amplitude harmonic waves traversing through materials. As materials undergo deformation, wave behaviour is impacted (Thurston and Brugger 1964). One option to understand the relationship between large-scale deformation and small-scale wave behaviour is perturbation analysis. After covering notation 2.

Kinematics, harmonic wave motion and notation
Here, we briefly review the basic kinematics in large deformation mechanics in order to later describe the macroscale deformation of phantom and tissue material subjected to MRE. For more complete dispositions, see for example (Malvern 1969;Wang and Truesdell 1973;Graff 1991;Holzapfel 2000;Bonet and Wood 2008). The motion of a solid body, denoted by the region Ω 0 ⊂ R d , is characterised by the displacement field U : Ω 0 × [0, T ] → R d . In this case, any reference point X ∈ Ω 0 may be deformed to its physical position by the diffeomorphic Lagrangian mapping L : Ω 0 ×[0, T ] → R d so that its position at a time t ∈ [0, T ] is given by Due to the Lagrangian mapping, note that any mdimensional function v : Ω 0 ×[0, T ] → R m can be written as a function of the deformed domain Ω(t), e.g. v : Ω → R m , with the understanding that v(x, t) = v(L −1 (x, t), t). For ease of notation, the · is dropped where there is clear distinction.
The stresses resulting from material deformation can be separated into deviatoric and hydrostatic components. For the latter, the effect is quantified using the hydrostatic pressure P, which ensures some constraint on the volume or pressure-volume relation. The deviatoric components are, in turn, related to shape changes (or strains and strain history).
Important quantities for mechanics are the deformation gradient tensor F and its determinant J : Here, ∇ 0 refers to the gradient operator taken with respect to coordinates of Ω 0 , so that (∇ 0 ) k = ∂/∂ X k . For later use, When relating kinematics and kinetics, the right or left Cauchy Green strains are often considered. These may also be defined in terms of their isochoric variants (Ogden 1997;Bonet and Wood 2008) aŝ which are invariant to changes in volumetric properties of the material, i.e.Ĉ(a F) =Ĉ(F) for a real positive scalar a ∈ R + (Holzapfel 2000). A set of invariants are associated with these quantities. Here, the first and second invariants of any second-order tensor A ∈ R m×m are defined as in Bonet and Wood (2008) For later clarity, the double contraction ":" is used as defined in Bonet and Wood (2008), e.g. the contraction of secondorder square tensors A, B ∈ R m×m becomes a scalar 1 and the contraction of a fourth-order tensor A ∈ R m×m×m×m and a second-order tensor B yields a second-order tensor in R m×m 2 For further use throughout the paper, the multiplication of second-and fourth-order tensors AA or AB is defined to take place over the last index of the first quantity and first index of the second quantity, i.e.
Finally, the derivative of a second-order tensor with respect to another is defined as MRE relies on the propagation of shear waves through the body, assumed to reach a harmonic steady state. Being wave motion, these deformations are assumed to be linear and small. The large deformation U is perturbed by these small wave deformations where the scale of the perturbations is much smaller than the dominant length scale of the body or the large displacement it undergoes. When the perturbations are periodic-in-time and have reached a steady-state, the deformations and pressures can be written as where give the complex periodic steady-state behaviour with real and imaginary parts given by subscripts r and i, respectively.

Perturbation of Cauchy's equations of motion
For an incompressible body, at any time t ∈ [0, T ], motion of Ω 0 satisfies Cauchy's first law (conservation of momentum) and conservation of mass shown in Eqs. 3a and 3b, respectively (see Malvern 1969;Ogden 1997;Holzapfel 2000;Bonet and Wood 2008): Here, ρ is the material density, ∂ tt is the second time derivative, and S e = S(P, C) is the second Piola Kirchhoff (PK2) stress tensor depending on the hydrostatic pressure, P, and the right Cauchy Green tensor, C. The PK2 stress can be characterised as elastic, hyperelastic, viscoelastic, etc., depending on the constitutive relation that describes the material. A hyperelastic material, for example, can be written as a function of C (or other strain metrics) and space (X). In contrast, a viscoelastic material additionally depends on time (t). We recall briefly that, for hyperelastic materials-such as those commonly applied in biomechanics-S e is given as the derivative of a governing hyperelastic strain energy function, W = W e (C) (Bonet and Wood 2008), with respect to the right Cauchy Green tensor C, i.e.
This model is extensively applied in nonlinear biomechanics-for example in simulations of the heart (Wang et al. 2009, Nordsletten et al. 2011a, McCormick et al. 2013, 2014, Hadjicharalambous et al. 2014, breast (Rajagopal et al. 2010, Reynolds et al. 2011, Gamage et al. 2011, arterial wall (Holzapfel 2000, Gasser et al. 2006, Hariton et al. 2007, etc.-where specific material response is usually defined through biorheological experiments. The simplest material models are purely elastic, based on spring rheological elements. Complementary, pure viscosity is based solely on dashpot elements. Biological tissues usually exhibit viscoelasticity; hence, a suitable tissue model should consider a combination of these springs and dashpots. For instance, the Maxwell model is formed by considering a spring and dashpot connected in series, whereas the Kelvin-Voigt model is formed by considering a spring and dashpot connected in parallel. These types of models have been analysed in literature (Liu and Bilston 2000, Bilston et al. 2001, Banks et al. 2011), but have shortcomings in predicting creep (Maxwell model) or stress relaxation (Kelvin-Voigt model) (Liu et al. 2006, Banks et al. 2011. While the Maxwell constitutive equation was found to be better suited for modelling fluids (Houten et al. 2000;Sinkus et al. 2005a), the Kelvin-Voigt constitutive equation was employed for purposes similar to ours, modelling soft tissues subjected to elastography testing (Sinkus et al. 2005a, Huwart et al. 2006, Huwart et al. 2008) (albeit not considering the effects of large deformations on the wave behaviour). However, the Kelvin-Voigt model is not able to accurately capture the power-law dependence on frequency that is usually observed in tissues (Chui et al. 2004, Nicolle et al. 2010, Nicolle 2015. Thus, an adaptation of this model is commonly used, which replaces the dashpot with a springpot. This form was observed to describe tissue better (Kiss et al. 2004).
To capture the viscoelastic behaviour exhibited generally by tissues and polymers, we proceed by assuming that the PK2 stress can be decomposed into an additive sum of elastic, viscoelastic and hydrostatic components, i.e.
The elastic part S e is defined through a hyperelastic strain energy function W e (C) as in Eq. 4, the hydrostatic part as S p = J PC −1 , and the viscoelastic part is defined using the Caputo formulation (Caputo 1967) of the fractional-order derivative: Here, S v depends on a viscoelastic strain energy function such as S v = 2∇ C W v . Considering the coefficient of derivation α to be 0, the viscoelastic stress S v acts as a hyperelastic term, while α = 1 indicates the first-order derivative ∂ t S v .
Intermediate values of α would lead to non-trivial transitional results between S v and ∂ t S v . Therefore, this formulation of the viscoelastic term leads to a fractional nonlinear Kelvin-Voigt type of model and the PK2 tensor for our viscoelastic material becomes In order to investigate the effect of deformation on wave behaviour, we consider the case where high-frequency/lowamplitude harmonic waves (see Eq. 1) are imposed onto our body that satisfies the set of Eq. 3. For our purposes, we assume that the introduction of these micro-deformations u ε does not disrupt the natural state of macro-deformation U, such that the total observed state (U ε , P ε ) can be characterised by with (U, P) satisfying the original Eq. 3 (without perturbation conditions) and with (U ε , P ε ) satisfying a perturbed version of Eq. 3, i.e.
where superscript ε indicates quantities dependent on the perturbed state variables. A simplification of the set of Eq. 8 can be achieved by expanding about the state variable U and then linearising with respect to the small perturbations u ε , p ε (Thurston and Brugger 1964). For completeness, details about the expansion and linearisation processes can be found in Appendix 1. Considering the case of harmonic waves, a set of periodic nonlinear viscoelastic wave equations can be derived in the reference domain Ω 0 . By recasting the reference frame gradient and divergence operators into their physical domain counterparts, the set of Eq. 8 can be transformed into the physical domain as: where ∇ = ∂/∂ x on Ω. Here, G and G are the real and imaginary stiffness moduli influencing the apparent wave dynamics of u c and taking the form of fourth-order tensors: where F, S e and S v depend on the unperturbed macrodeformation U and hydrostatic pressure and P. In the absence of deformation, under the assumption of an isotropic material defined by a standard law (e.g. a Neo-Hookean model for S e and S v ), these stiffness moduli are simply given by G + iG = G * I, where G * is the true complex shear modulus characterising the material and I i jkl = δ ik δ jl . The wave dynamics become symmetric, since I : ∇u c = ∇u c +∇u T c . However, under deformation, the components of ∇u c are scaled in a non-trivial way. By understanding the form of the scaling influencing the G and G moduli, the bias introduced by deformation on the wave behaviour can be predicted. Furthermore, the deformation can be undone, that is, the moduli that reflect the undeformed state can be recovered.

Loading bias of planar shear waves in pure compression
For an intuitive illustration, let us consider the case of a simple isotropic Neo-Hookean material described by the strain energy function where μ e is the elastic material parameter. The elastic stress is derived to be The viscous term is assumed to be zero, i.e. S v = 0. Under these assumptions, the real modulus G becomes To better understand the changes introduced by large deformations U on the wave dynamics, we analyse the case of plane waves through a material under pure compression. Thus, considering a cuboid that is compressed along the height (e 3 direction), imposing incompressibility, the deformation U is described by where λ is the ratio of the final height over the initial height (i.e. λ < 1). In this case, the hydrostatic pressure is a constant given by P = μ e (1/λ − λ 2 )/3 and the kinematic tensors describing this deformation are with J = 1. A pure shear wave which will perturb the macrodeformation U can be written as This form ensures that ∇u c has only off-diagonal components, i.e. ∇u c : I = 0, thus satisfying the incompressibility assumption of the material. Since the deformation considered employs no shearing (i.e. no off-diagonal components), then also B : ∇u c = 0 and therefore this wave form allows for simplifications when being contracted, as needed in Eq. 9a, with the stiffness modulus by removing the terms containing δ ml andB ml from Eq. 11: Additionally, under the divergence operator, ∇u T c vanishes (together with IB, P and J , which are constant in space), and thus, the elastic contribution can be further simplified to Let us now analyse particular examples of shear wavesnamely planar waves-travelling through the compressed cuboid. In the first case, consider that the front face of the cuboid is vibrated along the e 1 direction (see Fig. 1), an idealised wave of the form is created, with the only nonzero entry in the wave displacement gradient being (∇u c ) 12 . Therefore, the deformation contribution in Eq. 14 simply becomes a scaling by 1/λ: As such, the loading bias is influencing the apparent stiffness of the cuboid, which appears to be μ e /λ, whereas the intrinsic stiffness is μ e . Thus, under this setup of the wave, the material stiffness appears to be higher, since λ < 1.
In the second case, consider a planar wave created by vibrating the top face of the cuboid, this time along the e 2 direction (see Fig. 1). Then, the idealised wave displacements are probing the cuboid in the direction in which it was originally compressed, since the only nonzero component of the wave gradient is (∇u c ) 23 . This results in a stiffness scaling by λ 2 , i.e.
With λ < 1, the biased stiffness λ 2 μ e appears to be softer than the material's true stiffness μ e . From analysing these two simple instances of planar waves (Eqs. 15 and 16), we gain an understanding of how the waves probing a material under compression experience different loading bias-stiffer in the expanded direction e 1 and softer in the compressed direction e 3 . Had the waves been more complex, like in Eq. 13, then the effect of the loading would yield Thus, even under simplified conditions like an idealised shear wave travelling through a Neo-Hookean material, we get a feeling for how biased the wave propagation becomes under pure compression. It can be presumed that, under less than ideal conditions (more complex materials and deformations), the harmonic wave motion becomes increasingly intricate.

Nonlinear viscoelastic characterisation of PVA
Experiments done in controllable media, where shapes and deformations are simple, constitute the first step in evaluating the theory explained in the previous section. A first key component is using a material with a known rheological behaviour, as this governs the effective loading bias (e.g. Eq. 6). The following sections describe the fabrication process of PVA material and its testing in a rheological setup. A viscoelastic material law is also formulated. This law will then be later (Sect. 2.3.2) integrated with deformation and MRE data.

PVA phantom preparation
A material suitable for our purposes should be able to withstand large deformations without rupturing and to have low wave attenuation. For this, polyvinyl alcohol cryogel (PVA-C) was selected, which is suitable for mimicking soft tissue properties and has a high MR signal (Surry et al. 2004, Sinkus et al. 2005a. Following a protocol similar to the one described in Xia et al. (2011), phantoms were created by mixing polyvinyl alcohol powder (P1763 Sigma-Aldrich Company Ltd., UK) with deionised water in a concentration of 7%. A magnetic stirrer was used for 2 h, in order to fully dissolve the powder into the water, while heating the concoction to 90 • C. The mixture was covered throughout the process, to avoid evaporation, left to cool down at room temperature and poured into cuboid moulds of dimensions 64 × 48 × 42 mm. The moulds were then placed in a freezer at − 20 • C and underwent three cycles of freezing (14 h) and thawing (10 h) (F-T). These specifications ensured that the material was suitable for testing under large deformations. Laboratory experimentation with higher PVA concentrations (e.g. 10%) led to phantoms that were not easily deformable, whereas reducing the number of F-T cycles led to unstable phantoms that would leak water under loading. During the freezing process, volumetric expansion occurred. Therefore, the phantoms' top was cut, in order to achieve cuboid shapes. This led to differences in the phantoms' height, which ranged between 32 and 42 mm. For a meaningful variance testing, 14 homogeneous phantoms were created, out of which 7 were used for rheological experiments and 7 for MRE experiments. However, one phantom was excluded in the rheological tests due to rupture. The time between phantom fabrication and testing in either rheology or MRE was within two weeks, to avoid potential material degradation.

Rheological testing of PVA phantoms
The phantoms allow for suitable examination methods that can test the theory developed in Sect. 2.1. It was seen that, in order to account for the nonlinear complex stiffness moduli G and G (Eq. 10), understanding of the nonlinear behaviour and frequency dependence of the material must be inferred. This can be done through rheological testing, which is used to characterise material behaviour by quantifying its response to different types of applied stress. As such, each of the six phantoms created for rheology experiments was tested in a Bose Electroforce 5500 test instrument. The setup, seen in Fig. 2, consisted of a supporting fixed platen and an adjustable upper platen, which was movable in the vertical direction. A load cell of capacity 225N measured the force required to ensure a user defined deformation.
For our aims, three types of tests were employed, with the purpose of assessing PVA's nonlinear behaviour. The first test was designed to investigate the material's response to dynamic micro-oscillations around a state of large static load. This design replicates closely the in vivo MRE scenario, where low-amplitude harmonic waves are propagating in an organ which may be subjected to large deformation. For reproducing the nonlinear behaviour, four different static loading states were investigated. Thus, the moving platen was programmed to compress the phantom by 2, 4, 6 and 8 mm, respectively (∼5-20% uniaxial compression), and then to oscillate with an amplitude of 0.15 mm (∼0.4%). This test will be referred to as the micro-oscillatory test.
The second test was designed to investigate the nonlinear large deformation response in time. As such, large oscillations were imposed around a state of large deformation: the platen compressed the phantom by 5.8 mm (∼15% uniaxial compression) and then oscillations of 4 mm (∼10%) amplitude were imposed. This test will be referred to as the macro-oscillatory test.
In both tests, oscillatory tests were carried under a frequency sweep, from 0.1 Hz to 10 Hz (the instrument's limit): 0.1, 0.5, 1, 2, 5, 7.5 and 10 Hz. A clear depiction of the testing protocol can be seen in Fig. 2. This was done to investigate the frequency response of PVA, with the aim of extrapolating it to the higher frequencies used in MRE. In all tests, but clearly noticeable in the macro-oscillatory one, the displacements' amplitude is decreasing with increasing frequency. This is probably due to the instrument's limitation which, in the faster frequency regime, cannot reach the instructed displacements.
At the higher frequencies (≥ 5 Hz), it was observed that the moving platen gains acceleration and introduces a bias on the data. As such, calibration data were acquired, where the same protocol was followed, but with no sample in between the platens. Subtracting the calibration data from the phantom data eliminated the platen's momentum bias.
A final stress relaxation test was conducted by subjecting the material to a constant large deformation, in order to capture the specific viscoelastic effect. Hence, the phantoms were held under a compression of 11 mm (∼28.5%) for 5 min.
Overall, each phantom underwent six tests: four microoscillations at different loads (with seven different oscillatory states acquired), one macro-oscillation (with seven different oscillatory states acquired) and one relaxation. In between tests, the phantoms were left to rest in water for 15-20 min, since PVA is best stored in water (Surry et al. 2004), for hydration reasons. This resting time allowed for an efficient testing protocol of the phantoms, and its duration was sufficient to observe reproducibility.

Constitutive modelling of PVA
PVA and other hydrogels have been modelled, previously, as hyperelastic materials, employing laws that are generally used to describe rubber-like materials. Often, a Mooney-Rivlin type of model, employing two parameters, was found to be sufficient to capture the elastic behaviour of PVA (Anseth et al. 1996, Pazos et al. 2009). Extension to viscoelasticity was sometimes accomplished by considering Maxwellian elements (Anseth et al. 1996, King et al. 2011). As such, here the aim is to model the PVA material as viscoelastic, yet using springpot elements. Following previous investigations, we start the modelling process by considering the elastic part to be defined by a Mooney-Rivlin type of material law: where C 1 and C 2 (Pa) are material parameters. Considering the definition of the elastic part of the PK2 stress given, generally, as S e = 2∇ C W e , it follows that S e has two elastic parts scaled by the material parameters C 1 and C 2 : The elastic terms are derived to be As previously mentioned in Sect. 2.1.2, we consider the viscous part of PK2 to depend on a fractional-order derivative of a viscoelastic stress S v , which is defined here as the sum δ 1 S e1 + δ 2 S e2 . Hence, the PK2 tensor becomes Here, α is the fractional derivative coefficient, with α = 0 preserving the viscoelastic stress (D α t S e = S e ), which turns into a purely elastic contribution, and α = 1 indicating the first-order derivative, which becomes a purely viscous contribution. The hydrostatic stress S p = J PC −1 completes the definition of our PK2 tensor.

Rheology data analysis and model fitting
The viscoelastic model derived in the previous section (Eq. 19) needs to be tailored to the PVA material used in the rheological experiments by finding appropriate parameters C 1 , C 2 , δ 1 , δ 2 (Pa) and α.
As previously mentioned, the acquired rheological data consist of force measurements required to ensure a predefined displacement, each corresponding to a point in time. This was transformed into traction data in the e 3 direction by dividing the force readings by the cuboid's top area adjusted to the deformation, i.e.
Here, tr d (t) is the traction obtained from the data, F(t) is the force reading, and A(t) is the area at time t, which depends on the initial top area A 0 and the ratio of initial cuboid height Likewise, the traction in the e 3 direction arising from the model (tr m ) was computed using where σ is the Cauchy stress and n is the normal to the surface (in our case, we are interested in the normal to the top surface, i.e. n = [0, 0, 1] T ). Since the deformation at any point in time is given by the elastic and viscoelastic tensors S e1,2 and D α t S e1,2 can be determined, together with their counterparts in the Cauchy stress. They take the form of diagonal second-order tensors. The hydrostatic pressure P can also be determined at each time point by noticing that the traction in the e 1,2 directions is 0. Hence, P counteracts the elastic and viscoelastic contributions in both the e 1,2 directions and can be used as such in computing the traction in the e 3 direction. Thus, the only unknowns in the traction model are parameters C 1 , C 2 , δ 1 , δ 2 , α. In broad terms, since α is the only nonlinear parameter of the model, the parameterisation process is done by iterating over fixed values of α and then solving a minimisation problem which yields a best fit of the remaining linear parameters. It is desired that all tests (micro-, macro-oscillations and relaxation) carry the same importance in the fitting process, and that the error is not dominated by certain tests (e.g. those done at higher compression levels, which employ larger tractions). As such, for each of the six tests (four micro-oscillation, one macro-oscillation and one relaxation), the traction data and model were standardised according to the maximum traction value in the test, i.e.
where tr d (t) and tr m (t) are kept for simplicity, to avoid addition of further notation.
Three characteristics were considered to be important when fitting the model to the data: the compression, relaxation and oscillatory response. Thus, the error to be minimised was designed to comprise two parts, the first one dealing with the compression and relaxation behaviour, and the second one with the oscillatory behaviour. For the first part, in order to ensure that the error is not dominated by compression only, but also by relaxation, it is important to consider the data broken down into cycles, i.e. one full oscillation starting from the lowest load state, equivalent to one period (see Fig. 2). Each of these cycles is confined within a generic time interval [t 1 , t 2 ]. A mean traction value (t r d ,t r m for data and model, respectively) is defined as the integral of the traction over the cycle, e.g.
over generic cycle k. This mean value captures the compression level and, by considering locally every cycle in a test, captures the relaxation behaviour. Hence, the first part of the error deals with the sum, over all cycles, of the mean traction value difference between the data and the model, normalised by the data: The second part of the error is designed to measure the difference in the oscillatory behaviour between data and model. Thus, the traction amplitude corresponding to a generic time point t j in a cycle k is found as the difference between the traction value and the traction mean over the cycle, e.g.
The difference in amplitudes between the model and the data is investigated across all time points in a cycle, over all cycles, and is subsequently normalised by the data, as We note that the standardisation of the data was done over each of the six tests separately, to ensure that they carry equal importance within the error, while the normalisation was done over the tests altogether, to facilitate the interpretation of the error, which is finally outlined as err = err 1 + err 2 2 .
Having defined the error to be minimised, a more detailed parameter fitting process can be described. The nonlinear parameter α is iterated between 0.05 and 1 (close to the elastic and viscous limits, respectively), with an incremental step of 0.05. Since the PVA material used has a low viscoelastic response (and α = 0.05 was observed to be the most suitable), the incremental step was refined to be 0.01 between 0.01 and 0.1. For each fixed α, the linear parameters C 1 , C 2 , δ 1 , δ 2 were searched by solving the least squares problem corresponding to error Eq. 20, using the inbuilt function lsqnonlin of MATLAB and Statistics Toolbox Release 2015a, The MathWorks, Inc., Natick, Massachusetts, USA. The nonlinear solver was used in order to allow for a positivity constraint on the parameters. For each phantom, a full set of parameters was determined.
To ensure parameter uniqueness, once the best α value was determined, a new parameter search was performed. Thus, each individual linear parameter (C 1 , C 2 , δ 1 and δ 2 , respectively) was fit using the process described above, but while imposing the remaining three linear parameters to be 0 Pa. This was done in order to gauge the value of each model component and to identify potential parameter coupling.

Harmonic wave motion in deformed PVA
After characterising the PVA material through a viscoelastic model, the next essential element in estimating the intrinsic material parameters and bypassing the loading bias is merging the information on deformation and rheology with MRE data. Thus, in what follows, the MRE experimental features used to obtain simple deformations are described. The core component of integrating material behaviour and deformation knowledge into the data analysis process is outlined.

MRE experimental setup
When designing the MRE experiment, it was important to obtain shear waves that would propagate through the whole body of the phantom, while maintaining a fixed position of the phantom. The setup used is illustrated in Fig. 3. A transducer indenting the phantom transmits small amplitude compressional waves. Since the PVA material is aqueous, it could easily slide on a flat surface. For this reason, the phantom's back side was in contact with the support, ensuring no penetration, hence converting the compressional waves into A coil vibrating to the frequency established by the MRE sequence is connected to a flexible lamella, which causes an attached rod to vibrate longitudinally. A piston fit at the end of the rod is indenting the phantom, thus generating compressional waves. The phantom is resting on a smooth support and is in contact with the back support plate, which prevents the phantom from slipping and thus helps converting the compressional waves into shear waves. An upper plate compressing the phantom is kept in place by bolts fixed to the side plates. Reception coils are placed around the phantom, for signal enhancement. (Bottom) 2D top view of the setup shear waves. A large uniaxial deformation could be obtained by simply compressing the phantom with a smooth top plate, as seen in Fig. 3. Due to the sliding feature of the PVA material, the phantoms did not bulge under the uniaxial compression just described.
Using this setup inside the MR scanner, coronal wave data were recorded using MRE in undeformed and deformed configurations. Each phantom was scanned in the reference state (uncompressed), first deformation state consisting of uniaxial compression of ∼5 mm (∼14%) and second deformation state consisting of uniaxial compression of ∼12 mm (∼35%).
MR acquisition consisted of two primary scans. First, a geometric scan was acquired providing 55 cross-sectional slices with a field of view (FOV) of 96×96 mm (resolution of 96 × 96 pixels and slice thickness of 1 mm). Scan parameters were adjusted to enhance signal of the PVA material (TE/TR= 120.00 /5000.00 ms).
The second scan was a MRE eXpresso Gradient Echo sequence (Garteiser et al. 2013) providing 10 cross-sectional slices with a FOV of 96 × 96 mm (resolution 96 × 96 pixels and slice thickness of 2.0 mm). Driver frequency and motion encoding gradients were set, by turn, at 120, 130 and 140 Hz, and TE/TR at 6.91 /151.82 ms. The sequence encoded each of the three displacement component direction-e 1 , e 2 , e 3 and a reference non-motion encoded image removing background phase shifts due to the MRI gradients. Examples of the geometric and wave images can be seen in Fig. 4, for the undeformed and deformed configurations.

MRE data analysis
The geometric data acquired were used to provide detailed quantitative information of the geometrical configuration of the material in reference and deformed states. Assuming that the cuboid phantoms were compressed uniformly, then the deformation gradient can be written as in Eq. 12, where λ is the compression defined as the ratio of the compressed and uncompressed heights. Hence, due to the design of the MRE experiment, this form of the deformation gradient can be employed with the stiffness moduli given by Eq. 10.
From the recorded small amplitude displacements, two reconstruction methods were used. Firstly, a divergence free finite element reconstruction solving the set of Eq. 9 (Fovargue et al. 2018a) was used to retrieve the stiffness of the PVA phantom, merging the datasets acquired at 120, 130 and 140 Hz, for a better signal to noise ratio. This method takes as input the material density ρ, frequency ω and wave displace- The piston indentation can be seen, as well as the expansion of the phan-tom in the e 1 − e 2 directions under compression. Wave displacements from the MRE imaging protocol can be seen in the e 1 direction (third column), in the e 2 direction (fourth column) and in the e 3 direction (fifth column). It can be observed in the wave displacements, between the three rows, that the wavelength increases under compression. The wave images in the third, fourth and fifth columns have been cropped around the phantom area, to exclude the surrounding noise ments u c and reconstructs the complex stiffness modulus G * by assuming This is usually separated into the (real) storage modulus G and (imaginary) loss modulus G (where G * = G + i G ).
As explained in Sect. 2.1.2, the form of G * presented in Eq. 21 holds true in the absence of deformation and under an isotropic standard material law, and the wave gradient contribution becomes symmetric, as I : ∇u c = ∇u c +∇u T c . Hence, this method is suitable for reconstructing the stiffness modulus when no deformation is employed, but will lead to a biased result under deformation. Throughout the rest of the paper, this method will be referred to as uncorrected reconstruction (UR).
Secondly, the same reconstruction was adapted here, by integrating the formal definition of the real and imaginary stiffness moduli G and G (Eq. 10) in order to account for the bias in apparent stiffness introduced by the deformations. Conceptually, instead of relying solely on the wave displacements and its gradient ∇u c , the reconstruction integrates the scaling due to large deformations into the wave displacements gradient. Hence, instead of solving for G * , a modified form of the equations was developed, given below: where a solution is sought for M +i N. This second approach will be referred to as the corrected reconstruction (CR). As previously mentioned, within the UR the wave gradient is contracted with the fourth-order identity tensor scaled by the complex stiffness modulus, as (G * I) : ∇u c . By contrast, the CR allows the wave to be contracted with two distinct fourth-order tensors, defined by the user. One result is to be integrated with the real part and one with the imaginary part of the complex stiffness modulus, as M(G : ∇u c )+i N(G : ∇u c ), where M +i N is to be found. This is a novel way of analysing MRE data and is of utmost importance, since we saw, in Eq. 10, that the real and imaginary moduli G and G are non-trivial and change differently with deformation. These changes depend on metrics derived from the known deformation gradient F (B, I C , etc.), and, for our PVA material, are scaled by the parameters determining the viscoelastic Eq. 19. The specific form can be found in Appendix 2, which is applicable to the PVA material used here. Thus, assuming that the viscoelastic model describes the data well, then G and G already incorporate the true real and imaginary stiffness moduli. Hence, it is predicted that M and N will be reconstructed as unity.
With this understanding of the newly developed CR, the PVA material model was parameterised, this time employing the MRE data. The reference and second compression states were used, and the left hand side of Eq. 22a was minimised over all pixels. That is, the linear model parameters were sought such that M and N were reconstructed as closely as possible to unity for each pixel, over all pixels. The nonlinear parameter α was fixed to be the one indicated by rheology tests. The set of parameters thus obtained from the reference and second compression states was tested against the first compression state, for prediction value. The error was measured as the pixel-wise difference between the reconstructed values M and unity in uncompressed and second compression case, normalised by the number of pixels, as where n is the total number of pixels considered and subscript k iterates over those pixels.

Results and discussion
The first essential step in the experimental work was tailoring the viscoelastic law 19 to the PVA phantoms used in the rheology testing. The results of this process are presented in Sect. 3.1 and help us understand the rheological behaviour of the PVA material used, which directly influences the loading bias. The deformation bias is shown in Sect. 3.2, by investigating uniaxial compression in the MRE experimental setup. In Sect. 3.3 it is shown that, by incorporating the information on rheology and deformation into the CR, the intrinsic stiffness of the PVA material can be retrieved from the loaded states of the phantoms.

Nonlinear viscoelastic model for PVA
The aim of the modelling process was to fit the viscoelastic model 19 to the PVA data acquired in the rheological setup. Hence, for each of the six phantoms used in the rheology experiment, the parameters were fit accounting for all oscillatory and relaxation tests simultaneously. It will be seen that the parameters cannot be uniquely identified, which stems from the fact that the PVA has a small viscoelastic response, as it will be explained shortly. As such, a model with a reduced number of parameters is sought, which ensures unique parameter estimation. Initially, the 5-parameter model described by Eq. 19 was fit to the data, considering α values between 0 and 1. Figure 5 (left) depicts the errors for each of the six phantoms used in rheology. All phantoms yield a minimal error around an α value of 0.03-0.09. The minimal error increases only slightly at higher α values. The reason for the shallow changes in minimal error at high α values is that the 2 viscoelastic parameters δ 1 and δ 2 become very small; hence, the data are fit mostly with the elastic parameters C 1 and C 2 . Due to the PVA material having a small viscoelastic response, employing only an elastic model is not highly detrimental to the model fit, hence the small errors even at α = 1. The best parameter fit for each phantom is summarised in Table 1. The errors are small (less than ∼7%), which indicates that a good fit is ensured for each test. An example of the model fit can be seen in Fig. 6 (left). It can be observed that the compression level and relaxation behaviour are properly captured, while microoscillatory amplitudes are slightly underestimated. Notably,  Table 1 (5-parameter model) and Table 3 (3-parameter model), is enhanced  Table 1, that the phantoms for which a high δ 1 was yielded had a small C 1 and vice versa.
Due to the small viscoelastic response of the PVA phantoms, parameters {C 1 , δ 1 } and {C 2 , δ 2 } tend to exhibit very similar features. As such, in order to identify each parameter contribution to the model, α was set to 0.03 (the best choice according to Fig. 5 (right)) and each of the linear parameters was fit individually. The effects of each individual parameter on the model are illustrated in Fig. 7. When considering only C 1 , the model component considered is S e1 , which is equivalent to the Neo-Hookean elastic law. This parameter captures reasonably well the compression level of each test, indicating that the PVA material used has a dominant linearly elastic response. However, it lacks the ability to capture the changes in oscillatory amplitudes-in the four repetitions of the micro-oscillatory tests, the data indicates an increased amplitude at higher compression levels, whereas the model shows the same amplitude. The decrease in amplitude seen in the macro-oscillatory test is a direct result of the decrease in the actual displacement amplitude data, as seen in Fig. 7. This amplitude behaviour is entirely expected from the Neo-Hookean model, and the mismatch between (Right) The model fit to the data using a reduced number of parameters: C 2 , δ 1 , and α (with an error of 1.78%). The parameter values can be seen in Tables 1 and 3 Fig. 7 Contribution of each parameter to the model (black) compared to the data (red) (here, in phantom 4). Each of the four linear parameters (C 1 (top left), C 2 (top right), δ 1 (bottom left), δ 2 (bottom right)) was fit to the data, while setting the other 3 to be 0 (α was kept constant at 0.03). In each quadrant, the traction is depicted along the y-axis and the time along the x-axis (waiting times between tests not plotted for convenience). The parameters and errors can be seen in Table 2  Table 2 Best fit parameter for each PVA phantom; the error, computed as per Eq. 20, can be seen below each parameter; each parameter was fit considering the remaining three parameters to be 0 the model and data indicates that the PVA displays nonlinear characteristics. Lastly, no relaxation behaviour is displayed, which is also expected due to the model component being purely elastic. Conversely, the fractional derivative of the S e1 component, scaled by δ 1 , is the viscoelastic counterpart of the Neo-Hookean model. As such, investigating the value of δ 1 leads to similar observations as for C 1 , with the obvious difference that it can capture the relaxation behaviour. It is stressed, once again, that the similarity occurs due to α being very small. When investigating C 2 , the nonlinearity of the S e2 component becomes obvious-a linear increase in the compression level (2 mm between the micro-oscillatory tests) leads to a nonlinear response in the traction force. Hence, this parameter alone provides an unsuitable fit to the data; however, it adds value to the full model by being able to adjust for the increased amplitude in the four micro-oscillation tests. As before, the viscoelastic counterpart δ 2 captures the relaxation behaviour. Table 2 presents the best fit parameter for each PVA phantom. By investigating Fig. 7, it is expected that the C 1 and δ 1 parameters provide a better fit, and hence yield smaller errors than the C 2 and δ 2 parameters.
Having looked at each individual component, it is reasonable to conclude that two out of the four linear parametersone elastic and one viscoelastic non-counterpart-are sufficient for describing the PVA material behaviour. By keeping either C 1 or δ 1 , the compression level is ensured. By adding either δ 2 or C 2 , the amplitude fit can be improved for both the micro-and macro-oscillations. The relaxation is ensured by considering either viscoelastic component-δ 1 or δ 2 , although from the individual parameter fit it looks like δ 1 is more suitable for this aspect. These observations indicate that the best parameters to fit are C 2 and δ 1 , while setting C 1 and δ 2 to 0 Pa. Indeed, when investigating all pairwise tests (i.e. fitting either {C 1 , δ 1 }, {C 1 , δ 2 }, {C 2 , δ 1 } or {C 2 , δ 2 }, with α = 0.03), the smallest errors were obtained when considering C 2 and δ 1 . Hence, the initial model was reduced to 3 parameters-C 2 , δ 1 and α-, with δ 1 modelling the Neo-Hookean viscoelastic component, and C 2 spanning the nonlinear behaviour. The simplified PK2 tensor that characterises the PVA material becomes The new set of parameters obtained for the 3-parameter model can be seen in Table 3. It is notable that the error increase is very small; thus, the new parametrisation does not significantly deteriorate the quality of the fit. An example of the parameterised 3-parameter model fit to the rheology data can be seen in Fig. 6 (right). The error plot of the 3-parameter model can be seen in Fig. 5 (right). Each phantom indicated a best fit for a very small α (0.02 or 0.03). At low α values, the errors are generally small, since the PVA is only slightly viscoelastic. However, at high α values, there is a sharp increase in the minimal error, which eventually plateaus. This happens when the contribution of the viscoelastic parameter becomes insignificant (δ 1 becomes very small, to counteract the effect of α) and the only remaining parameter is thus C 2 . In this case, as seen in Table 2 as well, the errors rise up to ∼30-35%. In the low α regime, employing the 3-parameter model increases the errors only slightly compared to the 5-parameter model, which is preferred due to reduced complexity.

Loading bias in harmonic wave motion under pure compression
In this study, the PVA phantom material is described by the moderately complex viscoelastic Eq. 24. Harmonic waves With this experimental design, when reconstructing with UR the MRE phantom stiffness in the uncompressed and two compression states, it is observed that the material appears increasingly stiffer with increasing compression. This behaviour is observed in all seven phantoms, as shown in Fig. 8 (column 2). There it can be seen that, in the uncompressed case, the intrinsic stiffness is ∼ 4.5 kPa. At the first compression level, the loading bias leads to an apparent stiffness of ∼ 6 kPa, and at the second compression level to an apparent stiffness of ∼ 7.5 kPa. A summary of the average stiffness of every phantom in each deformation state can be seen in Fig. 9. The statistics are based on a region of interest (as seen in Fig. 8) which excludes the phantom's margins, in order to avoid peculiar boundary effects. A nonlinear stiffening trend can be observed with increased compression. The standard deviations increase as well, from ∼10 to ∼21%, suggesting a higher stiffness heterogeneity under compression, which could be a result of the complex wave behaviour. Since the phantom PVA material becomes apparently anisotropic under deformation, a shear wave (as in Eq. 13) probing, locally, different directions would lead to higher variability in the stiffness measurements.
Recall that in Sect. 2.1.3 we discussed the motion of a planar shear wave through a compressed body in order to gain an intuition of the biasing effects of loading. It was seen that, for an idealised planar wave (like in Eq. 15) moving through a compressed Neo-Hookean body, the material appears to be stiffer. Although the MRE experiment employed in this study carries additional complexity compared to the idealised case (particularly in the type of material used and wave structure), the stiffening trend observed in PVA is in accordance with the simplified theoretical outcomes just summarised. The same stiffening effect was observed in an ultrasound study as well (Gennisson et al. 2007), under similar conditions (idealised planar wave and uniaxially compressed material).

Intrinsic stiffness estimation
By understanding how loading biases stiffness estimation, it is possible to estimate intrinsic stiffness values when using CR. This requires knowledge of the deformation employed and the material law. Since this study deals with pure compression, the deformation was estimated using Eq. 12. Rheological testing of the PVA led to a 3-parameter material law (Eq. 24). It was seen, in Sect. 3.1, that there is variation in estimating the linear parameters C 2 and δ 1 , whereas the nonlinear parameter α is stable around 0.02-0.03. The optimisation performed in rheology was not repeated for MRE, as this was observed to lead to a parameter coupling between α and δ 1 . Thus, α was fixed at 0.03, as indicated by the rheological tests, and C 2 , δ 1 were estimated as described in the data analysis Sect. 2.3.2. The resulting parameters are presented in Table 4. It is notable that the C 2 parameter is generally lower in MRE (Table 4) than in rheology (Table 3), with the exception of phantom 13. This could be related to the degree of relative compression achieved. In rheology, the maximum relative compression is 28.5%. In MRE, however, relative compressions reach up to 35%, with the exception of phantom 13 (only 18%), as seen in Fig. 9. At lower deformation levels, the viscoelastic component is less employed; thus, the purely elastic part becomes more dominant. As C 2 is scaling the elastic component, its estimate is higher when less uniaxial compression is achieved. Within each MRE phantom, experimental imperfections lead to variability between pixels. Looking at the apparent stiffness maps in Fig. 8, excluding the boundary area (as shown in column 2), a small degree of heterogeneity is observed, which accentuates under higher compressions. One reason behind this is speculated to lie in the crystallisation process during the F-T cycles of PVA. Assuming that the PVA freezes first at the mould edges and then in the centre, this could affect the local structures formed. Another reason could be that the top side of the phan- toms is not perfectly flat, as it was cut manually. As such, the compression could have not been ideally uniform. Other arguments could revolve around the compression that the piston is exerting on the phantom, but an analysis in this sense is intricate and beyond the purpose of this study. Regardless of the underlying reasons, the small local variations in each phantom lead to the errors presented in Table 4, as slightly different pixels are attempted to be fit to the same measure. Although the errors are higher than in rheology (Table 3), the error metrics are not directly comparable. In rheology, a transient behaviour is investigated, that is, a data point is followed in time (Eq. 20), whereas in MRE a large spatial sample size is examined-all the pixels across several slices of each phantom, in uncompressed and second compressed state (Eq. 23). The parameters obtained in Table 4 were sought by employing the CR on each phantom in uncompressed and second compression states, as previously described. Subsequently, they were used directly in reconstructing the corrected maps M for the first compression case, corre-spondingly for each phantom. The corrected maps M are illustrated, in a phantom, in Fig. 8 (column 3). Since the parameters are fit to best describe the data in uncompressed and second compressed states, it was indeed expected that the M map's average for these cases is close to unity. Additionally, the corrected map for the first compression case, whose average also lies close to unity, demonstrates the prediction value of the parameters.
The quantified average of map M for all phantoms in all deformation states (uncompressed, first compression and second compression) is shown in Fig. 9 (right). Same as for the apparent stiffness G , only a ROI was considered inside each phantom, to avoid boundary effects. Depending on each phantom's height, a variable compression λ describes the first or second compression. For instance, the data points on the graph corresponding to λ between 0.75 and 0.8 depict the first compression for taller phantoms and second compression for shorter phantoms. It is observed that the standard deviations are sizeable (10-19%) and tend to be higher at the larger compression levels. In the stiffness measurements, this observation was attributed to the apparent anisotropy effect and waves' local probing directions. In the corrected maps M, the apparent anisotropy effect is undone by incorporating the knowledge on the deformation gradient, given by Eq. 12. Since the variance does not decrease in the M maps compared to G maps, it suggests that the predicted apparent anisotropy in compression is not the main cause of the increasing standard deviations. Instead, it could be due to experimental imperfections which are difficult to control. For instance, the slight unevenness of the phantoms' top or the compression induced by the piston indentation are unavoidable, yet for simplicity their effects are not accounted for in the deformation gradient, which employs an idealised form. Despite these, in all but one cases, the average of corrected map M lies within 10% of the ideal value of 1. This indicates that the model, parameterised with the values presented in Table 4, captures the bias induced by loading on the wave motion. This shows that the parameters found can describe the material at different deformation states and are suitable for estimating the intrinsic stiffness of the PVA material.
Integrating rheology and MRE is a difficult problem to tackle. While rheology can be employed to directly measure a material's stress-strain response, MRE relies on investigating the wave propagation behaviour in order to estimate material properties. Here, we tried to integrate the two methods by validating a viscoelastic model using rheological data and then using it to derive the G and G moduli to be used within CR of MRE data. The model was shown to capture the loading bias, thus proving the successful integration of the methods. Furthermore, model parameterisation was done separately for rheology and MRE, yet it is notable that the standard deviations around the mean of the C 2 and δ 1 parameters are overlapping across the two experimental methods. This, together with the undoing of the loading bias, demonstrates the improvement that the newly developed CR brings in analysing MRE data.

Extension and impact in tissue
The theoretical framework presented in Sect. 2.1.2 showed that the stiffness moduli are directly related to the material law and deformation. The work done in phantoms supports the theoretical concept. It was shown that, by incorporating the knowledge on material behaviour and deformation into the MRE reconstruction, the intrinsic stiffness can be estimated. This approach could be extended and applied in tissues. As a first step, a biomechanical model describing the tissue could be inferred and further adapted from animal rheological work. Additionally, by employing registration tools on anatomical images of an organ in undeformed and deformed states, deformations can be quantified. Together, a better characterisation of specific tissue rheology and large deformation recordings using MRI could enable proper interpretation of MRE data.
The most immediate extension of our approach could be breast MRE. This imaging technique relies on the breast being fixed in between two plates, which leads to tissue compression. In the literature, the loading bias was observed in cases where the breast was fixed tightly, compared to cases where the breast was fixed more loosely (Sinkus et al. 2005a), but was not accounted for. Undoing this bias using subject specific deformation would help in correcting the observed stiffness for both healthy and diseased areas. For instance, provided that the risk of breast cancer is associated with tissue stiffness (Boyd et al. 2014), accounting for the loading bias would ensure that the risk assessment is not underestimated or overestimated. Clinically, this could improve the prognostic value and disease management.
Another extension of our study could be in the liver, which is subjected to motion due to respiration. While this results mainly in translational motion, parts of the liver can be strained (Kang et al. 2012). This can be advantageous, since some liver pathologies were shown to become more apparent, compared to healthy tissue, under large strains (in ex vivo) (Yeh et al. 2002). Measuring stiffness in strained tissue could increase diagnostic accuracy, but would require techniques to eliminate loading bias.
The heart is a particularly complex organ which presents structural anisotropy, has dynamic stiffness, and undergoes large motion during the cardiac cycle. The latter constitutes one of the many challenges in assessing the unbiased stiffness of the myocardium (Kolipaka et al. 2010, Glaser et al. 2012, since the effects of deformation non-trivially interfere with the intrinsic tissue properties. The framework presented in this study would constitute a step forward in assessing the tissue characteristics, by eliminating the deformation bias. In order to follow this paper's proposed pathway of retrieving intrinsic material parameters from MRE data, one would have to consider Eq. 9 in the material frame as a starting point. Since the moduli G and G depend on constitutive behaviour and deformation (e.g. Eq. 10), then knowledge on the PK2 tensor and deformation gradient is necessary. The material behaviour can be modelled by considering a constitutive equation for the strain energy function W , like in Eq. 17, from which tensor S can be derived, or by proposing an adaption for S, e.g. Eq. 19. A formulation for the deformation gradient F is required, like the one presented in Eq. 12, which needs to describe the unperturbed macro-deformation that the body is undergoing. From these, a form of the moduli G and G can be inferred, which can be employed within CR (Eq. 22) to eliminate the loading bias.

Study limitations
A limitation encountered during the analysis process of the data presented in this study was the estimation of the viscous modulus from MRE data. The rheological tests performed on PVA phantoms identified a very low α, indicating that there is little viscous response of the material. This is a consequence of choosing a highly elastic material, with low wave attenuation. As such, the characteristics of the viscous component were difficult to identify under MRE testing, leading to an overestimation of the reconstructed viscous modulus. This effect was already identified with the employed MRE reconstruction (Fovargue et al. 2018a), where the overestimation of the viscous modulus was observed in a purely elastic phantom, particularly with increased noise. However, this drawback is predicted to disappear when using in vivo data, as tissue has a stronger viscous component. Nevertheless, in this phantom study, the bias due to the spectrum of stiffness prevented the characterisation of the full complex modulus, thus shifting the focus only on the real stiffness modulus.
This study provides the means through which a material's intrinsic characteristics can be retrieved, despite it experiencing loading. However, a reasonable quantification of intrinsic material parameters relies on a good estimation of the deformation state and on an appropriate constitutive material law. In more complex scenarios, like in vivo measurements, these prerequisites are more challenging to attain. The deformation metrics would presumably need to rely on precise image registration techniques. As medical image registration is a greatly developed field, identification of suitable methods is surely attainable. However, the material constitutive law could be difficult to model, particularly in tissues which exhibit anisotropy. Nonetheless, much work has been done in this sense, as tissue models have been developed in many forms-polynomial, exponential, logarithmic, power laws, etc., and can be further improved and adapted. Undeniably, a model with increased complexity would lead to convo-luted G and G moduli, yet they can be attained numerically (by contrast with the analytic form presented here, for PVA, in Appendix 2). Hence, despite an increased computational cost, the use of a complex model would not be hindered.

Conclusions
This paper presents a framework for integrating knowledge on material constitutive behaviour and deformation into wave dynamics in order to retrieve intrinsic material stiffness properties. First off, a theoretical foundation was established by perturbing Cauchy's equations of motion and transcribing the resulting equations into material frame. This determined a form of the stiffness moduli which was seen to depend on the material model and deformation metrics. To support the theory, experiments on PVA phantoms were employed. The PVA material was tested in a rheological setup and its examined behaviour was modelled using a viscoelastic law. The same material was subjected to MRE experiments, where a known deformation (uniaxial compression) was employed. It was seen that, under the designed testing conditions, using standard MRE analysis, the material appeared to be stiffer with increasing compression. Thus, the known material model and deformation were integrated into the stiffness moduli and subsequently into the MRE reconstruction. By doing this, it was shown that intrinsic material stiffness parameters can be estimated, thus undoing the loading bias observed using standard analysis. Hence, this study provides a framework demonstrated to work in phantoms, which can be adapted and applied to MRE in more complex instances.
J ε may be clustered as 3 Since ∂ J ∂ F = J F −T (Bonet and Wood 2008), it follows that This approximation, considered to hold strongly when truncated to the first-order, reduces the mass balance in Eq. 8b to Here, note that J is assumed to satisfy the original balance law in Eq. 3b, which leads Eq. 27 to become F −T : ∇ 0 u ε = 0.
Using Eqs. 25 and 28, it is now clear that J ε = J +O(ε 2 ), which simplifies the momentum term in Eq. 8a to To further reduce the momentum balance, we want to expand the stress term F ε S ε in Eq. 8a using Taylor expansion in (P, F). As a result, letting S = S e (C(F)) + D α t S v (C(F)) + S p (C(F), P), using ∇ F to denote the derivative with respect to F, i.e.
(∇ F ) i j = ∂/∂ F i j , and the deltas for the difference between the perturbed and unperturbed states, i.e.
3 Here, we denote by the directional derivative of f at a point x in the direction of increment u.
we may focus on the deformation gradient expansion the hydrostatic part expansion and the PK2 viscoelastic part expansion where it is assumed that S v and ∇ F S v are constant in time.
Hence, the components of the stress term F ε S ε are given by where the square bracketed terms denote fourth-order tensors resulting from differentials of the PK2 stress tensor with respect to F. For clarity, in indicial notation, Thus, considering Eqs. 29 and 35, the conservation of momentum can be split between macro-deformation (in square brackets) and micro-deformation components, i.e. ρ J ∂ tt U − ∇ 0 · (F S) + ρ J ∂ tt u ε − ∇ 0 · F∇ F (S e + S p ) + I S : ∇ 0 u ε Noting that the macro-deformation U satisfies Eq. 3a, the first square bracketed term on the left hand side is zero. Hence, a set of linearised elastic wave equations in the reference frame can be formulated: While the set of Eq. 38 deals generically with small scale perturbations, in the case where the frequency is sufficiently high or the time interval sufficiently long, harmonic waves may be considered. Thus, inserting the wave form given by Eq. 2 into the set of Eq. 38, and noting the Laplace transform as D α t e iωt ≈ (iω) α e iωt ∀t >> 0, it follows that +I S] : ∇ 0 u c + p c J F −T = 0, When comparing with real data, it is desirable to map the momentum equations into the physical frame. Using the property (Ogden 1997) that ∇ y = (∇ 0 y)F −1 , the gradients with respect to the reference frame can also be transformed into the current domain. Hence, the terms in the set of Eq. 39 become In indicial notation, Moreover, by noticing the property linking the divergence of a tensor between the physical and reference domains (Gurtin et al. 2010) the divergence of the stress terms can be recast in terms of the physical gradient operator ∇, i.e.
where the fourth-order tensor C encompasses Combining Eqs. 40-42 and noting the property linking the divergence of a scalar in the reference domain to the physical frame (Gurtin et al. 2010) the equations analogous to the classical set of Eq. 3, may be written as, with the understanding that the real and imaginary scaling moduli G and G are given by G = Re(C), G = Im(C), i.e.

Viscoelastic moduli estimates for PVA material
In this section, the aim is to understand the effect of macrodeformations on the wave behaviour in a viscoelastic material described by Eq. 19. So far, general elastic and viscoelastic stresses S e and S v have been used. Recall that, for this work, they were considered to comprise two parts each, based on a Mooney-Rivlin type of law, as follows: Hence, the fourth-order tensor C can be taken one step further and written as while the G and G moduli become The moduli G and G are influencing the wave dynamics ∇u c , leading to an apparent wave dynamic. Now, considering the analytic formulation of the contraction (G +iG ) : ∇u c , consider that the wave ∇u c satisfies the harmonic wave motion described in Sect. 2.1.1. After a cumbersome derivation and simplification process, the moduli defining the PVA material used can be formulated as