1.

Introduction

Optical coherence elastography (OCE) uses optical coherence tomography (OCT) to detect local tissue displacement following an applied force, enabling the quantification of biomechanical properties within a region of interest. OCE can potentially identify pathologies opaque to classical OCT using local maps of tissue mechanical properties (namely Young’s modulus), which can vary greatly between healthy and disease states. Detecting spatial heterogeneities in a mechanical structure can help monitor and treat disease progression in numerous fields, including ophthalmology,^1–5 dermatology,⁶ cardiology,⁷ gastroenterology,⁸ and oncology.^9–11 While it is a general belief that dynamic OCE is capable of high-resolution (on the order of the imaging modality, OCT) mapping of tissue mechanical properties, the spatial resolution in wave-based optical elastography has not yet been characterized.

Elastography, like all medical imaging modalities, can be characterized in terms of spatial resolution (how small a feature can be seen), contrast resolution (how small a difference in the imaging parameter, shear modulus for the case of OCE, can be detected), and artifact level (how well are spatially resolved, high contrast features related to true tissue properties). In general, spatial resolution in an image is defined as the minimum distance between two objects at which they can be considered individuals.^12–16 The conventional definition of spatial resolution based on the Rayleigh criterion is often adopted in OCE and characterized in terms of the system point spread function (PSF) (i.e., the apparent size of a point source in an image). This definition was recently applied in the field of compression OCE.¹⁷ Using this approach, the propagating elastic wave used to infer modulus in dynamic OCE is typically assumed to behave as a step function at a material interface perpendicular to the propagation direction. Thus, the spatial resolution in dynamic OCE is usually defined by the lateral (perpendicular to light propagation direction) resolution of the optical measurement system. As dynamic OCE does not image isolated sources (propagating mechanical waves are tracked to determine local speed), the image depends on factors such as the size of the object being imaged and the characteristics of the mechanical load. Consequently, the image cannot be considered simply as the convolution of an input with a PSF, and the general definition of spatial resolution must be modified.

A more appropriate definition for spatial resolution in dynamic OCE is thus the minimum size of an inclusion, in a homogeneous background, for which propagation speed (and thereby, elasticity) can be differentiated from the background. Although more appropriate for OCE, this definition is still incomplete because it is independent of image contrast. That is, an inclusion may be differentiable from the background, but an elasticity image based on the local wave speed for dynamic OCE in most cases does not present the true elastic modulus contrast between inclusion and background. Indeed, in many cases where bounded geometries are imaged, even relative wave speed differences (i.e., relative contrast) are influenced not only by shear elastic modulus but also by geometric factors such as medium thickness (e.g., in the cornea).^18,19 Indeed, it is not clear whether high spatial resolution maps of the mechanical wave speed represent faithful images of the true tissue elastic modulus.

In this study, we identify and resolve propagating waves from both sides of an interface to define spatial resolution for the simplest case encountered in dynamic OCE. The resolving power of the dynamic OCE system is quantified for the case of highly accurate contrast resolution (i.e., image contrast is quantitatively close to elastic modulus contrast) with minimal artifacts (i.e., all image features are closely related to true heterogeneity in the elastic modulus) around a single interface normal to a propagating plane wave. In application to dynamic OCE, we define resolution as the smallest possible inclusion that can be accurately detected without sacrificing quantitative contrast, which considers the probing mechanism (elastic wave propagation), detection mechanism (OCT), and the mechanical model.

In OCE, the probing mechanism determines the appropriate model used to estimate mechanical properties. In dynamic OCE, the typical assumption of an infinite, homogeneous, nearly incompressible elastic material would suggest that the elastographic resolution can be described by the OCT image resolution alone, yet, spatial resolution has been reported to range from the micron scale to multiple millimeters for similar reconstruction methods based on elastic wave propagation.^20–28 Recently, it was suggested that the mechanical model of the medium is directly related to resolution in the closely related field of compression OCE,¹⁷ thus indicating that spatial resolution is not always defined by the resolution of the imaging method (OCT).

Here, we raise the main question of this paper—how do elastic wave characteristics and the detection capabilities of OCT influence the overall resolution in dynamic OCE? It has been previously demonstrated that the lateral resolution of acoustic radiation force-based ultrasound (US) elastography is defined by both excitation and detection transducers, i.e., the convolution of excitation and detection beam width profiles.²⁹ In dynamic OCE, the PSF of a typical OCT system is much narrower (usually tens of micrometers) than the PSF of the excitation source (from hundreds of micrometers to a few millimeters). Thus, the convolution between excitation and detection beam width profiles will largely resemble that of the excitation beam. If the excitation beam width truly determines spatial resolution in OCE, OCT cannot provide an advantage in resolution over conventional US elastography if we imagine that the same pushing source is used for elastic wave excitation. As we will show, spatial resolution in dynamic OCE is determined primarily by the properties of propagating mechanical waves and not by limitations of the typical OCT imaging system, indicating that high-resolution dynamic OCE can only be achieved by spatial and temporal shortening of the excitation beam width.

To explore OCE resolution under highly controlled conditions, we used both detailed simulations and measurements on a specific model system. Specifically, propagating elastic waves with various spatial pulse widths (PW) were both numerically simulated and experimentally measured. The numerical simulation was designed to mimic experimental conditions in terms of wave excitation and model geometry, highlighting wave behavior close to an interface between two regions with different shear moduli. We primarily focused on Rayleigh waves propagating in the subsurface layer along an air–medium interface to explore scattering (reflection and transmission) effects on estimated local wave speeds near a vertical boundary (i.e., boundary perpendicular to the free surface) between the two regions. Experimentally, Rayleigh wave signals were generated within a tissue-mimicking phantom using acoustic microtapping ($\mathrm{A}\mu \mathrm{T}$) with a home-built, air-coupled ultrasound source³⁰ and detected by a home-built M-B mode, phase-sensitive, spectral-domain OCT (SD-OCT) system. Group velocity maps within multipart phantoms based on simulated and measured displacement fields were reconstructed and used to quantify lateral (direction perpendicular to the vertical boundary) resolution based on wave speeds near a boundary for varying OCE operating conditions.

As will be demonstrated below, spatial resolution is fully defined by the spatial PW of the mechanical wave, the spatial sampling used to track the mechanical wave (OCT-limited), and the assumptions of the mechanical model used to infer modulus. For sampling typically used in OCT, we show that the spatial resolution in dynamic OCE is primarily determined by a characteristic wavelength of the mechanical wave pulse and further complicated by complex wave-interface interactions.

2.

Background

When a temporally short axial load is applied within a spatially infinite, isotropic, purely elastic solid, local energy is converted into a pulsed elastic wave (i.e., wave packet) propagating within the medium. The goal of all wave-based elastography systems in the elastic limit is to determine the bulk shear wave speed based on the displacement of a propagating elastic wave.

Starting with an infinite, homogeneous, nearly incompressible elastic material, consider a plane longitudinal ($l$-) or bulk shear ($s$-) wave propagating across the boundary (at plane $x=0$) between two isotropic elastic media characterized with Lamé constants ${\lambda }_{1,2}$ and ${\mu }_{1,2}$, correspondingly [Fig. 1(a)]. The propagation velocities of longitudinal $[C_l^{(1,2)}]$ and shear $[C_s^{(1,2)}]$ waves from both sides of the interface are³¹

Fig. 1

(a) Bulk media void of surface boundary conditions. Incident and reflected plane bulk waves consist of both longitudinal ($l$-) and shear ($s$-) wave packets propagating normal to the interface. (b) Boundary conditions alter propagating Rayleigh wave shape. Tractional forces on the vertical boundary scatter energy in the form of bulk and Stoneley wave components in addition to the reflected and transmitted Rayleigh wave energy.

Now consider plane waves propagating along the $x$-axis (normal to the $z$-axis), located within the bulk material. At the plane interface ($x=0$), propagating waves are both reflected and transmitted according to the reflection and transmission coefficients primarily determined by the ratio $\left(\frac{C_{l,s}^{(2)}}{C_{l,s}^{(1)}}\right)$,³² and the propagation speed is discontinuous at the boundary.

Only plane waves propagating in the $\pm x$ direction are allowed under these conditions, and no other wave modes will be generated at the interface. The amplitude of counterpropagating wave packets is simply given by the transmission and reflection coefficients. These waves will not change their temporal shape; however, superposition of the counterpropagating plane waves in the $x<0$ region close to the interface will affect the apparent shape of the incident wave packet. Preservation of temporal wave shape is exploited in traditional time-of-flight reconstruction to locally estimate propagation speed and, thereby, infer the elastic modulus.³³ In the spatial domain, the original waveform can stretch or compress as a function of local propagation speed. This scaling has been exploited to locally estimate propagation speed based on the wavenumber, thereby inferring elastic modulus.^24,34,35

To accurately reconstruct the propagating wave speed near the boundary for this simple case, the two counterpropagating modes must be differentiated to locally estimate the propagation speed of the wave packet traveling in the transmission direction. For shear-wave imaging in dynamic elastography systems, a directional filter³⁶ is usually applied to wavefield data gathered by the imaging system to extract one of the wave packets. If an ideal directional filter can be designed that does not distort the wave packet, then the spatial resolution of the elasticity map can approach the resolution of the primary imaging system (i.e., OCT-scale resolution for dynamic OCE). However, practical directional filters applied to data with finite bandwidth (BW) and finite signal-to-noise ratio (SNR) cannot preserve the spatial resolution of the primary imaging system even for this highly simplified case of plane wave propagation in an infinite medium. While this degradation is only amplified for more practical conditions encountered in dynamic OCE, it is not difficult to design a directional filter that improves reconstruction accuracy without limiting resolution to the finite BW of the filter. In other words, directional filtering will inevitably limit spatial resolution in the ideal case of wave packets traveling only in the $\pm x$ direction but does not serve as the limiting factor in typical dynamic OCE.

Typical dynamic OCE systems are confined to imaging depths of about 1 to 2 mm due to the penetration limits of OCT. Consequently, mechanical waves tracked in OCE are launched and propagate in highly bounded media, typically with a free surface (air) defining one of the boundaries. Many mode types are possible in tissue with such boundary conditions, namely, Rayleigh, Love, Stoneley, Scholte, or other types of guided waves.^5,32 Here, we consider the simplest case of boundary effects typically encountered in dynamic OCE [Fig. 1(b)], where a temporally short axial load applied to the free surface ($z=0$) between air and the propagation medium (i.e., soft tissue) is used to generate a bulk $l$- and $s$-wave that travels into the medium, as well as a Rayleigh wave that travels along the free surface.

First, consider Rayleigh waves propagating along the $x$-axis far from the excitation source where they can be easily isolated from bulk waves. For a nearly incompressible material such as soft tissue, the Rayleigh wave speed is simply related to the shear-wave speed:³²

To understand the effects of mode conversion in the neighborhood of the interface, the mechanical boundary conditions on strain ($\sigma$) must first be specified at ($x=0$, $z>0$) and ($x$, $z=0$) interfaces

The boundary conditions must be applied to solutions of the general mechanical equations of motion in the medium

For the given boundary conditions, not only are transmitted and reflected Rayleigh waves produced, but also additional modes such as a Stoneley wave propagating along the $z$-axis and bulk shear waves traveling at an angle to the $z$-axis are launched.³⁷ The specific energies converted from the pure incident Rayleigh wave to secondary modes (longitudinal, shear, Stoneley, etc.) can be difficult to quantify and will vary greatly for different geometries and stiffnesses.³⁸

Because the values of each converted mode differ depending on factors such as the ratio ${\mu }_2/{\mu }_1$, it is difficult to generalize how each mode will specifically affect the reconstruction with depth for the wide range of conditions typical of OCE. Scattered body waves emerging in a particular direction create a finite region where, even for different mixed-mode amplitudes, the signal will be distorted. This region of undefined wave behavior, i.e., a mixture of multiple propagating modes, will corrupt Rayleigh wave tracking and, hence, its speed reconstruction. We also note that the Rayleigh wave waveform itself changes when passing the interface by transferring energy to other modes. As we show below, it is hard to separate the Rayleigh wave from other modes near the interface.

Finding an analytical solution for all modes near the $x=0$ interface is not a trivial problem, especially for broadband Rayleigh-wave pulses. Thus, to explore OCE spatial resolution under these conditions, we used both detailed simulations and experimental measurements on model systems. The specific geometry is designed to highlight wave-behavior close to an interface between two lateral regions with different shear moduli. Specifically, Rayleigh waves with varying pulse width were used to determine how close single-mode surface waves could approach a vertical interface without corruption due to scattering and mode-conversion. The ability to uniquely identify Rayleigh waves from both sides of the interface enabled quantitatively accurate reconstruction of the elastic modulus.

We present a detailed analysis based on the time-of-flight method to reconstruct changes in wave velocity around a vertical boundary and demonstrate the effects of mixed-mode propagation. Here we demonstrate that in a noiseless scenario, the Rayleigh wave pulse width primarily determines the spatial resolution of elasticity maps preserving the true elastic modulus contrast between media. Based on this analysis, we quantify lateral ($x$-direction for this model) resolution as the width of the transition zone around a vertical interface ($x=0$, $z>0$) where the Rayleigh wave velocity cannot be accurately measured. This lateral resolution is related to the minimal size of an inclusion that can be detected with accurate elastic modulus contrast.

3.

Numerical Simulation

A finite element method (FEM) simulation was designed to provide an ideal, noise-free comparison between theoretical and experimental data. The elastic wave equation [Eq. (4)] was solved using OnScale (OnScale 2018, onscale.com, Redwood City, California). The domain consisted of two linearly elastic regions separated by a vertical plane (Fig. 2). Simulations were run with input bulk $s$-wave speeds ($C_s$) on either side of the vertical plane as $C_{s^{(1)}}=4.4\mathrm{m}/\mathrm{s}$ (19.36 kPa) and $C_{s^{(2)}}=2.5\mathrm{m}/\mathrm{s}$ (6.25 kPa), typical values for biological tissue. Density ($\rho =1000\mathrm{kg}/{\mathrm{m}}^3$) and $l$-wave speed ($C_l=158.2\mathrm{m}/\mathrm{s}$) were uniform throughout the domain. These properties correspond roughly to a Poisson’s ratio ($v$) of 0.4995. Numerical tests showed that above this longitudinal wave speed, solutions do not change appreciably (See Sec. 8, Appendix A).

Fig. 2

Geometry and boundary conditions used in the two-part elastic wave propagation model.

Along the top boundary of the domain ($z=0$), we applied a spatially and temporally varying pressure load to simulate the acoustic push applied to the surface by the $\mathrm{A}\mu \mathrm{T}$ transducer in our experiments. To define the spatial shape of the pressure load, we acquired a raster scan of the acoustic field [Fig. 3(a)] near the focus of the air-coupled $\mathrm{A}\mu \mathrm{T}$ transducer used for the experiments shown below with a needle hydrophone (HNC-1000, Onda, Sunnyvale, California). Because the transducer was tilted by ($45\pm 2$) deg to the surface in our experiments, we sampled the field along a 45-deg line passing through the focus of the transducer [dashed line in Fig. 3(a) and black line in Fig. 3(b)]. The acoustic intensity $I(x)$ along this line is approximately Gaussian with a full-width-at-half max (FWHM) of $\sim 600\mu \mathrm{m}$ [Fig. 3(b)]. We used this Gaussian model to define the spatial push profile:

Fig. 3

(a) Acoustic field intensity measured near the focus of the $\mathrm{A}\mu \mathrm{T}$ transducer and 45-deg tilted cut line (dashed, white). (b) Measured acoustic intensity along the 45-deg tilted cut line (black) and Gaussian approximation with $600\mu \mathrm{m}$ full-width-at-half-max (red).

The temporal push profile is given by a super-Gaussian function

For a transient excitation that strongly satisfies the stress confinement condition (i.e., infinitesimally short push), the mechanical wave BW is determined by relaxation of the stress across the acoustic excitation beam, which occurs at the speed, $C_s$:

If the push intensity is not confined during the time of stress relaxation [push duration ($T$) exceeds stress-relaxation time], BW is constrained by $T$.

The spatial pulse width (PW) of the Rayleigh wave will obey the inverse relationship in Eq. (7), expanding in space as a function of the temporal pulse duration. Limiting the BW of the pulsed wave will increase the spatial pulse width of the generated mechanical wave, regardless of what pushing method is chosen. For our experiment, the spatial profile $I(x)$ was set to have an equivalent $d=600\mu \mathrm{m}$, and the push duration ($T$) was varied (250, 300, 450, and $545\mu \mathrm{s}$) to reduce spectral content in the simulated mechanical wave.

The bottom and right boundaries were set to be absorbing layers to minimize reflections. The domain was discretized with linear FEMs on a regular rectangular grid with a minimum of 40 elements per wavelength (equivalent to a grid spacing of $2.75\mu \mathrm{m}$). Because underintegrated linear elements are prone to spurious solutions, we applied Belytchko–Bindeman strain hourglass suppression.³⁹ The equations were integrated in time using an explicit time-stepping method to generate the full vibration velocity vector (time derivative of the dynamic displacement vector) at each space and time point. Only the vertical component of the vibration velocity was analyzed for direct comparison to experimental results (Fig. 4, Video 1).

Fig. 4

Snap-shots of the simulated 2-D vibration-speed field at various time-steps for the model geometry shown in Fig. 2. (a) The incident elastic wave energy in the $-x$ region is due solely to the Rayleigh mode. Fast wave energy due to a secondary Rayleigh wave solution (seen in the $+x$ region) separates from the incident Rayleigh wave and does not affect the incident pulse shape. (b)–(f) Upon interaction with the $x=0$ boundary, mode conversion can be clearly seen in each panel. Vertical and horizontal dashed white lines denote the $x=0$ and $z=0$ boundary, respectively. The dynamic range is scaled to 30% of the maximum value and an uneven aspect ratio is used to highlight mode conversion in depth (see Video 1). Experimental (top) and simulated (bottom) wave propagation ($C_{s^{(1)}}=4.4\mathrm{m}/\mathrm{s}$ and $C_{s^{(2)}}=2.5\mathrm{m}/\mathrm{s}$). The dotted black line (bottom) highlights the spatial reconstruction region used in the results. Mode-conversion, in addition to reflection, is clearly seen. Fast wave energy due to a secondary Rayleigh wave solution (seen in the $+x$ region) separates from the incident Rayleigh wave and does not affect the incident pulse shape. Vertical dashed white line denotes the $x=0$ boundary. The dynamic range is scaled to 30% of the maximum value and an uneven aspect ratio is used to highlight mode conversion in depth. (Video 1, MP4, 5772 KB [URL: https://doi.org/10.1117/1.JBO.24.9.096006.1]).

It is worth noting that the numerical solution to Eq. (4) for the imposed boundary conditions includes wave energy propagating ahead of the Rayleigh-mode. This disturbance is evidence that the Rayleigh wave equation has additional solutions besides the primary Rayleigh root. One possible solution is a super-shear evanescent wave coupled to a plane shear wave in the medium^40,41 and has been observed experimentally in multiple studies.^42–44 For this analysis, we treat the primary Rayleigh root solution as the incident wave-mode and ignore incident wave energy leading the primary disturbance. We also note that this “fast-mode” is generated due to mode-conversion at the vertical interface and serves as a further disruptor when tracking the Rayleigh wave speed in the near-field region close to a boundary.

4.

Phantom Experiment

4.1.

Tissue Mimicking Phantom

Boundary effects were explored experimentally using an elastic tissue-mimicking material [polyvinyl alcohol (PVA)-based] with controllable mechanical properties. PVA-based phantoms were created using protocols adapted from Ref. 45. They were made by first mixing a $4:1$ ratio of dimethyl sulfoxide (DMSO, CAS: 67-68-5, EMD Millipore Corp. Billerica, Massachusetts) to water using a stir plate. The optical properties were tuned by introducing titanium dioxide microparticles. A concentrated stock solution of 0.3 wt. % titanium fine powder 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 the microparticles in solution. The concentrated microparticle solution was then added to the DMSO and water solution to achieve a concentration of 0.025 wt. %. 4 or 6 wt. % PVA (146 to 186 kDa, $>99\%$ hydrolyzed, CAS: 9002-89-5, Sigma-Aldrich Corp., St. Louis, Missouri) was then added to the solution. The solution 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, it was degassed using a vacuum chamber to remove any bubbles before casting in molds and stored at $-20°\mathrm{C}$ until the phantom hardened. The entire mold was $10\mathrm{cm}\times 10\mathrm{cm}\times 10\mathrm{cm}$.

Phantoms with dissimilar mechanical properties were prepared using this protocol but with varying concentrations of PVA (4% and 6%). The 4% PVA solution was allowed to cool down to 35°C in a mold with a central barrier prior to casting the 6% PVA solution. The microparticle concentration was slightly increased in the 4% PVA solution to visualize separate regions. The mold barrier was removed and the 6% PVA solution was cast into the pre-existing phantom, resulting in a $10\mathrm{cm}\times 20\mathrm{cm}\times 10\mathrm{cm}$ welded-surface, two-part material (Fig. 5). Once cast, the mold was stored at $-20°\mathrm{C}$ for up to 12 h to solidify. Finally, the casted, hardened phantom was placed in a water bath. During this process, the DMSO slowly diffused out of the phantom in exchange for water. The water bath was regularly changed for a minimum of 48 h to remove all DMSO. The phantoms were stored in deionized water to prevent phantom dehydration.

Fig. 5

(a) Two-part PVA phantom block with welded boundary. Red dotted line denotes OCT B-scan line. (b) OCT structural image of two-part welded phantom. Hard 6% PVA material on the left and the soft 4% PVA on the right. A nonuniform aspect ratio is used for visualization.

5.

Results

5.1.

Medium with Vertical Interface

The spatial pulse width of the Rayleigh wave in experimental and simulated data is fully defined in the hard (${\mathrm{PW}}_1$) and soft (${\mathrm{PW}}_2$) regions as the full-width-at-half-maximum of the signal envelope [Figs. 6(a) and 6(b)]. The spectral content is defined by the center frequency (normalized first moment of the spectral amplitude bound by $-20\text{-}\mathrm{dB}$),⁵⁵ $F_c$, and the maximum frequency (upper frequency at $-6\text{-}\mathrm{dB}$), $F_{\text{max}}$, of the power spectrum following an FFT of the wave disturbance far from any source or boundary [Fig. 6(c)]. As wave speed is constrained by the material properties of the medium, wavelengths are long for low-frequency excitation. Conversely, the spatial pulse width is shorter for broadband (high-frequency) excitation.

Fig. 6

(a) Normalized vibration speed of example experimental Rayleigh wave fully in the hard and soft region (${\mathrm{PW}}_{\text{mean}}=1.53\mathrm{mm}$). Full-width-at-half-maximum of the envelope (dotted red line) is used to quantify pulse width on both sides of the $x=0$ boundary, noted with a vertical dotted line. (b) Normalized vibration speed of example simulated Rayleigh wave quantified fully in the hard and soft region (${\mathrm{PW}}_{\text{mean}}=1.54\mathrm{mm}$), respectively. Both (a) and (b) are at time points $t_1=2.43\mathrm{ms}$ and $t_2=4.62\mathrm{ms}$ relative to the wave-front entering the field of view (see Video 1). (c) Spectral power of experimental and simulated Rayleigh wave just below the surface one pulse width from the $x=0$ boundary in the $-x$ direction. The center frequency ($F_c$) and upper limit ($F_{\text{max}}$) are used to define the spectral characteristics of the pulsed wave. (d) Measured mean pulse width ${[(\mathrm{PW}}_{\text{mean}}=\frac{1}{2}({\mathrm{PW}}_1+{\mathrm{PW}}_2)]$ and the corresponding spectral characteristics for all simulated and experimental data sets.

In a linear elastic material, the frequency of a harmonic wave must be conserved as it passes through the interface between two media. As a pulsed wave is simply the superposition of signals over a wide range of frequencies, the pulsed signal spectra will also be conserved, and the pulse width will scale in space as a function of the wave speed. The measured pulse width in the soft portion approximately (due to mode conversions at the interface) scales linearly with the change in wave speed (${\mathrm{PW}}_2={\mathrm{PW}}_1[\frac{C_{{\mathrm{s}}^{(2)}}}{C_{{\mathrm{s}}^{(1)}}}]\cong (\frac{2.5\mathrm{m}/\mathrm{s}}{4.4\mathrm{m}/\mathrm{s}})·{\mathrm{PW}}_1$) for all data sets.

To simplify the analysis of resolution at the interface, here we introduce a mean pulse width $({\mathrm{PW}}_{\text{mean}})$ as ${\mathrm{PW}}_{\text{mean}}=\frac{1}{2}({\mathrm{PW}}_1+{\mathrm{PW}}_2)$, defining a characteristic wavelength of the Rayleigh wave in the region close to the interface. As we will show below, this characteristic wavelength better captures the mechanical resolution of OCE than the optical resolution of the OCT imaging system.

Two-dimensional group velocity maps reconstructed from simulated and experimental data with different pulse widths are shown in Fig. 7 {${\mathrm{PW}}_{\text{mean}}=1.54\mathrm{mm}$ [Fig. 7(a)], ${\mathrm{PW}}_{\text{mean}}=4.40\mathrm{mm}$ [Fig. 7(b)] for simulation, and ${\mathrm{PW}}_{\text{mean}}=1.53\mathrm{mm}$ [Fig. 7(c)] and ${\mathrm{PW}}_{\text{mean}}=4.24\mathrm{mm}$ [Fig. 7(d)] for experiment}. As demonstrated, mode-converted waves on the order of the pulse width affect the estimated wave speed close to the boundary, resulting in a finite transition size before stabilizing to the true Rayleigh wave-speed. The mode-converted bulk shear and Stoneley waves traveling into the medium degrade reconstruction accuracy, manifesting as features in the group velocity map that are not reflective of the actual geometry (i.e., artifact). The directional filter can mitigate this effect in the negative-$x$ direction ($-x$) but is still insufficient to remove scattered wave-packets that travel at an angle normal to the $z=0$ surface, particularly in the positive $x$ direction ($+x$). In addition, the temporal wave-shape change due to boundary traction forces will limit the NCC from accurately reconstructing a propagation speed directly on either side of the boundary. The structural OCT image confirmed that the vertical boundary in the phantom was normal to the surface ($<1\mathrm{deg}$) to ensure that no unexpected wave modes were generated.

Fig. 7

Group velocity reconstruction from simulation data with (a) ${\mathrm{PW}}_{\text{mean}}=1.54\mathrm{mm}$ and (b) ${\mathrm{PW}}_{\text{mean}}=4.40\mathrm{mm}$. Experimental group velocity reconstructed with (c) ${\mathrm{PW}}_{\text{mean}}=1.53\mathrm{mm}$ and (d) ${\mathrm{PW}}_{\text{mean}}=4.24\mathrm{mm}$.

The depth-averaged group velocity measurements and the corresponding sigmoid-fit profiles are normalized and shown in Fig. 8. The transition size was defined as the full-width-at-half-maximum of the first derivative of the sigmoid fit, with respect to the lateral direction of the fitting function. Note that the small “bump” appearing on the $+x$ side of the boundary right after the transition zone [see Fig. 8(a)] is an artifact resulting from averaging the 2-D velocity reconstruction at a depth relative to the wavelength. Because the reconstructed group velocities at each depth were centered before averaging, and depth averaging was performed relative to the elastic pulse width, any artifact in the reconstruction will affect each data set in the same way. Aligning the boundary prior to depth averaging also ensures that transition zone “blurring” was not a result of depth averaging.

Fig. 8

Depth-averaged transition edge profiles with superimposed antisymmetric sigmoid fittings (solid blue line) and first-order fitting derivatives (dashed line; green for hard part, red for soft part). FEM simulation examples are shown for (a) ${\mathrm{PW}}_{\text{mean}}=1.54\mathrm{mm}$, and (b) ${\mathrm{PW}}_{\text{mean}}=4.40\mathrm{mm}$. Experimental examples for (c) ${\mathrm{PW}}_{\text{mean}}=1.53\mathrm{mm}$ and (d) ${\mathrm{PW}}_{\text{mean}}=4.24\mathrm{mm}$.

Results show that the total size of the transition zone region scaled linearly as a function of incident and transmitted spatial pulse widths (Fig. 9). The spatial width of the pulsed wave changed as a function of the local velocity, thus the quantification based on the first derivative of the asymmetric fitting function determined the transition size relative to the pulse width uniquely in the hard ($-x$) and soft ($+x$) portion of the model. When a single Rayleigh wave-mode was incident on a vertical boundary, experimental and simulated results demonstrated a linear relationship between the pulse width and the size of the blurred region, where the transition out of, or into, the other material occurs in space at approximately $\frac{1}{4}$ the size of the spatial pulse width in each material. This result indicates that to accurately differentiate between two regions with different moduli, a finite region equal to at least $\frac{1}{2}$ the mean pulse width, ${[(\mathrm{PW}}_{\text{mean}}=\frac{1}{2}({\mathrm{PW}}_1+{\mathrm{PW}}_2)]$ must be considered. Thus, the ultimate resolution in dynamic OCE may be approximated by $\frac{1}{2}{\mathrm{PW}}_{\text{mean}}$.

Fig. 9

Transition zone as a function of spatial pulse width on either side of the $x=0$ interface (labeled as $-x$ and $+x$, respectively, for the left and right sides of the $x=0$ interface). The full transition zone (labeled “$\pm x$”) relative to the mean pulse width in both media (${\mathrm{PW}}_{\text{mean}}=\frac{{\mathrm{PW}}_1+{\mathrm{PW}}_2}{2}$) is also shown for both experiment and simulation.

5.2.

More Complex Media Boundaries

To examine the relationship between resolution and contrast for more complex internal boundaries, the wave field generated by a Rayleigh wave passing through a three-part, horizontally layered material was simulated. The two outer layers were assigned shear wave speeds of $C_s=2.5\mathrm{m}/\mathrm{s}$, whereas the inner layer, representing a stiff inclusion, was assigned a shear wave speed of $C_s=5\mathrm{m}/\mathrm{s}$. The simulated push excitation was applied in the first outer layer, generating a Rayleigh wave with a spatial pulse width of 0.5 mm. This corresponds to a maximum potential spatial pulse width of 1 mm in the hard inclusion. For the wave speeds in this arrangement (2.5 and $5\mathrm{m}/\mathrm{s}$), considering that frequency is conserved, the mean pulse width is ${\mathrm{PW}}_{\text{mean}}=0.75\mathrm{mm}$. We varied the width of this middle layer (1, 0.5, 0.25, and 0.1 mm) to examine whether group velocity reconstruction can accurately resolve the inclusion.

The reconstruction presented the inclusion at nearly full contrast with minimal artifact [Fig. 10(a) and 10(b)] when the inclusion width is greater than the mean pulse width. The contrast begins to deviate proportionally when the inclusion size reduced relative to ${\mathrm{PW}}_{\text{mean}}$ {as seen for the cases of 0.5, 0.25, and 0.1 mm inclusion widths [Figs. 10(c)–10(h)]}. At an inclusion size equal to one-tenth the ${\mathrm{PW}}_{\text{mean}}$, [Figs. 10(g) and 10(h)], the reconstruction detected the vertical inclusion but with very poor contrast. We conclude that the minimum inclusion size that can be reconstructed with full contrast is approximately equal to the mean pulse width ${\mathrm{PW}}_{\text{mean}}$. Note here that the minimum inclusion size is twice the transition zone between two media because the inclusion consists of two vertical interfaces, with a transition region of $\frac{1}{2}{\mathrm{PW}}_{\text{mean}}$ at each. A corollary of the previous analysis is that full resolution at maximum contrast may be achieved for smaller inclusion sizes using shorter pulse width (i.e., higher BW) waves.

Fig. 10

Wave-speed maps in soft ($2.5\mathrm{m}/\mathrm{s}$) media with a hard ($5\mathrm{m}/\mathrm{s}$) inclusion. The incident spatial pulse width was measured as 0.5 mm, and the maximum potential pulse width was 1 mm, resulting in ${\mathrm{PW}}_{\text{mean}}=0.75\mathrm{mm}$. The inclusion width was (a), (b) 1 mm; (c), (d) 0.5 mm; (e), (f) 0.25 mm; and (g), (h) 0.1 mm. The theoretical velocity value (dotted gray-line), and the depth averaged reconstruction curve (black-line) for the corresponding group velocity reconstruction is illustrated in (b), (d), (f), and (h).

As more complex geometries are introduced, boundary conditions become even more disruptive to reconstruction accuracy. Due to the elliptical nature of particle motion in a Rayleigh wave, local stresses transfer energy along both lateral and axial directions. These multidimensional traction forces can generate complicated wave fields at all interfaces. Once a single-mode wave interacts with an inclusion boundary, wave components are guided along the boundary, mode-convert, diffract, and interfere with the original wave shape, reducing correlation and negatively affecting spatial resolution in complex boundary conditions.

To illustrate this point, a $5\text{-}\mathrm{m}/\mathrm{s}$ semicircular inclusion (6 mm wide and 3 mm deep) was simulated at the center of a $2.5\text{-}\mathrm{m}/\mathrm{s}$ media [Fig. 11(a)] and probed using a Rayleigh wave with mean pulse width, ${\mathrm{PW}}_{\mathrm{mean}}$, of 0.75 and 4.5 mm, respectively (see Video 2). The boundaries of the semicircular inclusion are better defined when reconstructed using ${\mathrm{PW}}_{\text{mean}}=0.75\mathrm{mm}$ [Fig. 11(b)] compared with ${\mathrm{PW}}_{\text{mean}}=4.5\mathrm{mm}$ [Fig. 11(c)]. Note that boundary regions are corrupted in the reconstruction even for ${\mathrm{PW}}_{\text{mean}}=0.75\mathrm{mm}$ due to mode conversions. For longer pulse widths (lower frequency) waves, the reconstructed image identified the inclusion but could not accurately quantify modulus, highlighting the importance in defining both image contrast and resolution.

Fig. 11

A hard ($5\mathrm{m}/\mathrm{s}$) semicircle (3 mm radius) inclusion imbedded in a soft (2.5 m/s) media was simulated. Prescribed model wave speeds are shown in (a). Two-dimensional group velocity reconstruction from incident wave with (b) ${\mathrm{PW}}_{\text{mean}}=0.75\mathrm{mm}$ and (c) ${\mathrm{PW}}_{\text{mean}}=4.5\mathrm{mm}$ (see Video 2 highlighting the complex wave field used to reconstruct group velocity). Simulated wave propagation in a soft ($2.5\mathrm{m}/\mathrm{s}$) media with a hard ($5\mathrm{m}/\mathrm{s}$) semicircle (6 mm diameter) imbedded inclusion. (b) 0.5-mm incident pulse width and (c) 3-mm incident pulse width. Vibration speed magnitude is scaled to 30% of the maximum amplitude to highlight mixed modes. (Video 2, MP4, 5843 KB [URL: https://doi.org/10.1117/1.JBO.24.9.096006.2]).

6.

Discussion

The typical dynamic OCE procedure is to: (1) launch a broadband mechanical wave within a sample using a temporally short, spatially sharp excitation force (akin to shear wave elasticity imaging,^56,57 acoustic radiation force imaging,^54,58 and transient MR elastography⁵⁹); (2) track the propagating elastic wave in real time (B-M mode)^30,60 or with M-B mode^56,61 OCT scanning; (3) use local displacement (vibration velocity) fields to estimate mechanical properties.^5,21,62 Under highly controlled conditions, the estimated wave velocity can define the shear modulus (stiffness) within a small region.^{5,20,56,61,63,64} Even with typical limitations imposed by a practical measurement system (i.e., measurement noise will introduce additional imperfections), we have shown using numerical simulations and experimental studies in tissue-mimicking phantoms that OCE spatial resolution is fundamentally limited by the properties of propagating mechanical waves.

As shown above, the fundamental resolution limit in dynamic OCE is determined by the pulse width of the elastic wave used to probe a material. As complicated boundaries can change wave shape and generate interfering wave-modes, a quantitative, artifact-free time-of-flight reconstruction cannot be done until the parent mode has separated from the complex displacement field. It follows that the local material properties will play a role in the achievable resolution. Thus, resolution cannot be defined simply by the detection mechanism of OCT. It is interesting that, although Qian et al.²⁹ reached similar conclusions and correctly defined resolution in conventional US elastography, they incorrectly reported that the effective lateral resolution in dynamic OCE is defined by the capabilities of the OCT imaging system. As shown in this work (as well as by Hepburn et al.¹⁷), the resolution in both wave-based and compression OCE is not defined solely by the detection mechanism of OCT.

If a Rayleigh wave encounters an inclusion with a modulus different from that of the original propagation medium, the pulse width of the incident wave, and any generated wave-modes, will scale as a function of the relative differences in $C_s$, and the resolution will be defined by the local pulse width. As we have shown, the minimum resolvable feature (e.g., any interface with elastic modulus differential) in a system will be approximately equal to ½ the mean spatial pulse width at all locations within the imaging field. The minimum inclusion size reconstructed at full contrast is twice this limit and corresponds approximately to the mean pulse width at the inclusion interface.

Note that when the inclusion modulus is unknown, the resolution in OCE can be defined only approximately because the mean pulse width depends on both the inclusion and surrounding medium moduli. Defining the resolution by the probe pulse width is not accurate and can either underestimate or overestimate it depending on inclusion properties.

Even with noise-free conditions near a boundary, mode-conversion will reduce correlation between reference and delayed wave forms and lead to errors in wave speed estimates. For an elastic plane wave traveling in a known direction within a nearly incompressible homogeneous material bound by a single surface, the maximum of the NCC between any reference and delayed signal will be very close to the theoretical limit of one at all locations, enabling direct quantification of the shear modulus.⁶⁵ When this wave approaches a heterogeneity, the discontinuity imposes boundary conditions [Eq. (3)] producing multiple mode conversions near the interface (Fig. 4, Video 1) changing the shape of the transmitted Rayleigh mode. When multiple modes contribute to the displacement field, errors due to decorrelation will shift the maximum, producing an inaccurate measurement of local time-delay.

As bulk shear and Stoneley waves launched at an interface travel in different directions from the parent Rayleigh wave packet, mode mixing will be spatially dependent, producing an artifact close to the interface in the wave-speed reconstruction (Figs. 7, 10, and 11). This effect has been recognized as an “edge-effect” in MRE,⁶⁶ and modulus quantification is often only considered accurate when viewed approximately ½ of a pulse width from any boundary or source.^67,68 Thus, image details with spatial variation less than about a half of the elastic pulse width cannot be used to directly determine local elastic modulus assuming a single-propagating wave-mode.

System noise, relative to wave amplitude (SNR), will further reduce reconstruction resolution and accuracy. In practice, phase-stability of the OCT system determines the minimum detectable displacement (vibration velocity) based on the time-delay between reference and delayed signals. For a homogeneous material, noise will shift the maximum of the cross-correlation function away from the actual peak, referred to as jitter error. The predicted jitter magnitude is the minimum error achievable by any unbiased time delay estimation algorithm. As SNR decreases, jitter will increase. It has been shown that time-delay estimation error is worse for very short window sizes and can be improved with higher SNR.^51,52

Assuming single-mode propagation, detector spacing must be wide relative to the mechanical wavelength to minimize jitter using a single kernel window size. For an inclusion, the effective detector spacing to accurately reconstruct time delays determines OCE resolution. The lateral resolution of OCT (determined by the numerical aperture of the objective lens in the imaging system, typically tens of micrometers) can provide spatial sampling much finer than typical acoustic wavelengths,⁶⁹ allowing an overdetermined set of time-delay terms to greatly improve the accuracy of local wave-speed measurements,⁷⁰ a feature used here to improve reconstruction accuracy. While current systems can detect nm-scale displacements,^71–73 both OCE resolution and reconstruction accuracy may be optimized by sub-nm-scale displacement sensitivity, allowing accurate time-delay reconstruction of low amplitude waves with very small detector spacing.

Another way to improve resolution is to utilize more complicated reconstruction methods to quantify material stiffness. To accurately recover spatial resolution to a fraction of the pulse width for an arbitrary imaging environment, full wave field reconstruction should be performed at a fine spatial scale similar to current coarse spatial scale methods in magnetic resonance elastography (MRE).^66,74–77 The general case of wave field reconstruction often assumes local homogeneity, a valid assumption when probing regions spaced by approximately ½ the wavelength of a harmonic elastic wave. Consequently, extrinsic variables can significantly influence the spatial resolution as defined in this study.

Directional filtering can reduce the influence of reflected and scattered wave-packets to improve reconstruction accuracy by attempting to isolate a single wave packet traveling in an assumed direction.³⁶ However, wave scattering can make it difficult to isolate a Rayleigh wave, particularly when there are mixed-mode projections traveling in the same direction as the parent wave. Practical directional filters applied to data with finite BW and finite SNR cannot preserve the spatial resolution of the primary imaging system. However, practical directional filters can be designed such that the BW and SNR preserve the parent wave pulse width, reducing artifact without limiting spatial resolution. As a thought experiment, if a filtering approach could be designed to isolate and remove independent wave-modes, it may be possible for OCE to approach OCT resolution. Of course, this requires that all possible modes are known.

7.

Conclusions

In conclusion, the spatial resolution in dynamic OCE when tracking mechanical wave propagation near an interface is ultimately limited by the characteristic wavelength of the mechanical wave or by the pulse width for broadband signals. The ability to accurately reconstruct an inclusion with contrast nearing that of the physical modulus breaks down when the mean pulse width is greater than the size of the physical inclusion. While the ultimate resolution in dynamic OCE struggles to match that of its parent imaging modality, OCT, an advantage that cannot be understated is that high-resolution structural information can be acquired simultaneously. Thus, OCE can provide unique information on both mechanical modulus and tissue light scattering at a scale and speed difficult to achieve with other imaging modalities. It follows that OCE and OCT can be combined with optical microangiography⁷⁸ to potentially provide quantitative image contrast based on tissue elasticity, optical properties, and microvasculature in a single acquisition that can be done within seconds. The combination of relatively coarse quantitative elastic modulus mapping with the high-resolution capabilities of OCT and OCT angiography may provide valuable clinical information that will likely continue to drive the future development of this technique.

8.

Appendix A. Effect of Poisson’s Ratio on the Finite Element Model for a Semi-Infinite Medium

Numerical simulations were performed with the FEM in OnScale (OnScale, Redwood City, California—formerly PZFlex). We performed extensive testing of the OnScale solver for the propagation of mechanical waves in nearly incompressible linear elastic media to ensure that numerical solutions for Rayleigh wave excitation and propagation were not corrupted by numerical dissipation or dispersion and were accurate for the input Poisson’s ratio ($v=0.4995$).

Soft biological tissues typically exhibit a longitudinal wave speed ($\approx 1540\mathrm{m}/\mathrm{s}$) orders of magnitude higher than the shear wave speed (1 to $10\mathrm{m}/\mathrm{s}$). This property presents a number of computational challenges including long simulation times, increased floating point error, and the formation of spurious solutions due to hourglassing³⁹ (i.e., artificial, zero-energy strains). In practice, it is not necessary to directly model the true longitudinal wave speed, as numerical solutions will tend to converge in the limit of $C_s\ll C_l$ when the longitudinal wave speed is sufficiently high. To verify this, we compared numerical solutions for different longitudinal wave speeds from 12.2 to $353.6\mathrm{m}/\mathrm{s}$ with the theoretical solution predicted by the elastodynamic Green’s function convolved with the spatial and temporal push profiles.^79,80 In all cases, the shear wave speed was set as $2.5\mathrm{m}/\mathrm{s}$. Far-field waveforms show close agreement to the theoretical solution, particularly when the longitudinal wave speed exceeds $158.2\mathrm{m}/\mathrm{s}$ (Poisson’s ratio $>0.4995$) (Fig. 12). The difference between numerical and exact solutions is quantified using the ${\ell }_2$ error. Results indicate that it is not necessary to set the true value of $C_l=1540\mathrm{m}/\mathrm{s}$ (typical for soft tissues) in the simulations. Any $C_l$ 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, spurious solutions due to hour glassing can be virtually eliminated, and floating point error in the numerical solution can be minimized.

Fig. 12

(a) Rayleigh wave profiles with various input longitudinal wave speed ($v/C_l$) at a lateral distance of 10 mm from the excitation. (b) Signal shapes converge to the analytic solution (defined by the Green’s function) for Poisson’s ratios higher than 0.4995 ($C_l=158.2\mathrm{m}/\mathrm{s}$). Signal excitation parameters for pressure load were $d=600\mu \mathrm{m}$ and $T=240\mu \mathrm{s}$. The shear wave speed was $2.5\mathrm{m}/\mathrm{s}$.

9.

Appendix B. Multiple-Kernel-Weighted-Averaging

Due to jitter error in time-of-flight reconstruction, it is practical to increase detector spacing to improve the accuracy of the measured wave speed. However, increasing the kernel size will reduce spatial resolution. To accurately reconstruct wave speed in experimental data, while minimizing spatial resolution loss, the MKWA method described in Sec. 4.3 was employed. The trade-offs between resolution and accuracy in MKWA and single kernel methods are explored below.

For direct comparison between MKWA and the single kernel method, an equivalent kernel size based on the weighting function [Eq. (10)] was defined as follows:

Reconstruction accuracy as a function of increasing weighted averages [increasing $n$ in Eqs. (9) and (10)] was compared to the single kernel time-of-flight reconstruction using the experimental data set corresponding to $P{W}_{\text{mean}}\cong 1.54\mathrm{mm}$ [Fig. 6(a), Video 1]. Assuming a quasi-semi-infinite homogeneous phantom (i.e., forward propagating Rayleigh wave only), with wave speed values generated using an unbiased estimator, variation in the local wave velocity, $C_R$, can be used to define reconstruction accuracy as

Table 1

Comparison of measured SNRCR and spatial resolution loss between single kernel and MKWA methods.

Increasing the kernel size to compute spatial gradients along the transverse dimension improves accuracy at the expense of OCE spatial resolution. The loss (%) in OCE resolution as a function of increasing kernel size for both single kernel and MKWA methods was quantified using the antisymmetric sigmoid fitting function described in Sec. 4.3 for the equivalent noise-free simulated data set. The minimum kernel size, equal to the grid spacing ($2.75\mu \mathrm{m}$), is representative of super-resolution OCT.^81,82 The minimum kernel size will not reduce resolution and was used as a baseline for determining any relative reduction in resolution with increasing $k$. Spatial resolution losses for all equivalent kernel sizes are summarized in Table 1.

The tradeoff between ${\mathrm{SNR}}_{C_R}$ gain and spatial resolution loss using MKWA versus single kernel processing is summarized in Table 1.

For a kernel size $\frac{1}{10}$th the minimum pulse width, MKWA provides an 8.1 dB increase in reconstruction accuracy and a negligible increase in the measured transition zone. Conversely, the single kernel method reduced spatial resolution at the equivalent kernel size and is more affected by noise-induced jitter error. Based on this analysis, all wave speed images constructed from experimental data used MKWA processing with an equivalent kernel size of $\frac{1}{10}$th the minimum pulse width.