Spectroscopic optical coherence elastography

Steven G. Adie; Xing Liang; Brendan F. Kennedy; Renu John; David D. Sampson; Stephen A. Boppart

doi:10.1364/OE.18.025519

1. Introduction

Optical coherence tomography (OCT) is a high-resolution medical and biological imaging modality analogous to ultrasound imaging except reflections of low-coherence light are detected rather than those of sound [1]. Optical coherence elastography (OCE) is an extension of OCT that can image sample mechanical properties with resolution and penetration depth governed by its parent modality [2]. Early studies utilized speckle tracking methods to resolve micro-scale displacements of tissue resulting from the application of quasi-static compression [2–7]. A drawback of speckle tracking is that the amplitude and shape of OCT speckles varies with sample deformation, making it difficult to perform micro-scale cross-correlation of speckles. Doppler OCE, introduced to overcome the decorrelation problem of speckle tracking, is based on detecting phase shifts between A-scans to compute maps of velocity, displacement, and strain rate under quasi-static mechanical compression [7,8]. The method has also been used to track phase changes between B-scans during slow compression [9].

Since its inception, developments in OCE have migrated from the use of static compression to dynamic mechanical excitation, in parallel with developments in ultrasound and MRI elastography [10–12]. As with OCE, these methods provide spatial resolutions governed by their respective parent modalities. In contrast to quasi-static OCE techniques [6], dynamic OCE methods use mechanical waves in the audio frequency range for excitation [13–19]. These studies presented OCE images at a fixed frequency, or at most two separate mechanical vibration frequencies.

The development of dynamic excitation in ultrasound has naturally led to mechanical spectroscopy of tissue [12,20,21]. There have been similar developments using low-coherence spectral-domain optical detection to measure the natural vibration frequencies of a sample, based on step internal excitation with magnetically-actuated nanoprobes [22]. A recent extension of this work demonstrated measurement of the frequency-dependent complex mechanical response of liver tissue undergoing fixation [23]. In addition to the elastic response, these methods provide access to the viscous damping properties of a sample that describe the energy dissipation of the applied mechanical excitation. They target a need to obtain better correlations between tissue viscoelastic properties and human disease. These studies utilized internal excitation using magnetically-actuated nanoprobes, and were based on M-mode OCT, i.e. did not present B-mode images of the mechanical response.

In this paper, we present a method to image the frequency-dependent complex mechanical response of a sample undergoing externally-applied audio frequency excitation. We demonstrate, for the first time, mechanical spectroscopy based on B-mode OCT imaging. Sinusoidal mechanical excitation was performed in the audio-frequency range using a piezoelectric (PZT) actuator, and a spectral-domain (SD)-OCT system was used to track the resulting phase shifts (between A-scans) in B-mode images using a method previously described [8]. Displacement maps were computed from the phase shifts, which were used to calculate the local displacement amplitude and (mechanical) phase response. In Section 2, we describe the theory of mechanical resonance, and the B-mode OCT signal during mechanical excitation. Section 3 presents details of the experimental setup and data processing steps to compute the complex mechanical response. In Section 4, we characterize the frequency response of the mechanical hardware and present mechanical spectroscopic measurements of silicone phantoms, a rat mammary tumor margin, and rat muscle tissue. In Section 5, we discuss the limitations of the technique and outline further work, and conclude in Section 6.

2. Theory

2.1 Mechanical resonance

In this section we review the basis of resonant behavior in a viscoelastic system with a single degree of freedom, consisting of a mass, m, coupled to a spring with spring constant k, and a viscous damper with viscous damping constant c. This review is based on derivations previously reported [24]. Under harmonic excitation at angular frequency Ω and force amplitude $F_{0}$ , the equation of motion can be written as

m Δ z^{″} (t) + c Δ z^{'} (t) + k Δ z = F_{0} e^{i Ω t},

where the time-dependent displacement

Δ z (t)

of the mass is constrained to the z-axis, and

Δ z^{'} (t)

and

Δ z^{″} (t)

are the velocity and acceleration, respectively. Assuming a linear system, a particular solution to Eq. (1) is expected to be harmonic and the steady-state solution can be written as

Δ z (t) = \frac{F_{0}}{\sqrt{{(k - m Ω^{2})}^{2} + {(c Ω)}^{2}}} e^{i (Ω t - Φ)},

where the phase angle (mechanical phase) of the response is

Φ = \tan^{- 1} (\frac{c Ω}{k - m Ω^{2}}) .

The undamped natural (or resonant) frequency of the system is given by $Ω_{n} = \sqrt{k / m}$ . This indicates that a higher spring constant is thus associated with a higher resonance frequency. Three regimes of the mechanical response can be identified. At resonance, $Ω / Ω_{n} = 1$ , and the phase angle $Φ = π / 2$ . Below resonance, $Ω / Ω_{n} < 1$ , the phase angle $Φ \to 0$ with decreasing Ω, which occurs more rapidly for smaller damping coefficients. Above resonance, $Ω / Ω_{n} > 1$ , the phase angle $Φ \to π$ with increasing Ω, which occurs more rapidly for smaller damping coefficients. The behaviors in these three regimes result in a resonant peak with a broader linewidth for higher values of the viscous damping coefficient.

2.2 OCT signal with applied vibration

During B-mode OCT imaging the time-dependent phase of the applied vibration can be expressed in the lateral spatial coordinate of the image, x, using $x = v_{x} t$ , where $v_{x}$ is the transverse scan velocity. Under sinusoidal mechanical excitation at angular frequency Ω and for amplitude less than the source coherence length, the signal at depth z and lateral scan position xcan be written as

\tilde{S} (x, z, Ω) = {\tilde{S}}_{0} (x, z) e^{i ϕ (x, z, Ω)},

where

{\tilde{S}}_{0} (x, z)

is the complex OCT signal in the absence of applied vibration. The optical phase,

ϕ (x, z, Ω)

, is related to the local displacement

d (x, z, Ω)

according to

ϕ (x, z, Ω) = (4 π n / \bar{λ}) d (x, z, Ω)

[8], where

\bar{λ}

is the mean optical wavelength, and n is the background refractive index of the sample. We express

d (x, z, Ω)

in terms of the PZT excitation and mechanical phase response

Φ (x, z, Ω)

, and write the optical phase

ϕ (x, z, Ω)

as

ϕ (x, z, Ω) = (4 π n / \bar{λ}) Re {| \tilde{h} (x, z, Ω) | e^{- i Φ (x, z, Ω)} Δ z_{Drive} e^{i (Ω x / v_{x} + ϕ_{0})}} .

Here we have expressed the optical phase in terms of the local complex mechanical response of the sample

\tilde{h} (x, z, Ω) = | \tilde{h} (x, z, Ω) | e^{- i Φ (x, z, Ω)}

, and the complex analytic form of the PZT displacement,

Δ z_{Drive} e^{i (Ω x / v_{x} + ϕ_{0})}

. The amplitude

| \tilde{h} (x, z, Ω) | = d_{Amplitude} (x, z, Ω) / Δ z_{Drive}

represents the ratio of the local displacement amplitude to that of the applied excitation, and the factor

e^{- i Φ (x, z, Ω)}

encodes the mechanical phase from Eqs. (2) and 3. The constant

ϕ_{0}

is the initial phase of the PZT drive function at

x = 0

.

3. Methods

3.1 Experimental setup

The OCE system, as shown in Fig. 1 , was based on a spectral-domain OCT system which used a Nd:YVO₄-pumped titanium:sapphire laser as a broadband source, providing a center wavelength of 800 nm and a bandwidth of 100 nm. The full-width half-maximum axial and transverse resolutions of the OCT system were approximately 3 µm and 13 µm, respectively. The average power incident on the samples was 10 mW.

Fig. 1 Schematic diagram of the OCE system.

Download Full Size | PDF

The sample arm of the OCE system employed a PZT stack (AE0505D16F, Thorlabs, Newton, NJ) to sinusoidally excite the sample (distal to the incident optical beam) in the axial direction. The PZT was driven with a maximum displacement of 4.5 µm at vibration frequencies within the range DC to 1 kHz. The sample was bounded from below by a coverslip, with approximate thickness of 125 μm, that was fixed with epoxy to the PZT rod (~5 mm × 5 mm), and from above by a round wedge prism (PS810-B, Thorlabs, Newton, NJ) with a 2° angle fixed to form a semi-rigid upper boundary to the sample. The glass wedge was mounted on a translation stage, allowing movement against its spring. During alignment, the wedge was lowered to make contact with the sample and then translated a further 10 μm to apply a constant preload. The axial depth scans in the OCE images (depth × lateral pixel dimensions of 1024 × 4000), were detected using a CCD line-scan camera at acquisition line rates of 1 kHz, 5 kHz and 10 kHz. The camera acquisition was synchronized with a transverse scanning galvanometer and the PZT excitation signal.

3.2 Data processing

The data processing steps, summarized in Fig. 2 , extract the complex local displacement $\tilde{d} (x, z, Ω) = \tilde{h} (x, z, Ω) Δ z_{Drive} e^{i (Ω x / v_{x} + ϕ_{0})}$ from the optical phase term given in Eq. (5). First, the complex OCT image $\tilde{S} (x, z, Ω)$ is obtained by Fourier transforming the spectral-domain data after background subtraction, dispersion correction and resampling to k-space. The lateral phase derivative is computed from Eq. (4), and encoded in complex form as

\tilde{S} (x + d x, z) {\tilde{S}}^{*} (x, z) = e^{i [ϕ (x + d x, z, Ω) - ϕ (x, z, Ω)]},

where it is assumed that the separation between A-scans,

d x

, is small relative to the lateral resolution, to ensure that the OCT signal in the absence of vibration

{\tilde{S}}_{0} (x + d x, z) \approx {\tilde{S}}_{0} (x, z)

. (Note that the spatially-dependent amplitude of the product in Eq. (6) has been omitted). For data with small phase shifts between adjacent A-scans, the phase difference between every m ^th line can be computed and divided by m, provided that their separation is significantly less than the lateral resolution. Values of

m > 1

has the benefit that in low SNR regions, the random phase noise contributions are reduced by a factor of m.

Fig. 2 Flowchart of the processing steps used to compute the spatially-localized complex mechanical response of the sample.

Download Full Size | PDF

Signal conditioning of complex interferometric data before evaluating the phase has shown to be an effective method to filter phase noise [25–28]. The measurement of interferometric phase in electronic speckle pattern interferometry is notorious for a high level of noise due to the statistical nature of speckle [25]. Additionally, it is known that the error in estimating interferometric phase is proportional to $1 / \sqrt{S N R}$ [29], leading to significant underestimation of the velocity in Doppler OCT near the limits of the velocity range [28]. We have found that median and low-pass filtering (along the lateral dimension) of the real and imaginary parts of Eq. (6) provides a robust method to minimize phase noise before taking the angle to obtain the real-valued phase difference map

Δ ϕ (x, z, Ω) = ϕ (x + d x, z, Ω) - ϕ (x, z, Ω) .

This signal conditioning is particularly useful when the SNR is low, or the sample velocities are near the limits of the velocity range ( $| Δ ϕ | \approx π^{-}$ ) when the addition of phase noise results in phase wrapping. Phase wrapping manifests as a multiplicative factor (the factor is approximately 1 for $| Δ ϕ | \approx π^{-}$ , and approximately −1 for $| Δ ϕ | \approx π^{+}$ ). This multiplication, or amplitude modulation of the real-valued phase difference map, modifies the effective frequency content of the vibration signal. Noting the convolution theorem, phase noise can thus shift the frequency content of signal (corresponding to the vibration frequency) outside the passband of the subsequently applied bandpass filter. If the magnitude of phase differences due to the underlying physical motion is large (i.e. $| Δ ϕ | > π$ ), then these signal conditioning steps cannot minimize the effects of wrapping, and phase unwrapping techniques should be applied.

Residual phase noise is minimized by bandpass filtering $Δ ϕ (x, z, Ω)$ along the lateral dimension, and the displacement map is obtained by integration and multiplication by $\bar{λ} / 4 π n$ [8]:

d (x, z, Ω) = (\bar{λ} / 4 π n) \int_{0}^{x} Δ ϕ (x^{'}, z, Ω) d x^{'},

where

x^{'}

is a dummy variable for integration.

The complex displacement, i.e. the argument of Re{} in Eq. (5), can be calculated by projecting $d (x, z, Ω)$ onto the basis given by the (unit-amplitude) complex PZT drive function:

\tilde{d} (x, z, Ω) = \int_{x - Δ x / 2}^{x + Δ x / 2} d (x^{'}, z, Ω) e^{i (Ω x^{'} / v_{x} + ϕ_{0})} w (x^{'}) d x^{'} .

This is analogous to performing a short time Fourier transform (STFT) of the displacement amplitude over

Δ t = Δ x / v_{x}

. A window function,

w (x)

, such as a Hanning window, can be used to shape the spatial point-spread function in the complex mechanical response images. The displacement amplitude is given by

| \tilde{d} (x, z, Ω) |

, and the displacement phase is given by

Φ (x, z, Ω) = \arg {\tilde{d} (x, z, Ω)}

.

The results in Sections 4.2 to 4.4 were processed in Matlab using the following parameters, unless otherwise noted. The lateral phase derivative was computed from adjacent A-scans (i.e. $m = 1$ ). Median filtering was performed separately on the real and imaginary parts of the complex phase difference map using a $5 \times 5$ kernel for vibration periods exceeding 20 transverse pixels, and with a $5 \times 3$ kernel for vibration periods less than 20 transverse pixels (the lower lateral kernel width was to maintain a kernel size that is small relative to the vibration period). Raw and median filtered phase difference maps were examined for evidence of large sample motions, and data sets corresponding to sample motions with $| Δ ϕ | \geq π$ were excluded from the analysis. Low-pass Butterworth filtering (third order) of the complex phase difference map was performed with a cutoff at $2.4 f_{vib}$ , where $f_{vib} = Ω / 2 π$ is the vibration frequency. The band-pass filter (third order) had a 20% bandwidth, i.e. the pass band was given by $[0.8 f_{vib} 1.2 f_{vib}]$ . The digital filtering was performed using the ‘filtfilt’ command, so as to minimize phase distortion. The PZT was driven with a sinusoidal waveform given by $V_{P Z T} (t) = V_{0} \cos (Ω t + ϕ_{0})$ , with complex analytic representation ${\tilde{V}}_{P Z T} (t) = V_{P Z T} (t) + i H {V_{P Z T} (t)}$ . The integration in Eq. (9) was performed over 3 vibration cycles, a (unit area) Hanning window was utilized for $w (x)$ , and a factor of two, i.e. $2 w (x)$ , was applied for amplitude correction. For all results, sample displacement was calculated using a sample background refractive index of 1.4. OCT images are displayed on a decibel scale that was calculated using a reference noise floor obtained from an optical depth of approximately 725 μm. A ramp SNR mask was used to attenuate the amplitude of the computed complex mechanical response in regions of low OCT signal. This mask was used since the phase difference map has random phase in the absence of OCT signal, which has an in-band component that passes the band-pass filter. When calculating image SNR for the mask, an averaging kernel with dimensions 45 μm (depth) × 12 μm (lateral) was first applied to smooth speckle. The multiplicative mask was unity in regions above a threshold of SNR = 3, and for pixels below the threshold, the mask values were equal to the linear scale SNR/3.

3.3 Data acquisition and processing constraints

The B-mode data acquisition parameters simultaneously determine the spatial sampling of the OCT image as well as the temporal sampling of the vibration waveform. In order to resolve dynamic displacements at the vibration frequency, the data acquisition parameters should provide respective spatial and temporal sampling that exceeds the Nyquist criterion. The data processing outlined in Section 3.2 assumes that the phase is correlated between every m ^th A-scan used to compute the transverse phase derivative, in order to minimize the contribution of sample structure ${\tilde{S}}_{0} (x, z)$ to the final displacement amplitude. (A lateral band-pass filter applied to ${\tilde{S}}_{0} (x, z)$ , i.e. in the absence of applied vibration, would isolate in-band contributions resulting from sample structure, or from phase noise). This condition is closely tied to the requirement of phase correlation of the vibration signal in order to calculate a valid Doppler-based OCE signal [8].

Another practical constraint on data acquisition is to maintain a phase shift between A-scans within the range $(- π, π)$ . Phase shifts that exceed this range require phase unwrapping to accurately reconstruct the waveform representing sample motion. Maintaining phase shifts within the range $(- π, π)$ requires that

- π < \frac{m v_{z} (x, z, Ω)}{f_{A-scan} (\bar{λ} / 4 π n)} < π,

where

v_{z} (x, z, Ω)

is the axial velocity of the sample, and the phase difference is calculated between every m ^th A-scan recorded at an A-scan frequency

f_{A-scan}

. Values of

m > 1

can be beneficial when the phase differences between adjacent A-scans are comparable to or smaller than the phase noise. Equation (10), along with the spatial and temporal sampling constraints discussed above, provide the constraints on the data acquisition parameters, such as the A-scan rate, the lateral field of view and number of A-scans per image, and also on the choice of the lateral A-scan spacing over which to compute the phase difference. As can be seen from Eq. (10), fixed values of

f_{A-scan}

and m will cover a given range of velocities within the sample.

Tuning the vibration frequency over a wide range, while keeping the PZT amplitude constant, requires a corresponding coverage of a wide range of velocities $v_{z}$ (the velocity $v_{z}$ is proportional to the PZT frequency). Coverage of a wider range of sample velocities than possible at a single A-scan frequency was achieved by acquiring separate hyperspectral data sets at three different A-scan frequencies, similar to an earlier study [23]. This ensures adequate temporal sampling at all vibration frequencies. Additionally, for a given vibration frequency, the A-scan rate also determines the image acquisition time and, therefore, the number of vibration cycles within each image, which limits the OCE spatial resolution [17].

3.4 Silicone phantom fabrication

Tissue phantoms with similar optical scattering and mechanical properties to tissue were used in the experiments. Silicone-based tissue phantoms were chosen due to their permanence and the ability to control their elastic moduli. Phantoms were fabricated from pure polydimethylsiloxane (PDMS) fluid (50 cSt viscosity, ClearCo, Inc.), a room-temperature vulcanizing silicone, and its associated curing agent (General Electric RTV-615 A and B, respectively, Circuit Specialists, Inc.). To obtain different elastic moduli, a series of phantoms were fabricated with different concentration ratios of these three components. Polystyrene microspheres (mean diameter 0.53 µm, Bangs Laboratories, Inc.) were embedded in the tissue phantoms with a concentration of 1 mg/g to function as optical scatterers for OCE imaging. Due to smaller refractive index variations, phantoms with polystyrene microspheres enable relatively deeper penetration depth for OCT imaging [30]. The phantom solutions were mixed thoroughly in an ultrasonicator for 30-60 minutes at room temperature and then poured into 35 mm plastic Petri dishes, with a thickness of approximately 2 mm for single-layer phantoms. The phantoms were then cured at 80°C for approximately 12 hours, and subsequently at room temperature for 24 hours. Tissue phantoms with inclusions were also fabricated. A homogeneous silicone phantom with elastic modulus of 195 kPa was first ground into particles with average dimensions ~1 mm. These particles were mixed with silicone fluid to create a background elastic modulus of 66 kPa, and then cured. All silicone tissue phantoms were sectioned into cylindrical shapes with diameter of about 12 mm for experiments.

4. Results

4.1 Characterization of the mechanical hardware

These measurements were conducted to facilitate the separation of the sample contributions from the frequency-dependent response of the overall coupled system, by determining the natural vibration modes (or resonances) of the hardware. The mechanical frequency response of the hardware was characterized with and without single layer silicone samples coupled between the PZT actuator and glass wedge. Each measurement comprised a vibration frequency sweep from DC to 1 kHz during acquisition of an M-mode image, which was acquired with a 10 kHz A-scan frequency and an exposure time of 95 μs. The data were analyzed by computing the phase difference between A-scans, which is proportional to the Doppler velocity. Spectra were calculated from the phase difference map, at a depth of approximately 160 μm, and are shown in Fig. 3 .

Fig. 3 Frequency response spectra of the mechanical hardware calculated from M-mode chirp measurements. (a) Piezo actuator with coverslip affixed, and (b) Coupled mechanical response of PZT, mounting hardware and sample, for four single-layer silicone samples with elastic moduli of 36, 66, 101 and 141 kPa. Note: For clarity, the spectra for 101 kPa and 141 kPa are shown with 5 dB and 10 dB subtracted respectively.

Download Full Size | PDF

Figure 3a shows the frequency response of the PZT actuator with attached coverslip. The overall increase in signal power with frequency is expected from constant displacement amplitude, due to the higher velocities associated with higher vibration frequency. There is a minor but sharp peak observed at ~44 Hz, a (negative) peak at 121 Hz, and a minor peak at ~324 Hz.

Figure 3b shows the results from measurements of single layer silicone phantoms of varying stiffness that were sandwiched between the PZT coverslip and the glass wedge. It can be seen that there are several peaks that are common to two or more samples. We conclude that these peaks represent the natural vibration modes (or resonances) of the mechanical hardware. For all samples, a dominant peak that can be attributed to the mounting hardware is observed at ~122 Hz. As sample stiffness increases, it can also be observed that the signal power in the 500-1000 Hz range, relative to the peak at 122 Hz, is increasing. This can be attributed to the mechanical transfer function of stiffer samples favoring higher vibration frequencies. Minor peaks can be observed at frequencies in the range 44-970 Hz. Of interest for the subsequent results in this paper, is a minor peak at approximately 315 Hz, followed by a minimum at 326 Hz, both of which become more prominent with stiffer samples. Additionally, the peak at ~365 Hz, followed by a minimum at ~380 Hz, is relevant. The 380 Hz minimum also becomes more prominent with stiffer samples.

4.2 Silicone phantom with inclusions

A hyperspectral data cube was acquired over a frequency range of 25 Hz to 1 kHz, and over an OCT field-of-view of 2 mm, centered on a single inclusion (~1 mm). To satisfy Eq. (10), the data cube was assembled from three subsets – the first was acquired with a 1 kHz A-scan rate over a 20 Hz to 125 Hz excitation range, the second was acquired with a 5 kHz A-scan rate over a 100 Hz to 625 Hz excitation range, and the third was acquired with a 10 kHz A-scan rate over a 200 Hz to 1 kHz excitation frequency range. Data above 83 Hz from the first subset, and data above 417 Hz from the second subset was not utilized since the phase difference magnitude between adjacent A-scans was close to or exceeded π.

Figure 4 shows the results of imaging an inclusion laterally centered within the OCT field of view. The outline of the inclusion can faintly be distinguished from the OCT image – the inclusion-to-background scattering contrast is low due to the use of similar polystyrene concentrations. A lateral dependence can be observed in the 71 Hz displacement amplitude image, with higher displacement amplitude in the central region corresponding to the spatial location of the inclusion. The displacement amplitude image shows contrast along both axial and lateral dimensions that outlines the overall shape of the inclusion. Complex mechanical spectra computed from the background silicone matrix (blue ROI), and the inclusion (red ROI) are shown. An SNR threshold of 10 dB was applied to select the points within each ROI that were used to compute the average complex mechanical response, resulting in a similar number of points, 14,792 and 15,394 for the background and inclusion ROIs respectively. The mechanical spectra are consistent with the overall higher displacement amplitude of the inclusion. Additionally, the displacement amplitude spectrum of the sample inclusion for the frequency range above 120 Hz has a higher weighting relative to the 122 Hz hardware resonant mode. In particular, the sample response about the 313 Hz peak is more prominent in the ROI corresponding to the stiffer inclusion.

Fig. 4 Complex mechanical spectra of a silicone phantom containing an inclusion with closely matched optical scattering, but with a different elastic modulus of 195 kPa compared to the background matrix (66 kPa). (a) Displacement amplitude (top) and phase (bottom) of the background matrix, (b) OCT image (top) and displacement amplitude image corresponding to vibration frequency of 71 Hz (bottom), and (c) displacement amplitude (top) and phase (bottom) of the inclusion. Image dimensions are 750 μm (optical depth) by 2 mm (lateral). The solid lines in the displacement amplitude and phase plots represent a linear interpolation between the experimentally measured points.

Download Full Size | PDF

B-mode images of the complex mechanical response are shown in Fig. 5 . The displacement amplitude is found to increase between 109 Hz and 119 Hz, and the mechanical phase changes from an in-phase response towards $Φ = π / 2$ , as the vibration frequency gets closer to the dominant hardware resonant mode at 122 Hz. Above this hardware vibration mode, the displacement phase shows a layer with $Φ \approx π$ near the surface of the sample that is associated with reduced displacement amplitude. This region represents the interaction between the sample and the glass wedge located above the sample, which represents a vibration source driving the top of the sample with a mechanical phase of $Φ \approx π$ . A plausible interpretation of the reduced surface displacement is destructive interference between two counter-propagating pressure waves resulting from the motions with $Φ \approx π$ . These results suggest that both amplitude and phase images can be useful to detect the presence of spatially-varying mechanical properties, particularly around frequencies corresponding to resonance.

Fig. 5 Images of the complex mechanical response of a silicone phantom containing an inclusion at vibration frequencies of (a) 109 Hz, (b) 119 Hz, and (c) 125 Hz. Displacement amplitude images are shown in the top row, and phase images are shown in the bottom row. Image dimensions are 750 μm (optical depth) by 2 mm (lateral).

Download Full Size | PDF

4.3 Tumor margin

A hyperspectral data cube was acquired over the frequency range of 25 Hz to 1 kHz from a tumor margin in rat mammary tissue. The data cube comprised two subsets – the first was acquired with a 1 kHz A-scan rate over a 25 Hz to 167 Hz excitation range, and the second was acquired with a 10 kHz A-scan rate over a 313 Hz to 1 kHz excitation frequency range. Strain rate amplitude images calculated from this data set were previously presented at 45 Hz and 313 Hz, along with the same corresponding histology image [17]. Here we present complex mechanical spectra over the range of 25 Hz to 1 kHz, and show displacement amplitude and phase images computed from the same hyperspectral data cube.

Figure 6 shows complex mechanical spectra from ex vivo rat mammary tissue consisting of adjacent adipose and tumor tissue. The complex mechanical response was averaged over selected regions of interest (ROIs) corresponding separately to adipose and tumor. An SNR threshold of 10 dB was applied to select the points used to compute the average complex mechanical response, providing a similar number of above threshold points - 26,134 and 26,045 for the adipose and tumor regions, respectively. The standard error of the mean was calculated for the points comprising the amplitude spectra, but the error bars are not shown because they are on the order of 10⁻⁴ μm. The amplitude spectrum for the adipose ROI shows a distinctive peak at approximately 83 Hz that is not present in the tumor region. Both tumor and adipose ROIs show a prominent peak at 385 Hz and a minor peak at 313 Hz, but with lower amplitude in the adipose ROI. The minor peak at 313 Hz can be attributed to excitation of a mounting hardware vibration mode, but the major peak at 385 Hz approximately coincides with a minimum in the hardware frequency response (see Section 4.1). Additionally, the hardware vibration modes in the 300-500 Hz range are minor with respect to the peak at approximately 122 Hz, whereas the 385 Hz peak is prominent. Consequently, we attribute the 385 Hz spectral peak to a peak in the mechanical response $| \tilde{h} (x, z, Ω) |$ of tumor tissue.

Fig. 6 Complex mechanical spectra of a rat mammary tumor margin. (a) Adipose displacement and phase spectra averaged over the blue ROI box in the OCT image, (b) OCT image with nearby histology (adipose on the left and tumor on the right), and (c) tumor displacement and phase spectra averaged over the red ROI box in the OCT image. Image dimensions are 750 μm (optical depth) by 2 mm (lateral). The scale bar in the histology images is 200 μm. The solid lines in the displacement amplitude and phase plots represent a linear interpolation between the experimentally measured points.

Download Full Size | PDF

Images of displacement amplitude and phase are shown in Fig. 7 corresponding to spectral peaks in Fig. 6. At 83 Hz excitation, the displacement amplitude image highlights the adipose region of the sample and a central higher scattering region. The phase image shows that the majority of the sample is responding in phase with the applied vibration except for an oval-shaped region with sharply delineated boundaries (red arrow) that has a mechanical phase of $Φ \approx + π$ . The amplitude image also outlines the corresponding oval structure, but with lower contrast. Structures with similar dimensions can be observed in the histology image (Fig. 6b). This oval feature was present in the phase difference maps at vibration frequencies of 125 Hz and 167 Hz as well as below 83 Hz, down to the lowest frequency of 25 Hz (data at these frequencies not shown). This suggests that it has significantly different and sharply delineated mechanical properties from its surroundings, and is possibly a fluid-dense follicle or vacuole.

Fig. 7 Complex mechanical response images of rat mammary tumor adjacent to adipose tissue for vibration frequencies of (a) 83 Hz, (b) 313 Hz and (c) 385 Hz. Displacement amplitude images are shown in the top row, and phase images are shown in the bottom row. Image dimensions are 750 μm (optical depth) by 2 mm (lateral). Note that the maximum value of the displacement amplitude scale in (a) is ~0.12 μm, but is 0.20 μm in (b) and (c).

Download Full Size | PDF

Figure 7 also shows the displacement amplitude and phase at 313 Hz and 385 Hz, i.e. frequencies attributed to hardware resonance and tumor respectively. The images for 313 Hz and 385 Hz were computed using an A-scan spacing $m = 2$ . Apart from the central high-scattering region at the adipose-tumor margin, the amplitude images at 313 Hz and 385 Hz highlight complementary regions to that at 83 Hz. The amplitude image at 385 Hz faintly outlines the oval feature seen in the displacement phase and amplitude images at 83 Hz. At 313 Hz, the mechanical phase of $Φ \approx π / 2$ over the imaged depth, indicates a bulk resonant response. This coincides approximately with a minor vibration mode of the mounting hardware found at 315 Hz (see Section 4.1). The phase image at 385 Hz shows that the tissue bulk is responding in phase.

4.4 Muscle tissue

A three-dimensional data set was acquired from rat muscle tissue, corresponding to vibration frequencies in the range of 25 Hz to 167 Hz using an A-scan frequency of 1 kHz. Results are shown in Fig. 8 . The OCT image shows a heterogeneous, spatially localized tissue structure that corresponds to similar regions in the displacement amplitude image at 83 Hz. This indicates that spatially localized variations of biomechanical properties can be differentiated from the hardware resonances. Additionally, there are also distinct regions (orange, red and green arrows) in the 83 Hz displacement amplitude image that demonstrate different contrast compared to the OCT image. These regions have relatively high OCT scattering, but relatively low displacement amplitude.

Fig. 8 Complex mechanical response images for rat muscle tissue. (a) OCT image, and displacement amplitude (top) and phase images (bottom) at vibration frequencies of (b) 83 Hz and (c) 125 Hz. Image dimensions are 750 μm (optical depth) by 2 mm (lateral). Arrows in the displacement image at 83 Hz indicate regions with different OCE and OCT contrast.

Download Full Size | PDF

The complex mechanical response at 125 Hz clearly demonstrates excitation of the dominant hardware resonant mode at approximately 122 Hz. This is characterized by a large displacement amplitude and mechanical phase of $Φ \approx π / 2$ . Evidence of interaction between the region of the sample with a $Φ = π / 2$ phase shift and the underlying regions that are being driven in phase with the PZT can be observed in the 125 Hz phase image, at a depth of approximately 600 μm. Although the SNR decreases with increasing depth, a (dark blue) band with mechanical phase $Φ \approx - π$ can be discerned. (Note that a uniform in-phase response can be discerned at similar depth from the phase image at 83 Hz). As postulated for the 125 Hz silicone phantom images, the interaction of mechanical pressure waves with different phases may produce spatially-resolved structure in the mechanical phase response.

5. Discussion and future work

These results demonstrate spectroscopic contrast in displacement amplitude and phase data sets from silicone phantoms and biological tissue. In particular, the 83 Hz peak assigned to the adipose-rich region of the rat tumor and the higher 385 Hz peak assigned to the tumor tissue (Fig. 6) are consistent with softer adipose tissue and a stiffer tumor. This contrast was distinguishable from the natural vibration modes of the mechanical hardware based on the hardware characterization of Section 4.1 (Fig. 3). For the tumor margin case, the 385 Hz peak was attributed to tumor tissue on the basis that it corresponded to a local minimum in the hardware frequency response. This highlights the importance of characterizing the frequency-dependent response of the mechanical hardware. In general, it would be preferable to utilize a vibration-stable setup with a flat mechanical frequency response. The current setup was assembled from off-the-shelf optomechanical components, but a better alternative may be to use custom-designed plastic mounting hardware with a high damping coefficient.

For most cases presented here, the amplitude of the complex displacement at a given vibration frequency provided a greater level of contrast than the phase. Strong amplitude-based contrast was found in scattering silicone phantoms containing stiffer inclusions, which had scattering properties closely matched to the background matrix (Figs. 4 and 5). Higher displacement amplitude was also found to highlight complimentary regions of a rat adipose-tumor boundary, in particular at 83 Hz and 385 Hz (Fig. 6). In the muscle data at 83 Hz, distinct regions of lower displacement amplitude were found to correlate to optical scattering-based contrast present in the OCT image (Fig. 8).

Phase-based contrast at a given vibration frequency was generally associated with frequencies near the natural resonant modes of the mounting hardware. A system with N degrees of freedom will, in general, have N natural vibration modes [24]. Immediately above the 122 Hz hardware resonance mode, the silicon phantom with stiffer inclusions showed a top layer with $Φ \approx π$ relative to the underlying sample (Fig. 5). Phase-based contrast was also seen in the muscle tissue data at 125 Hz (Fig. 8). A notable exception (to mechanical phase contrast occurring near hardware resonant modes) is the oval structure in the phase image of the tumor sample at 83 Hz (Fig. 7), with boundaries that are most sharply delineated in the phase image. Further work is required to definitively correlate these findings to tissue pathology.

A biological sample may in general be described by a mechanically-coupled three-dimensional system with multiple degrees of freedom. Additionally, the local frequency-dependent amplitude or phase contrast may be influenced by the sample geometry and boundary conditions. Future work on spectroscopic OCE would benefit from theoretical modeling or FEM simulations [31] to quantitatively link the experimental findings to underlying elastic and viscous properties of the sample, and investigate the extent to which internal structure versus boundary conditions affect the local mechanical response of a sample. Simulations, based on an accurate mechanical model for tissue, provide a method to investigate the three-dimensional (both within and out of the imaging plane) dynamic interaction of pressure waves and mechanical coupling within a heterogeneous sample. For example, the presence of the 385 Hz peak in the adipose ROI (Fig. 6a) may be enhanced by mechanical coupling to tumor regions of the sample. The level of mechanical coupling between regions of a sample is also expected to reduce the spatial resolution with which local mechanical properties could be measured, and inverting the mechanical model may improve this resolution. Another interesting candidate for theoretical modeling or simulation is the interaction of coupled vibration modes that are closely spaced in frequency.

The transverse resolution limit of the displacement amplitude and phase images is inversely proportional to the band-pass filter bandwidth utilized to isolate the vibration frequency. When using a constant fractional bandwidth (20% in our case), this limit is proportional to the vibration frequency. Intuitively, as vibration frequency decreases, there are fewer vibration cycles within each image. The OCE resolution limit can be obtained from calculations in a previous study [Eq. (2), Ref. 17] by substituting the number of cycles used in the STFT calculation in place of the downsampling factor.

In this paper we have presented images based on displacement amplitude and phase. Another method is to display relative displacement amplitude by dividing the displacement amplitude images by a corresponding image acquired with the PZT turned off, permitting images to be displayed on a decibel scale [17]. This has the advantage of reducing the effects of the in-band contribution of phase noise in regions with low OCT SNR.

Further work is required to extend the vibration frequency coverage, to characterize the mechanical response of tissue with finer frequency resolution and over a wider frequency range. Extension of the frequency range to 10 kHz would require A-scan rates of several tens of kilohertz to 100 kHz, so as to exceed Nyquist sampling criteria for resolving the vibration waveform. In vivo spectroscopic OCE may be achieved through the use of a ring actuator design that performs co-incident mechanical excitation and optical imaging [16]. For imaging in this configuration, the mechanical rigidity of the underlying tissue or underlying bone provides a semi-rigid boundary against which to compress the imaged tissue.

6. Conclusions

We have presented a method to image the complex mechanical response of a viscoelastic sample, and perform mechanical vibration spectroscopy based on externally-applied excitation. This work has demonstrated, for the first time, mechanical spectroscopy of silicone phantoms and biological tissue using B-mode OCT imaging. Hyperspectral data cubes, comprising phase-resolved OCT images acquired for a range of discrete vibration frequencies, provide a rich source of information about the mechanical response of heterogeneous samples. Spectroscopic contrast based on displacement amplitude was observed between adipose and tumor regions of an ex vivo rat tumor margin. Displacement phase provided spatially-resolved contrast within the tumor region that was not evident in the OCT image, but could be associated with the dimensions of structures observed in the associated histology. The complex mechanical response also demonstrated spatially-dependent contrast in silicone phantoms with stiffer inclusions, and in ex vivo rat muscle. Further work is required to determine the frequency range that is best suited for applications such as tumor detection or sensing of atherosclerotic plaque vulnerability.

Acknowledgements

The authors gratefully acknowledge the laboratory assistance provided by Eric Chaney and Dr. Haohua Tu from the Biophotonics Imaging Laboratory, and discussions with Marko Orescanin from the Ultrasonic Imaging Laboratory, at the University of Illinois at Urbana-Champaign. We also gratefully acknowledge discussions with Dr. Amy Oldenburg, Dr. Daniel Marks and Adeel Ahmad. Brendan Kennedy acknowledges funding support from a UWA Research Collaboration Award. This research was supported in part by the National Institutes of Health (Roadmap Initiative, R21 EB005321; R01 EB005221; and RC1 CA147096 ARRA). Additional information can be found at http://biophotonics.illinois.edu.

References and links

1. W. Drexler, and J. G. Fujimoto, Optical Coherence Tomography: Technology and Applications. (Springer, New York, 2009).

2. J. M. Schmitt, “OCT elastography: imaging microscopic deformation and strain of tissue,” Opt. Express 3(6), 199–211 (1998). [CrossRef] [PubMed]

3. A. S. Khalil, R. C. Chan, A. H. Chau, B. E. Bouma, and M. R. Mofrad, “Tissue elasticity estimation with optical coherence elastography: toward mechanical characterization of in vivo soft tissue,” Ann. Biomed. Eng. 33(11), 1631–1639 (2005). [CrossRef] [PubMed]

4. G. van Soest, F. Mastik, N. de Jong, and A. F. W. van der Steen, “Robust intravascular optical coherence elastography by line correlations,” Phys. Med. Biol. 52(9), 2445–2458 (2007). [CrossRef] [PubMed]

5. J. Rogowska, N. A. Patel, J. G. Fujimoto, and M. E. Brezinski, “Optical coherence tomographic elastography technique for measuring deformation and strain of atherosclerotic tissues,” Heart 90(5), 556–562 (2004). [CrossRef] [PubMed]

6. H. J. Ko, W. Tan, R. Stack, and S. A. Boppart, “Optical coherence elastography of engineered and developing tissue,” Tissue Eng. 12(1), 63–73 (2006). [CrossRef] [PubMed]

7. S. J. Kirkpatrick, R. K. Wang, and D. D. Duncan, “OCT-based elastography for large and small deformations,” Opt. Express 14(24), 11585–11597 (2006). [CrossRef] [PubMed]

8. R. K. Wang, Z. H. Ma, and S. J. Kirkpatrick, “Tissue Doppler optical coherence elastography for real time strain rate and strain mapping of soft tissue,” Appl. Phys. Lett. 89(14), 144103 (2006). [CrossRef]

9. R. K. Wang, S. J. Kirkpatrick, and M. Hinds, “Phase-sensitive optical coherence elastography for mapping tissue microstrains in real time,” Appl. Phys. Lett. 90(16), 164105 (2007). [CrossRef]

10. J. F. Greenleaf, M. Fatemi, and M. Insana, “Selected methods for imaging elastic properties of biological tissues,” Annu. Rev. Biomed. Eng. 5(1), 57–78 (2003). [CrossRef] [PubMed]

11. M. Fatemi, A. Manduca, and J. F. Greenleaf, “Imaging elastic properties of biological tissues by low- frequency harmonic vibration,” Proc. IEEE 91(10), 1503–1519 (2003). [CrossRef]

12. M. Orescanin, K. S. Toohey, and M. F. Insana, “Material properties from acoustic radiation force step response,” J. Acoust. Soc. Am. 125(5), 2928–2936 (2009). [CrossRef] [PubMed]

13. X. Liang, and S. A. Boppart, “Dynamic optical coherence elastography and applications,” in Asia Communications and Photonics Conference and Exhibition, Technical Digest (CD) (Optical Society of America, 2009), paper TuG2.

14. X. Liang, A. L. Oldenburg, V. Crecea, E. J. Chaney, and S. A. Boppart, “Optical micro-scale mapping of dynamic biomechanical tissue properties,” Opt. Express 16(15), 11052–11065 (2008). [CrossRef] [PubMed]

15. S. G. Adie, B. F. Kennedy, J. J. Armstrong, S. A. Alexandrov, and D. D. Sampson, “Audio frequency in vivo optical coherence elastography,” Phys. Med. Biol. 54(10), 3129–3139 (2009). [CrossRef] [PubMed]

16. B. F. Kennedy, T. R. Hillman, R. A. McLaughlin, B. C. Quirk, and D. D. Sampson, “In vivo dynamic optical coherence elastography using a ring actuator,” Opt. Express 17(24), 21762–21772 (2009). [CrossRef] [PubMed]

17. X. Liang, S. G. Adie, R. John, and S. A. Boppart, “Dynamic spectral-domain optical coherence elastography for tissue characterization,” Opt. Express 18(13), 14183–14190 (2010). [CrossRef] [PubMed]

18. X. Liang, M. Orescanin, K. S. Toohey, M. F. Insana, and S. A. Boppart, “Acoustomotive optical coherence elastography for measuring material mechanical properties,” Opt. Lett. 34(19), 2894–2896 (2009). [CrossRef] [PubMed]

19. X. Liang and S. A. Boppart, “Biomechanical properties of in vivo human skin from dynamic optical coherence elastography,” IEEE Trans. Biomed. Eng. 57(4), 953–959 (2010). [CrossRef]

20. M. Fatemi and J. F. Greenleaf, “Ultrasound-stimulated vibro-acoustic spectrography,” Science 280(5360), 82–85 (1998). [CrossRef] [PubMed]

21. M. Fatemi and J. F. Greenleaf, “Vibro-acoustography: an imaging modality based on ultrasound-stimulated acoustic emission,” Proc. Natl. Acad. Sci. U.S.A. 96(12), 6603–6608 (1999). [CrossRef] [PubMed]

22. V. Crecea, A. L. Oldenburg, X. Liang, T. S. Ralston, and S. A. Boppart, “Magnetomotive nanoparticle transducers for optical rheology of viscoelastic materials,” Opt. Express 17(25), 23114–23122 (2009). [CrossRef]

23. A. L. Oldenburg and S. A. Boppart, “Resonant acoustic spectroscopy of soft tissues using embedded magnetomotive nanotransducers and optical coherence tomography,” Phys. Med. Biol. 55(4), 1189–1201 (2010). [CrossRef] [PubMed]

24. S. S. Rao, Mechanical Vibrations, (Addison-Wesley, Reading Massachusetts, 1986).

25. B. Ströbel, “Processing of interferometric phase maps as complex-valued phasor images,” Appl. Opt. 35(13), 2192–2198 (1996). [CrossRef] [PubMed]

26. C. J. Tay, C. Quan, and W. Chen, “Dynamic measurement by digital holographic interferometry based on complex phasor method,” Opt. Laser Technol. 41(2), 172–180 (2009). [CrossRef]

27. C. Quan, C. J. Tay, and W. Chen, “Determination of displacement derivative in digital holographic interferometry,” Opt. Commun. 282(5), 809–815 (2009). [CrossRef]

28. A. Szkulmowska, M. Szkulmowski, A. Kowalczyk, and M. Wojtkowski, “Phase-resolved Doppler optical coherence tomography--limitations and improvements,” Opt. Lett. 33(13), 1425–1427 (2008). [CrossRef] [PubMed]

29. R. A. Leitgeb, and M. Wojtkowski, “Complex and coherence noise free Fourier domain optical coherence tomography,” in Optical Coherence Tomography: Technology and Applications, W. Drexler and J. G. Fujimoto, eds., (Springer, New York, 2008).

30. B. W. Pogue and M. S. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,” J. Biomed. Opt. 11(4), 041102 (2006). [CrossRef] [PubMed]

31. P. J. Prendergast, “Finite element models in tissue mechanics and orthopaedic implant design,” Clin. Biomech. (Bristol, Avon) 12(6), 343–366 (1997). [CrossRef]

Spectroscopic optical coherence elastography

Abstract

1. Introduction

2. Theory

2.1 Mechanical resonance

2.2 OCT signal with applied vibration

3. Methods

3.1 Experimental setup

3.2 Data processing

3.3 Data acquisition and processing constraints

3.4 Silicone phantom fabrication

4. Results

4.1 Characterization of the mechanical hardware

4.2 Silicone phantom with inclusions

4.3 Tumor margin

4.4 Muscle tissue

5. Discussion and future work

6. Conclusions

Acknowledgements

References and links

Cited By

Figures (8)

Equations (10)

Optics Express