Shear wave pulse compression for dynamic elastography using phase-sensitive optical coherence tomography

Abstract. Assessing the biomechanical properties of soft tissue provides clinically valuable information to supplement conventional structural imaging. In the previous studies, we introduced a dynamic elastography technique based on phase-sensitive optical coherence tomography (PhS-OCT) to characterize submillimetric structures such as skin layers or ocular tissues. Here, we propose to implement a pulse compression technique for shear wave elastography. We performed shear wave pulse compression in tissue-mimicking phantoms. Using a mechanical actuator to generate broadband frequency-modulated vibrations (1 to 5 kHz), induced displacements were detected at an equivalent frame rate of 47 kHz using a PhS-OCT. The recorded signal was digitally compressed to a broadband pulse. Stiffness maps were then reconstructed from spatially localized estimates of the local shear wave speed. We demonstrate that a simple pulse compression scheme can increase shear wave detection signal-to-noise ratio (>12  dB gain) and reduce artifacts in reconstructing stiffness maps of heterogeneous media.

Shear wave pulse compression for dynamic elastography using phase-sensitive optical coherence tomography Thu-Mai Nguyen, a, * , † Shaozhen Song, a,b, † Bastien Arnal, a Emily Y. Wong, a Zhihong Huang, b Ruikang K. Wang, a,c and Matthew O'Donnell a 1 Introduction Elastography is a recently developed field to evaluate the biomechanical properties of soft tissue for medical diagnosis. Stiffness variations have indeed proven to be good indicators of pathological state, providing valuable additional information to structural images in many applications: breast tumor detection, 1 liver fibrosis diagnosis, 2,3 and vascular pathologies characterization. 4,5 Several techniques have been developed in the past 20 years, most of them involving magnetic resonance imaging (MRI) 6 or ultrasound 7,8 modalities. More recently, optical coherence elastography (OCE) has emerged 9 based on the optical coherence tomography (OCT) imaging. OCT provides high spatial resolution (micron scale) and is, therefore, an ideal technology to examine small organs such as skin layers, ocular tissues (cornea, intraocular lens), and peripheral vessels. The earliest OCE techniques used imaging of tissue displacements following a static deformation. 10,11 Static elastography, however, suffers from the following limitations: operator dependence (variability of the applied deformation) and nonquantitative estimates of the stiffness due to high sensitivity to boundary conditions (uneven stress distribution).
Phase-sensitive OCT (PhS-OCT) [12][13][14][15] is one of the latest improvements enabling displacement detection with much higher sensitivity (nanometer scale) than the previous optical intensity-based techniques. Dynamic elastography (i.e., based on dynamic excitations [16][17][18] and later on transient elastography methods [19][20][21] ) have been implemented using PhS-OCT. Transient elastography has been introduced to address the abovementioned limitations of static elastography. In the transient regime, propagating shear waves are induced in the sample and their propagation is tracked with PhS-OCT. The stiffness of the sample is then characterized by retrieving the shear modulus from the shear wave propagation speed. PhS-OCT has also been combined with multifrequency excitations to investigate the frequency dependence of the mechanical response of viscoelastic media. 22 In this work, we propose to combine dynamic elastography with a digital pulse compression method to improve shear modulus estimation. Pulse compression was first introduced for radar detection in the early 1960s 23,24 and later extended to ultrasound applications. 25 The main idea is to spread the instantaneous peak energy level over a longer period while simultaneously improving the axial resolution and increasing the signal-to-noise ratio (SNR) through pulse compression. Here, we propose to implement pulse compression with shear waves to: (1) improve the SNR of detected displacements induced in the sample and (2) compress spatially and temporally the shear wave for better reconstruction of shear modulus maps. For that purpose, we propose to generate a broadband frequency-modulated (chirp) mechanical excitation and apply digital pulse compression to the resulting displacement field detected with PhS-OCT. We evaluate the performance of this pulse compression scheme by comparing the results to those from monochromatic excitation as employed in our previous studies. 26

Materials and Methods
We performed shear wave elastography experiments on agar phantoms. Our setup consists of a PhS-OCT imaging system combined with a mechanical shear source.

Sample Preparation
Agar phantoms were made to mimic the elastic properties of biological soft tissues. Two types of tissue-mimicking phantoms were constructed using the same methods detailed in the previous studies: 26 • A homogeneous phantom was made from a 0.5%-agar (w/v) solution. This simple medium was used to assess the performance of our pulse compression algorithm.
• A heterogeneous phantom was made from a 0.5%-agar background and a 1%-agar, 1-mm diameter cylindrical inclusion. This phantom constitutes a crude model of a diseased tissue (e.g., tumor) that we used to demonstrate the benefits of pulse compression for elastic modulus reconstruction in complex media.
A few drops of 10% latex microspheres suspension (0.3-μm diameter, Duke Scientific Corporation, California) were added to both phantoms as optical scatterers. Such scatterers were used in the previous studies 27 and found to provide adequate backscattering to reach imaging ranges up to 2-mm deep.

Shear Wave Generation
Shear waves were induced in the sample using a piezoelectric actuator (AE0505, Thorlabs, New Jersey) driven by a function generator (Tektronix, Oregon) amplified by a power amplifier with controllable gain (AE Techron, Indiana). Typical voltages of 80 V were used to induce ∼8 μm amplitude displacements of the actuator. For the purpose of this study, we used a linear frequency modulation to generate an 8-ms long chirp with a frequency content ranging from 1 to 5 kHz as the driving signal for the actuator. In addition, to compare with the previous techniques, we also repeated all the measurements using an 8-ms toneburst at a single frequency of 3 kHz. This excitation method is similar to that used in our previous work 26 and serves as a reference to evaluate the performances of the pulse compression method.
The coupling between the actuator and the sample was done through a polished stainless steel tip attached to the actuator. The tip was wedge-like in shape, with a minimal width of 1.7 mm. The waves generated by this tip can therefore be considered approximately planar. The tip was positioned in such a manner that it touches the sample side without deforming it (Fig. 1). Accurate positioning was performed using a precision translation stage and monitored using real-time OCT B-mode images. The vertical vibration of the tip generates displacements that are mainly polarized axially (i.e., along the OCT probing beam axis) and propagate laterally across the sample, creating a propagating shear wave. The tip was placed on the side of the sample so that significant displacements were induced across the whole sample depth, generating bulk shear waves.
Such a shear wave has a propagation speed that depends on the sample elastic moduli. For a locally homogeneous, incompressible, and isotropic medium that does not exhibit strong shear viscosity, the shear wave speed c S is 28 where μ is the shear modulus of the sample and ρ is its density. . Light in the reference arm is delivered to a stationary mirror while that in the sample arm is focused into a tissuemimicking phantom. Backscattered light from both arms is recombined through an optical coupler. The resulting spectral interferogram is detected by a line-scan camera (14 bit, 1024 pixel, InGaAs detector, Sensors Ltd, New Jersey) capable of 47,000 A-scans per second. The system has a measured dynamic range of ∼105 dB at 0.5-mm depth and a measured phase noise of 3.0 mrad. Shear waves were detected by operating the PhS-OCT system in M-B mode, as detailed in our previous publications. 26 In short, each M-scan is a set of 512 successive acquisitions performed at one given lateral location at a 47 kHz repetition rate, yielding a ∼10 ms recording time. The incident light beam is then swept across the sample to perform a B-scan while the mechanical excitation is repeated. Each B-scan consists of 128 A-lines with a 23.4 μm spacing, covering a lateral imaging range of 2.5 mm. The axial pixel size is ∼5 μm. A full M-B scan therefore consists of 512 × 128 A-lines and is acquired in ∼1.4 s. The actuator excitation was triggered at the beginning of each M-scan using custom LabVIEW code.

PhS-OCT Imaging System
Each M-scan was used to retrieve the displacements induced by shear wave propagation. Phase differences Δφ between two successive A-lines yield the local axial displacement u z ðx; z; tÞ according to u z ðx; z; tÞ ¼ Δφðx; z; tÞ 4πn λ; (2) where ðx; zÞ are the spatial coordinates of a given pixel, t is a given time, λ ¼ 1310 nm is the central wavelength of the light source, and n≈1.33 is the refractive index of the sample. Thus, the axial displacement field was computed as a function of time for each pixel of the region of interest, providing a "movie" of shear wave propagation (512 frames at 47 kHz).

Pulse compression
Digital pulse compression is theoretically achieved by performing the autocorrelation of the displacement field. In practice, we used an inverse filter method to optimize the autocorrelation, mainly by reducing side lobes of the autocorrelation function. This approach closely parallels that discussed detailed in Ref. 25 for the conventional ultrasound imaging, where the filter was designed using the following steps: tÞ is computed by averaging along the depth of the axial displacements at x ¼ 0 (nearest point to the shear source).
• An ideal output pulse pulseðτÞ is derived from the autocorrelation of the reference signal where wðτÞ is a weighting function designed to suppress side lobes of the autocorrelation function, ⊗ and * denote respectively the convolution and the multiplication.
• An inverse filter fðτÞ is then designed by inverting where U ref z ðτ; τ 0 Þ is the convolution matrix of the reference signal. The inverse filter is then applied (i.e., convolved) to the autocorrelation of the axial displacement field at each location in the region of interest.

Shear modulus reconstruction
The local shear wave speed was determined by applying a custom time-of-flight algorithm to the axial displacements after pulse compression. The details on the principle of time-of-flight estimation for shear wave elastography can be found in Refs. 1 and 29. Briefly, correlation along the time of two adjacent Alines of the axial displacement field provides the travel time of the shear wave between these two adjacent locations, yielding the local shear wave speed (group velocity) at each location within the region of interest. Maps of the local shear modulus are then retrieved using Eq. (1).
A "time-gating" of the displacement field can be applied prior to the time-of-flight algorithm to reduce the effects of reflected shear waves on the reconstruction in heterogeneous media. Time-gating selects only the earliest incident shear wavefront that is less likely to be affected by reflections (on elastic heterogeneities) occurring at later times (i.e., "near-ballistic propagation"). For that purpose, we developed an algorithm to detect the earliest arriving shear wave (at each location of the imaging plane) and applied a Gaussian temporal window to eliminate later wavefronts. Figure 2 shows the broadband displacement field detected in a 0.5%-agar phantom resulting from an 8-ms long chirp excitation. As illustrated in Fig. 2(a), the shear wave propagates laterally (from left to right) with lower frequencies (i.e., larger wavelengths) at early times and higher frequencies (i.e., smaller wavelengths) at later times. The SNR was computed from the spectrum of the displacements at one location [ Fig. 2(b)] as the ratio between the energy contained in the [1 to 5] kHz range and the energy out of that range. For this example, the SNR was estimated to be 10.2 dB.

Homogeneous Phantom
The distortion of the displacement field at the phantom surface is probably due to interferences between the surface and bulk shear waves, which propagate with slightly different speeds. Thus, the surface region was not considered for the pulse compression.   Fig. 2. The result of that inverse filter is shown in Fig. 4. Figure 4 shows the result of pulse compression obtained by applying the filter shown in Fig. 3(b) to the data shown in Fig. 2. The initial displacement field has been spatially and temporally compressed to a short broadband pulse of 1-ms duration, using the procedure described in Sec. 2.4.1. Figure 4(b) is an example of the axial displacements at one given location. The SNR was estimated to be 23.0 dB from the spectrum shown in Fig. 4(b), which corresponds to a 12.8 dB increase compared to the initial chirp.

Figures 3(a) and 3(b) show, respectively, the ideal output pulse [as defined in Eq. (3)] and the inverse filter [as defined in Eq. (4)] derived from the displacement field shown in
The compressed displacement field at a given depth can also be represented in the ðt; xÞ domain as displayed in Fig. 5(a). In the ðt; xÞ domain, the local slope of the signal represents the local shear wave propagation speed. The local shear wave speed values are converted to the local shear modulus values, leading to the shear modulus map shown in Fig. 5(b). For this phantom, the median shear modulus was estimated to be 4.14 AE 0.63 kPa.

Heterogeneous Phantom
We performed experiments on a heterogeneous phantom to investigate the advantages of pulse compression in a more complex medium. Figure 6 shows a B-mode image of the phantom, made from a 0.5%-agar background and a 1%-agar inclusion. This structural image exhibits very low contrast between the inclusion and the background.
We first performed measurements using a 3-kHz single-frequency excitation for the comparison with our previous work. 26 Figure 7(a) shows the displacement field detected at a given depth in the ðt; xÞ domain. The distortions of the shear wavefronts can be observed on most cycles. They are due to interference between the incident shear wave and reflections of the shear wave at the boundaries of the inclusion. Multiwave superposition due to reflections can significantly affect the reconstruction, as shown in Fig. 7(c): the corresponding shear modulus map exhibits strong artifacts in the inclusion (unexpected soft region inside the inclusion and soft stripe near the right edge) as well as downstream of the inclusion (soft stripes on the right side of the region of interest). The median shear modulus of the expected inclusion region and the background estimated from this map are 6.86 AE 4.26 kPa and 4.86 AE 1.53 kPa, respectively.
The effects of the reflected shear waves can be reduced using a time-gating method, as described in Sec. 2.4.2. The resulting time-gated displacement field is shown in Fig. 7(b). The corresponding shear modulus map is shown in Fig. 7(d). Only minor artifacts have been suppressed and the primary ones remain. This suggests that the first wavefront has too large a period (0.33 ms) to avoid interfering with the reflected waves. The median shear modulus of the expected inclusion region and    the background estimated from this map are 7.20 AE 5.11 kPa and 4.51 AE 1.30 kPa, respectively. Figure 8 summarizes the results obtained on the heterogeneous phantom using the pulse compression method. Figure 8(a) shows the displacements detected at a given depth resulting from the chirp excitation. Note that only the first 6 ms of the 8-ms long signal is displayed here. Similar to the single-frequency excitation, interferences between the incident waves and the reflected waves can be observed on the later cycles of the chirp. Figure 8(b) shows the displacement field after pulse compression. The initial chirp has been compressed to a main broadband pulse of about 0.21-ms width. Again, side lobes and distortions appear on the later cycles as a consequence of shear wave reflections at the boundaries of the inclusion. However, time-gating can be applied to isolate the earliest wavefront [ Fig. 8(c)]. The corresponding shear modulus map is displayed in Fig. 8(d) with the exact same color scale as that of Figs. 7(c) and 7(d). Artifacts within the inclusion have been eliminated, allowing clear identification of the inclusion as a stiffer region than the background with high elastic contrast. The median shear modulus of the inclusion and the background are 14.81 AE 5.12 kPa and 5.00 AE 0.79 kPa, respectively.

Discussion and Conclusion
We have implemented digital pulse compression with an OCTbased shear wave elastography technique. We have performed the experiments on tissue-mimicking phantoms. A broadband ([1 to 5] kHz) frequency-modulated mechanical excitation was applied to phantoms through a piezoelectric actuator. The subsequent displacement field was detected and recorded using a PhS-OCT system. The inverse filtering was performed on the recorded displacements to spatio-temporally compress the initial propagating chirp into a broadband propagating pulse. The pulse propagation was then used to either reconstruct maps of the local shear modulus or compute dispersion curves (see below).
We have demonstrated that pulse compression significantly improves the SNR of detected displacements: a 12.8 dB increase has been obtained in a homogeneous agar phantom. PhS-OCT already has a great sensitivity (nanometer scale) to detect axial displacements compared to other imaging modalities such as ultrasound or MRI. Increasing the SNR will still be valuable for in vivo applications and particularly for ocular tissues. For instance, the cornea is both a stiff and fragile tissue that cannot undergo large displacements.
We have also demonstrated that pulse compression is an efficient method to simplify the reconstruction of shear modulus maps in heterogeneous media. The incident shear wave can indeed be reflected by the boundaries of elastic heterogeneities. For the case of nonlocalized (spatial and temporal) excitations, reflected waves have a high probability of interfering with the incident wave, leading to artifacts in straightforward time-offlight reconstruction. In contrast, pulse compression enables time-gating on a localized (spatial and temporal) signal to isolate the incident waves from reflected waves. For instance, for the heterogeneous phantom presented in this study, the expected shear modulus values are, respectively, about 4 and 15 kPa for the 0.5%-agar background and the 1%-agar inclusion according to the values reported in the literature. 30 The shear modulus of the inclusion estimated from the nonlocalized excitation ( Fig. 7) is strongly underestimated (7.20 kPa) and has a high spatial relative variance (71% of the shear modulus), whereas that estimated from the pulse compression approach (Fig. 8) is much more accurate (14.81 kPa) and has a lower spatial relative variance (34% of the shear modulus). The shear modulus contrast between the inclusion and the background is much higher using the pulse compression method (50%) than using single-frequency excitation (23%).
Other refinements in postprocessing can be added to timegating to improve the reconstruction. The shear modulus map shown in Fig. 8 (reconstructed after pulse compression and time-gating) still exhibits heterogeneities downstream of the inclusion. These apparent heterogeneities are likely caused by diffraction of the shear wave by the inclusion. In the previous work on elastographic reconstruction for shear waves, directional filters have been applied to displacement fields to reduce the effects of multidirectional propagation due to diffraction. 31,32 Using a two-dimensional Fourier representation of displacement fields, plane waves propagating in different directions occupy different regions of Fourier space. Consequently, directional filters can be applied to eliminate propagation directions that are highly angled compared to the incident one. Even without these more sophisticated reconstruction approaches, we have proven the impact of pulse compression. Optimizing reconstructions with more sophisticated algorithms will be considered in future work.
An additional feature of pulse compression is the possibility to perform shear wave spectroscopy with single-shot acquisition. Frequency-dependent information such as the phase velocity (i.e., propagation speed of each frequency contained in the wave) and the attenuation can indeed be extracted from a broadband shear wave. The details can be found in Refs. 33 and 34. Although the samples used in this study did not exhibit any significant dispersion over the frequency range of the applied chirp, this feature will be useful for investigating biological tissue. For instance, viscosity can be extracted from the dispersion curves. Besides, thin-layered tissues (e.g., the cornea, vessels walls, skin layers) constrain shear wave propagation to guided modes that induce strong dispersive effects. 5,34 Accessing the dispersion will thus be highly valuable, especially for OCT-based elastography, since OCT tends to examine small tissues due to its high spatial resolution. Finally, we used a setup combining a mechanical actuator and PhS-OCT similar to our previous studies. 26 The feasibility of performing elastography measurements in vivo with such a setup has been demonstrated in our previous studies on human skin. 20 Nevertheless, pulse compression can be extended to a wide range of dynamic elastography methods using different shear sources and different detection modalities.