Existing and Potential Applications of Elastography for Measuring the Viscoelasticity of Biological Tissues In Vivo

Mechanical tissue properties contribute to tissue shape change during development. Emerging evidence suggests that gradients of viscoelasticity correspond to cell movement and gene expression patterns. To accurately define mechanisms of morphogenesis, a combination of precise empirical measurements and theoretical approaches are required. Here, we review elastography as a method to characterize viscoelastic properties of tissue in vivo. We discuss its current clinical applications in mature tissues and its potential for characterizing embryonic tissues.


INTRODUCTION
Viscoelasticity, among other mechanical properties, is intrinsic to biological tissue. The term implies that tissue exhibits time-dependent responses to an applied force [1][2][3]. Characterization of the viscoelastic behavior of biological tissue has been performed in vitro [4,5], ex vivo [6], and in vivo [7,8] to gain insight into tissue stiffness and fluidity. In vivo assessment is preferable to determine tissue properties in their native environment, and elastography has the advantage of minimally disturbing that environment. In brief, elastography introduces a disturbance to displace specific regions within a tissue, which is subsequently imaged and analyzed to determine the local viscoelastic response.
It is worth noting that elastography was not originally developed for the purpose of measuring viscoelastic properties. As the name indicates, in its earlier development, only elasticity (in terms of Young's modulus) was at the center of attention, lacking viscosity information (coefficient of viscosity) [9,10]. However, as the properties of biological tissue were progressively understood, a variety of theoretical frameworks were developed and integrated with earlier elastography techniques. In recent years, although the name of elastography remains largely unchanged (some studies have adopted the term "viscoelastography" [11,12]), it has become more common that quantitative values of various moduli, as we will discuss in detail in later sections, include information of material viscoelasticity.
In clinical settings, elastography techniques for measuring the viscoelastic properties of mature tissues have been well established in vivo to diagnose and distinguish among different pathologies [7,[13][14][15][16]. The motivation of shifting from measuring elasticity information to viscoelasticity information lies within the fact that as mature tissues exhibit varying degrees of viscoelastic behavior, the sufficiency and validity of solely relying on elasticity information have been questioned [17]. For instance, elastography-based studies on assessing hepatic tumor malignancy [15,18], benign/normal breast tissues [19], and the efficacy of cornea disease treatment [20] have reported a more prominent effect by using viscosity as the differentiator rather than elasticity. Thus, assessing the viscoelastic behavior of mature tissues is necessary to provide more accurate in vivo measurements and characterization, subsequently improving tumor diagnostic accuracy.
In studies of embryonic tissues, it is increasingly understood how forces generated by cells can coordinate morphogenesis [21][22][23]. In response to forces, embryonic tissues, which are both liquid-and solid-like, exhibit viscoelastic behavior [3,24]. Cells have receptors that can sense physical forces as well as chemical cues. For example, cadherin molecules have both mechanical and sensory properties which are applied for adhesion and signaling, respectively [24]. It is currently challenging to perform loss and gain of function experiments of mechanical cues to define their roles in vivo as one can with chemical cues. Although increasingly precise tools are being developed to characterize the mechanical and viscoelastic properties of embryonic tissues [25], data supporting their efficacy are currently sparse. Albeit less explicitly discussed, the motivation to study viscoelasticity in the context of developmental biology is twofold. At the single-cell level, various membrane-enclosed and membrane-less organelles exhibit viscoelastic behavior to different degrees. The latter can take on the form of liquid droplets that undergo controlled dissolution and condensation via phase separation [26,27]. The occurrence of irreversible aggregation promotes further transition of some condensates from liquid-like to solid-like [28], underlying the pathologies in many neurodegenerative diseases [29]. Thus, probing viscoelasticity at the single-cell level would potentially facilitate our understanding of how changes in compartmental viscoelasticities correlate to the phase change in condensates and overall cell behavior. At the multicell level within an embryonic tissue, gradients of stiffness measured by the Young's modulus were discovered as a cue that potentially guides cell migration by a process called durotaxis [30]. However, due to the omission of the viscous properties of ECM, durotaxis may need to be reexamined in the context of viscoelasticity to incorporate the potential role of viscosity in impeding the migratory speed. Viscoelasticity and Developmental Biology is dedicated to the in-depth discussion of existing and potential roles of viscoelasticity in the context of developmental biology.
For the purpose of this review, the common assumption made in elastography techniques, that is, biological tissues are isotropic and homogenous, is also assumed. The definitions of elastic stress (σ), whether compressive or tensile, and shear stress (τ) are illustrated in Figures 1A,B, respectively. A glossary of the

ESTIMATION OF MATERIAL MODULI
The classical approach to estimating moduli that reflect the viscoelastic properties of a material is implemented by performing one of the following three canonical tests: creep, stress relaxation, and oscillatory loading. Once the data obtained from the tests, whether time-dependent or frequency-dependent, is fitted with the constitutive equation of a selected rheological model, the coefficient of each parameter within the model can be determined as an approximate representation of the viscoelastic properties of the material. Alternatively, if the constitutive equation is solved in the time domain, the Fourier transform can be taken to derive the solution in the frequency domain [31].

Classical Rheological Models
A common approach to the study of material viscoelasticity is to derive the expression of creep compliance, relaxation modulus, or their complex forms, from the constitutive equation of a linear theoretical framework. Depending on the properties of the material, differently structured rheological models can generate different levels of fit. Expressions of the creep compliance (J) and relaxation modulus (F) are commonly derived in the time domain via creep and stress relaxation tests. If their frequency-dependent complex forms, (J p ) and (F p ), need to be determined, oscillatory tests are conducted.

Elements of Rheological Models
Classical linear rheological models consist of various combinations of different numbers of Hookean springs, which model elastic behavior, and Newtonian dashpots, which model viscous behavior. Schematics of a Hookean spring and a Newtonian dashpot are shown in Figures 1C,D, respectively. Upon an applied elastic stress (σ), a Hookean spring exhibits linear elastic behavior and immediately produces an elastic strain (ε). The elasticity can be quantitatively represented by the Young's modulus (E), which determines the stiffness of the material (with units of Pa), via Hooke's law as follows: A Newtonian dashpot exhibits a viscous behavior that is proportional to the strain rate (_ ε), which represents the change in strain with respect to time. The viscous behavior is representative of the material fluidity and is expressed (with units of Pa s) via the coefficient of viscosity (η) as follows: Common two-element rheological models include the Maxwell model and the Kelvin-Voigt model. The Maxwell model consists of a Hookean spring of Young's modulus E in series with a Newtonian dashpot of coefficient of viscosity η. The Kelvin-Voigt model consists of a Hookean spring E and a Newtonian dashpot η in parallel. Models with three elements were also developed to more realistically characterize the creeprecovery and stress relaxation behavior of a viscoelastic material. These include notably the standard linear solid model, also known as the Zener model, which consists of two Hookean springs and one Newtonian dashpot arranged in one of the two equivalent configurations. The schematics of the above rheological models are illustrated in Figures 2A-C.

Creep and Stress Relaxation
There are two one-dimensional tests that utilize the application of a step input to examine material viscoelasticity, namely, creep and stress relaxation. The creep test assesses the time-dependent change in material strain, ε(t), upon the introduction of a step stress. A stress relaxation (or relaxation) response is the timedependent change in material stress, σ(t), after a step strain is introduced. Data from the creep and relaxation tests are often acquired and analyzed in the time domain. For a linear viscoelastic material, the stress response to a step strain input of ε 0 is defined as follows: Here, F(t) is a monotonically decreasing function of time defined as the relaxation modulus. Similarly, if a step input of stress, σ 0 , is introduced, the corresponding creep response is given as follows: Here, J(t) is a monotonically increasing function of time defined as the creep compliance. The mathematical representation of creep compliance and relaxation modulus of the three classical rheological models is summarized in Table 2.

Oscillatory Loading
In an oscillatory loading test, an oscillatory excitation, instead of a step input, is used to disturb the material. The material can be loaded in a time-varying manner such that the response can be examined over a range of frequencies. The complex modulus of the material, Y p , can be determined using dynamic mechanical analysis (DMA) on the data collected with dynamic mechanical analyzers and rheometers [4,32]. The complex modulus [33] consists of a real component, the storage modulus Y ' , which is the slope of the loading curve that characterizes the elastic behavior and measures the stored energy, and an imaginary component, the loss modulus Y '' , which is the area bounded by the loading and unloading curves that characterizes the viscous behavior and measures the energy loss. The repetitive depiction of the stress-strain relationship allows for individual estimation of the frequency-dependent storage and loss moduli.
If an oscillatory strain input with an amplitude of ε 0 is applied to a linear material, following the Euler's formula, its timedependent form can be written as an exponential: The corresponding stress response will also be oscillatory with an amplitude of σ 0 , but out of phase by δ with the strain input as follows: This phase delay, δ, is referred to as the loss angle of the material and reflects to what degree the material is viscoelastic. If the stress and strain curves are completely in phase, the loss angle is at its minimum value of 0 and the material is considered purely elastic. If, however, completely out of phase, the loss angle maximizes at π/2 and the material is purely viscous. The material is considered viscoelastic when the phase delay is away from the boundary limits.
The complex modulus, Y p , is calculated from the ratio of oscillatory stress to strain, in a similar fashion as how the creep compliance and the relaxation modulus were defined. The storage and loss moduli are directly equated to the real and imaginary components of the ratio as follows: The tangent of the loss angle, tan δ, as well as the magnitude of the complex modulus, |Y|, are alternative ways to characterize the complex modulus as follows: Similarly, if an oscillatory stress is applied such that σ(t) σ 0 e iωt , the ratio of oscillatory strain to stress is defined as the complex compliance J p that is also composed of a real storage compliance (J ′ ) and an imaginary loss compliance (J ″ ) [33]. The complex compliance is used less often due to the practical difficulty to apply and control an oscillatory stress in comparison to an oscillatory strain.
As the oscillatory loading test is often conducted in the frequency domain, the expression of the complex modulus is more often a function of frequency (Y p (ω)). Once fitted to a rheological model, the coefficients of each parameter within the model can be determined at a given frequency. The expressions of

Model name
Maxwell Kelvin-Voigt Zener Frontiers in Physics | www.frontiersin.org June 2021 | Volume 9 | Article 670571 the complex moduli derived from the constitutive equations of classical rheological models are summarized in Table 2.

Limitation of Classical Rheological Models
Although classical rheological models serve as a foundational theoretical basis for the characterization of material viscoelastic behavior, they are limited by modeling accuracy. The creep and stress relaxation behaviors of the three classical rheological models are shown in Figures 3A-F. When an assumed Maxwell material is subjected to a step stress input of σ 0 , the creep response follows a linear relationship as a function of time, which fails to realistically represent the "creeping" behavior as should be observed in a viscoelastic material ( Figure 3A). The Kelvin-Voigt model, in comparison, can predict the creep response of a viscoelastic material more realistically as an assumed Kelvin-Voigt material creeps following an increasing exponential decay ( Figure 3B). However, it is limited in the ability to model the relaxation response of a viscoelastic material since an impulse at t 0 is only idealistic ( Figure 3E). As models of more elements were developed, including the three-element Zener model and the more generalized Maxwell and Kelvin models, the modeling accuracy was consequently improved. However, as the constitutive equations of multielement models became more complex, the analysis of model parameters became more computationally expensive. In fact, studies using generalized models have revealed that the viscoelastic response of several biological tissues [34][35][36][37], such as the epithelial tissue, the kidney, and the liver, follows a distinctive power law behavior such that the stress-time, strain-time, and complex modulus-frequency relationships are approximately linear on a log-log scale [31,38]. For the purpose of modeling the viscoelastic behavior of biological tissues in the common elastography frequency range of 40-1,000 Hz [31], studies have shown that the incorporation of fractional calculus into classical rheology bears modeling advantages [31,38,39].

Fractional Models
In fractional models, a fractional derivative element, the fractional springpot ( Figure 1E), is introduced. The springpot is defined by its coefficient, c α , similar to E for a Hookean spring and η for a Newtonian dashpot. The stress-strain relation as defined by a springpot [38] is given as follows: Conceptually, a springpot is a generalization of the classical viscoelastic components. At the two limiting conditions, one such that the fractional exponent α 0, the springpot reduces to a Hookean spring and its coefficient c α becomes E, whereas when α 1, the springpot reduces to a Newtonian dashpot and its coefficient becomes η. Naturally, a springpot of varying α can exhaustively represent any intermediate viscoelastic behavior bounded by the limiting conditions.
The creep and relaxation behaviors as defined in Eqs 3, 4 can now be characterized by the fractional creep compliance and the relaxation modulus, respectively [38], as follows: The creep and relaxation behaviors of the springpot have been adapted to assess the viscoelastic response of a single cell alongside its subcompartments [40][41][42]. Under an oscillatory input, the complex modulus of a springpot, following the definition from Eq. 7, can be expressed as follows: Fractional models can be constructed from the classical models by replacing the Newtonian dashpot(s) within the original models with springpot(s). Common fractional analogs of the classical models include the fractional Maxwell model, the fractional Kelvin-Voigt model, and the fractional Zener model. Schematics of the fractional analogs are shown in Figures  2D-F. Expressions [38,43] of the moduli derived based on these fractional models are summarized in Table 3.

WAVE PROPERTIES AND PARAMETERS
This section includes the relevant wave equations that are necessary in facilitating our understanding of how a sample medium can be excited by the propagating waves and how wave parameters can be derived and used to quantify viscoelasticity of the medium. In the context of wave-based elastography techniques, the assumption of one-dimensional, time-harmonic waves that propagate rightward along the +x direction in an unbounded biological medium is sufficient. We, first, restrict our discussion to pure elastic waves that propagate with no energy loss, before expanding to consider wave attenuation in a viscoelastic medium.

Wave Propagation in a Purely Elastic Medium
For a one-dimensional, time-harmonic wave that propagates rightward along the +x direction in an unbounded biological medium, several equations (Eqs 12)-b(16b) have been adapted from [31,33,44] to demonstrate the derivation of relevant wave parameters. The displacement of particles along x at time t can be defined as follows: Here, A is the wave amplitude, and k is the wave number which is given as follows: In an elastic medium with known values of Poisson's ratio (v) and the Young's modulus (E), the propagation speed of compressional wave (c c ) is as follows: Step input Creep compliance Mittag-Leffler function [173].
The propagation speed of shear wave (c s ) in an elastic medium is related to the shear modulus (G) of the medium, medium density (ρ), angular frequency (ω), phase delay (δ), and displacement between two points that the wave propagated through (Δx), which is given as follows: The shear modulus (G) can be inferred from the angular frequency (ω) and wave number (k), or known values of Poisson's ratio (v) and the Young's modulus (E), which is given as follows: Another important parameter is the propagation speed of surface acoustic waves (c sf ), which can be related to c s if evaluated in an elastic medium [45,46], and is given as follows: For nearly incompressible medium, the Poisson's ratio is approximately 0.5, allowing Eqs 14, 16b, 17 to be further simplified.

Wave Propagation in a Viscoelastic Medium
When waves propagate in a viscoelastic medium, hysteresis occurs due to its viscous nature such that the waves attenuate as they propagate, dissipating energy [31,47]. In comparison to a purely elastic medium, the stress-strain relationship of a viscoelastic medium is no longer linear. Instead of an instantaneous response to a step input, a complex response is characterized by the complex modulus Y p , as mentioned in Estimation of Material Moduli. The storage modulus is determined by the restoration of energy due to the elastic property of the medium, while the loss modulus is related to its ability to dissipate energy [33]. The fraction of stored-todissipated energy determines whether the medium behaves more like a viscoelastic solid or viscoelastic liquid. Depending on whether compressional (longitudinal) waves or shear (transverse) waves are propagated, the notation of the generalized complex modulus Y p takes on E p or G p , respectively. The equations of the propagation speed of various wave types subsequently reflect the change in medium property. If shear wave propagation is used as an example, Eq. 15b becomes as follows: The wave number also takes a complex form of k p . Eq. 16a becomes as follows: Equating the real and imaginary components of G p , expressions of the storage and loss moduli can be individually obtained [47], which is given as follows: Equation 12 can be rewritten to reflect the added influence of attenuation as follows: The real component of the complex wave number, k ′ , remains in the same form as Eq. 13. Combining with Eq. 15a, k ′ can be related with the 1D gradient of the phase delay as follows: The term e k″x describes an exponential decay along the +x direction, of which the attenuation coefficient k ″ has a leading term relating to the first power of frequency [31] as follows: For biological tissues that exhibit a power law viscoelastic behavior, the wave attenuation may also follow a power law behavior such that ω is replaced with ω α . The significance and inclusion of these equations will become clear as we move into the details of how elastography techniques are built upon various actuation methods.

VISCOELASTICITY MEASUREMENT WITH ELASTOGRAPHY
Palpation has long been a principal method to externally examine tissue stiffness. Although it is still commonly used as a preliminary assessment method to detect abnormal tissue, it is unable to provide quantifiable data. With an increasing interest in probing properties deeper within the tissue, the concept of elastography was brought to attention by Ophir et al. in 1991 [9]. Elastography techniques first require a method of actuation that introduces disturbance to the tissue, which is then assessed with an imaging tool, most commonly including ultrasound, magnetic resonance imaging (MRI), optical coherence tomography (OCT), alongside photoacoustic (PA) imaging, and Brillouin spectroscopy (BS). In combination with an actuation method, these form the measuring basis for the earliest ultrasound elastography (USE) [9], magnetic resonance elastography (MRE) [48], optical coherence elastography (OCE) [49], photoacoustic elastography (PAE) [50], and Brillouin Microscopy (BM) [51], respectively. Elastography techniques can often output quantitative values of the storage and loss moduli and/ or parameters of a fitted rheological model. In this section, we provide a classification of the types of elastography that have been used on biological tissues and how viscoelasticity can be measured through each method of actuation.
The actuation methods are conceptually illustrated in Figure 5. Mechanical, acoustic, optical, and magnetic means of actuation have been principally developed. They could be quasi-static, transient, or oscillatory. Data gathered through elastography techniques based on oscillatory input are usually analyzed in the frequency domain, and the others in the time domain. In the case of optical actuation, the disturbance is caused by photon absorption or scattering within the tissue, whereas the other methods result in mechanical deformation. A summary of the elastography techniques discussed in this section is provided in Table 4.

Mechanical Actuation
The mechanical load can be quasi-static ( Figure 5A) or dynamic ( Figure 5B). Elastography was first developed on the basis of quasistatic mechanical actuation to assess only the elasticity of a sample [9]. In the experimental setup, a compression plate was placed onto the sample surface to alter the local strain within the sample. An ultrasonic transducer was used to send echo signals into the sample. By cross-correlating the pre-and post-compression curves of the echo amplitude, the time delay between two segments of A-lines was determined at the point where the maximum correlation value occurred. From a series of time delays estimated from the crosscorrelation of multiple A-line pairs, the axial displacements, a strain profile ε(x), and a Young's modulus profile E(x), if the applied compressional stress is known, as a function of depth (x) can be estimated. In combination with ultrasonic imaging, this technique is referred to as quasi-static elasticity imaging or strain elastography. Cross-correlation became a fundamental analytical method in strain elastography and was later used to evaluate the viscoelastic behavior of a sample. Once the data gathered through crosscorrelation were fitted to a rheological model [52], the Young's modulus and coefficient of viscosity of the model elements could then be calculated using the constitutive equation of the fitted model. The ease of implementation and cost efficiency in computation have allowed quasi-static compression to be combined with other imaging modalities besides ultrasound. For instance, OCT and MRI have been used in the development of compression OCE [53,54] and compression MRE [55], respectively. Similarly, upon a quasi-static load or strain, a timedependent creep or relaxation profile can be obtained via the corresponding imaging modality and used to derive the parameter coefficients once fitted to a rheological model.
Dynamic mechanical actuation requires the placement of a mechanical vibrator onto the sample surface ( Figure 5B). The vibrator can produce transient impulses or oscillatory waves that propagate deeper into the sample along the +x direction. If the actuation is induced transiently, a single cycle of sinusoidal wave at low frequency (∼50 Hz) is typically applied to the sample surface, generating both compressional and shear waves that propagate spherically into the sample with distinguishable wave properties [56]. In particular, the propagation speed (c s ) and attenuation (k '' ) of the shear wave as functions of depth (x) Frontiers in Physics | www.frontiersin.org June 2021 | Volume 9 | Article 670571 8 can be deduced if combined with an appropriate imaging modality of a frame rate in the kHz range [57]. Once fitted to a rheological model, the coefficients of model parameters can be determined [58][59][60]. Combined with ultrasound-based imaging techniques, this technique is referred to as compression transient elastography [56]. The earlier OCE was also coupled with a transient stepwise mechanical pulse [49].
Sonoelastography is an ultrasound-based elastography technique that applies continuous mechanical vibration to produce low amplitude (less than 0.1 mm) and low frequency (less than 1 kHz) harmonic shear waves [61][62][63]. Parameters from the vibration patterns at various input frequencies, including the propagation speed of shear wave (c s ) and phase delay (δ), are analyzed from the Doppler shift and fitted into a rheological model to obtain coefficients of viscoelasticity [64]. Using MRI, one of the first viscoelasticity studies was performed by Muthupillai et al. in 1995 in conjunction with dynamic mechanical actuation [48]. This combination was termed dynamic compression MRE, where harmonic mechanical waves of frequency on the lower spectrum (less than 1 kHz) were propagated to induce shear stress. The threedimensional displacement fields, including u(x,t), and the phase delay extracted from MRI with harmonic motion-sensitizing gradient waveforms were then used to either reconstruct the viscoelastic parameters by inversion of the Helmholtz equation [7,13], or directly calculate the frequency-dependent complex shear storage G'(ω) and loss moduli G''(ω) of the sample [13,15]. In the latter case, a rheological model is needed to further determine the shear modulus and coefficient of viscosity of the model elements. Some OCE studies have also been coupled with dynamic compression, including a branch of shear wave OCE that is actuated by direct contact [65][66][67][68]. The complex wave number and shear wave speed can be extracted and used to derive the complex shear modulus. A rheological model is necessary if coefficients of viscoelasticity need to be quantified. Spectroscopic OCE (S-OCE) [69] is another representative of dynamic compression OCE. A frequency sweep (0-1,000 Hz) allows the frequency-dependent viscoelastic behavior of the sample to be examined within a range of frequencies. The raw OCT data often undergo several steps of processing to eventually arrive at relationships that can depict the complex viscoelastic response of the sample. A complex OCT signal is first obtained by sampling in the k-space and filtered to minimize phase noise, from which the phase difference, frequency-dependent complex displacement, modulus, and strain rate can be estimated [69,70]. If the complex modulus is further fitted with a rheological model, individual coefficient of the components within the model can be determined.

Acoustic Actuation
Acoustic actuation relies on focused ultrasound beams produced by an ultrasound transducer to propagate acoustic waves within the tissue sample ( Figure 5C). The displacement fields and properties of the propagating acoustic waves can be obtained through the coupled imaging tool to further characterize the viscoelastic behavior of the sample.
A technique referred to as acoustic radiation force impulse (ARFI) imaging developed by Nightingale et al. combined contactless transient acoustic actuation with ultrasound-based imaging [71,72]. ARFI relies on focused ultrasound beams to deliver a high-intensity burst that generates acoustic radiation force, deforming the sample within a specific region of interest (ROI). The spatiotemporal relaxation behavior of the sample at the focal point is then recorded by the same transducer, and the resulting time-dependent tissue displacements are mapped. By fitting the relaxation behavior to a rheological model, the viscous and elastic parameters of the model may be determined [73].
Another type of acoustic actuation used is based on the properties of acoustic shear wave. A technique that utilizes oscillatory shear waves induced by and propagating in the orthogonal direction of the acoustic actuation is shear wave elasticity imaging (SWEI) [74]. SWEI offers a quantitative assessment of local material viscoelasticity in terms of shear storage modulus (G ' ) and shear loss modulus (G '' ) derived from the propagation speed of shear wave (c s ) and phase delay (δ). The quantification of sample viscoelasticity is typically carried out by fitting to appropriate rheological models [75]; however, model-independent methods have also been developed  [12,76,77]. To allow for an extended imaging region and an increased data acquisition speed (∼5,000 fps), supersonic shear imaging (SSI) [78] was developed with the added ability to focus the impulses at multiple focal depths such that conical shear waves can be generated [14,79]. Another popular coupling modality is OCT, upon which surface acoustic wave OCE (SAW-OCE), acoustic radiation force OCE (ARF-OCE), and shear wave OCE (SW-OCE) have been developed. In SAW-OCE, the phase velocity of the propagating surface wave (c sf ) laterally across the surface of the sample and its dispersion curve can be determined from the phase delay via OCT. Values of the parameters can then be fitted into the Rayleigh wave dispersion equation [80] or a rheological model [81] to extract the viscoelastic parameters. A subset of SAW-OCE termed aircoupled OCE utilizes contactless air puffs to generate SAW, under the influence of which the sample's time-dependent deformation can be mapped with OCT [45,80,82]. The slope and area bounded by the hysteresis curve can be calculated to derive the loss in energy that corresponds to the viscous behavior and storage modulus to the elastic behavior [82]. A rheological model can be used to derive numerical values of the Young's modulus and coefficient of viscosity. In ARF-OCE, similar to ARFI, an ultrasound transducer is used to propagate pulsed ARF to remotely initiate local sample displacements [83]. The phase delay between adjacent A-lines can be obtained under OCT and used to estimate time-dependent axial displacements and strain. When model-dependent, further estimation of the parameter coefficients can be achieved. ARF-OCE can also be modelindependent where the complex shear modulus can be estimated from a direct measurement of the propagation speed of surface wave [83]. Acoustically actuated SW-OCE can sometimes be induced by ARF [84], where the propagation speed of bulk shear waves is used to derive model parameters once fitted.

Optical Actuation
Optical actuation is enabled by the optical properties of biological tissues including absorption and scattering. For example, photoacoustic imaging ( Figure 5D) was developed based on the photoacoustic effect by which a sample of biological tissue absorbs the energy from optical beams and releases it in the form of acoustic waves. A technique referred to as photoacoustic viscoelasticity (PAVE) imaging was developed by Gao et al., using intensitymodulated laser beams emitted toward the tissue sample [85]. Developed based on the photoacoustic effect as well, the tissue absorbs incoming light waves and undergoes thermal expansion. The viscoelastic nature of biological tissues introduces a phase delay (δ) in the detected photoacoustic signal, which relates to the viscosity-elasticity ratio (η/E) of the tissue and the modulation frequency (ω) [85], which is given as follows: thus providing a qualitative comparison between the viscous and elastic behavior of the tissue sample. Recent advances have demonstrated the feasibility of quantifying the Young's modulus by establishing a photoacoustic shear wave model [86]. The viscosity parameter can be subsequently determined via Eq. 24.
Using a combination of optical beams, an external vibration source, and a measurement technique based on photoacoustic computed tomography (PACT), photoacoustic elastography (PAE) has emerged in recent years, attempting to measure strain concurrently with the functional parameters of a tissue sample [50,87,88]. A strain profile can be obtained from crosscorrelating A-line pairs using photoacoustic imaging. If the stress applied by the external vibrator is known, a model can be fitted to quantify the elasticity of the sample. Although only quantitative values of the Young's modulus have been demonstrated, an extension to viscosity may be achieved if the strain profile can be obtained as a function of time, in a way that is similar to how strain elastography can be used to quantify viscoelastic parameters.
Another optically actuated elastography technique, Brillouin microscopy, uses inelastic Brillouin light scattering as a contrast mechanism to measure properties of biological tissues ( Figure 5E). The incident light interacts with the inherent acoustic phonons within the tissue sample. A frequency shift, referred to as the Brillouin shift, is then introduced to the outgoing light. This change in frequency can be directly related to the mechanical properties of the sample by calculating the complex longitudinal modulus containing a real and an imaginary part [51]. The real part yields a measurable frequency shift (v B ), which can be used to derive the longitudinal modulus of the tissue that characterizes its elastic properties [89], whereas the imaginary part is related to the spectral width (Γ B ) that can be further used to determine the longitudinal coefficient of viscosity of the material [51,90].
In combination with OCT, several optically actuated OCE techniques are worth mentioning, namely, photonic force OCE (PH-OCE) and pulse laser OCE. PH-OCE offers a contactless method that induces a harmonically modulated ultra-low radiation pressure force generated by a low-numerical aperture beam [91,92]. The radiation pressure force subsequently causes sub-nanometer oscillations of predeposited microbeads in a sample. A linear model is used to decouple the mechanical and photothermal responses of the sample to the incoming radiation pressure such that the mechanical response can be isolated. The oscillation amplitude of a microbead can then be captured under OCT and estimated, as it is a function of the input radiation pressure force and radius of the microbead. The oscillation amplitude is also directly related to the complex shear modulus of the sample from which further extraction of viscoelastic parameters can be achieved depending on the choice of a suitable rheological model. In pulse laser OCE [93], wave generation is achieved by focusing pulsed laser irradiation toward the sample. By absorption of the light and localized thermal expansion within the sample, the irradiation is converted into compressional and shear waves that further propagate within the sample. The propagation of the laser-induced elastic waves can be profiled and the local time-dependent displacements can be mapped. A derivation of pulse laser OCE uses dye-loaded perfluorocarbon nanodroplets [94]. When excited by the laser pulse, the nanodroplets undergo rapid liquid-gas phase transition, inducing elastic waves.

Magnetic Actuation
Magnetic actuation introduces a disturbance to the sample using a magnetomotive force generated by an external magnetic field. First, magnetic nanoparticles (MNPs) are predeposited into the sample that is placed within the external magnetic field ( Figure 5F). The magnitude and direction of the magnetomotive force can be controlled by adjusting gradients of the magnetic field. The viscoelastic behavior of the surrounding microenvironment of the MNPs in response to the disturbance can be inferred by analyzing the time-dependent displacement of the MNPs. One common imaging modality that magnetic actuation, whether static or dynamic, can be implemented with is OCT. Collectively, this technique is referred to as magnetomotive OCE (MM-OCE) [95][96][97]. MM-OCE allows for contactless manipulation, through which local MNP displacement in the subnanometer range can be detected. Static MM-OCE typically uses OCT to map the time-dependent phase variation in response to a step stress, from which the nanoscale timedependent local displacements, amplitude with decay, and resonant frequency can be deduced. A rheological model can be fitted to determine the elastic modulus and coefficient of viscosity of the model parameters. Swept-frequency loading techniques have also been implemented in dynamic MM-OCE studies to determine the frequency-dependent viscoelastic behavior of the tissue [97].

VISCOELASTICITY AND DEVELOPMENTAL BIOLOGY
Physical forces, whether externally applied or internally generated, drive tissue shape changes during embryogenesis [24]. For instance, anisotropic tensional stress orients ectodermal remodeling in the early mouse limb bud [98], anteroposterior (AP) tensile force mediates Drosophila germ band extension [22,99], contractile force of an actin cable initiates Drosophila dorsal closure [100], tensile convergence force drives Xenopus blastopore closure [23], and the formation of head fold in a chick embryo is driven by a mechanical force that is likely associated with neurulation [101].
Since embryonic tissues exhibit both solid-and liquid-like characteristics, they undergo viscoelastic changes under an applied force [3]. Substrates with various storage-to-loss moduli ratios can be artificially manufactured to mimic properties of the extracellular environment. It has been demonstrated that cells are sensitive to both elasticity and viscosity, resulting in morphological changes, proliferation, spreading, stiffening, softening, migration, and differentiation of cells [102][103][104]. Changes in the viscoelastic properties of extracellular matrix (ECM) have been shown to affect the functional and migratory behaviors of cells [105]. In an in vivo study of the Drosophila embryo during morphogenesis, myosin II pulses were used to assess how viscous dissipation in response to transient forces can stabilize local cell deformation [106]. Thus, to understand how physical forces regulate morphogenetic movements, characterization of the viscoelastic response within an embryo at both the cellular and tissue scales is necessary. In this section, emphasis is placed on the phase separation and transition of the intracellular membrane-less organelles, durotaxis and viscotaxis, respectively, along with the capability of existing techniques to examine viscoelastic properties at cellular and tissue scales.

Liquid-Liquid Phase Separation and Liquid-Solid Phase Transition
Intracellular mechanics have been shown to directly correlate with intracellular rheology and mechanotransduction [107]. Many intracellular compartments such as subcompartments within the nucleus, stress granules, germ cell granules, actin bodies, and other ribonucleoprotein (RNP) bodies are membrane-less and exhibit liquid-like behavior [108][109][110]. At the subcellular level, thermodynamic force drives intracellular compartmentalization toward a more energy-favorable state, sequestering molecules that assemble into phase-separated liquid droplets [26,111]. The process by which condensates form is referred to as liquid-liquid phase separation (LLPS) [112]. During the dynamic dissolution and condensation of intracellular condensates, a cell will likely exhibit any combination of elastic and viscous behaviors as its properties reflect the collective properties of membrane-less and membranebound intracellular compartments. In vivo studies have suggested that some condensates, such as the fibrillarin (FIB1) protein, exhibit viscoelastic behavior [27,108]. As many of the condensates reside within the nucleus and are constituted by RNP bodies [113], the assembly and disassembly of intracellular condensates will likely have an effect on the protein concentration [27,108,114] and subsequently the rate at which the condensates contribute to genetic activities that regulate the functional behavior of a cell [113]. For instance, it has been demonstrated in vitro that an increase in the concentration of the RNA LAF-1 within RNP bodies can decrease their viscosity [114]. However, there have been rather limited studies on the mechanical properties of intracellular condensates. Thus, a deeper look into the viscoelastic properties of intracellular and, especially, intranucleolar condensates will likely refine our understanding of the role viscoelasticity plays in nuclear and cellular activities.
During LLPS, liquid droplets can, often undesirably and irreversibly, undergo further liquid-solid phase transition (LSPT) during which compartments progress to a more viscoelastic or even solid-like state [26,27,112]. In studies of neurodegenerative diseases, LSPT has been proposed as an explanation for undesirable stiffening and inhibited molecular dynamics [26] due to protein aggregation. For example, LSPT involving a mutant FUS protein is believed to contribute to the progression of ALS [115]. Hence, a shift on the viscoelastic spectrum exhibited by a condensate that is normally fluid-like toward a more solid-like state can suggest a pathological tendency. The protein aggregation driven by LSPT is often irreversible, suggesting it is of value to investigate the exact phase boundary where liquid droplets become solid-like in vivo [116]. The long-term motivation is to detect the onset and monitor the progression of neurodegenerative diseases to assess the adequacy of medical intervention.

Durotaxis and Viscotaxis
Durotaxis was proposed in the early 2000s by Lo et al. as a mechanism by which cells migrate collectively toward greater substrate stiffness [117]. Despite the fact that biological tissues are intrinsically viscoelastic, pioneering in vitro studies in which stiffness gradients were artificially generated to mimic the ECM environment were often conducted on purely elastic substrates [118,119]. Later studies identified stiffness gradients in vivo during embryogenesis, where the stiffness is commonly quantified using elastic (Young's or shear) modulus. For example, the apparent elastic modulus of Xenopus cranial mesoderm was measured using atomic force microscopy (AFM), revealing that the collective migration of neural crest cells is triggered by mesodermal stiffening [120]. Also using AFM, it was shown that axonal growth, while not strictly durotaxis, is guided by a stiffness gradient in the brain of the Xenopus embryo [121]. Magnetic tweezers were used to uncover a mesodermal stiffness gradient along which cells migrate within the early mouse limb bud [122].
As the stiffness of ECM is sometimes attributed to the abundance of collagen [123], which can be realistically modeled as a viscoelastic material with quantifiable storage and loss moduli [124], cell-matrix interactions and resulting cell migration thus exhibit dynamic time-dependent mechanical responses [105]. Paths of migration and areas that cells migrate toward are likely of different viscosity in addition to higher elasticity, thereby impacting the timescale of the responses of migrating cells to stiffness gradients [125]. To elucidate the effect of ECM viscosity on cell movements, several in vitro studies have either isolated the effects of the viscosity of substrates or implemented tunable storage-to-loss moduli ratios of substrates. An increase in the substrate loss modulus can lead to viscous drag and cell-matrix energy dissipation that impedes the migratory speed of cells [104,105,[126][127][128][129][130]. In a separate study using viscoelastic substrates, it was proposed that within a viscous environment, greater migratory speed corresponds to greater apparent modulus [131]. Therefore, the empirically measured Young's moduli for characterizing local tissue stiffness may differ from the theoretical values in the absence of ECM viscosity considerations. A more recent study constructed a substrate with gradients of viscosity while keeping the elasticity constant, and demonstrated that human mesenchymal stem cells migrate along a viscosity gradient [132]. Viscositydependent cell migration has been given the name "viscotaxis" in some studies [132,133].
To further test the validity of these in vitro observations and to more physiologically evaluate how mechanical cues correlate with cell migration patterns, we believe there is sufficient motivation to examine the role of viscoelasticity gradients in vivo. Results obtained from studies of stiffness, viscosity, and viscoelasticity gradients should be compared to more rigorously delineate the collective migratory behavior of cells in response to mechanical cues.

Existing Techniques for Measuring Embryonic Tissue Viscoelasticity
A handful of methods have been applied to measure the viscoelastic properties of embryonic tissues, including, predominantly, AFM ( Figure 6A), magnetic tweezers using magnetic particles or ferrofluid droplets ( Figure 6B), and micropipette aspiration ( Figure 6C). AFM has been used to quantify the viscoelastic properties of embryonic tissues by indenting the tissue surface and analyzing the approach-retraction hysteresis using rheological models [134,135]. For instance, AFM has been applied to quantify the elasticity and viscosity of the mouse mandibular arch in vivo [136]. Magnetic tweezers have been used to probe tissue viscoelastic behavior in the early mouse limb bud [122], mouse mandibular arch [137], and Drosophila embryo during cellularization [138]. Other magnetically actuated systems to measure viscoelasticity include ferrofluid oil droplets deployed in the zebrafish tailbud [139] and the Drosophila embryo [140,141]. Confocal imaging has been combined with computer-aided cell tracking to calculate local tissue strain rate (a measurement of viscosity) during Drosophila germ band extension [22,99]. Micropipette aspiration, which applies a known suction force to deform a region of interest, has been applied to determine the local viscoelastic response of a tissue or a single cell [142,143]. For the purpose of measuring the viscoelastic properties of intracellular condensates, microrheology techniques have also been used to quantify the viscosity of the condensates through the Stokes-Einstein relation by assuming the condensates behave as equilibrium Newtonian liquids [112,114]. For other intracellular structures such as the nuclear actin network, microrheology has also been used to probe the viscoelastic creep response [144]. Magnetic micron-size wires in combination with rotational magnetic spectroscopy have been adapted to measure the shear viscoelasticity of cytoplasm and to quantify the shear modulus and coefficient of viscosity [145]. At the single-cell level, a rheometer has been used in vitro to estimate the timedependent power law creep function [40], relaxation behavior, and frequency-dependent complex modulus [146].
However, existing in vivo measurement techniques either require direct contact, are invasive, lack cellular spatial resolution, lack adequate throughput, or are unable to retrieve depth information. A limitation of AFM is that only twodimensional surface measurements can be reliably taken. Deeper tissue assessment by AFM either requires dissection for exposure or model-based derivation [135,136], both of which are suboptimal. Magnetic devices allow measurements within a bulk tissue sample but require the injection of magnetic particles or fluids by skilled individuals. Although the use of magnetic particles allows for the simultaneous assessment of multiple locations within a 3D sample, existing magnetic techniques lack broad spatial coverage within an embryo. To define mechanisms of cell migration in 3D, ECM properties and cellular behaviors need to be considered, and techniques that provide noninvasive, continuous, volumetric, and spatially resolved measurements without sacrificing acquisition speed would be ideal.

Adult Tissue
Transient elastography has been used in vivo for the detection and mapping of aggregated breast tumors within surrounding soft tissues [57], the assessment of stages of liver fibrosis [60], and the quantification of properties of blood clots [59]. Dynamic mechanical actuation has been adapted in compression MRE on assessing a variety of adult tissues such as liver fibrosis [15], breast lesions [7,147], gray and white matter within the brain [13], and glioblastoma [148]. In the case of liver fibrosis, for instance, the purpose is usually to diagnose and classify stages of fibrosis. Among breast lesions, it is mainly for detecting and distinguishing between malignant and benign tumors. In glioblastoma, viscoelastic measurements are for the purpose of presurgical evaluation. More recently, dynamic compression MRE has been used to assess the differences of subcortical gray matter among adults in different age-groups [149].
SWEI-based elastography has been used to probe the viscoelastic properties of a healthy liver [150], a posttransplant liver [77], liver fibrosis [75], and normal breast tissue [151]. Sonoelastography has been used to assess the shear modulus and viscosity of healthy skeletal muscle [16]. During contraction, an increase in both shear modulus and viscosity has been observed in comparison to the relaxed state [16]. In vivo applications of supersonic shear imaging (SSI) include assessment of the complex shear modulus of breast lesions [79], measuring the viscoelastic and anisotropic properties of muscle tissue [14], and liver tissue [152].
OCE (compression, ARF based, and shear wave.) is extensively used to assess the viscoelastic properties of the human cornea in vivo as the superior submicron scale resolution of OCT can capture subtle changes in wave properties [153] and tissue properties, making it suitable for the detection of early-onset disease such as keratoconus [154]. In addition, in vivo dynamic OCE has been performed on human skin to assess its mechanical properties [155]. BM has, in recent years, been used in noncontact assessment of corneal biomechanics with and without keratoconus [156]. Although the current BM application focuses only on extracting the elastic information in terms of the longitudinal modulus, it is capable of decoupling the viscous component if needed [157].

Embryonic Tissue
Elastography has emerged to measure the viscoelastic properties of embryonic tissues. Most applications have been optically actuated due to the intricate nature of light's ability to achieve subcellular spatial resolution. The attenuation of light in embryonic tissues is also less of a concern since the tissue can be considered as homogenous and isotropic, such that there is a minimal difference in the refractive indices of different parts of the embryo. Recent advances in OCE have acknowledged its potential for facilitating our understanding of the viscoelastic properties of embryonic tissues [53,158,159]. For instance, transient-compression OCE was used for the quantification of viscoelasticity of a living chicken embryo in vivo [160]. BM was used to quantify the longitudinal modulus and viscosity of spinal cord tissue [161], and ECM [162] in vivo in a zebrafish larva. For the purpose of examining the viscoelastic properties of intracellular condensates within a developing cell, a form of magnetic actuation has been used in combination with optical imaging [144].

CHALLENGE AND OUTLOOK
Elastography techniques for assessing the viscoelasticity of mature tissues have been established, and future considerations may include tissue anisotropy and heterogeneity to improve measurement accuracy. Substantial instrument development is still necessary to facilitate further application of the aforementioned elastography techniques to embryonic tissues in vivo. For elastography actuated by quasi-static compression, quantification of the elastic modulus cannot be achieved without known values of the applied stress. Moreover, coupling of the stress applied to the embryo surface requires additional attention and may vary case by case. USE-and MREbased techniques, regardless of the coupled actuation method, have spatial resolutions of 100-200 μm at best, and therefore generally cannot achieve a spatiotemporal resolution sufficient for capturing the viscoelastic behavior of embryonic tissues [163][164][165]. Acoustically actuated elastography often requires additional time gain compensation as it suffers from the tradeoff between spatial resolution, which increases with increasing frequency of the incoming wave, and attenuation of the wave intensity as the wave travels deeper into the tissue, which also increases with higher frequency. Actuating magnetically, as previously mentioned, requires the deposition of foreign particles into the tissue sample and can only provide measurements at a limited number of locations within a volume. Optical actuation methods are thermally induced. Overabsorption of light energy by the tissue may lead to phototoxicity and tissue damage. In addition, the advantageous cellular resolution offered by OCT is achieved at the sacrifice of depth penetration, although this may be insignificant as the depth of embryonic tissues may be up to several millimeters at most. A unified consideration, therefore, leaves the potential advancements of BS-and OCE-based elastography techniques in the spotlight.
However, there are some limitations. A pitfall of implementing BM on embryonic tissues is that a standardized systematic approach is not yet formulated. The biophysical basis on which the measurements and relationships were derived as well as the appropriateness of applying high frequencies (in GHz range) on biological tissues remain contentious [51,166]. Furthermore, to be able to quantify the longitudinal modulus, known numerical values of the optical properties of the tissue, including the local refractive index and tissue density, are prerequisites [51]. The complex longitudinal modulus, a parameter also used to characterize viscoelastic behavior, differs from the complex modulus (Young's or shear) used in existing elastography techniques. The lack of a universal correlation between complex moduli may require additional calibration.
A crucial issue hampering the assessment of mechanisms of morphogenesis is the lack of an appropriate approach for mapping viscoelastic properties in vivo. Despite current drawbacks, BM and OCE are among the few techniques that have the ability to generate 3D in vivo quantitative viscoelasticity data for an embryonic sample. BM can achieve submicron-scale spatial resolution [25,167], which is sufficient for visualizing subcellular structures. A recent study used BM to reveal a liquid-solid phase transition in intracellular stress granules by measuring the elastic longitudinal modulus [168]. Although no direct measurement was made on the viscoelastic properties, it demonstrated the feasibility of probing intracellular compartments with BM. OCT enables rapid 3D imaging within seconds with millimeter-scale depth penetration [169], can achieve a spatial resolution of 1-10 μm [170,171], and is able to resolve local displacements at the nanometer scale [172]. Therefore, it may be a fit for more embryonic animal models and possibly human embryos.
To advance the applicability of elastography techniques for assessing embryonic tissue properties, an ideal derivative of current methods would entail an actuation approach that is contactless, or at the minimum, easy to couple, noninvasive, and nondestructive. The technique should neither disturb the surrounding environment of the embryo nor interfere with its natural development. In addition, nontoxicity, which is a property that most optically actuated techniques lack, needs to be considered and carefully calibrated. The imaging method should provide adequate contrast, signal-to-noise ratio, data acquisition speed, and spatial resolution (submicron scale) to capture time-dependent movement of individual cells during a substantive developmental interval. During the next few years, we hope to witness the emergence of such elastography techniques and their applications to developmental biology.

AUTHOR CONTRIBUTIONS
KZ drafted the manuscript. MZ, ET, SH, and YS edited the manuscript.

FUNDING
This work was supported by a project grant from the Canadian Institutes of Health Research (#168992).