Viscoelasticity Imaging of Biological Tissues and Single Cells Using Shear Wave Propagation

Changes in biomechanical properties of biological soft tissues are often associated with physiological dysfunctions. Since biological soft tissues are hydrated, viscoelasticity is likely suitable to represent its solid-like behavior using elasticity and fluid-like behavior using viscosity. Shear wave elastography is a non-invasive imaging technology invented for clinical applications that has shown promise to characterize various tissue viscoelasticity. It is based on measuring and analyzing velocities and attenuations of propagated shear waves. In this review, principles and technical developments of shear wave elastography for viscoelasticity characterization from organ to cellular levels are presented, and different imaging modalities used to track shear wave propagation are described. At a macroscopic scale, techniques for inducing shear waves using an external mechanical vibration, an acoustic radiation pressure or a Lorentz force are reviewed along with imaging approaches proposed to track shear wave propagation, namely ultrasound, magnetic resonance, optical, and photoacoustic means. Then, approaches for theoretical modeling and tracking of shear waves are detailed. Following it, some examples of applications to characterize the viscoelasticity of various organs are given. At a microscopic scale, a novel cellular shear wave elastography method using an external vibration and optical microscopy is illustrated. Finally, current limitations and future directions in shear wave elastography are presented.


INTRODUCTION
Changes in mechanical properties of biological soft tissues are often associated with physiological dysfunctions. Viscoelasticity is an important mechanical biomarker to characterize structural changes and/or constituents of tissues. However, the assessment of elasticity through imaging has been more often exploited than viscosity, and an historical perspective of development has been to replace manual palpation by physicians and to answer an ultimate and natural question: 'is the region hard or soft?' The elasticity, represented by the Young's modulus E, is able to characterize tissue deformation by using a linear relationship between stress σ and strain ε as E σ/ε. The rationale behind the elasticity assessment of biological soft tissues is that Young's moduli of different types of human tissues differ by a few orders of magnitude [1], and are affected by the presence of a pathology.
Since biological soft tissues are hydrated, they are not only represented by their solid-like behavior using elasticity, but also by their fluid-like behavior using viscosity. Viscosity represents the hysteretic effect between stress and strain applied on a tissue. It is becoming an important biomarker of pathological changes in biological tissues. Mechanical test is the most basic method to measure the viscosity of ex vivo soft tissues. To this end, a constant strain is applied to a test specimen. The stress relaxation with time is used to characterize the viscosity of the tissue. Another way of mechanical test is dynamic mechanical analysis. Periodic strain or stress is imposed on a specimen. The dynamic responses in terms of different incentive frequencies are associated with tissue viscoelasticity. Relevant works on mechanical testing of biological soft tissues can be found in [2][3][4][5][6][7][8][9]. Alternatively, imaging-based approaches are becoming more popular for in vivo viscoelasticity assessment. These approaches exploit tissue deformations in acoustic, magnetic or optical fields to characterize viscoelasticity in a non-invasively manner. Soft tissues can maintain their function during the measurement avoiding destructive testing [10]. Based on a clinical perspective, the in situ and localized assessment of tissue viscoelasticity through imaging had major impacts on diagnosis (e.g., cancers, liver fibrosis, musculoskeletal disorders, cardiovascular diseases, etc. . .). Since biological tissues consist of cells, extracellular matrices, and structural proteins, a recent field of development has been to study sub-cellular biomechanical properties associated with pathological processes through imaging. This finding encouraged researchers to impel bioelasticity research further into a microscopic scale.
This review aims to provide a state-of-the-art summary of developments made in the field of shear wave elastography, which concerns elasticity and viscosity imaging through mechanical shear wave analysis. This technology requires a shear wave source, the tracking of shear wave propagation through imaging, and the processing of the shear wave propagation characteristics through physical models or image processing algorithms. Shear waves can be generated by external or internal vibrating sources. An external mechanical actuator in physical contact with an organ or cell is a common way to induce shear wave propagation from the surface to the core, whereas acoustic radiation or Lorentz forces can be used as internal in situ localized shear wave sources. The detection of the shear wave propagation is usually performed using ultrasound (US), magnetic resonance (MR), optical or photoacoustic imaging methods. Elasticity and viscosity can be obtained from estimations of shear storage and loss moduli, which require the determination of the shear wave velocity and attenuation into the interrogated medium. In the following sections, the generation and tracking of shear waves are described. Determining elasticity and viscosity maps through the solution of an inverse problem based on elastic wave propagation equations and underlying assumptions are also addressed throughout the text. Note that this review is not intended to detail artifacts and confounders of shear wave imaging because these are organ, tissue structure, and tissue pathology specific.
Nevertheless, such information is presented briefly in some sections.

BIOMECHANICAL PRINCIPLES OF SHEAR WAVE ELASTOGRAPHY General Concepts in Shear Wave Elastography
The major approach since the early steps of elastography imaging has been to approximate the tissue as isotropic and purely elastic. This situation has been widely described so the main constitutive relations essential to shear wave elastography (SWE) understanding are recalled here. The relationship between the applied stress and the strain response of the solicited tissue is given by Hook's law: where T and S are the stress and strain tensors, respectively, and c ijkl contain the elastic parameters of interest. Under the assumption of small deformations, the strain tensor is given by: where u is the 3-dimensional motion field in unit of (m) induced by stressing the tissue. The time domain wave equation describing the propagation of local displacements is obtained using Newton's second law: where ρ is the material's density in [kg/m 3 ] and f [Nm −3 ] is the source term. After Fourier transform into the frequency domain, Eqs 1-3 are expressed as: T ij (r, ω) c p ijkl (r, ω)S kl (r, ω) S kl (r, ω) −ρω 2ũ (r, ω) ∇ ·T(r, ω) +f (r, ω) where ω 2πf with f the frequency, and the tilde (∼) and star (*) notations refer to complex numbers. Full development of Eq. 3 (available in Ref. [11]) leads to the governing equation of motion propagation in elastic solids known as the Navier's equation: where μ is the shear modulus in (Pa) reflecting the amount of energy the tissue can store as elastic deformation, and λ is the first Lamé coefficient in (Pa) reflecting the tissue's compressibility. The Navier's equation does not rely on rheological models and conveys an approximation of the material's natural properties based on physical assumptions. The same development applied to Eq. 6 leads to: where G p G′ + jG″ is the general notation of the complex shear modulus in the frequency domain of Navier's equation, G′ is the shear storage modulus that reflects the amount of mechanical energy stored as shear deformation in the solid, and G″ is the shear loss modulus reflecting the amount of mechanical energy dissipated due to shear viscosity. Similarly, λ p λ′ + λ″ is the complex Lamé coefficient, where λ′ and λ″ are the compression storage and loss moduli, respectively, reflecting the amount of energy stored and lost in the solid due to compression deformation and compressional viscosity. In purely elastic solids, the wave field does not dissipate and shows an instantaneous response to load. In such cases, the loss moduli equal zero, and G p G ′ μ and λ p λ. The Young's modulus E, characterizing the solid's resistance to deformation under loading, and the Poisson's ratio ] characterizing the tissue's compressibility, are defined by: and ] λ 2 λ + μ The motion field u propagates as compression and shear waves of which velocities v c and v s are, respectively, given by (see Ref. [11] for details): v c λ + 2μ ρ (11) and v s μ ρ (12) The large difference between the two velocities observed in biological tissues (typically v c is around 1,540 ms −1 and v s is found around 1-10 ms −1 ) suggests that λ is much greater than μ, thus allowing for the following approximation often used to report SWE measurements: This observation is directly linked to the tissue incompressibility assumption, which causes λ to approach infinity [12]. Additionally, λ was shown to vary slightly as opposed to μ, which spans a few orders of magnitude in biological tissues [13]. Consequently, the shear modulus or the Young's modulus is the mechanical parameter considered in elastography reconstruction processes. Another approach much less often used to reconstruct material mechanical properties in SWE is to model the solid as isotropic and viscoelastic instead of purely elastic. Here, the response of the loaded tissue shows a delay with respect to actuation and elastic waves attenuate due to energy dissipation in the solid, which is specific to viscous materials. Attenuation may be accounted for in the time-domain Navier's equation by introducing a damping term, and linking it to viscosity using rheological models such as the generalized Maxwell, standard linear solid, or Kelvin-Voigt, which are further discussed in the next section. The following Table 1 presents the stress-strain relationships used to derive wave equations from Newton's second law for the three aforementioned rheological models [14].
A different option to integrate viscoelastic properties into the description of the material is to formulate the problem in the frequency domain, provided a harmonic actuation. In this case, the complex shear modulus, as described in Eq. 8, has a non-zero imaginary part accounting for shear wave dissipation due to the material's shear viscosity. Relating G″ values to actual viscosity values depends on rheological modeling, as discussed in the next section. For instance, in the case of a solid described by the Kelvin-Voigt model, the complex shear modulus is given by G(r) p G(r)′ + jG(r)′ μ(r) + jη(r)ω, where η(r) is the local shear viscosity.
For G′ and G″ estimation, the displacement fieldũ may be used as the solution to the Navier's equation in direct or iterative inversion, or the shear wave velocity at frequency ω 2π may be measured. Here, the complex wave number is noted as , in analogy with the complex shear modulus notation, and the dispersion relation is given by k ′ ω v p , where v p is the phase velocity of the shear wave at the frequency ω 2π . Considering a plane wave decomposition of the wave field, i r i . The imaginary number k″ is often noted α and is the shear wave attenuation coefficient (m −1 ). Thus, a linear system of two equations may be raised and independent experimental evaluations of v p and α from the displacement field allows assessing G′ and G″. Also, the Young's modulus becomes complex-valued in viscoelastic models. However, most viscoelasticity reconstruction processes stick to the evaluation of G′ and G″. Finally, it is to be noted that the equivalence between longitudinal and compression waves on one hand, and transverse and shear waves on the other hand, is true for plane waves only.
To date, isotropic elastic and viscoelastic characterization of soft matters have mostly been considered in shear wave elastography, owing to the availability of various inversion schemes. However, the anisotropic and poroelastic nature of certain biological tissues, such as the brain, has long been acknowledged. In poroelasticity, the medium is modeled as a porous solid matrix crossed by flowing fluid, and thus contains two separate phases, as opposed to more common models containing one phase. Consequently, the motion field measured in imaging protocols is not only due to the solid tissue deformation but also to the pressure gradient in fluid pores. Although poroelasticity was early studied using quasi-static deformations [15], implementation in shear wave elastography imaging remains in its infancy. Oscillatory deformations in poroelasticity have first been described by Biot [16,17], and later by [18]. Assuming a viscous fluid flow, fluid saturation in pores, and a compressible linearly elastic solid, poroelasticity equations of propagation are given by: where µ is the shear modulus, λ the first Lamé's parameter, u the complex time harmonic displacement field, β the effective stress coefficient (dimensionless), p the complex time harmonic pressure field, ω the actuation frequency, ρ the bulk density, ρ f the pore fluid density, ϕ p the material porosity, κ the hydraulic conductivity, and ρ a the apparent mass density. Finally, anisotropy has been considered for biomechanical modeling in the context of shear wave elastography. In such cases, the material's mechanical response (strain) is dependent on the direction in which it is solicited (stress). Full derivation of relevant mechanical parameters under different symmetry assumptions is beyond the scope of this review, and the interested reader is referred to the excellent pedagogical development in [19]. Briefly, Hooke's law describes the relationship between applied stress and material strain: where σ and ϵ are the stress and strain tensors, respectively. In isotropic materials, the tensor C ijkl is fully described by two parameters, E and ]. In anisotropic materials, more constants are needed to account for the direction dependance of ϵ. The most used anisotropic model is transverse isotropy, which is particularly used to characterize fibrous tissues (e.g., muscles). Transversely isotropic materials are organized in layers where inplane mechanical properties are isotropic, and out-of-plane ones are anisotropic. In such cases, 5 parameters are necessary to describe C ijkl . Two of them, μ 13 and μ 12 , characterize the shear motion along and perpendicular to the fiber axis, respectively. The other three, E 1 , E 2 , and E 3 characterize the compression motion along, perpendicular in-plane, and perpendicular out-ofplane to the fiber axis, respectively. Injection of the Hooke's law into Newton's second equation allows to derive equations of shear wave propagation along directions of dependency the same way as to derive the Navier's equation of elasticity. Other anisotropic models exist, such as the orthotropic one which shows a lower level of symmetry (three perpendicular planes) than the transverse isotropy model. The orthotropic model has been used to develop waveguide elastography, which describes the propagation of the different polarizations of shear waves along separate directions. The orthotropic tensor along with equations of polarized wave propagation are described in detail in [20].
Frontiers in Physics | www.frontiersin.org June 2021 | Volume 9 | Article 666192 4 the complex shear modulus by tracing and analyzing a propagated SW in the specimen.
During the past years, studies in the field of SWE measurements were often targeting absolute values of elasticity |µ| and viscosity η. To do so, rheological models of the material are needed to derive those parameters. For instance, when the Kelvin-Voigt model or Maxwell model is considered, the complex shear modulus can be written as follows [28,34,35]: where G KV represents G* of the Kelvin-Voigt model, and G M is that for the Maxwell model, both satisfy G* G ′ + jG″. In Eq. 18, ω S is the angular frequency of the SW. Alternatively, by solving the wave equation in Eq. 8, a general solution of G′ and G″ without considering a rheological model can be obtained, as in [28]: Here, as synergized with Eq. 18, one can see that both |µ| and η are functions of v S , α S , and ω S . Since v S and α S can be experimentally measured at certain ω S , |µ| and η, in a specified rheological model, can be thereafter calculated.
In common practice, when the viscosity η is taken into account for tissue characterization, it can be determined either directly using both v S and α S with knowing the corresponding ω S , or alternatively by evaluating the dispersion of v S with respect to ω S without determining the value of α S , i.e., by knowing multiple pairs of v S and ω S [27,36,37]. However, one should also notice that the viscoelastic property of biological tissues are rather complex, depend on the tissue type, and on the presence of a pathological condition, so that there is not simply a best, or a most appropriate material model for all tissues. Furthermore, pathological changes with time of a tissue could also lead to a major change of its viscoelastic property, thus a certain material model may no longer be suitable for the tissue when it becomes abnormal and progresses toward a more severe pathological state. Meanwhile, the derivation of elasticity |µ|, and viscosity η, are rather different among different material models. That means, using different models with same measures (such as v S and α S ) would lead to different results [28], and hence would be meaningless to clinical studies. Therefore, nowadays it is always suggested that no rheological model is assumed, and instead of that, directly access the shear storage G′ and loss G″ moduli would be not only more rigorous and appropriate, but also mathematically convenient to describe tissue viscoelastic properties [28,37].

Generation and Detection of Shear Waves
A shear wave, also called a transverse wave, is a moving mechanical wave that consists of particle oscillations occurring perpendicular to the direction of the energy transfer [38]. As briefly introduced, SWs in ultrasound imaging can be generated either from an external vibration source (such as a mechanical vibrator/shaker) [39][40][41][42], or internally by an acoustic radiation force (ARF) [13,[44][45][46][47][48][49][50][51], as illustrated in Figure 1. The control of the SW amplitude and frequency of the ARF is considered in [13]. In terms of waveforms, the SW can also be generated as continuous waves [39,40,42,[44][45][46] or impulse waves [41,[47][48][49][50][51], as can be seen in the examples of Figure 2. In ultrasound imaging, the probe fires longitudinal pressures and detects particle displacements along the axial direction, therefore only SWs that propagate along the lateral direction of the ultrasound beam, or SW components whose displacements occurred on the axial direction, can be detected.
Many remarkable techniques were invented over the past 30 years based upon different combinations of external or internal SW sources, and continuous or impulse SWs. In 1988, Lerner et al. proposed a method to map the propagation of low frequency SWs with a Doppler ultrasound displacement detection technique to assess tissue stiffness [39]. Later in 1990, Yamakoshi et al. proposed a dynamic measurement method to determine the SW speed v S using an external mechanical vibration source [40], as seen in Figure 2A. The parameter v S was determined by analyzing the wavelength of a continuously propagated SW using the Doppler ultrasound technique. Catheline et al. developed in 1999 an impulse SW measurement method [41]. In this method, an ultrasonic probe was located at one side of the specimen to capture the propagation of the impulse SW generated by a mechanical vibrator located at the other side of the specimen. A plane wave ultrasound system was used, enabling a high frame rate in detection mode, then the parameter v S was determined through the time-of-flight (TOF) technique applied on successive images. See below for more information on the TOF method. Since such a method used an impulse SW, it is also called transient SW imaging, or transient elastography (TE). Nowadays, the largely used clinical device Fibroscan [52,53] is based on TE. In 2004, a method using two external SW sources to generate continuous SWs toward each other with slightly different SW frequencies was proposed [42]. Due to the frequency difference, interfered SW patterns termed as "crawling waves" moved with a much slower speed than the expected v S , This allowed the observation of propagated crawling waves with a conventional low frame rate B-mode imaging system. Once the speed of the crawling wave is obtained, v S could be derived. In 1998, Sarvazyan et al. developed a SW measurement method, termed shear wave elasticity imaging (SWEI), by using a SW remotely generated by an ARF of a focused ultrasound beam [13]. In this method, a transient SW pulse was firstly produced at the focus of the ultrasound beam and propagated along sideways. Then, imaging transducers were used to trace the moving of SW fronts and viscoelastic parameters were derived thereafter. The same year, Fatemi and Greenleaf developed a method to produce an oscillatory ARF by mixing two ultrasound beams with different frequencies [44]. A short period of harmonic (continuous) or tone-burst SWs was generated and propagated along sideways, which made the SW narrow-band (while the SW generated by Sarvazyan's method was broadband) and then v S could be determined by finding the phase difference of the SW at two apart locations [45]. By repeating the measurement with continuous SWs at different frequencies, or retrieving different frequency components of a tone-burst SW, both G′ and η could be derived. This method is termed as SW dispersion ultrasound vibrometry (SDUV). Later in 2004, based on the combination of Catheline's impulse SW method [41] and the ARF technique [13], Bercoff et al. developed an advanced SW measurement method known as supersonic shear imaging (SSI) [47]. With this method, an ultra-high-speed scanner is used, then multiple ARF impulses are triggered consecutively and very quickly at different depths. Each impulse produces a SW point-like source then all these SWs are interfering constructively and result in two SW planes propagating in opposite directions, as can be seen in Figure 2B. A two-dimensional v S image is obtained with this method. Moreover, since this technique is creating broadband SWs, tissue viscosity can also be estimated through the v S dispersion method (using the same principle as SDUV). In 2012, Song et al. developed a SW method, which also used multiple lateral ARFs as in [48]; the method was termed as comb-push ultrasound shear elastography (CUSE) [49]. It firstly generates multiple ARF excitations at different spatial locations to produce multiple impulse SWs by using the push mode of the ultrasound probe, and then quickly switches to the scanning mode of the probe to detect the SW propagation. Therefore, v S could be measured through the TOF technique by tracing the movement of the SW front from each SW source. The use of the comb-push excitation provided multiple SW sources in the specimen so that such method is effectively compensating for the worse signal-to-noise ratio (SNR) due to the SW attenuation at the location far away from a given SW source.
Knowing that v S λ S × f S , one may like to measure v S through the TOF technique, or instead to determine λ S in the spatial domain. For most biological tissues, v S travels at a few m/s, which means it only takes a couple of ms for a SW to travel through the entire field of view (FOV) of a common ultrasound probe. Physically, a focused ultrasound system triggers transducer elements sequentially from one edge to another to complete a B-mode scan, as a result the frame rate is typically less than 100 frames/sec in such a system, which is not fast enough to measure the TOF without a particular modification of the experimental setup [42]. Therefore, measuring λ S becomes the realistic option. In this scenario, the spatial resolution is determined and limited by λ S . Although increasing f S is reducing λ S , and so is improved the resolution, one should also notice that a higher frequency would cause a quick attenuation of the SW propagation. Thus, empirically f S is usually adjusted to a few hundreds of Hz, which leads the lateral resolution of elastography images to a subcentimeter level with using standard focused ultrasound beamforming. On the other hand, a plane wave system can trigger all transducer elements of the probe at the same time to emit a plane compression wave, enabling it to have a very high frame rate in B-mode (up to 10,000 frames/sec) [47]. Thus, in this scenario, the TOF technique is applicable, and theoretically the distance that a SW travels within two consecutive frames could be as short as a few-tenth of mm. Since the distance is comparable to the physical interval of two adjacent transducer elements in a common ultrasound array probe, the lateral spatial resolution of SWE with using a plane wave system is approximately the same as that of B-mode imaging [54]. It is also worth noting that, when tissue boundaries/layers exist under the FOV, and a propagating SW passes through those interfaces, physical phenomena such as reflection, refraction, diffraction, and mode conversion could occur at the interfaces and cause artifacts. Although directional filters are usually applied to mitigate those effects [55,56], practically it is still difficult to remove all the unwanted waves. Therefore, the v S measured within approximately one λ S from the interfaces are usually considered unreliable, which to some extent would downgrade the spatial resolution at those areas [57].

Viscoelasticity Reconstruction
As introduced earlier, ultrasound SWE contributed to the noninvasive assessment of mechanical properties of soft tissues [58][59][60][61][62]. One method largely used clinically is transient elastography (TE) [63][64][65][66][67], which utilizes a dynamic compression generated by the vibration of the transducer on the skin to produce shear waves. No structural imaging is provided with this method to guide the measure. Moreover, in some patients with morbid obesity and ascites, the attenuation of shear waves travelling from the surface of the body to the organ of interest (typically the liver) may avoid reliable measurements [68,69]. The acoustic radiation force impulse (ARFI) [60,70] and supersonic shear imaging (SSI) [47,71,72] methods use a radiation pressure to locally induce shear waves within the organ of interest. TE, ARFI and SSI are assessing tissue elasticity (no viscosity) based on the measurement of the shear wave speed [73]. The Young's modulus is estimated and displayed as an image using E 3 |µ| ρv 2 s . Alternatively, different approaches have been developed to retrieve and display the viscous component of a tissue [50,51,74,75]. Such approaches have not yet been validated on large clinical cohorts, nor implemented on clinical scanners. Details on technologies proposed to determine tissue viscoelasticity are given next.
Most studies utilized a rheological model, which has been introduced earlier, e.g., the Kelvin-Voigt model, to find the viscosity after the reconstruction of elasticity [45,[76][77][78]. The complex shear modulus G* G′+ jG″ was also estimated using different approaches, such as measuring the acoustic radiation force-induced creep [74], solving the Navier's wave equation numerically [79], inverting analytically the solution of the shear wave scattering from a mechanical inclusion [80][81][82], and using a finite-element based method [83]. Note that a torsional SW source was used in the latter method, as originally proposed in [84].
Kazemirad et al. [50] developed a method for the quantitative measurement of viscoelastic parameters G′ and G″ at various frequencies, based on the assumption of a cylindrical shear wave front produced by a radiation pressure, allowing to avoid wave diffraction effects. Other studies have also used the same geometrical assumption for quantitative viscoelastic measurements [75,85,86]. Notice that the cylindrical wave front assumption would not necessarily hold when considering inhomogeneous media, such as a tissue embedding a tumor. A recent method for estimating tissue viscosity without geometrical assumption on the wave front was proposed by [51], which utilized the shear wave velocity v S and attenuation α S computed by the frequency shift method [87]. Recently, [88] performed a study to characterize viscoelastic properties of oil-ingel viscoelastic phantoms and in vivo human livers. They found that the shear wave dispersion and attenuation were linked together and related to the tissue viscosity. As reviewed above, the shear wave speed and attenuation are widely used for reconstructing viscoelastic properties. Experimental methods to obtain those shear wave properties are separately described below.

Shear Wave Speed
One of the widely used methods implemented on clinical scanners is the group velocity [89][90][91][92] that assumes the tissue as elastic, homogeneous, isotropic, linear, and of infinite dimension with respect to the wavelength. The group velocity is estimated using time-of-flight (TOF)-based algorithms for particle displacement or particle velocity assessments in the time domain [72,93]. TOF-based algorithms are usually based on cross-correlation (CC) [94] and time-to-peak (TTP) methods [90]. Basically, the CC provides a moving average estimate of the shear wave speed using all sample points, and performs multiple crosscorrelations along the direction of the wave propagation, which may result in artifacts for periodic shear wave patterns, whereas TTP estimates the velocity based on the tracking of the movement of one point on the waveform [90]. Although group velocity estimation methods are considered robust [91], they are theoretically applicable to strictly elastic materials thus requiring resorting to other techniques for the evaluation of the viscous behavior [95].
The variation of the shear wave velocity with frequency refers to the wave dispersion happening in a viscoelastic medium [45]. Some methods have used this phenomenon to evaluate viscoelastic properties of tissues [45,96]. Measuring shear wave velocities at specific frequencies is known as phase velocity estimation [45]. Beside the viscoelastic property of a tissue, its finite thickness can also affect the dispersion due to reflections during propagation, which may result in wave mode conversion [97][98][99]. The phase velocity and the group velocity are not equal in the presence of dispersion. It was shown that the phase velocity has a lower value by a factor of 8-9% compared with the group velocity in soft tissues [100].
One technique to measure the phase velocity is the phase gradient approach, which estimates the velocity using the phase difference evaluated at different spatial locations for specific frequencies [45,85,96]. An alternative method to estimate the shear wave phase velocity is performed by two-dimensional Fourier transform (2D-F) analysis, which converts spatiotemporal data to a wavenumber in the frequency domain, and uses the peak magnitude distribution to estimate the phase velocity [101,102]. The dispersion either from the phase gradient or 2D-F can be fitted to rheological models to quantify viscoelastic parameters of the medium [45,103]. The attenuating nature of a tissue is the cause of the dispersion of the phase velocity. The shear wave dispersion and attenuation can be estimated by the computation of a power law coefficient, with the assumption of a power law rheological model for the tissue [100,104].
Local phase velocity imaging (LPVI) [105,106] is another method that can produce a phase velocity map. The LPVI requires applying bandpass filters to obtain the maximal frequency range for the phase velocity. Although this method demonstrates good reconstructions of 2D shear wave phase velocity, results are sensitive to the frequency range selected, and they may change when using different transducers, focal configurations, and focal depths [78].

Shear Wave Attenuation
The dependency of the wave amplitude with distance is attributed to geometrical spreading of the wave energy and to viscoelastic attenuation. Wave diffraction by geometrical spreading can be reduced by using the method of [50] that is considering cylindrical shear waves produced by a supersonic radiation pressure source. Other methods including this assumption were based on a 2-D Fourier transform and the computation of the spectral width to assess the frequency dependent attenuation [75,107,108]. A robust method assuming a cylindrical wavefront and no rheological model is the attenuation-measuring ultrasound shear wave elastography (AMUSE) algorithm [107]. This method, however, does not provide any attenuation map since the computation requires all datasets within the selected region-of-interest. Since biological tissues such as the kidney, muscles, and tendons are anisotropic, and because the wave produced by a linear SW front may no longer be cylindrical in those media; then, abovementioned algorithms may lead to inaccurate results [96,[109][110][111][112].
Frequency-shift methods used for compression and seismic wave analyses [113][114][115] inspired the field of shear wave elastography to assess tissue viscosity. Frequency-shift methods are not based on wave amplitude, so the dependency of these methods to geometrical wave spreading is released [113]. Bernard et al. developed such a frequency-shift method for shear wave attenuation by model-fitting of the amplitude spectrum [87]. This method made a few assumptions, which may not hold in all viscoelastic media such as fatty liver. A two-point frequency-shift method was later proposed by Kijanka and Urban to soften assumptions made by Bernard et al.; in their report, they considered a varying shape parameter of the gamma distribution used to fit the shear wave amplitude spectrum [116]. This technique used only two spatial points instead of all points along the propagation path, as in [87], to estimate the attenuation coefficient [116]. Using only two spatial points reduces the computation time but may affect robustness in cases of noisy shear wave displacement maps. Viscosity maps based on the cylindrical wavefront assumption of [50] or frequency-shift method of [87] can be found in [51].

Applications
A few examples of ultrasound shear wave viscoelasticity imaging applications are presented next. The reader may refer to recent review papers on this subject for other examples [117][118][119]. The focus below is on the liver and breast as those organs were largely investigated in clinical studies using SWE.

Liver
Liver fibrosis occurs when an abnormal large amount of liver tissue becomes scarred. It can lead to cirrhosis, its long-term sequel, and further evolve as hepatocellular carcinoma [120]. Liver fibrosis can be differentiated into 5 categories, from F0 for a normal liver to F4 for cirrhosis; these categories have been obtained by biopsy intervention, which is the gold standard for liver classification. However, liver biopsy is invasive, could lead to bleeding or worse outcomes, and even death [121], and because a small amount of tissues is taken, it is not always representative of the full liver due to sampling errors [118,122]. Fibrosis is one pathology known to increase liver stiffness [69,[123][124][125][126] along with inflammation, edema, congestion and extra hepatic cholestasis [64,127,128]. Shear wave elastography was mainly used to classify fibrosis based on liver elasticity using different cutoff values. This imaging method is accurate to assess liver fibrosis of stage 2 and higher [69,123,[129][130][131], has a good repeatability [132], and may allow to diminish the number of biopsy [133]. Yet the impact of steatosis on liver stiffness is uncertain [88,[134][135][136]. To overcome this, some teams proposed investigating viscous properties. If no clear consensus is reached yet, a few studies showed promising results based on shear wave dispersion and attenuation to assess steatosis stages [88,137] or necroinflammation [138,139]. Shear wave elastography presents some limitations for liver imaging, such as difficult measurements in obese patients, and confounding impact of factors such as inflammation, which can increase liver stiffness or change the liver stiffness threshold for classification. In addition to fibrosis assessment and classification, SWE was also proven useful to follow patients with chronic liver disease [65]. An example of a liver SWE image is given in Figure 3.

Breast
Shear wave elastography is used to help identify breast cancers, since it has been shown that malignant tissues appear stiffer than its healthy counterpart [140][141][142]. X-ray mammography, MR imaging and ultrasonography are used to detect tissue lesions or to classify suspicious masses into different categories, typically classification 0 for incomplete data to 6 for histologically proven malignancy. Nonetheless, excluding the expensive MR imaging method, these approaches have poor specificity and mammography often find false negative results in dense breasts [143]. Category 4, which corresponds to suspicion for malignant tissues, has a degree of certainty varying from 2 to 95% to assess malignancy proven by biopsy, and has a cancer detection rate of 10-30% [144]. Shear wave elastography allowed improving breast lesion characterization [145][146][147][148], and reducing the number of unnecessary invasive biopsy due to the improvement in specificity [148][149][150]. Elasticity parameters such as the maximum or mean Young modulus E within the lesion, and in surrounding tissues, are used to separate benign from malignant masses. Recent works investigated viscosity behavior using the shear viscosity [151], linear dispersion slope [152], and storage and loss moduli [153] to differentiate malignant from benign tissues. Ultrasound data on the viscous behavior of breast lesions are scarce but may prove to be of clinical value in the future. Some studies investigated the use of SWE as a tool to monitor cancer treatment performance [154], with a decrease in malignant mass elasticity during treatment, or for early prediction of therapy successes [155,156], with a better treatment response for softer tumors. Figure 4 gives examples of Young's modulus elasticity maps of breast lesions.

Other Applications
Although SWE has targeted mainly the liver and breast, other organs and techniques have been developed. Prostate cancers [158], thyroid cancer nodules [159], and blood clot characterization [160] have been investigated, to name a few examples, with the Young's modulus as the descriptive mechanical parameter. If the assumption of an isotropic medium is generally accepted for most organs, it is not the case for muscles and tendons. Anisotropic and transversely isotropic models using shear waves have been recently investigated [24,161,162], some other teams explored viscoelastic properties using different probe orientations [109,[163][164][165]. A non-exhaustive list of SWE clinical applications can be found in Table 2.

MAGNETIC RESONANCE ELASTOGRAPHY
Magnetic resonance elastography (MRE) is another non-invasive imaging technology for assessment of mechanical properties of soft tissues. Since its first description by Muthupillai et al. in 1995 [166], MRE has been integrated into clinical routines for liver disease detection, and has shown great potential for other organs, notably the brain, of which only MRE can assess the in-vivo viscoelastic components without surgical intervention. Principles of MRE investigation are similar to those of any SWE method ( Figure 5). A major feature of MRE resides in its ability to measure 3D displacement fields by simply changing the axes of encoding gradients, which is an advantage over other imaging devices operating elastography. The main drawback may be found in the longer scan times relative to ultrasound elastography for instance. MRE has a poor temporal resolution and relies on a stroboscopic-like recording arrangement to generate time resolved images, as opposed to ultrasound SWE where burst measurements are performed at a high acquisition rate. Typically, MRE data contain 4 to 8 images per harmonic actuation cycle. Spatially, MRE is sometimes referred to as a super resolution imaging modality as measured displacement amplitudes are much smaller than the image pixel size (tens of microns versus one to 3 mm). We review in this section the main three steps in MRE investigation, namely motion generation strategies, motion encoding techniques, and inversion methods. Finally, applications to the liver and brain are discussed. These organs were subjectively chosen as liver disease diagnosis is the only MRE protocol clinically established, and non-invasive in-vivo brain mechanics assessment is not enabled by any elastography techniques other than MRE.

Generation of Acoustic Waves in MRE
In MRE, most applications involve the generation of timeharmonic wave fields using external surface actuators. These actuators must meet the requirements imposed by magnetic resonance safety rules, in other words they must be made of non-magnetic materials and be adaptable to fit into the experimental or clinical magnet bore of the scanner. Design of actuators has been shown to be application dependent; we review in this section the main techniques to induce motion in soft tissues in the context of MRE. Loudspeakers have been widely used to transmit motion to tissues and may be divided into pneumatic and rigid categories. In the pneumatic one, air pulses are transmitted from the loudspeaker pulsing  [167][168][169][170][171][172][173][174][175][176], and the electricity-free transmission system (pneumatic) requiring no electrical current inside the magnetic resonance imaging (MRI) room. The main drawback may reside in the limitation to simple mono-frequency waveforms only. In the rigid category, loudspeakers transmit motion to tissues via a rigid rod attached to the membrane on one end, and to the patient on the other end. This configuration is also versatile [177][178][179][180][181][182][183][184] and handles arbitrary waveforms but requires the loudspeaker to be inside the MRI room. Lorentz-coil actuators have also been used to generate motion in various organs [185][186][187][188][189][190][191][192], and rely on the coupling of the MRI magnetic field B 0 with an electrical current injected into the coil. These actuators must be designed according to the targeted organ as they are placed relative to B 0 . Piezoelectric drivers allow to deliver arbitrary waveform pulses to tissues while avoiding the constraint of positioning relative to B 0 . Significant displacement fields could be obtained using this technique in various conditions (human abdomen [193,194], human brain [195,196], human breast [197], and mouse brain [198,199]). A major bottleneck in generating sufficiently high motion deflection in tissues is the decrease of the motion amplitude with increasing excitation frequencies. Using a wide range of mechanical excitation finds its application in the analysis of frequency dependent mechanical behavior of soft tissues [8,178,191,198,[200][201][202][203]. The more frequencies, the more valuable is the information. Centrifugal force based MRE drivers have been proposed to circumvent this limitation at high frequencies. The centrifugal force allows to maintain the displacement amplitude high regardless of the frequency [204]. The "air-ball" actuator and "gravitational actuator" are the first examples of the centrifugal force implementation. The "air-ball" actuator [205] consists in a ball circulating in a circular chamber under injection of compressed air. The revolution speed of the ball is imposed by the pressured air of which the pressure determines the actuator vibration frequency. The "gravitational transducer" is made of a mass attached to an axis and rotating around this axis. The rotation speed is driven by the rotating axis connected to a motor [206]. All the aforementioned techniques consist in shaking the surface of the probed tissue, which implies that elastic waves propagate to the region of interest with sufficient amplitude. This can be an issue if the imaged domain is deep under the surface thus increasing the risk of high attenuation. Producing a wave field in situ may be an alternative way to ensure that a sufficient amount of displacement remains in the region of interest. In that regard, focused ultrasounds have been used to generate shear waves along with an MR scanner for motion detection [207]. This technique requires a heavy experimental setup compatible with the MR environment and has not been, to date, more than a proof of concept. Instead of using external devices to produce motion at chosen locations, the concept of intrinsic actuation taking advantage of natural internal vibrations has gained interest. This method consists in encoding motion induced in organs by the natural pulsation of the heart and arteries, and presents the significant advantage of not requiring any additional equipment. Tailored MRE protocols must be adopted to adapt to the low frequency characteristics of natural pulsations (around 1 Hz). For now, intrinsic actuation has been applied to the brain [208][209][210][211].

Acoustic Wave Detection
Whereas most magnetic resonance imaging protocols attempt to reduce or compensate for motion, MR elastography seeks to take advantage of small vibrations in the scanned tissue. Numerous MRE specific pulse sequences [chronologically sorted application of radiofrequency (RF)-pulses and magnetic gradients to generate and manipulate the MR signal] have been designed to acquire driven or natural motion in biological tissues, while maintaining reasonable scan times and image quality. We review, in this section, the main concepts of MR elastography pulse sequences allowing for detecting acoustic wave propagation. More specific details and theory, along with fast acquisition strategies are available elsewhere [19,212]. Motion encoding principles in MRI were first introduced by measuring sea-water velocity [213], and further applied in the context of angiography to measure blood flow [214]. The proof consisted in relying spin velocities to the phase shift spins experienced when space and time varying bipolar magnetic fields are applied. A similar concept, leading to MRE, was developed in which motion of spins around their position at rest is encoded in the phase of the complex MR signal using magnetic-field gradients, named motion encoding gradients (MEGs) [166]. The accumulated net phase of moving spins varies according to their trajectory while a time dependent MEG is applied. This net phase thus allows to track local motion of tissue eventually providing the displacement maps required for retrieving mechanical parameters. MRE sequences are generally based on existing MR encoding outfitted with MEGs. Not only are their setting key to be compatible with the characteristics of the induced motion but so are other inherent MR imaging parameters, leading to a broad variety of MRE sequences.
The timing of the chosen MR based-sequence rises a certain amount of constraints regarding scan time, motion sensitivity, and image quality. The MR sequence design and underlying physics are beyond the scope of the present review, thus only the main concepts relevant to overall MRE understanding is briefly discussed. For in-depth details, see [19,215]. When a tissue is placed in the strong static magnetic field B 0 of an MR scanner, a net magnetization aligned with the magnetic field is produced from the contribution of each uncoupled individual nuclear spin (those of hydrogen nuclei in clinical scanners) [215]. The key concept in MR signal generation consists in FIGURE 5 | Workflow in MRE investigation. Small time harmonic deformations are generated in the scanned tissue using an external actuator or natural body pulsation. These small deformations are recorded by the MR scanner using tailored acquisition sequences to produce time series propagation images. The complex MR signal is then processed to extract the motion field data, which is finally used as an entry to physical modeling allowing for calculation of mechanical parameters.
Frontiers in Physics | www.frontiersin.org June 2021 | Volume 9 | Article 666192 tipping the net magnetization of the desired portion of the scanned tissue out of its resting state using some excitation radiofrequency pulses tuned at the Larmor frequency of the spins in the scanner [216]. Magnetization enters a precession motion about the B 0 axis under the effect of this RF-pulse. Once excited, the magnetization is no longer aligned with B 0 and tends to realign and reach its resting state back again. This process called relaxation occurs at a certain time rate dictated by the interactions between spins themselves and with their environment. The time constant T 1 characterizes the exponential regrowth of the magnetization parallel to B 0 (longitudinal component of the magnetization) due to spinlattice interactions. Time constants T 2 and T 2* characterize the exponential decay of the transverse magnetization component (perpendicular to the longitudinal one) due to spin-spin interactions and magnetic field inhomogeneity (combined with the spin-spin interaction), respectively. Timewise, T 2* < T 2 < T 1 and decay rates are given by R 2* 1/ T 2* and R 2 1/T 2 , where R 2* > R 2 . The time evolution of the MR signal immediately following the RF-pulse excitation is named free induction decay (FID) and is governed by T 2* effects. Receive coils are used to monitor relaxation by measuring the voltage induced by the precessing magnetization according to Faraday's law. Manipulation of the magnetization allows to generate MR signal peaks at adjustable delays after the application of the excitation RFpulse, i.e., during the FID, or later when the signal appears to have vanished. MR sequences may be divided into two main categories of mechanisms leading to different timing for data acquisition. Spin echo (SE) sequences employ a second RF-pulse called refocusing RF-pulse, occurring after the FID, and allowing for compensation of T 2* effects (magnetic field inhomogeneity). Some of the MR signal can thus be recovered after the FID. In this case, the limiting time constant becomes T 2 > T 2* . Passed the exponential decay due to T 2 effects, the MR signal can no longer be recovered. The peak of the recovered MR signal, the echo, occurs at the "echo time" TE after the application of the RF-pulse excitation. Gradient recalled echo (GRE) sequences, however, typically operate within the FID (occurring immediately after the RF-pulse excitation) and do not allow for magnetic field inhomogeneity effects compensation. The operating window in such sequences is thus limited by the T 2* weighted decay, and thus leads to much faster acquisition protocols. Magnetic field gradients are used, instead of a second RF-pulse, to manipulate the magnetization and generate a signal echo at TE. As aforementioned, MRE sequences usually consist in incorporating motion encoding gradients into an MR basedsequence. This modification is consequently subjected to timing limits of the base-sequence. The short timing of gradient echo type sequences presents a narrower time slot for the MEGs to operate than that of spin echo type sequences. The impact of such inherent characteristics is discussed below.
The first descriptions of the motion encoding mechanism in MRE were reported in Refs. [166,217]. As aforementioned, a spin moving in the presence of a magnetic-field gradient G experiences a phase shift ϕ: where c is the gyromagnetic ratio of the material [rad s −1 T −1 ] and r is the time-dependent position vector of the spin. From this equation appears that the phase shift depends on both the spin trajectory r and the applied G. Consequently, a given arbitrary spin motion results in different accumulated phases depending on the magnetic gradient waveform. Hence, the remaining definition of G sets the type of motion the encoding process is sensitive to. Since the inherent function of magnetic field gradients is to add a controlled spacedependency to the static and homogeneous magnetic field B 0 , even static spins experience a space dependent phase accumulation while MEGs are switched on. In order to cancel this unwanted phase accumulation, G can be set to oscillate in time allowing the phase accumulated during the first half of the gradient oscillation period to be compensated during the second half. This technique is called zeroth moment nulling [19]. Non-oscillating MEGs are called unbalanced gradients and are thus rarely used in conventional MRE sequences. Additionally, the effect of constant velocity and constant acceleration background components in moving spins may also need to be cancelled. This can be achieved by applying first and second moment nulling, respectively [19]. Both consist in adjusting the MEGs oscillation profile so that the accumulated phase in Eq. 21 goes to zero for unwanted spin motion. Many MRE applications have resorted to full wave encoding (MEGs tuned to the same frequency as that of the motion oscillation) with zeroth [188,197,198,[218][219][220][221] and first moment nulling [166,172,173,175,181,203]. For all types of oscillating gradients, the area under the curve of MEGs over the operation time must equal zero for proper motion encoding. Early in MRE, several examples of motion encoding strategies were derived in the case of timeharmonic excitation, leading to time-harmonic spin trajectories and involving full-wave encoding [217]. Spin trajectories r can then be written as: where r 0 is the position vector of spins at rest, ξ 0 is the spins displacement amplitude, k the wave vector, r the position vector, ω the oscillation frequency, and θ some initial phase offset. Solving Eq. 21 using Eq. 22 and the MEG temporal profile allows to quantify the encoding efficiency, which is defined as the amount of phase shift in the signal per displacement unit. Encoding efficiency formulas for common MEG waveforms are available in Refs. [19,212]. Full wave encoding scheme ensures a good motion sensitivity but compels the minimum achievable value of the echo time dependent on the driver actuation period. GRE sequence short timing due to T 2* effects is well suited to high actuation frequencies (short time slot for MEGs to operate); however, it limits the applicable number of MEG cycles [166,188,197,220,222,223]. The optimal setting resides at the trade-off between the SNR increase permitted by short TEs and the higher motion sensitivity permitted by multiple MEG cycles. A similar conclusion can be drawn with regards to SE sequences, which present a more flexible timing enabled by their inherent longer echo and repetition times [215,224]. Multiple MEG cycles can thus be incorporated into the sequence while maintaining the echo intensity sufficiently high. SE sequences are often implemented with fast readout strategies, for instance echo planar imaging (EPI) necessitating only one or few combinations of excitation-and refocusing-RF pulses to generate a whole image, which allows to circumvent the use of many long repetition times (TRs) [172,175,181,203].
So far, single actuation frequency cases have been presented. Full-wave encoding can also be used to extract a frequency of interest from a multi-frequency oscillating wave field by selecting proper MEG profile (frequency and number of cycles) [177,222,225]. This configuration presents little interest in cases where the actuation frequency is chosen by the user. The only way of performing multi-frequency acquisitions using full-wave encoding is to repeat the encoding for each frequency separately [188,191,198,201,219,223,226], which has an impact on the acquisition time. However, multiple frequency components can be simultaneously encoded in a single acquisition using wide band MEGs, and manually selected using a temporal Fourier transform [177,178,180,182,210,[227][228][229][230]. The main advantage of simultaneous multi-frequency encoding is the time saving making them more suited to in vivo studies compared with repeated single frequency acquisitions over a given frequency range. The main drawback is the overall lower motion amplitude at each frequency of the multi-frequency actuation compared with the repeated acquisition scheme, due to total energy deposition divided into the total number of frequencies.
The wideband property of MEGs has been further extended to fractional encoding where the frequency of the mechanical oscillation is smaller than that of MEGs, and the mechanical time period is larger or equal to the repetition time [186]. With shorter MEG time periods, scan duration can be reduced and higher SNR can be obtained by shortening the echo time accordingly [190]. A major advantage of this approach is found in measurement of low frequency induced motion, such as heart pulsation driven actuation (around 1 Hz), where full-wave encoding would lead to unpractical echo times [208,209,211,231]. Despite the lower motion sensitivity in fractional encoding, this method has proven successful using fast acquisition protocols (spoiled GRE and GRE/SE equipped with EPI readout strategy) [169, 170, 179, 190-192, 201, 209, 232-234]. Besides multi-frequency acquisitions, reduced TE and TR permitted by fractional encoding have also been exploited in balanced steady state free precession MRE [186,235,236], despite the original development circumventing the use of MEGs [237]. Although high phase-to-noise ratios were reported, this sequence type presents significant timing constraints (actuation frequency linked to TR), and non-linear phase accumulation between consecutive TRs leading to additional signal post-processing steps. It has consequently been used only sporadically [19,190,212].

Inverse Problem in MRE
The previous section reviewed some acquisition approaches to measure motion induced in the tissue of interest. The last essential step in elastography consists in relying these displacements to mechanical parameters using physical models. This section addresses the most reported inversion schemes employed in MRE. More specific information about processing times and modeling details can be found in Refs. [238,239]. A major strength of magnetic resonance is the capacity of encoding motion in the three directions of space, allowing for full 3D inversion of the Navier equation. This strength comes at the cost of overall longer scan times for which alternatives have been discussed above. Despite the availability of fast 3D MRE sequences, all mechanical parameter reconstruction methods do not make use of complete displacement data sets and take advantage of physical assumptions allowing for processing of reduced dimension displacement data. We propose to classify inversion schemes into two categories. The direct approach consists in formulating the inverse problem with the mechanical parameters as unknowns. Experimentally obtained displacement data are inserted into the equations of elasticity, and quantities of interest are extracted through direct inversion. The iterative approach consists in iteratively solving the forward problem for displacements starting from an initial set of guessed mechanical parameters. These mechanical parameters are iteratively updated to minimize the difference between experimental displacement data and computed displacement solution. The final solution is the set of mechanical parameters that makes that difference converge to a global minimum.

Direct Methods
The first reported inversion method in the context of MRE, assuming isotropy, local homogeneity, no attenuation, and incompressibility, consisted in estimating the local wavelength of the measured wave field. This technique is termed LFE (local frequency estimation). Briefly, pairs of filters centred on spatial frequencies usually separated by one octave are applied to the wave field. The ratio of displacements filtered by each filter of one pair equals the local wavelength [240]. To ensure that local spatial frequency is included in the bandwidth of the filter pair, the process is repeated over a certain range of frequencies. From the evaluated wavelength (inverse of spatial frequency), the magnitude of the shear modulus |μ| is retrieved using μ ρv 2 s ρ(λ s f ) 2 , where λs is the local wavelength and f is the temporal actuation frequency. We recall here that ρ is the tissue density, vs the shear wave speed, and that such assessment assumes a purely elastic tissue (no viscosity). The original publication describing LFE [225] employed log-normal quadrature filters but other functions have been studied [241,242]. This method has been widely used in all types of study [166,169,171,174,175,221,223], as it is fast and only requires a single component of the displacement field. LFE has also proven a certain robustness against noise as it does not directly compute spatial derivative of the image, thus circumventing noise enhancement. To date, LFE is the only reconstruction method used and marketed for routine clinical practice. Although the LFE in itself provides no insight into the viscous behavior of the investigated tissue, this method has been combined with an attenuation model to estimate both G′ and G″, thus avoiding calculation of 2nd or 3rd order derivatives [223]. Phase gradient methods allow for simple estimation of the wave number k, similarly to LFE, which is used to quantify elasticity only. They have been sporadically used given their insensitivity to wave attenuation and their dependency on planar waves [243].
Early in MRE were also reported direct methods assuming viscoelastic materials, as described in General Concepts in Shear Wave Elastography, and using the strong formulation of the Helmholtz equation [198,220,244]. The underlying assumptions of isotropy, local homogeneity, and incompressibility are used to neglect the stiffness gradient across the tissue, and to decouple motion components in the equation system. From there, a single motion component can be used to retrieve the complex-valued shear modulus (G* G′ + jG″). Planar assumption allows to consider the 2D curvature of the wave field instead of its 3D, which further decreases the required amount of data to solve the inverse problem. This method has also been widely used [177, 178, 182, 219, 227-229, 243, 245] given its simplicity and low computational cost but strongly depends on data filtering and evaluation of second derivatives [244,246]. Using this scheme, stiffness reconstruction was shown to be altered by the neglected first Lamé parameter (Eq. 8) [247]. Applying the Helmholtz decomposition to the wave field allows to separate divergencefree (shear) from irrotational (compressional) components. Taking the curl of the Helmholtz equation increases the differentiation order but physically isolates the shear component of interest. This approach has become prominent when the Helmholtz equation is employed to retrieve storage and loss moduli [172,187,189,190,198,199,232,247,248]. To improve resolution and stabilize the direct inversion of the monofrequency Helmholtz equation, a multi-frequency approach named MDEV (multi-frequency dual elasto visco inversion) was introduced [200]. A multi-frequency wave field is built upon data sets of individual different frequencies, ignoring the dispersion of mechanical parameters with respect to frequency, reducing the risk of nodes due to standing waves and stabilizing the equation system by adding equations with same unknown. Various studies have resorted to inversion schemes based on this method [8,193,194]. As an alternative to the 2nd order derivative assessment required by MDEV, a phase-gradient based method termed k-MDEV was proposed. It consists in evaluating the complex wave number k of a plane wave (see General Concepts in Shear Wave Elastography), which can then be related to both phase velocity for elastic modulus estimation and attenuation for viscous behavior quantification [194].
A finite-element (FE) based inversion method was recently proposed, also assuming local homogeneity, for storage and loss moduli assessment. It takes advantage of the weak form of the equations of motion to reduce the differentiation order, and exploits divergence-free test functions to lower the impact of the compression field [249]. So far, most of the discussed approaches have in common the assumption of local homogeneity, neglecting the gradient of mechanical parameters across the tissue, and incompressibility, invoked to neglect terms involving the divergence of displacements (unless the curl operator is applied). The local homogeneity assumption was shown to alter the reconstructed mechanical parameters in regions where the latter are not constant [218], that is in most clinical cases and notably tumorous tissues. Direct methods, still employing the strong form of the Navier equation and neglecting the divergence of the wave field, were proposed to consider heterogeneity using single [250] and multifrequency (HMDI-heterogeneous multifrequency direct inversion) [251] approaches. To address the compression aspect in nearly incompressible materials, which was shown to lead to artefacts and inaccuracies [252,253] when disregarded [253] or processed using displacement formulation only [254,255], direct FE formulations of the inverse problem have been proposed using curl-based [256] and mixed displacement-pressure [253,256,257] schemes.

Iterative Methods
Overall, iterative methods make less restrictive assumptions on tissue mechanical properties than direct ones relying on the algebraic strong formulation of the elasticity equations, and have been reported to solve for more unknowns than FE-based direct methods by adjusting the number of parameters to update in the minimization process. From a computational standpoint, solving forward problems, that is, mapping data information (i.e., displacement field) from source information (i.e., elasticity distribution) is a smoothing process. On the contrary, inverse problem consisting in mapping source information (i.e., elasticity distribution) from data information (i.e., displacement field) is a noise-enhancing process [258]. Consequently, iterative approaches tend to be more robust against noise than direct ones. Following FE discretization approaches, similar to direct FE ones, iterative schemes have been developed [259][260][261]. Near incompressibility is often assumed and requires to modify the formulation of elasticity equations in order to solve for pressure in addition to displacements (aforementioned mixed "pressuredisplacement" formulation), and material heterogeneity is mostly considered [257,262]. In MRE, the subzone technique has gained significant interest amongst iterative processes [263]. It consists in dividing the imaged domain into overlapping subdomains termed subzones, and solving iteratively the forward problem for displacements in each subdomain parallelly [197]. Once the solution in each subzone has converged, subzones are randomly redistributed over the domain and the iterative solution calculation is performed again. Retrieved mechanical parameter distributions, corresponding to each subzone distribution, are finally averaged to form the final solution. This reconstruction method has been applied to phantoms, brain [176,203,264,265] and breast [266,267] data, and has proven its capacity to reconstruct multiple variables at various actuation frequencies using elastic and viscoelastic physical models (compressible elastic [268], compressible viscoelastic [269], and nearly incompressible viscoelastic [176, 184, 203, 231, 264-266, 268, 270-272]). Additionally, poroelastic models have been introduced for accurate consideration of the biphasic nature of entangled solid-liquid structures in biological tissues [168,184,211,270,271,273]).

Applications
Liver MRE can be applied to virtually any organ provided sufficient displacement data quality and suitable inversion scheme. We mainly restrict our discussion to liver applications due to the clinical availability of the technique, and to brain of which in-vivo mechanical properties have yet only been non-invasively accessible using this elastography method. MRE in clinics has so far been restricted to liver scanning for fibrosis, and diagnosis of chronic liver diseases using purely elastic models, i.e., assessing the shear modulus only and ignoring the tissue's viscous behavior. Under these assumptions, meta-analyses over the past few years on liver MRE have highlighted the high performance of the method in distinguishing liver fibrosis stages in non-alcoholic fatty liver disease, and cirrhosis considering the stiffness increase under this condition [274,275]. Viscoelastic parameters in human liver diseases have also been early investigated using multi-frequency MRE to assess the frequency dispersion of storage and loss moduli (G′ and G″) in healthy and fibrotic patients [180]. Results showed an increase in both G′ and G″ in fibrotic with respect to healthy livers. Evaluation of the elastic modulus and viscosity using a standard linear solid model in healthy and fibrotic human livers led to a similar conclusion, where both elasticity and viscosity increased in pathological liver tissues [178]. This trend was also reported in [276], and in a performance study of MRE in the detection of fibrotic livers [277]. Overall, stiffness only or both stiffness and viscosity increases have been observed against the fibrotic stage. Interestingly, a case of liver steatosis in rats where only viscosity varied while stiffness remained unchanged has been reported in [278]. Despite these findings, the storage modulus was found to correlate much better with stages of liver fibrosis than viscosity [279]. Additionally, measurement of wave damping for viscosity characterisation is influenced by reflections off boundaries and renders its measurement troublesome. To date, stiffness variations for estimation of liver fibrosis severity has been mostly investigated.
Brain MRE has also been proven successful and robust in the brain [280]. Its high water content and the observed shear wave attenuation in MRE acquisitions suggest that restriction to purely elastic models may lack of accuracy. The healthy brain's viscoelastic behavior has been highlighted by evaluating the dispersion of reconstructed mechanical parameters at varying frequencies in humans [8] and rats [198]. Both storage and loss moduli tended to increase with frequency. Additionally, cerebral viscoelasticity was shown to follow a frequency power law, where all reconstructed parameters vary independently [203]. These reconstructions suggested that the falx cerebri's viscous behavior is singular in comparison with other brain regions. Along similar lines in healthy brain characterization, viscoelasticity changes due to physiological aging have been considered [227,229,281]. From these studies appear that the brain softens but sees its relative viscous-to-elastic behavior unchanged over time. Such investigations along with MRE of neurological diseases have underlined the high potential of this technique in detecting neurodegenerative pathologies [282]. Non-invasive differentiation of natural structures of the brain based on their mechanical response to stimulus imparts MRE a significant advantage. For instance, the cerebellum has been shown to be softer and tends to be less viscous than the cerebrum [189].
High resolution mapping of stiffness and dispersion effects have suggested that cortical white matter is stiffer and more viscous than grey matter [191]. MRE has also been used to quantify the viscoelastic changes of altered brains. Notably, Alzheimer's disease was shown to reduce the brain's stiffness (elasticity only) [181]. Glioblastoma has been shown to take lower stiffness and viscosity values using multi-frequency MRE in humans [201], and mono-frequency MRE in a rat model [283]. A similar softening trend was observed in multiple sclerosis [182,228,284], where viscoelasticity was assessed using a global parameter. On the other hand, normal pressure hydrocephalus appears to trigger the opposite effect [175,211]. Recent research on brain viscoelasticity has opened new avenues in the understanding of connexions between cerebral functions and tissue mechanical behavior. For instance, joint investigation of hippocampus viscoelasticity (shear stiffness and damping ratio), and relational memory has allowed to correlate hippocampal viscoelastic variations to performance in completion of spatial reconstruction tasks [285]. Viscoelasticity was characterised using an adjusted damping ratio that indicates the dominant tissue behavior between elasticity and viscosity. Results showed that better relational memory performance correlated with a rather elastic mechanical behavior of the hippocampus. This constituted the first observation of the kind.
The same principle was applied to assess the correlation between cardiovascular health through aerobic fitness exercises, relational memory performance through spatial reconstruction tasks, and hippocampal viscoelasticity using MRE. The study showed that better memory performance was associated with higher values of the adjusted damping ratio, which was itself associated with better aerobic fitness performance [286]. Light fitness exercise has also been shown to have a potential impact on hippocampal viscoelasticity and associated cerebral functions in multi-sclerosis patients [287]. These investigations laid the first stone for the characterisation of relationships between physical and cerebral functional behaviors, and brain viscoelasticity [288][289][290][291].

Brain Anisotropy and Poroelasticity
Finally, most advanced improvements in viscoelasticity characterisation embed tissue anisotropy, which is particularly relevant in the brain given its fibrous structure. As a deviation from brain applications: the first use of anisotropy in MRE was proposed for breast tumor detection through the evaluation of an assumed symmetrical stiffness tensor [292]. Results suggested that carcinoma have an anisotropic structure revealing a preferred orientation, certainly due to vascularisation, and suggesting transverse isotropy. Breast cancer was then also characterized assuming a transversely isotropic model leading Frontiers in Physics | www.frontiersin.org June 2021 | Volume 9 | Article 666192 to a 5-parameter reconstruction [247]. Again, results showed a preferred orientation in the tumor structure. Transversely isotropic mechanical property recovery was experimentally validated in fibrous tissues using MRE and diffusion tensor imaging (DTI) [293]. In-vivo brain anisotropic stiffness assessment in humans assuming separately orthotropy and transverse isotropy suggested that white matter exhibits a transverse isotropy structure [20]. Shear wave speed analysis was performed from prior knowledge of fiber orientation using DTI. Two shear wave modes were then observed, a faster longitudinal mode relatively to a slower transverse mode. The transverse anisotropy of white matter was later nuanced in favor of a mild only anisotropy in a study on exvivo porcine brain, where a purely transverse shear wave mode was generated and used to estimate three isotropic parameters in the absence of longitudinal modes [294]. Human brain anisotropy was also highlighted using variations in reconstructed stiffness distribution depending on actuation direction [176]. This constitutes a direct measurement of the anisotropy impact in isotropic models. Such observation was quantified using a finite element formulation of a heterogeneous, nearly incompressible, and transverse isotropic model providing benchmark displacement fields for inversion testing [265]. In addition to the significant research effort in evaluating and understanding cerebral viscoelasticity, the high water content of the brain has motivated to consider it as poroelastic, i.e., made of two, solid and liquid, entangled phases. The impact of poroelasticity versus viscoelasticity on reconstruction has been shown to be relevant at low frequencies (a few Hertz) using the forward problem formulation in the harmonic regime, and the aforementioned subzone iterative scheme [184]. At higher frequencies, viscoelasticity seems to remain a more suitable model than poroelasticity. Overall, poroelasticity and low frequency intrinsic actuation thus constitute an interesting and original package in MRE investigation. This setup circumventing resorting to pulsing equipment has been used in a few studies, and holds promise for more accurate detection of brain pathologies [208,211]. Another approach to highlight brain poroelasticity was to solve for both shear and bulk moduli using the algebraic inversion technique, which resulted in a bulk modulus much lower than expected, confirming the poroelastic nature of the brain (compressible solid matrix and incompressible fluid channels) [209]. More recently was proposed an improvement in MR poroelastography acquisition processes allowing to separate solid and fluid contributions to the shear motion field using an inversion recovery sequence adapted to MRE, along with a tailored MR signal modeling [295]. To conclude, Figure 6 illustrates typical wave maps and elastograms from MRE acquisitions in the liver, breast, heart, and brain. Table 3 presents an overview of main components constituting MRE investigations, from motion generation techniques to inversion categories described in previous paragraphs.

OPTICAL SHEAR WAVE ELASTOGRAPHY
Elastography based on ultrasonography or MRI has found popular clinical applications facilitated by the implementation of those imaging technologies on clinical systems. These tools can provide images over centimeter to whole-body depth ranges.
However, many applications require millimeter-scale spatial resolution images, which can only be made possible using optical means. For example, on the cellular scale, the measurement of mechanical properties requires higher resolutions to focus on the understanding of how cells respond to physical forces. Thus, the use of optical elastography provides an opportunity for microscale imaging and for numerous applications in fundamental research [298].
Within the last 2 decades, developments in this new area of imaging led to multiple scientific advancements at the interface between optics and mechanics, which included biomedical applications in ophthalmology, oncology, and cell mechanics. The following subsections discuss recent developments in cellular and optical elastography, and their applications across biomedical and life sciences.

Cellular Shear Wave Elastography
Tissue elasticity at a microscopic scale is determined by the cell and the extracellular matrix elasticity. Main components of a cell are the membrane, cytoplasm, and cytoskeleton. The latter structure contributes to the cell mechanical stability and characteristics, and to its morphology. An imbalance in the mechanical homeostasis and defect in the cellular Brain [8, 172, 175-177, 181-184, 191, 195, 198, 201, 203, 209, 219, 227-229] [ 184,208,210,211,223] [ 8,177,182,184,198,203,208,210,211,219,[227][228][229] [ 8,172,175,176,181,183,191,195,201,209,223] [ 8,177,182,184,191,195,201,209,210,[227][228][229] [172, 175 [172, 175-177, 181-184, 209, 223, 227-229] [191] [ 8,195,198,201,203,219] [208-211] [8,172,175,177,181,182,191,195,198,201,209 Breast [185] [ 197] [ 185,197] [ 185,197] [185] [ 197] [ 185] [ 197] Spleen [193,194] [ 193,194] [ 193,194] [ 193,194] [ 193,194]  mechanotransduction can contribute to various human diseases. Therefore, cellular viscoelasticity can be viewed as a biomarker for determining the cellular state [299]. Determining mechanical properties of a cell during different stages of the disease progression could help to develop novel treatments by considering the role of mechanical factors into genetic and drug therapies [300]. As recently reviewed [301], there are several techniques for measuring cell mechanical properties, most of them requiring a mechanical stress (e.g., micropipette aspiration and atomic force microscopy). However, in this review, the focus is on rheological properties assessed using a mechanical stress based on acoustical shear wave propagation. Grasland-Mongrain et al. [302] developed a novel method called optical microelastography, also labeled as "cell quake elastography". This technique uses a high frequency shear wave excitation and an optical microscope to assess cell elasticity. High frequency shear waves inside the cell are produced by a vibrating micropipette at a wavelength comparable to the cell's size. The wave propagation is captured optically by a high frame rate camera coupled to the microscope. The sampling rate of the camera is selected to avoid shear wave frequency aliasing with sufficient samples per wavelength to allow efficient speckle tracking. The spatial resolution of captured images should also be sufficient to track the shear wave speed from displacement maps (knowing the time elapsed between images). The proof-of-concept in [302] was made by using an ultrasound speckle tracking method adapted to optical images for obtaining displacement maps [303]. A passive elastography algorithm was used as a reconstruction method to obtain shear modulus images [304]. The passive elastography method was inspired by the seismology field [210], so the name "cell quake elastography" for this method. The main advantage of this technology compared with other cell elasticity methods is the time resolution of a few microseconds to produce elasticity maps with a good spatial resolution, which may allow studying dynamic cellular processes. First experiments were performed on mouse oocytes, making a shear wave inside the cell by a 15 kHz vibrating micropipette, and capturing the traveling wave optically at a 200,000 acquisition frame rate [302]. The sensitivity and spatial resolution of the technique allowed to distinguish the shear modulus of different regions/zones of a cell (Figure 7). The technique was recently applied on mouse macrophage-like RAW 264.7 cell clusters, Abelson leukemia virus-transformed cell line derived from mice, using a 18 kHz stimulation and a 100,000 frames per second image capturing rate. Shear wave displacement maps at 10 µs intervals are given in Figure 8. Future developments should aim at assessing the viscous component of single cells, likely using finite-elements modeling (FEM) reconstruction methods.

Optical Coherence Elastography Imaging
Optical coherence elastography (OCE) is a technique that can non-invasively assess tissue mechanical properties by measuring the localized deformation, strain or shear wave propagation properties inside a sample [305]. In OCE, a stimulation technique is utilized to load the tissue and its response is recorded with an OCT based detection method [306][307][308]. The high resolution structural images of OCE (1-10 μm invivo) provides it an advantage over the ultrasound or MRI modality [309], that stretches its potential for micron and submicron imaging of elastic properties of biological tissues. Microstructures of biological tissues can be quantified based on optical scattering properties of the tissue under investigation. OCE holds great potential for diagnosis of many clinical conditions and pathologies, particularly for detection and monitoring of cancers [310], cardiovascular diseases [311], and eye diseases [312].
While optical contrast signals are detected based on differences in two or multiple optical scattering events, the mechanical contrast requires only one scattering event to obtain an OCT signal. Thus, structural inclusions that cannot be detected by OCT can be revealed by OCE if a mechanical contrast exists for the inclusion. The first few studies in OCE development focused on static mechanical contact loading (i.e., no shear wave involved) [313,314]. Later, the emergence of phase resolved OCT, which is detecting the interferometric phase information from complex OCT signals, enabled the assessment of tissue deformation with a high accuracy for tissue elasticity reconstruction [315][316][317][318]. A shear wave stimulus was involved in studies of [316] and [318].
The latest developments include OCE resolution to improve over the range from several microns to hundreds of microns [308,319,320]. The lowest range of OCE spatial resolution is similar to the cell quake elastography imaging method described earlier. In comparison, the spatial resolution of ultrasound or MRI elasticity imaging methods remain at a macroscopic level with a typical resolution of hundreds of micrometers to several millimeters, respectively [54,191]. OCE is a great alternative to traditional elastography methods in terms of spatial resolution, acquisition speed, sub nanometer mechanical displacement sensitivity, but at the cost of a lower penetration depth into the probed tissue than ultrasound or MRI [308]. Additionally, shear wave OCE as a 3D imaging modality may enable its clinical applications in many areas, such as ophthalmology and cardiology using intravascular devices [321][322][323]. Shear wave based OCE has shown potential for measuring local elasticity changes of mouse brains [324,325]. Details on these methods are given next.

Systems and Methods
An OCE system has two main components: a loading system that can deform the biological tissue, and an OCT imaging system for detection. Shear wave methods in OCE are relatively in the very early stages of development. Shear waves-based OCE utilize an excitation from a noncontact air-puff or air-coupled ultrasonic probe [326][327][328], or piezo-transducers (PZT) [320]. In addition, an OCT mechanism is then employed to detect the displacement field of generated shear waves. By monitoring the shear wave propagation in the sample, elasticity, shear wave speed, or the shear modulus can be quantified. Shear wave visualization was performed in tissue mimicking phantoms with phase sensitive optical coherence elastography [329]. Razani et al. [318] were one of the first to measure the shear wave speed and its associated properties with OCT phase maps. They utilized an external acoustic radiation force mechanism for excitation and a sweptsource OCT system to acquire phase images. The central wavelength of the laser was 1,310 nm and the bandwidth was ∼110 nm. The system could register a lateral resolution of 13 μm in gelatin mixed with titanium dioxide phantoms. Images could be acquired at a depth of 3 mm. Song et al. [320] used a piezoelectric point loading to generate shear waves within samples. More recently, Zhu et al. [330,331] developed a PZT-based system to induce longitudinal shear waves and they visualized the signals using OCT for the quantified mapping of shear moduli. A brief detail of their technique is presented next.
The OCE system included an OCT imaging unit and a PZT excitation unit, as shown in Figure 9A. Elastic waves were induced by a ring PZT actuator driven by a PZT amplifier. The vibrating mechanism of the PZT system could excite three types of waves in the sample under investigation: 1) Rayleigh waves, 2) compressional waves travelling from the top surface to the deep region, and 3) transverse and longitudinal shear waves traveling through the interior of the sample, as shown in Figure 9B. Rayleigh waves propagate at the surface of the sample. Compressional waves propagate parallel to the oscillation direction of the vibrator. Transverse shear waves propagate perpendicular to the displacement direction.
Additionally, in the near field of the planar vibration source, which contained multiple sub-sources, a longitudinal shear wave much slower than the compressional wave also propagate along the displacement direction. This longitudinal shear wave is present due to the sum contributions of diffracted transverse shear waves [331]. These longitudinal shear waves could be visualized with the attached OCT imaging unit. The OCT system was based on a swept source at a central wavelength of 1,310 nm, and a wavelength tuning range of 141 nm. Axial and lateral resolutions of the employed OCT unit were 7.6 and 17.7 μm, respectively. The PZT unit utilized for excitation was driven by a function generator producing a sine wave cycle with a frequency of 1 kHz. The displacement observed in the near field was close to 10 μm.
As introduced above, noncontact shear wave imaging optical coherence tomography (SWI-OCT) system has been developed using a focused air-puff device for localized tissue deformation [333]. The non-contact mechanical excitation in a sample could be performed with a PZT transducer that was specially designed to launch an US beam through air that was focused onto the air-medium interface. The reflection of the beam at this interface could produce significant acoustic radiation force toward the sample medium. This induced a transient displacement at the surface, including shear waves. The large difference in acoustic impedances of air and soft tissues could increase the efficiency of the acoustic energy conversion even to the extent of one hundred percent. Ambrozinski et al. [332] developed a non-invasive system that needed transient displacements to be only about 1 μm, and the acoustic pressure only a few kPa, which was within safety limits for clinical applications. This acoustic micro-tapping method had enabled 4D imaging of tissue stiffness by employing a focused air-coupled US to induce mechanical deformations at the boundary of a tissue [332]. A schematic representation of their system is shown in Figure 9C. In here, the cornea surface was aligned at the transducer focus and the US radiation push was sent through with a repetition period of 3 ms. The driving signal was having a bandwidth range of 0.95-1.05 MHz. The measured pressure amplitude at the transducer focus was about 7 kPa. Wang et al. performed a quantitative biomechanical characterization of cardiac muscles and corneas using a noncontact SWI-OCT system [334,335]. Shear waves had a frequency range of 0-2.5 kHz. This method employed a multiwave imaging technique, where shear wave measurements in the tissue enabled mapping of the mechanical contrast in elastograms, and the OCT unit enabled improving the imaging resolution from a millimeter scale to a micron scale [334]. The system was capable of simultaneously providing structural images with depth wise maps of the tissue stiffness [335]. Recently, a confocal air-coupled US probe could also be co-focused with a phase-sensitive OCT system to generate elastic waves up to a 4 kHz frequency for quantitative elastography [336]. These noncontact excitation methods have found wide applications in ophthalmology and dermatology [337][338][339].
Spatial resolution in dynamic shear wave based OCE is governed by temporal and spatial characteristics of mechanical waves rather than optical waves. Hence, the mechanical resolution in dynamic OCE is different from the usual optical resolution of OCT systems [57]. Spatial resolution ideally should match the spatial resolution of the detection system, however, propagating mechanical waves undergo mode conversions at tissue interfaces causing artifacts in the elasticity image. The geometry of the tissue interface and its elasticity contrast can produce complex propagating fields near the tissue boundary affecting both the spatial resolution and contrast of the final reconstructed image [57].
Recent dynamic OCE systems provided elasticity information from local group velocity measurements [321,330,339], however, the complex geometry of bounded tissues like the cornea may not reflect a simple relationship between group velocity and elasticity [338]. Dynamic OCE has been successfully utilized in elasticity mapping of the cornea using noncontact excitation methods based on air-puffs and acoustic micro tapping [332,333,[340][341][342]. Inversion of moduli from experimental data, especially in the case of bounded and anisotropic tissues such as cornea, is a challenging and complicated process in dynamic elastography. Recently, a nearly-incompressible transverse isotropic (NITI) model addressed this challenge and characterized corneal biomechanics while accounting for corneal microstructure and anisotropy, and presented a more accurate model for cornea shear moduli computation [337]. Viscosity assessment in shear wave OCE is in its early phase of development. Proposed methods used shear wave frequency dispersion [343][344][345][346], storage and loss moduli using a rheological model [347], and the elastic wave attenuation [345].
A trade-off in OCE is its reduced depth of field while evolving for higher resolution measurements due to the requirement of higher numerical aperture for such systems. On the other hand, the ability to measure and record depth scans with a single spectral acquisition can be used as an advantageous feature to enable phase-sensitive displacement measurements. Of course, the tissue penetration attained with OCE, although sufficient for numerous applications, is not comparable to ultrasound or MRI elastography methods. Song et al. implemented a beam-steering US as a wave source for shear wave optical coherence elastography of retinal and choroidal tissues within a porcine eyes ball ex vivo. Shear wave propagation imaged on a porcine retina by their system is shown in Figure 10

Photoacoustic Elastography
Photoacoustic elastography (PAE) research is rapidly growing due to its potential and promising features of clinical interest [348][349][350]. PAE can exhibit a mechanical contrast in biological tissues while also providing high spatial resolution images and an excellent penetration depth compared to commercially available optical imaging modalities [351]. It has the promise to provide great scalability, ranging from cellular levels to entire body with multiple resolution levels. Recent studies have demonstrated recovery of mechanical properties of biological tissues using PAE [352][353][354][355]. Several studies demonstrated computation of elastic properties of soft tissues [354,[356][357][358][359]. Nevertheless, clinical translation of PAE is still far way for research studies to accomplish, the development of the PAE technology has shown the potential to be used in life threatening diseases, such as breast and prostate cancers, and brain tumors [350,360]. Photoacoustic elastography can be used for mapping elastic properties of diseased tissues with highly vascularized structures, such as carcinoma and glioblastoma [351]. Most PAE studies have focused on qualitative imaging and quantitative PAE is still a challenge. Moreover, PAE using propagating shear waves still need to be clearly addressed. A recent study by Wang et al. did develop a PA viscoelasticity technique for quantitative imaging of liver cirrhosis based on a PA shear wave model [359]. This viscoelasticity imaging model was inspired by the acoustic radiation force impulse (ARFI) technique (see Ultrasound Shear Wave Elastography). In this model, a laser beam was focused into a tissue that resulted in the tissue thermal expansion and a PA pressure field was generated. The pressure field induced a localized ultrasound impulse similar to ARFI, and subsequently a tissue displacement field could be observed. The study assumed that these forward propagating PA waves could be modeled using shear wave equations.
As a summary of methods addressed in this review,

Ultrasound Shear Wave Elastography
Many notable ultrasound elastography methods have been translated into clinical applications, and adopted by clinicians for diagnosis of several organs, as introduced in Applications. A limitation lay in the depth of SW penetration due to attenuation, especially for the diagnosis of liver fibrosis and steatosis, which could result in unsuccessful measurements with large patients or patients with ascites [31]. Note that SW attenuation is a concern for any shear wave elastography method. Consequently, measurements on superficial regions showed a higher success rate, such as the diagnosis of breast lesions and tumors. Another limitation lies in the assumption often used in shear wave elastography; notably considering the tissue as isotropic and homogeneous. Certain tissues such as muscles or tendons do not respect the isotropy hypothesis and are rather considered as anisotropic or transverse isotropic media. To answer this problem, teams have developed stiffness tensors for assessing shear wave propagation and for evaluating mechanical properties in several directions [161,361], even in three dimensions [162]. Bones, brains, or lungs are parts of the body that can be considered porous and for which the assumption of homogeneity is limited. Poroelasticity based on the estimation of the temporal response of tissues to compression [15,[362][363][364] is a technique derived from strain elastography. Although a little off topic because it does not use shear waves, its development in the characterization of tumors is promising [364,365]. Other applications, such as the characterization of muscles, Achilles tendons, the cardiovascular system, and lymph nodes [110,[366][367][368], have shown good results that reflected the difference between normal and abnormal tissues. At present, the measurement of the tissue elasticity has dominated the field, and technologies, such as the transient elastography, SSI, ARFI, and comb push ARFs are available on clinical scanners [67,110,144,[369][370][371][372][373][374][375][376]. In fact, most manufacturers have today a shear wave elastography package for clinical use, and the spatial resolution of those elastography systems are as good as ultrasound B-mode imaging. However for measurements on more complex tissues, such as anisotropic,

Magnetic Resonance Shear Wave Elastography
Magnetic resonance elastography has been proven successful and robust in a broad range of applications, from clinical diagnosis of liver diseases [386] to brain [280], and breast [387] pathologies. However, generalization of MRE in clinics is impacted by the scan time required for the acquisition of complete data sets necessary for accurate reconstruction. Clinical sequences currently acquire one component of the motion field at one frequency, and rely on a 1D direct inversion assuming isotropy, incompressibility, local homogeneity, and pure elasticity. This package is fast and guaranties results. Mechanical solicitation of the tissue under multi-frequency loads may be a first step into the characterization of tissue viscoelastic responses by providing valuable information on the frequency dependent biomechanics in pathological cases [201]. Motion encoding may also reach a limit at high frequencies in the case of oscillating gradients due to peripheral nerve stimulation [388]. The actuation regime also determines the physics to consider for inferring mechanical parameters from experimental datasets [184]. An elegant avenue circumventing the use of external actuators is intrinsic actuation from low frequency heart beats along with poroelastic modeling of the tissue [211]. Along the same lines, natural vibrations in the brain have been exploited using the novel passive elastography technique based on time reversal concepts to quantify the vibration wavelengths, assumed to be related to brain stiffness [210]. A potential solution to shear wave attenuation in soft viscoelastic tissues may be to place motion sources closer to the region of interest. This may be achieved using ultrasound transducers generating an acoustic radiation force impulse (ARFI) along with MRE acquisition [207,389]. Alternatively, a promising approach that may prove feasible with clinical MRI scanner is Lorentz force elastography for in situ actuation at different frequencies [390,391]. An elasticity reconstruction obtained using a Lorentz force and a clinical MRI scanner is displayed in Figure 11, in the case of a gel phantom. On a similar note, localized motion generation using ARFI has been employed with MRE to measure elasticity changes during high intensity focused ultrasound ablations in ex-vivo porcine muscle samples [392]. Such monitoring requires sufficient displacement amplitude [393]. Assessment of stiffness changes due to ablation or percutaneous procedures has been performed both during [394,395], and separately before and after ablation [396], all cases reporting a stiffness increase after the intervention.
On the acquisition side, significant amount of effort has been put into MR sequence developments to reduce scan time while preserving 3D motion encoding and signal amplitude. Although equipping an MR sequence with bipolar magnetic field gradients prevails, recent implementations took advantage of MR sequence inherent gradients to encode motion, thus keeping the timing shorter than conventional use of MEG [397]. As in any MRI scan, artefacts may occur due to patient motion. Sequence dependent artefacts include and are not restricted to signal loss in GRE sequences due to irregular geometries, and associated magnetic field inhomogeneity and distortion in echo-planar sequences. Specific to MRE, phase wrapping occurs when motion cannot be encoded in the [-π, π] range leading to phase jumps within this range. Three solutions appear and consist in either decreasing the gradient sensitivity, decreasing the motion generator strength, or using a phase unwrapping algorithm. Whilst a weakness of MRE may be viewed as a lack of universal protocol applicable to any organ, its strength resides in the capacity of providing usable data in multiple cases owing to various hardware and MR sequences, and in the availability of physical models to process produced experimental data for stiffness estimation. This technology is still mainly used in the context of clinical research, and additional validations might be required for robust viscosity, porosity, and anisotropy assessments.

Cellular Shear Wave Elastography
Although the optical microelastography technique has an unprecedented high temporal resolution with the capability of producing elasticity images, its spatial resolution is not as good as other rheology methods, such as atomic force microscopy. The spatial resolution of this technique is currently limited to 10 μm, approximately, but it could be improved to resolutions close to optical microscopy (≈1 µm) by utilizing a higher stimulus frequency and by improvements in the reconstruction process. This technique is currently limited to elasticity measurements but viscosity might become available through FEM modeling, or by considering shear wave attenuation [50].
More improvements need to be done beyond solving limitations mentioned above. The mechanical behavior of the cell should be further investigated in a range of frequencies. One remaining difficulty is the absence of good standard rheology method to validate the microelastography technique at frequencies in the kHz range. This makes it difficult to compare results from different techniques on a similar cell type. Other reasons for the lack of concordance of results, obtained with different cell elasticity technologies, were recently addressed [398] and apply to the reported microelastography method. Recent developments might allow using mechanically stable microgel bead to compare different cell elasticity methods [399], which may reduce variability when performing such comparisons at different (non-overlapping) frequencies.
Also, to make these techniques more applicable and practical for biologists, and to promote using cell mechanics as a biomarker, improvements are required for the technique to be automatic, high throughput while being robust, accurate, and sensitive with high time and space resolutions. This might be done by coupling the optical microelastography technique with other methods, such as microfluidics with a high throughput [400,401].

Optical Coherence Elastography
Higher resolution OCE may face computational challenges due to the fact that the speckle decorrelation length scales with the speckle size [402,403]. This would reduce the maximum displacement that can be measured between frames. Many studies have made progress to further demonstrate substantial improvements in resolution [404,405]. High resolution OCE systems can be used to assess mechanical properties of cells and a few preliminary studies have showed this potential [404,406]. These are relatively new developments and the hope remains that OCE would be able to characterize cell aggregates [305], with penetration depth going up to several hundred microns, whilst maintaining a sub-cellular scale resolution.
There has been several studies on elastogram image reconstruction in OCE by inverse problem approaches [407]. Sridhar et al. [408] used an inverse problem approach to understand how stromal tissues affect the broad spectrum of the viscoelastic response [409], by minimizing the mean squared error between computed and measured displacements. Different methods to constrain the optimization algorithm has been summarized by [410], in the context of ultrasound strain elastography. Basic principles are also applicable to OCE. However, one challenge that often prevails in these scenarios is the optimization of the regularization parameter for efficient reconstruction, especially in the context of in vivo experiments.
Another area of interest representing some challenges is the quantitative assessment of tissue viscoelastic properties with OCE. This research is still in its early years acknowledging the fact that the viscosity is not accounted for in the simple approach [411], but hopefully with the development of new models, OCE would be able to convert elastic wave speed and attenuation into quantitative values for clinical diagnosis based on tissue viscoelasticity.
Despite several advancements, very few studies have been done in the area of validation of performance. This would require phantoms that are developed for optics rather than mechanics [305]. Rigorous assessment of sensitivity and specificity for diagnostic applications would be required for translating the method to the clinics.

Photoacoustic Elastography
Photoacoustic elastography imaging is relatively a new development and it still needs to overcome many challenges. The PA signal contrast detected by ultrasound transducers is low due to lower variation in the tissue elasticity distribution in comparison to the optical absorption coefficient. This can potentially be overcome by additionally employing an ultrasound modulation of the laser pulse to provide external mechanical simulation of the tissue [412]. In addition, elasticity can also be estimated from the resonance frequency of the tissue material observed in the measurement of PA signal strength against the operating frequency of the external (ultrasound) mechanical stimulation [412,413]. Another challenge is the development of quantitative PAE imaging systems, as the first few studies in the field reported only qualitative assessment of elastic properties. However, Hai et al. were the first to develop a quantitative PAE system [414]. It would also be interesting to detect the contrast in the PA signal due to elastic property variation separately from that of other parameters (including the optical absorption coefficient). Grasland-Mongrain et al. generated shear waves in soft tissues in ablative and thermoelastic regimes with a 532 nm Nd:YAG laser [415]. However, it remained a challenge to keep the laser beam energy within safety limits for use in biomedical applications. This can potentially be overcome by use of other types of laser or by emission of the high energy laser beam onto a protective absorbing layer, such as a black sheet, that can cover the tissue externally. In conclusion, PAE is still in the beginning phase of its development compared to ultrasound, MRI or OCT elastography, and there definitely remains scope for many promising improvements to increase its potential for various imaging applications.

AUTHOR CONTRIBUTIONS
HL supervised and structured the review, wrote the introduction and the conclusion, contributed to the section on the generation and detection of shear waves, reviewed, and edited the final document. GF co-supervised and structured the review, wrote the magnetic resonance elastography and the general concepts on shear wave elastography sections, reviewed, and edited the final document. MB wrote the optical coherence elastography imaging and the photoacoustic elastography sections. ZQ wrote the characterization of tissue viscoelasticity section, and contributed to the generation and detection of shear waves section. SG wrote the cellular shear wave elastography section. LY wrote the viscoelasticity reconstruction section and applications in ultrasound. GB wrote the clinical applications of ultrasound shear wave elastography section. IR contributed to the section on ultrasound shear wave elastography. GC received the invitation to provide this review, contributed to the general supervision of all co-authors, did reviewing of all sections of this document for final approval, and provided editing and funding.