1.

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–3

Recently, high-resolution optical coherence tomography (OCT) systems have been used for elastography with both static⁴ 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–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¹⁴ and, especially, in ophthalmology.^{5,6,8,15–19} Of particular note is dynamic OCE to map corneal elasticity using noncontact mechanical loads based on air-puffs^{6,8,11,20–22} or acoustic microtapping ($\mathrm{A}\mu \mathrm{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–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 ($\mu$) using the simple relation

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.

2.

Method

In this study, we performed both numerical simulations and experimental studies of mechanical wave propagation in bulk and bounded tissue-like phantoms with mechanical characteristics mimicking soft biological tissue in the purely elastic limit (i.e., the effects of viscosity are negligible).

2.1.

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.

Fig. 1

Mechanical waves in the subsurface region of a nearly incompressible linear elastic medium (mimicking soft tissue) simulated in OnScale. Spatial and temporal parameters of the push are $d=250\mu \mathrm{m}$ and $T=100\mu \mathrm{s}$, respectively; shear wave speed is $C_s=5\mathrm{m}/\mathrm{s}$ ($\mu =25\mathrm{kPa}$). Left column: semi-infinite material with air at the upper boundary; and the right column: a 1-mm thick tissue layer bounded with a liquid from the bottom. (a), (b) Diagrams for wave excitation and propagation (taken from Ref. 5); (c), (d) wave fields (XT plots) along the medium surface displayed over an arbitrary dynamic range scaled to the signal maximum; (e), (f) temporal profiles of Rayleigh and guided waves at a 10-mm distance from the excitation point (shown with a white dashed line), respectively; and (g), (h) phase velocity dispersion plots for Rayleigh and guided waves, respectively: numerical simulation (red crosses) superimposed on the analytical solution (solid lines).⁵ Dynamic variations of the waveforms are presented in Video 1 (MP4, 498 MB [URL: https://doi.org/10.1117/1.JBO.24.7.076003.1]).

The shear wave speed was chosen in simulations to be $C_s=5\mathrm{m}/\mathrm{s}$ ($\mu =25\mathrm{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\mu \mathrm{m}$ Eq. (3) in Appendix A]; and a super Gaussian temporal profile with a characteristic time constant $T=100\mu \mathrm{s}$ [see Eq. (4) in Appendix A]. These parameters are similar to those used in $\mathrm{A}\mu \mathrm{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\mathrm{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 cross-correlation-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).⁵ 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–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.

2.2.

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 $\sim 1\mathrm{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°\mathrm{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\mathrm{m}/\mathrm{s}$.

Fig. 2

(a) Experimental setup showing a 0.5-mm-thick PVA phantom suspended on water boundary. An air-coupled $\mathrm{A}\mu \mathrm{T}$ transducer can be seen in the upper left. A red dotted line denotes an OCT B-scan line. (b) OCT structural image of thin-plate phantom. Note: uneven aspect ratio displayed to highlight phantom structure.

3.

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 $\mathrm{A}\mu \mathrm{T}$ technology can produce mechanical waves with significant components up to 8 kHz.¹⁵ Note that the high frequency cutoff in $\mathrm{A}\mu \mathrm{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.

Fig. 3

Mechanical wave propagation in a layer of nearly incompressible linear elastic medium bounded with water on the bottom [see Fig. 1(b)] simulated in OnScale. Spatial and temporal parameters of the push are $d=250\mu \mathrm{m}$ and $T=100\mu \mathrm{s}$, respectively; shear wave speed is $C_s=5\mathrm{m}/\mathrm{s}$ ($\mu =25\mathrm{kPa}$); layer thickness is $h=1\mathrm{mm}$. (a) XT plots (wave fields) calculated along the tissue surface and obtained for different signal bandwidths ($BW$) measured at $-6\mathrm{dB}$ cut-off so that ($BW*\frac{h}{C_s}=1,0.5,0.25,0.1$), respectively, for the rows from top to bottom. (b) The phase velocity dispersion computed from the wave fields of (a) red crosses, superimposed on the analytical solution⁵ (solid lines), respectively. (c) 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 $XZ$ plane with the correlation method for the same $BW*\frac{h}{C_s}$ values, respectively.

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\mathrm{mm}$, $C_s\cong 5\mathrm{m}/\mathrm{s}$) probed over a 800-Hz bandwidth (typical for air puff^20–22) would result in $f_d\cong 0.04÷0.08$; whereas a 5-kHz bandwidth for $\mathrm{A}\mu \mathrm{T}$^5,15,23 results in $f_d\cong 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\cong 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\cong 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\cong 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\cong 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. 4Fig. 5–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*\frac{h}{C_s}\cong 0.5$ and, thus, the upper row in Fig. 3 corresponding to $BW*\frac{h}{C_s}=1$ is absent from the experimental study. The shear wave speed was measured to be $C_s=3.8\mathrm{m}/\mathrm{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–6]. Note that the layer thickness measured with OCT structural imaging is $h\cong 0.5\mathrm{mm}$ (typical for human cornea), compared to 1 mm (close to that for porcine cornea) used to obtain Fig. 3.

Fig. 4

(d)–(f) Experimental study of mechanical wave propagation in a layer of PVA phantom (mimicking soft tissue) bounded with water on the bottom (see Fig. 2) is compared to numerical simulations in (a)–(c) OnScale for the same medium parameters. The dimensionless signal bandwidth $BW*\frac{h}{C_s}=0.5$ (maximum achievable in the experiment). Spatial and temporal parameters of the push are $d\cong 500\mu \mathrm{m}$ and $T=100\mu \mathrm{s}$, respectively; the elastic modulus measured for a thick phantom is $\mu =14.44\mathrm{kPa}$ (shear wave speed $C_s=3.8\mathrm{m}/\mathrm{s}$). The layer thickness was measured from the OCT structural image [Fig. 2(b)] to be $h\cong 0.5\mathrm{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⁵ (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 $XZ$ plane with the correlation method.

Fig. 5

(d)–(f) Experimental study of mechanical wave propagation in a layer of PVA phantom (mimicking soft tissue) bounded with water on the bottom (see Fig. 2) is compared to numerical simulations in (a)–(c) OnScale for the same medium parameters. The dimensionless signal bandwidth $BW*\frac{h}{C_s}=0.25$. Medium parameters and spatial and temporal parameters of the push are the same as for Fig. 4.

Fig. 6

(d)–(f) Experimental study of mechanical wave propagation in a layer of PVA phantom (mimicking soft tissue) bounded with water on the bottom (see Fig. 2) is compared to numerical simulations in (a)–(c) OnScale for the same medium parameters. The dimensionless signal bandwidth $BW*\frac{h}{C_s}=0.1$. Medium parameters and spatial and temporal parameters of the push are the same as for Fig. 4.

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 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–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*\frac{h}{C_s}\cong 0.5$) is related to high-frequency attenuation of mechanical waves during propagation [high-frequency wave field “modulation” disappears at distances larger than $\sim 5\mathrm{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*\frac{h}{C_s}\cong 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.

The 2-D group velocity images computed with experimental data [(f) in Figs. 4–6] and compared with numerically simulated [(c) in Figs. 4–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 high-frequency 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*\frac{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. The 2-D group velocity images computed for $BW*\frac{h}{C_s}=0.25$ [Fig. 5(f)] and for $BW*\frac{h}{C_s}=0.1$ [Fig. 6(f)] closely match simulated images [Figs. 5(c) and 6(c) respectively] in most details.

Figure 7 shows estimates of the group velocity averaged over the 2-D distributions (over $10\mathrm{mm}\times 1\mathrm{mm}$ area in case of a 1-mm thick medium and over $10\mathrm{mm}\times 0.45\mathrm{mm}$ for a 0.5-mm thick medium) shown in Figs. 3–6. The near field area (0.8 mm from the $\mathrm{A}\mu \mathrm{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\to 0$. We note that the high frequency asymptote for the lowest order mode is a Scholte wave having a propagation speed of $C_{\mathrm{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.

Fig. 7

Mean value of group velocity as a function of signal bandwidth. Blue (obtained for a 1-mm thick medium) and red (obtained for a 0.5-mm thick medium) circular dots correspond to averaging over numerically simulated 2-D group velocity distributions [(c) in Fig. 3 and (c) in Figs. 4–6, respectively]. Black triangles correspond to averaging over experimentally obtained distributions [(f) in Figs. 4–6] for a 0.5-mm thick PVA sample. Solid blue line presents the dispersion curve for a guided mode of the lowest order (same as in Figs. 4–6). Dashed black line corresponds to the Scholte wave velocity asymptote for this mode (derived from Ref. 5 for an air/tissue/water system). Dashed red line corresponds to the true shear wave velocity. Error bars show variance of group velocity estimates in the $ZX$ plane.

4.

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\ll 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.¹⁵ 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.⁴⁵ The authors reported not only the mean shear wave speed in the cell ($\sim 1\mathrm{m}/\mathrm{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 ($\text{thickness}\approx 10\mu \mathrm{m}$). For these experimental parameters, even an optimistic estimate of the dimensionless frequency gives $f_d\approx 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\to \infty$ can provide a careful estimate of tissue elasticity.^5,15 For more complex geometries, wave field inversion algorithms similar to those used in MRE⁴⁷ 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.

5.

Appendix A

5.1.

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, California – formerly 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²⁹).

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)].

Fig. 8

Finite element model in OnScale (formally PZFlex): (a) geometry and boundary conditions for a nearly incompressible semi-infinite medium and (b) a layer bounded on its bottom with water.

To simulate a spatiotemporally sharp push, we apply a time- and 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.

Eq. (3)

$$g(x)=e^{-4·\log 2·{\left(\frac{x}{d}\right)}^2},$$

Eq. (4)

$$s(t)=e^{-16·\log 2·{\left(\frac{t-t_0}{T}\right)}^4},$$

Eq. (5)

$$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\ll 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.⁴⁸

Fig. 9

(a) Temporal profiles of Rayleigh wave signals at different Poisson’s ratios $\nu$ (corresponding to different values of longitudinal wave speed $C_l$) and (b) convergence curve for the ${\Vert l\Vert }^2$ norm of vibration speed as a function of longitudinal wave speed $C_l$. Signal excitation parameters were $d=250\mu \mathrm{m}$ and $T=100\mu \mathrm{s}$ [see Eqs. (3)–(5)].

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⁴⁹) 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\mathrm{m}/\mathrm{s}$ and changes very little. The longitudinal wave velocity necessary for this convergence is $C_l=158.2\mathrm{m}/\mathrm{s}$. Thus, it is not necessary to set the true value of $C_l=1550\mathrm{m}/\mathrm{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 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⁵ shown in Fig. 3.

Fig. 10

(a) Temporal profiles of the Rayleigh wave signals at different grid discretization and (b) convergence curve for the Rayleigh signal vector ${\Vert l\Vert }^2$ norm. Signal excitation parameters were $d=250\mu \mathrm{m}$ and $T=100\mu \mathrm{s}$ [see Eqs. (3)–(5)].