Does group velocity always reflect elastic modulus in shear wave elastography?

Abstract. Dynamic elastography is an attractive method to evaluate tissue biomechanical properties. Recently, it was extended from US- and MR-based modalities to optical ones, such as optical coherence tomography for three-dimensional (3-D) imaging of propagating mechanical waves in subsurface regions of soft tissues, such as the eye. The measured group velocity is often used to convert wave speed maps into 3-D images of the elastic modulus distribution based on the assumption of bulk shear waves. However, the specific geometry of OCE measurements in bounded materials such as the cornea and skin calls into question elasticity reconstruction assuming a simple relationship between group velocity and shear modulus. We show that in layered media the bulk shear wave assumption results in highly underestimated shear modulus reconstructions and significant structural artifacts in modulus images. We urge the OCE community to be careful in using the group velocity to evaluate tissue elasticity and to focus on developing robust reconstruction methods to accurately reconstruct images of the shear elastic modulus in bounded media.


Introduction
Mapping tissue mechanical properties, or elastography, has become an important medical imaging modality. There is a large body of work using different imaging systems, such as MRI and ultrasound, to track internal displacements and strains resulting from either external or internal mechanical loads to infer tissue mechanical properties. [1][2][3] Recently, high-resolution optical coherence tomography (OCT) systems have been used for elastography with both static 4 and dynamic 5,6 loads. Although both load types can provide valuable information on tissue mechanical properties, dynamic elastography has a distinct advantage because quantitative elastic modulus maps can be obtained from local mechanical wave speed estimates for a wide range of practical operating conditions. 5,6 Due to the very high line rates of spectral-domain (tens of kHz [5][6][7][8][9] ) and swept-source OCT (a few MHz 10-13 ), dynamic elastography can operate at sub-mm spatial resolution using pulsed, temporally compact mechanical waves. In other words, dynamic optical coherence elastography (OCE) can acquire snapshots of shear wave temporal profiles in soft tissues propagating with speed many times smaller (a few m/s typically) than the speed of sound.
OCE has been used to evaluate tissue biomechanics 14 and, especially, in ophthalmology. 5,6,8,[15][16][17][18][19] Of particular note is dynamic OCE to map corneal elasticity using noncontact mechanical loads based on air-puffs 6,8,11,[20][21][22] or acoustic microtapping (AμT). 5,15,23 To date, dynamic OCE is primarily used to provide local elasticity information from estimates of the local group velocity 5,6,9,15,20,[24][25][26][27][28] of mechanical waves generated at the surface of the medium and propagating within that medium. Indeed, high quality two-dimensional (2-D) three-dimensional (3-D) maps of group velocity have been shown for the cornea using noncontact mechanical loads. 5,6,15,24 The group velocity is often interpreted as being a direct measure of bulk shear wave speed C s , which leads to an estimate of the shear elastic modulus (μ) using the simple relation E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 3 4 0 where ρ is the material density. This is often the case for bulk mechanical waves produced in magnetic resonance elastography (MRE) and in ultrasound-based shear wave elastography, where the medium viscosity is negligible but is not the case in highly heterogeneous materials or in media greatly influenced by mechanical boundary conditions (i.e., heterogeneous and bounded media).
Here, we investigate how well group velocity images represent shear modulus information for the conditions typically encountered in dynamic OCE of the cornea assuming a purely elastic system (i.e., in the absence of viscosity). As will be demonstrated below, not only do shear modulus estimates derived from a simple assumption of bulk-waves differ greatly from the actual modulus but significant spatial features in group velocity images often interpreted as shear modulus heterogeneities are actually artifacts produced by wave propagation in a bounded medium. characteristics mimicking soft biological tissue in the purely elastic limit (i.e., the effects of viscosity are negligible).

Numerical Simulation
Numerical simulations were performed with a finite element method in OnScale (OnScale, Redwood City, California-formerly PZFlex). We extensively tested the OnScale solver for the propagation of mechanical waves in nearly incompressible linear elastic media to ensure that numerical solution is not corrupted by numerical dissipation or dispersion, a common characteristic of commercially available software packages. We summarize the simulation method and the medium geometry in Appendix A and also provide Video 1 showing Rayleigh wave propagation in a bulk medium [see Fig. 1(a)] versus guided waves in a 1-mm thick layer [ Fig. 1(b)] with the same mechanical moduli but bounded with water on its bottom.
The shear wave speed was chosen in simulations to be C s ¼ 5 m∕s (μ ¼ 25 kPa), typical for many soft tissues. To excite mechanical waves, a spatially and temporally compact force was applied to the medium with a Gaussian profile in space [characteristic width of d ¼ 250 μm Eq. (3) in Appendix A]; and a super Gaussian temporal profile with a characteristic time constant T ¼ 100 μs [see Eq. (4) in Appendix A]. These parameters are similar to those used in AμT 5,15,23 experiments and define the shape and bandwidth of the Rayleigh wave. However, these simulation results are fully valid for any excitation method outside of the excitation zone.
In Fig. 1, the left column corresponds to the Rayleigh wave, and the column on the right corresponds to the guided wave. Figures 1(c) and 1(d) show XT plots of propagating Rayleigh versus guided waves along the surface of the medium. A horizontal slice of the XT plot (shown with a dashed white line) presents the temporal profile of the propagating signal for that spatial position at the medium surface (X ¼ 10 mm from the excitation point) in Figs. 1(e) and 1(f), respectively.
As seen in Fig. 1(c) and 1(e), the temporal profile of the Rayleigh wave has a simple shape that does not change with distance outside the excitation zone. Thus, the cross-correlation coefficient between two signal profiles measured at different spatial positions along the surface will equal one, independent of where measurement points are taken. This means that crosscorrelation-based calculation of the group velocity gives the same propagation speed over all X positions. Because medium viscosity is omitted in our analysis, the calculated group velocity is equivalent to the phase velocity of the Rayleigh wave and, therefore, can be directly used for shear modulus calculation with Eq. (1), noting that the Rayleigh wave speed, C R , is related to the shear wave speed C s by a constant multiplicative factor of 0.955 in nearly incompressible materials such as soft tissue. 5,25,29 Consequently, elasticity reconstruction based on the group velocity is accurate in this case.
If the medium is a layer bounded by water on its bottom [as shown in Fig. 1(b)], the signal shape is not preserved during propagation, showing dramatic perturbations of the signal temporal profile [see Fig. 1(d), 1(f), and Video 1] even if there is no viscosity. These effects result from highly dispersive guided modes [see Fig. 1(h)] that coexist and propagate along the medium surface instead of a simple nondispersive Rayleigh wave [ Fig. 1(g)]. Note that analytical solutions to wave propagation for this case were reported previously and can be found in Ref. 5.
We used a 2-D Fourier transform to process signal wave fields (XT plots) in wavenumber and frequency coordinates 30,31 and, finally, to compute phase velocity dispersion (phase velocity and frequency coordinates) from simulated data; results are superimposed on the analytical solution in Fig. 1(h). 5 In general, calculating the group velocity of multimode motion does not produce quantitative elasticity estimates. Nevertheless, this method has been used to evaluate tissue elasticity through Eq. (1) in multiple studies. 5,6,9,15,20,[24][25][26][27][28] Here, we show that this approach results in poor shear modulus maps, both in terms of quantitative accuracy and the presence of image artifacts.

Tissue Mimicking Phantoms
Thin-plate polyvinyl alcohol (PVA)-based phantoms were created using protocols adapted from Ref. 32. PVA phantoms provide nearly pure elasticity, can be tuned to closely mimic the mechanical properties of soft tissue due to their extracellular matrix, and can be easily fabricated into thin plate shapes that serve the purpose of this study. They were made by first mixing 4:1 ratio dimethyl sulfoxide (DMSO, CAS: 67-68-5, EMD Millipore Corp.) into water using a stir plate. The phantom's optical properties were tuned by introducing titanium dioxide nanoparticles. A concentrated stock solution of 0.3 wt. % titanium nanoparticles suspended in water was prepared, and a tip sonicator (Digital Sonifier 450, Branson, Danbury, Connecticut) was used for 3 min (30% duty cycle) at an amplitude of 30% to help disperse nanoparticles in solution.
The concentrated nanoparticle solution was then added to the DMSO solution to achieve a nanoparticle concentration of 0.025 wt. %. Two (2) wt. % PVA (146 to 186 kDa, >99% hydrolyzed, CAS: 9002-89-5, Sigma-Aldrich) was then added to the solution. It was covered and stirred on a hot plate maintaining a temperature of 120°C for ∼1 h to dissolve the PVA. Once fully dissolved, the solution was degassed using a vacuum chamber to remove any bubbles before casting in molds and storing at −20°C for up to 12 h to solidify.
Multiple thin-plate phantoms were cast in circular molds with a diameter of 5 cm. The thickness was controlled by adding just enough molten PVA solution to cover the bottom of the mold and then rotating the mold to evenly distribute the solution. Multiple thin slabs were created and then physically sized using OCT structural imaging. Finally, casted, hardened phantoms were placed in a water bath. During this process, DMSO slowly diffused out of solution in exchange for water. The bath was regularly changed for a minimum of 48 h to remove all DMSO. Phantoms were stored in deionized water to prevent dehydration.
During imaging, the PVA phantom was suspended on top of water to force asymmetric boundary conditions and match the medium diagram used in numerical simulations (see Fig. 2). To estimate the true shear wave speed of the phantoms, we performed group velocity estimates at the surface of a separate thick layer (15 mm thick) phantom made using the same protocol. Note that in the case of a thick layer, dispersive guided wave modes do not develop and the group velocity provides an accurate estimate of the shear wave speed, in this case C s ¼ 3.8 m∕s.

Experimental Setup
Elastic waves were tracked in the phantom using a phase-sensitive OCT (PhS-OCT) system operating in M-B mode, as detailed previously. 5,9,23,33 The A-line rate of the system was 46.5 kHz, and the optical resolution was ∼15 μm axially and ∼24 μm laterally. Each M-scan consisted of 512 A-scans in the same location repeated at 256 locations (B-scans) horizontally across the imaging plane (dx ¼ 58.6 μm), forming one complete M-B scan (1024 depth × 256 lateral locations × 512 frames) with an effective imaging range of 1.5 mm × 15 mm (axial × lateral). The total acquisition time for one M-B scan was 3.66 s.
To generate quasiplanar guided waves, AμT was used with a cylindrically focused air-coupled ultrasound (US) transducer. 5,15,23 It was fixed at an angle (20 deg to 30 deg to vertical) normal to the phantom to avoid blocking the laser light from the OCT system, providing a push beam with dimensions of d ≈ 0.5 mm lateral and L ≈ 9 mm elevational. The push was applied for T ¼ 100 μs (same as in numerical simulation) and focused to the edge of the OCT imaging range.
Note that the push width of d ≈ 0.5 mm in the experiment is close to the minimum width that can be applied with our current air-coupled transducer tilted to the surface normal by ∼30 deg. This value can be made smaller by positioning the transducer normal to the surface or by using a higher frequency ultrasound excitation. In one set of simulations, we intentionally used d ¼ 0.25 mm, i.e., twice narrower, to demonstrate the development of higher order guided modes.

Signal Processing
Similar to speckle tracking in ultrasound to detect shear waves, 34,35 PhS-OCT can be used to detect motion by analyzing the differential signal between successive A-scans at a single location following an applied load. 5,[36][37][38][39][40][41] The displacement term is defined as the change in scatterer position between scans and is more appropriately reported by the depth-resolved vertical component (along z) of vibration speed: 5,15,33 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 3 3 0 where f s is the scan rate, Δφ opt ðx; z; tÞ is the optical phase difference between successive scans, λ is the optical wavelength, and n ¼ 1. Group velocity maps in the XZ plane were obtained with the correlation method described in detail in Refs. 5, 42, and 43. Finally, the average group velocity over the 2-D image was calculated.

Results
Using the numerical model, we first explore the effects of intrinsic (shear modulus) and extrinsic (geometry and mechanical wave bandwidth) parameters on group/phase velocity measurements and their relation to elasticity maps. To simplify the analysis, we assume infinite signal-to-noise ratio (SNR) with the goal of identifying key experimental parameters that produce the best quantitation of shear modulus and minimize artifacts in elasticity maps. One of the important experimental parameters is the high frequency cutoff (or bandwidth) of propagating waves. Indeed, different excitation methods used in OCE produce different characteristic frequencies of generated mechanical waves. For example, air-puff techniques can generate broadband signals up to 400 Hz bandwidth. The newer AμT technology can produce mechanical waves with significant components up to 8 kHz. 15 Note that the high frequency cutoff in AμT 5,15,23 depends mainly on the operating frequency of the driving air-coupled US transducer. Increasing the transducer carrier frequency can produce a much higher cutoff frequency, extending even to tens of kHz.
The main question here is how shear modulus estimates using Eq. (1) differ from actual material properties when tracking guided mechanical waves of different bandwidths. Figure 3 provides an answer for a specific illustrative example.
To properly analyze guided waves, the frequency should be scaled to the medium thickness and shear wave speed, as in the central column of Fig. 3. Indeed, geometric dispersion of the guided-wave phase velocity has a universal character when plotted in dimensionless coordinates. For example, f d ¼ fh∕C s can be obtained in multiple ways with different shear wave speed and medium thickness, yet the character of wave propagation will be the same for the same dimensionless frequency. This means that our analysis can be applied to a wide range of experimental conditions encountered in OCE. Note that the large majority of prior experiments performed on the cornea, or other thin media such as the sclera and dermis, are over a range of f d smaller, or much smaller, than 1 (one). For example, a typical human cornea (thickness h ¼ 0.5 mm, C s ≅ 5 m∕s) probed over a 800-Hz bandwidth (typical for air puff [20][21][22] ) would result in f d ≅ 0.04 ÷ 0.08; whereas a 5-kHz bandwidth for AμT 5,15,23 results in f d ≅ 0.5.
To explore how the group velocity estimate depends on signal bandwidth (i.e., speed over a specific range on the dimensionless frequency scale), we applied a low-pass (LP) Gaussian filter to simulated data, artificially limiting the signal bandwidth. The filter cutoff frequency is indicated by a vertical dashed arrow in the central column of Fig. 3. For geometries without geometric dispersion (e.g., for Rayleigh wave propagation on a bulk material), the group velocity does not depend on signal bandwidth. For guided waves, it most certainly does.
The upper row of Fig. 3 presents the case where the signal spectrum is limited by f d ≅ 1. Here, a complicated wave field, shown in the left panel as the XT plot, produces (three propagating guided modes 5,29 shown by red crosses in the phase velocity dispersion panel (central panel). The elasticity map in the right panel is produced by calculating the group velocity at every position in the medium and normalizing it to the true shear wave speed. Clearly, this image is very heterogeneous, with severe deviations from the true shear wave speed. These heterogeneities have nothing to do with intrinsic material parameters (the shear elastic modulus is constant within the bounded medium) and are artifacts. Their details are dominated by extrinsic parameters such as the medium thickness and signal bandwidth.
If the signal spectrum is limited by f d ≅ 0.5 (see second row in Fig. 3), the XT plot dramatically changes and results in two guided modes and a totally different 2-D spatial distribution of estimated elasticity (again, not related to true medium elastic properties because the medium is homogeneous). This happens because limiting the signal bandwidth not only smooths the group velocity distribution but also affects the guided mode content that, in turn, creates a different group velocity structure.
Limiting the signal spectrum by f d ≅ 0.25 (see third row in Fig. 3) results in a single propagating guided mode and further changes the XT plot and the 2-D distribution of estimated elasticity. Further narrowing the signal spectrum by f d ≅ 0.1 (bottom row in Fig. 3) makes the 2-D group velocity distribution more homogeneous, but dramatically reduces the average group velocity relative to the true shear wave speed due to strong dispersion in the low f d range.
Experimental results shown in Figs. 4-6 show behavior similar to that of the simulations. To simplify the comparison, experimental results are split into three figures corresponding to a different range of dimensionless frequency f d and compared to numerical simulations in Onscale for the same medium parameters. The maximum value of the dimensionless bandwidth we could reach experimentally was BW Ã h C s ≅ 0.5 and, thus, the upper row in Fig. 3 corresponding to BW Ã h C s ¼ 1 is absent from the experimental study. The shear wave speed was measured to be C s ¼ 3.8 m∕s for the bulk 15-mm thick sample using the group velocity method. The same speed was used to compute theoretical curves for phase velocity dispersion [solid blue lines in (b) and (e) in Figs. [4][5][6]. Note that the layer thickness measured with OCT structural imaging is h ≅ 0.5 mm (typical for human cornea), compared to 1 mm (close to that for porcine cornea) used to obtain Fig. 3.
Wave fields (XT plots) of propagating mechanical waves near the sample surface from experimental data are represented in (d) in Fig. 4-6 for different signal bandwidths. To reduce the  Fig. 2) is compared to numerical simulations in (a)-(c) OnScale for the same medium parameters. The dimensionless signal bandwidth BW Ã h C s ¼ 0.5 (maximum achievable in the experiment). Spatial and temporal parameters of the push are d ≅ 500 μm and T ¼ 100 μs, respectively; the elastic modulus measured for a thick phantom is μ ¼ 14.44 kPa (shear wave speed C s ¼ 3.8 m∕s). The layer thickness was measured from the OCT structural image [ Fig. 2(b)] to be h ≅ 0.5 mm. (a) An XT plot (wave field) calculated along the phantom surface. (b) The phase velocity dispersion computed from the wave field of (a) red crosses, superimposed on the analytical solution 5 (solid lines), respectively. (c), (f) 2-D group velocity images normalized to the true shear wave speed (i.e., displayed value is 1 for a group velocity estimate equaling the true shear wave velocity) obtained for every point of the X Z plane with the correlation method.
Journal of Biomedical Optics 076003-5 July 2019 • Vol. 24 (7) signal bandwidth, LP filtering was applied to experimental data in the same way as described above. To improve SNR, measured wave fields were averaged in Z and X with a 2-D Gaussian filter having kernel size corresponding to one-tenth (1/10) of the characteristic signal wavelength C s ∕BW. They were then used to compute phase velocity dispersion via a 2-D Fourier transform [see red crosses in (e) in Fig. 4-6], which closely matches theoretical curves [see red crosses in (b) in Figs. [4][5][6].
To compare with experimental data, XT plots were also obtained from numerical simulations for the same medium parameters [(a) in Fig. 4-6]. Clearly, experimentally obtained wave fields match quite closely to numerically simulated ones. A little mismatch between (a) and (d) in Fig. 4 (BW Ã h C s ≅ 0.5) is related to high-frequency attenuation of mechanical waves during propagation [high-frequency wave field "modulation" disappears at distances larger than ∼5 mm from the source]. This behavior is clearly related with the medium viscosity not taken into account in simulations. On the other hand, for BW Ã h C s ≅ 0.1, i.e., for the narrowest bandwidth considered, the experimental wave field is affected by missing very low frequencies (below 100 Hz) due to high-pass filtration of experimental data to reduce environmental noise.   [4][5][6] show very strong variations in the XZ plane. They are mostly due to multiple dispersive modes propagating in the bounded layer. Image artifacts depend on multiple parameters, such as signal bandwidth, layer thickness, shear wave speed, and the temporal profile of the excited mechanical wave. Note that in addition to structural artifacts, there is some highfrequency noise present in group velocity images [especially in Fig. 4(f)]. Because the signal shape changes drastically from point to point, quite a small kernel size (on the order of the spatial signal wavelength, C s ∕BW) was used to calculate the group velocity. In addition, structural artifacts are smoothed at distances larger than 5 mm for BW Ã h C s ¼ 0.5 [ Fig. 4(f)] due to medium viscosity limiting the signal bandwidth. The cross-correlation coefficient between signals was always better than 0.9. Further increasing the kernel size leads to signal decorrelation and, thus, reduces the accuracy of group velocity estimates.  Figure 7 shows estimates of the group velocity averaged over the 2-D distributions (over 10 mm × 1 mm area in case of a 1-mm thick medium and over 10 mm × 0.45 mm for a 0.5-mm thick medium) shown in Figs. 3-6. The near field area (0.8 mm from the AμT source) was excluded from group velocity averaging. Clearly, the broader the bandwidth, the closer the average group velocity is to the true shear wave speed (shown by a red dashed line in Fig. 7). However, broader bandwidths also yield higher variance, as evidenced by the artifacts seen in elasticity maps shown in Figs. 3 and 4. For f d significantly lower than 1, the group velocity estimate generally follows the dispersion curve for the guided mode of the lowest order, approaching zero in the limit of f d → 0. We note that the high frequency asymptote for the lowest order mode is a Scholte wave having a propagation speed of C Sch ¼ 0.846C S . Both numerical simulations and experimental studies result in very similar average group velocity estimates (both being wrong) when the signal spectra are limited with the same f d even though the details of the measurements (e.g., medium thickness, excitation width, and shear wave speed) are different.

Discussion and Conclusions
In this paper, we have analyzed the most common situation in dynamic OCE in which a high frame rate OCT system tracks propagating mechanical waves in soft tissue. We have confirmed that for a homogeneous bulk material, and negligibly small viscosity (no frequency dispersion in the phase velocity), the group velocity of Rayleigh waves propagating along the air/tissue interface is equal to its phase velocity and can, therefore, be used in Eq. (1) (noting that C R ¼ 0.955C s 5,29 ) to directly evaluate tissue elasticity.
For a bounded material [ Figs. 1(b)], however, the situation is very different due to the tissue/liquid (or other) interfaces on wave propagation. This situation is exactly what OCE faces in evaluating corneal biomechanics and is very close to that for skin. In other words, when OCE is used to image layered tissue, multiple dispersive propagating modes should be expected.
As shown in Sec. 3, the group velocity estimate is quantitatively incorrect (due to multiple propagating guided modes), has strong artificial fluctuations in the XZ plane (see images in right columns of Figs. 3-6), and the average value is highly underestimated due to a nearly linear reduction in velocity for the lower order mode in the low frequency limit.
Artifacts in group velocity images make it difficult to probe heterogeneous layered tissues because artificial spatial fluctuations in the group velocity cannot be distinguished from real structural heterogeneities in the elastic modulus. Furthermore, both artifacts in group velocity images and inaccurate group velocity average values are hard to predict and estimate in real experiments. The problem is that the position of the measurement point on the dispersion diagram C∕C s versus fh∕C s (see central columns in Figs. 3 and 4) depends on C s itself, i.e., depends on the value to be measured! Thus, even if the medium thickness is known (can be evaluated with OCT structural image), it is still unclear how far the group velocity is from the true shear wave speed. For example, when f d ≪ 1 (which is typical for many OCE cases), the group velocity can be 10 to 100 times smaller than C s depending on which f d corresponds to the mechanical wave signal bandwidth, which is not known a priori. Note that the problems associated with using the group velocity for elasticity evaluation were briefly discussed in our previous study of ex vivo porcine corneas. 15 Replacing estimated group velocities with phase velocities was also mentioned in Ref. 44.
In conclusion, we urge the OCE community to be careful using the group velocity to evaluate tissue elasticity and be extremely careful in describing 2-D and 3-D group velocity images. In many practical cases, when the medium under study consists of one or more layers, wave propagation splits into multiple, highly dispersive, guided modes producing group velocity estimates greatly influenced by extrinsic parameters and very far from the true value of the shear wave speed. Inversion of group velocity images in such cases using Eq. (1) will produce very poor renditions of true tissue elasticity maps.
A similar analysis applies when the wavelength of propagating waves is of the same order as the object size as, for instance, when dynamic elastography is used to quantify elasticity within a single cell. For example, an impressive study was recently performed to evaluate the elasticity of a mouse oocyte with propagating mechanical waves generated at a frequency of 15 kHz. 45 The authors reported not only the mean shear wave speed in the cell (∼1 m∕s) but also images of shear modulus across the cell. Wave propagation in the oocyte is most likely severely influenced by guided wave modes in the zona pelludica (thickness ≈ 10 μm). For these experimental parameters, even an optimistic estimate of the dimensionless frequency gives f d ≈ 0.15. This means that the true shear wave speed is at least twice the measured estimate and, therefore, the shear modulus is at least four times higher than the reported value. In addition, the dimensions and boundary conditions of the cell clearly affect the shear modulus image of the cell, leading to artifacts that cannot be ignored or considered insignificant. Finally, the OCE community must develop correction methods to accurately convert the group velocity into the shear wave speed, which may work for some experimental cases, or create alternative methods to reconstruct the shear modulus (elasticity) from wave field data. For example, a model-based wave speed dispersion diagram can be used to fit experimental data and determine both viscosity and elasticity, as recently demonstrated in Ref. 46. Alternatively, high frequency asymptotes of guided waves approaching the Rayleigh and Scholte wave speeds when f d → ∞ can provide a careful estimate of tissue elasticity. 5,15 For more complex geometries, wave field inversion algorithms similar to those used in MRE 47 must be considered. In any event, there is a great need to develop robust reconstruction methods for dynamic OCE to convert wave field data into quantitative maps of tissue shear modulus.

Finite Element Model for a Semi-Infinite Medium
To model the propagation of elastic waves in soft tissue, we constructed a 2-D finite element model using OnScale (Onscale, Redwood City, Californiaformerly PZFlex). This approximation matches well to our experimental situation of a line excitation source, 15,23 avoids diffraction effects during wave propagation, and maximizes the cross-correlation coefficient between signal profiles at different spatial locations (under this approximation the Rayleigh wave signal shape does not change during propagation 29 ).
First, we wanted to ensure that the simulation provides accurate results for nearly incompressible media, demonstrating no significant artifacts (numerical dispersion or attenuation) during wave propagation. To do this, Rayleigh wave excitation and propagation in a bulk linear-elastic material were investigated [see Fig. 8(a)].
To simulate a spatiotemporally sharp push, we apply a timeand spatially dependent pressure load to the top (z ¼ 0) boundary of the domain. The spatial push profile is taken to be a Gaussian function [Eq. (3)] with full-width-at-half-max (FWHM) given by the push width d. The temporal push profile is given by a super-Gaussian function [Eq. (4)] with FWHM given by the push duration T. A time delay t 0 is introduced to avoid impulsive loading at time t ¼ 0. The full pressure load Pðx; tÞ is given by Eq. (5). In addition, we assume that the top boundary is free of shear stress.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 5 3 9 Pðx; tÞ ¼ P 0 · gðxÞ · sðtÞ: Boundary conditions are assigned as follows. We assume the solutions are symmetric about the x ¼ 0 boundary. The lower and right boundaries are set to be absorbing layers that minimize the reflection of incident waves back into the domain. Because absorbing boundaries do not always perfectly prevent reflection of incident shear waves, the right boundary was set sufficiently far away from the region of interest so any partially reflected waves did not affect the results.
The domain was discretized with linear finite elements on a regular rectangular grid with a minimum of 40 elements per wavelength. The equations were integrated in time using an explicit time-stepping method and the vertical component of the vibration velocity field, equivalent to data recorded in OCE experiments, was extracted for analysis.
Because the Rayleigh wave speed does not depend on the longitudinal wave speed C l in a nearly incompressible medium, 5,29 the numerical solution should converge in the limit C s ≪ C l . We checked this hypothesis by varying the Poisson ratio (see Fig. 9). For nearly incompressible media such as this, underintegrated linear finite elements are prone to spurious solutions. To prevent this, Belytchko-Bindeman strain hourglass suppression was applied to all models. 48 Figure 9(a) shows how Rayleigh wave profiles change with increasing Poisson's ratio. Clearly, signal shapes converge to the analytic solution (defined by the Green's function 49 ) for Poisson's ratios higher than 0.4995. This is also seen in the quantitative error estimate in Fig. 9(b), which shows that for sufficiently high longitudinal wave speed, the numerical solution approaches the analytic solution at C l ¼ 1550 m∕s and changes very little. The longitudinal wave velocity necessary for this convergence is C l ¼ 158.2 m∕s. Thus, it is not necessary to set the true value of C l ¼ 1550 m∕s (typical for soft tissues) in simulations; any C l values corresponding to Poisson's ratio larger than 0.4995 will produce the same wave propagation over the time scales used in these simulations. Using this fact, simulation times can be dramatically reduced (from 8 h to 50 min per full simulation on a standard, 128 GB RAM workstation), spurious solutions due to hourglassing can be virtually eliminated, and floating point error in the numerical solution can be minimized.
To ensure that the numerical solution provides the correct outcome, we varied the grid size (number of points per wavelength). The results are shown in Fig. 10. As seen, profiles Journal of Biomedical Optics 076003-8 July 2019 • Vol. 24 (7) converge at more than 40 points per wavelength. In addition, we see that the solution does not change with propagation distance, which verifies the absence of significant numerical attenuation and dispersion within the simulation volume for the parameters used. The validity of results is further supported by a nearly perfect match of simulated wave speed dispersion curves to the analytical solution for guided waves 5 shown in Fig. 3.

Finite Element Model for a Bounded Medium
The bounded tissue model consisted of two layered regions, both of which are assumed to be linear elastic materials under plane strain [ Fig. 8(b)]. The upper layer models a medium under study with a thickness h, shear wave speed C s , density ρ ¼ 1000 kg∕m 3 , and longitudinal wave speed C l ¼ 158.2 m∕s. These properties correspond to a Poisson's ratio of about 0.4995, depending on the shear wave speed C s . We set the longitudinal wave speed in the lower medium to that of water (closely approximating the anterior chamber of the eye). This region can support longitudinal waves but not shear waves and generates appropriate reflections back into the upper tissue layer. To achieve this, we match the density and longitudinal wave speeds to the upper region (ρ ¼ 1000 kg∕m 3 , C l ¼ 158.2 m∕s) but set the shear wave speed to C s ¼ 0 m∕s. The thickness of the lower layer is set to four times the estimated wavelength in the top layer to minimize interactions between the "leaky" longitudinal wave and the lower boundary.

Disclosures
The authors declare they have no competing financial interests.  Journal of Biomedical Optics 076003-9 July 2019 • Vol. 24 (7) imaging of the human posterior tibial tendon. As a senior fellow in Department of Bioengineering, University of Washington, his research focuses on medical ultrasound imaging, shearwave elastography, and optical coherence elastography.
John Pitre received his BSE in biomedical engineering from Tulane University and MS and PhD degrees in biomedical engineering from the University of Michigan. He is currently a postdoctoral researcher at the University of Washington in the Department of Bioengineering. His research interests include ultrasound and optical coherence elastography, tissue biomechanics, and elastic wave propagation.
Mitchell A. Kirby received his BSc degree in biomedical engineering from Michigan Technological University. He is currently a doctoral student in bioengineering at the University of Washington. His research interests include biophotonics, imaging, and the combination of acoustics and optics in medicine. Journal of Biomedical Optics 076003-11 July 2019 • Vol. 24 (7)