Elastography of Soft Materials and Tissues by Holographic Imaging of Surface Acoustic Waves References and Links

We use optical interferometry to capture coherent surface acoustic waves for elastographic imaging. An inverse method is employed to convert multi-frequency data into an elastic depth profile. Using this method, we image elastic properties over a 55 mm range with <5 mm resolution. For relevance to breast cancer detection, we employ a tissue phantom with a tumor-like inclusion. Holographic elastography is also shown to be well-behaved in ex vivo tissue, revealing the subsurface position of a bone. Because digital holography can assess waves over a wide surface area, this constitutes a flexible new platform for large volume and non-invasive elastography. On the inverse boundary value problem for linear isotropic elasticity, " Inverse Probl. Global teleseismic earthquake relocation with improved travel times and procedures for depth determination, " Bull. Elasticity and inelasticity of thermoplastic polyurethane elastomers: Sensitivity to chemical and physical structure, " Polymer (Guildf. Elastic moduli of normal and pathological human breast tissues: an inversion-technique-based investigation of 169 samples, " Phys. Prospective comparison of transient elastography, fibrotest, APRI, and liver biopsy for the assessment of fibrosis in chronic hepatitis C, " Gastroenterology 128(2), 343–350 (2005). Imaging of the elastic properties of tissue-a review, " Ultrasound Med. resonance elastography by direct visualization of propagating acoustic strain waves, Detection and discrimination of optical absorption and shear stiffness at depth in tissue-mimicking phantoms by transient optoelastography, " Appl. In vivo three-dimensional optical coherence elastography, " Opt. OCT-based elastography for large and small deformations, " Opt. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps, " Geophys. Surface acoustic wave depth profiling of elastically inhomogeneous materials, " J. Excitation and propagation of surface waves on a viscoelastic half-space with application to medical diagnosis, " J. Estimation of tissue's elasticity with surface wave speed, " J. Pulsed digital holography for deformation measurements on biological tissues, " Appl. Finding the position of tumor inhomogeneities in a gel-like model of a human breast using 3-D pulsed digital holography, " J. Imaging the mechanical stiffness of skin lesions by in vivo acousto-optical elastography, " Opt. Toward soft-tissue elastography using digital holography to monitor surface acoustic waves, " J. Quantitative evaluation of skin vibration induced by a bone-conduction device using holographic recording in the quasi-time-averaging regime, " in Digital Determination of physical property gradients from measured surface wave dispersion, " IEEE Trans. Algorithms and software for total variation image …


Introduction
In linear elastic theory, it has been posed that knowledge of an object's surface deformation in response to mechanical stress is sufficient to uniquely solve for the internal elastic properties [1].Solutions to this type of boundary value problem have been employed in seismology [2] and non-destructive testing [3] to map interior elastic properties of the earth and of structural materials, respectively.In many cases, elastography is limited by the inability to collect enough surface data to compute an accurate inversion.High frame rate holography is uniquely suited for capturing nanometer amplitude vibrations over a wide area, particularly in soft materials where surface acoustic waves have millimeter wavelengths at kilohertz frequencies.Here we demonstrate high resolution elastography by holographic imaging of surface acoustic (Rayleigh) waves at multiple interrogation frequencies.
The elastic properties of soft materials, such as polymers and human tissues, are of interest for their effects on electronic device performance [4], thermoplastics used in manufacturing [5], and medical diagnostics for breast cancer [6] and hepatic fibrosis [7].In breast cancer, malignant breast tissue exhibits an elastic modulus 3-10 × greater than its healthy counterpart [6], motivating the development of biomedical imaging modalities that contrast elasticity.Current methods of biomedical elastography [8][9][10][11][12][13] are based on tracking interior tissue motion during mechanical excitation, which requires (3 + 1)-D data acquisition (3 spatial + 1 temporal dimension).To reduce the dimensionality of the data acquisition, we propose an elastography method based on high frame rate holographic imaging of Rayleigh waves, a (2 + 1)-D imaging system.By computing the solution to the boundary value problem, one may in theory extract all 3 spatial dimensions of elasticity.
As known in seismology [14] and non-destructive testing [15], Rayleigh waves provide depth-resolved elastography because their sensitivity kernel (i.e., the material space over which propagation is affected by elasticity) extends in depth on a scale proportional to the wavelength, as illustrated in Fig. 1.In effect, lower frequency waves probe more deeply than higher frequency waves, so that the apparent wave dispersion provides depth-resolved elasticity information.The use of Rayleigh waves to quantify the elasticity of layers in soft materials has previously been performed by scanning a single detector away from a point source [16,17].Imaging the spatial pattern of surface waves on tissues has been accomplished by pulsed digital holography [18,19] or dynamic speckle imaging [20], the latter of which was used to generate 2-D surface strain maps.We recently demonstrated imaging of the real-time propagation of Rayleigh waves with high frame rate holography [21], and showed that they qualitatively exhibit the correct dispersive behavior in layered media.It was also recently shown that high frame rate holography can visualize vibrations induced in skin by a bone conduction device [22].To combine the 1-D depth resolution of Rayleigh waves with 2-D imaging, here we report, for the first time, the application of an inverse method for quantitative elastography by using a sensitivity kernel equal to the Rayleigh wave amplitude in a homogeneous half-space.From multi-frequency holography data we reconstructed depth profiles that represent 1-D axial elastograms, while scattering in the 2-D transverse plane is observed in response to a subsurface elastic inhomogeneity only when it is within the sensitivity kernel.We perform this method in an experimental context of soft biological tissues to demonstrate its validity for biomedical applications, by investigating phantoms matched to the elastic modulus of breast tissue.Our previous work [21] provides a technical basis for the holographic imaging method and the qualitative features of SAW propagation in layered and inhomogeneous phantoms.In this work we demonstrate, for the first time, quantitative elastography in layered phantoms, SAW propagation in a breast phantom with a 5mm tumor-like inclusion, and well-behaved SAW propagation on an ex vivo tissue sample.

Theoretical model and inverse method
In a homogeneous, elastic half-space, Rayleigh waves exhibit a phase velocity of [23]: where E is Young's modulus, ρ is mass density and ν is Poisson's ratio.For heterogeneous materials, where Young's modulus varies spatially, c R becomes frequency dependent.In the case of a material with a 1-D elasticity variation, E(z), as a function of depth, z, a common approach for determining c R is to divide the heterogeneous half-space into N homogeneous layers [15,24].Wave equations are then written for each layer, along with the conditions of continuity across boundaries between the layers, a free top surface, and finiteness at infinity.This leads to a system of 4N equations, and the determinant of this system contains an implicit relationship between c R and frequency.While this approach is exact within the constraints of the model, it requires tedious numerical calculation and is not easily amenable to an inverse method for reconstructing E(z).
We propose an alternate approach that uses a Rayleigh wave sensitivity kernel equal to the depth-dependent Rayleigh wave amplitude in a homogeneous medium, specifically, the z component of particle displacement, which takes the form of a double inverse exponential [23].We hypothesize that Rayleigh waves at a given frequency f in a medium with elastic depth profile, E(z), behave as they would in a homogeneous medium with an effective Young's modulus, E eff (f), given by the average of E(z) weighted by the sensitivity kernel for the effective homogeneous medium.Thus, our forward model has the form: where the sensitivity kernel k(f,z) is given by: α and β depend only on Poisson's ratio: 0.87 1.12 2 1 1 and c R (f)/f is substituted for the Rayleigh wavelength to explicitly show the frequency dependence of the kernel.For forward computation, Eq. ( 2) can then be solved iteratively by substituting E eff (f) into Eq.( 1) to obtain the dispersion curve c R (f) for a given E(z).This model shares similarities with other forward models for an effective c R in layered media [15,25,26], but these require perturbative variations in E(z).We chose to linearize in E instead of R c E ∝ because the conditions of continuity for elastic waves are linear in E, analogous to permittivity for electromagnetic waves.As will be shown in below, our forward model can predict the dispersion curve c R (f) even in highly non-perturbative tissue phantoms.
From this model, we can now straightforwardly solve the inverse problem of determining the elastic profile E(z) from an experimentally measured dispersion curve c R (f).We rewrite Eq. ( 2) in a discrete vector-matrix form, E eff = AE, with E being the elastic depth profile to be determined at each discrete depth z, and E eff a vector representing the measured effective elasticity obtained by inversion of Eq. (1) using the c R (f) measured at multiple frequencies f.The observation matrix A is the sensitivity kernel k(f,z) integrated over each discrete depth interval and normalized.Due to the Laplace transform nature of Eq. ( 2), the inverse problem is highly unstable, and regularization is needed.We find that Total Variation Regularization [27,28], which assumes a piecewise-continuous solution, works best with our tissue phantoms.That is, the vector E is obtained by minimizing: where ( ) is the total variation of E, and between adjacent elements of the Young's modulus vector.γ is a parameter determining the amount of regularization applied, the optimum value of which is chosen by the L-curve method [29].

Experimental setup
An in-line, image plane, digital holography system, illustrated in Fig. 2(a), was used to image the SAWs.Briefly, the system was composed of a Mach-Zehnder interferometer and lowpower, long coherence length laser.The beam from a 633 nm, 5mW Helium Neon laser (Thorlabs Inc, Newton, NJ) was passed through a 92:8 beam splitter, with 92% of the power directed along the sample path.The sample beam was then passed through a 2.5 × beam expander composed of two lenses, with the lens spacing adjusted so that the beam was slightly divergent.A downward reflecting mirror illuminated an approximately 30 mm diameter area of the sample, and the diffusely reflected light was directed through a 200 mm focal length, 2 in diameter imaging lens, forming an image onto the sensor of a high speed CMOS camera.A 1:1 magnification was used to image a 17 × 17 mm area of the sample onto the 1024 × 1024 pixel array.An aperture in front of the imaging lens was used to adjust the subjective speckle size to be within the spatial resolution of the camera.The beam along the reference arm was passed through a 14 × beam expander, and coupled through a spatial filter (a converging lens combined with a 15 µm pinhole).The reference and sample beams were combined on-axis, forming an interferogram on the camera sensor.A variable neutral density filter was used in the reference arm to adjust the intensity, and optimize the interference pattern within the dynamic range of the camera.The path lengths of the reference and sample arms were carefully adjusted to be well within the coherence length of the laser (~200 mm).
As waves propagated across the phantom surface, out-of-plane (z) surface displacement induced an optical path difference producing a phase shift in the interferogram, while in-plane motion was not observed as it was small compared to the subjective speckle size (typically 30-50 µm).
A piezo-electric transducer (Kinetic Ceramics Inc, Hayward, CA) with a maximum travel range of 20 µm was placed just outside the illumination area to excite SAWs on the sample.The contact area between the transducer and sample is considered a point, as its diameter was much smaller than the SAW wavelengths observed.The waves were allowed to reach a steady state (~2-3 seconds after excitation) before the interferograms were recorded.The piezo-electric transducer was sinusoidally driven at various frequencies between 70 Hz and 300 Hz.The framerate of the camera was chosen to be at least 10 × the frequency of the SAW, at values of 1 kHz, 2 kHz or 3 kHz, corresponding to respective array sizes of 1024 × 1024, 640 × 640 and 512 × 512 pixels.This resulted in 17 × 17 mm, 10.6 × 10.6 mm and 8.5 × 8.5 mm imaged areas of the sample, respectively.A typical interferogram obtained by the system is illustrated in Fig. 2(b).Note that the fringes in this image are attributed to a cover glass plate in front of the sensor; these fringes are stationary and do not provide information about the SAWs.Instead, the interference between reference and sample beams can be observed in Media 1 as the slowly changing speckle pattern in response to the SAWs moving across the surface.
The series of images recorded by the camera was then transferred to a computer, where a phase reconstruction program was used to first bandpass filter the time-series of interferograms, and then apply a temporal phase shifting algorithm detailed in Ref. 21.For sinusoidal excitation of the surface, the reconstructed images using this algorithm represent the phase of the surface acoustic wave, modulo π, as illustrated in Fig. 2(c).Note that we do not track SAW amplitude, as only the phase velocity c R is needed for elastographic mapping.First, two-layer phantoms were prepared to test the ability to extract E(z) by our inverse method.One phantom was composed of a "soft," 15 mm layer on top of a "stiff," 40mm layer, while a second phantom was identical to the first but with the layers switched.These were prepared in two steps.First, the 40 mm lower layer was mixed and allowed to cure in the container.The 15 mm upper layer was then poured over the cured lower layer, and the sample was then allowed to cure for another 24 hours.
Next, a phantom was prepared to simulate the mechanical properties of a tumor-like inclusion embedded inside otherwise healthy breast tissue, to test whether SAWs may be sensitive to tumors based on their mechanical properties.The inclusion was approximately cylindrical with a diameter of 5mm and height of 5mm, and embedded at a depth of ~9 mm below the surface of the phantom.This phantom was prepared in four steps.The first 40 mm, "soft" lower layer was cured in the phantom container.Simultaneously, a 35 mm diameter by 5 mm height "stiff" cylindrical sample was mixed and cured.A 5 mm height cylindrical piece was then cut from this sample to act as the tumor inclusion.This piece was then laid on top of the "soft" layer, surrounded by another 5 mm layer of "soft" mixture, and allowed to cure.Finally, the top-most "soft" layer of approximately 9 mm height was poured on top and allowed to cure.The final result was thus an almost homogeneous phantom of Young's modulus of 4.15 kPa, with an inclusion in the center at a depth spanning from 9 to 14 mm, with Young's modulus of 31.9 kPa.B-mode ultrasound images were then acquired of the phantom to verify the location of the inclusion before holographic imaging.The ultrasound system was a SonixTOUCH (Ultrasonix Medical Corporation, Richmond, BC, Canada).

Elastic depth profiling of two-layer phantoms
First, we investigated the ability to extract E(z) using our model in section 2.1 by imaging SAWs on two-layer phantoms with known Young's moduli.The experimentally obtained dispersion curves c R (f) for each two-layer phantom are illustrated in Fig. 3.As described in our earlier work [21], the slope of the dispersion curve depends on the relative stiffness change for increasing depth.From Eq. ( 1), it can be seen that stiffer materials (i.e. higher Young's modulus) exhibit a higher phase velocity than softer materials.Thus, the SAW speed at higher frequencies (i.e.lower wavelength, and therefore shallower probing) tends towards what would be expected for a homogeneous phantom composed purely of the material in the upper layer.Conversely, at lower frequencies, the speed of the SAW tends towards that of a phantom composed of the material in the lower layer.Thus, for the phantom with a "soft" upper layer and "stiff" lower layer, the dispersion curve has a negative slope, while the "stiff" upper layer and "soft" lower layer phantom dispersion curve has a positive slope.
The use of phantoms with known Young's moduli now allows us to test our inverse method for quantitatively extracting E(z).First, we compared the measured dispersion curves with those computed using our forward model (Eq.( 2)) and the known distribution of E(z).While the experimental results seem to produce consistently larger values of c R , the overall shape of the dispersion curves is in good agreement for both phantoms.Next, we performed the inverse method upon the measured dispersion curves to predict E(z), the results of which are illustrated in Figs.3(c) and 3(d) in comparison with the known E(z) according to the texture analyzer.For both phantoms, the Young's modulus distribution obtained by inversion matches the true distribution fairly well (< 2.5 kPa on average), and the transition depth between the layers was accurately predicted within one sampling point (∆z = 2.6 mm).For the phantom with the soft upper layer, the inversion shows a transition from the upper layer to lower layer between 15.7 and 18.3 mm (compared to the true value of 15 mm).The inversion result for the phantom with stiff upper layer shows a tighter match to the true distribution, with a transition between the two layers predicted between 13.1 and 15.7 mm.

SAWs on a homogeneous phantom with a tumor-like inclusion
Digital holographic imaging of the homogeneous phantom containing a tumor-like inclusion is illustrated in Fig. 4. To verify the location of the inclusion, a B-mode ultrasound image of a region containing the inclusion is shown in Fig. 4(a), in comparison to a region to one side of the inclusion (Fig. 4(e)).The horizontal lines in the ultrasound images correspond to the separately cross-linked layers of the phantom during the preparation process.The corresponding Rayleigh wave phase reconstructions, obtained by our digital holography method, are shown in the x-y plane next to each corresponding x-z ultrasound image in Fig. 4. Note that the frame rate was increased from 1 to 3 kHz as the excitation frequency was increased from 96 to 228 Hz, impacting the noise level on the phase reconstructions.
As can be seen, the Rayleigh wavefronts are increasingly perturbed at lower frequencies in the region containing the inclusion.The higher frequency waves and the waves not passing over the inclusion, on the other hand, remain unperturbed, appearing like waves in an elastically homogeneous medium.This is because the Rayleigh wavelengths in the background region decrease from 12 mm at 96 Hz to 5 mm at 228 Hz, so that the sensitivity kernel, which scales as the wavelength, does not extend significantly into the inclusion for the highest frequency waves.Conversely, the lower frequency waves, which have a deeper sensitivity kernel, exhibit scattering and refraction as they pass over the inclusion.

Imaging of an ex vivo tissue sample
To test the viability of this technique in real tissue, a raw chicken thigh was obtained for holographic imaging.This ex vivo tissue was chosen because of its geometry: the thigh had an ellipsoidal shape, with an approximate length of 105 mm, width of 68 mm and a depth of 40 mm, with a bone running through the middle along the length of tissue at a depth of 20 mm, measured using calipers.The bone, which is much stiffer than the surrounding tissue, is approximately a cylindrical rod of diameter 10 mm.An ultrasound image illustrating the location of the bone and corresponding digital holographic Rayleigh wave phase reconstructions at various frequencies are shown in Fig. 5.
Similar to the results obtained with the phantom containing an inclusion, the higher frequency Rayleigh waves are more circular and regularly spaced, with little apparent scattering, suggesting that the tissue at this depth scale (λ R ≈11 mm) is elastically homogeneous.The c R at 432 Hz was approximately 4.78 m/s, corresponding to an E eff of 72 kPa in the soft tissue region above the bone.The lower frequency Rayleigh waves, on the other hand, have a deeper sensitivity kernel and are highly deformed, consistent with the presence of the much stiffer, subsurface bone.Overall the Rayleigh waves appear to be wellbehaved in soft tissue, highlighting the feasibility of our approach for tissue elastography.

Conclusion
The results shown here demonstrate the potential for a novel elastography method using digital holography to image surface acoustic waves.In two-layer phantoms we were successfully able to reproduce the 1-D elastic depth profiles, up to a depth of 55 mm, using an ad hoc model of Rayleigh wave dispersion behavior.On average, the breast size in the U.S. population would require an imaging depth between 40 and 80 mm to reach the chest wall without pre-compression [31].To obtain this depth at average breast elasticity [6] would require sampling SAWs between 12 and 25 Hz, which could be obtained in our experimental setup by more onboard memory or by a lower available frame rate than 1 kHz.Importantly, this method provides greater penetration depth than other optical elastography methods reported to date [11][12][13], as it is not limited by the penetration depth of light into the tissue, but rather that of the Rayleigh wave.
It was also shown that real-time monitoring of Rayleigh wave propagation provided the ability to identify the approximate transverse position and depth of a tumor-like inclusion based on multi-frequency data.Thus, while the frequency-dependent behavior of SAWs contains information about the elastic properties as a function of depth, the transverse propagation contains information about the lateral profile.In future work, the application of a more complete elastic scattering model of Rayleigh wave behavior may allow one to obtain 3-D elastograms, through a similar inverse problem approach.
In this demonstration, importantly, the SAW propagation exhibited a well-behaved profile in both phantoms and real tissues that should be amenable to quantitative analysis.The ultimate goal is for the approach to be used to study in vivo human tissue.Before this can be achieved, some experimental and theoretical issues need to be addressed.For example, the effect of tissue motions other than that due to SAWs will need to be investigated.In this work, we found that by bandpass filtering at the excitation frequency, the system became more robust to motion artifacts.Chirped (frequency-swept) SAW excitation may also enable more rapid acquisition to minimize motion artifacts.We note that the nanometer displacement sensitivity of digital holography affords measurements using low strains and strain rates ( ε < 10 −5 and ɺ ε < 0.02 s −1 , respectively), which may respectively minimize artifacts due to nonlinear and viscoelastic mechanical properties of breast tissue.Another need will be the use of a larger imaging area and multiple transducers for whole breast imaging.The holography system is easily adaptable to larger field-of-view by modifying the imaging optics, as long as sufficient resolution to observe the smallest desired Rayleigh wavelengths (typically 3-4 mm) is maintained.Furthermore, wide field imaging provides a more complete set of data that leads to high resolution elastic reconstructions.
This elastography technique constitutes a novel, non-invasive approach that takes advantage of the flexibility, nanometer-scale displacement sensitivity, and wide area imaging provided by high frame rate digital holography.Furthermore, the use of SAWs provides depth-dependent elasticity information up to several centimeters into soft materials, which is of high relevance for investigating soft tissues and for cancer detection applications.

Fig. 1 .
Fig. 1.Cartoon diagram of depth-dependent out-of-plane (z) motion of Rayleigh waves as a function of wavelength (not to scale).Scanning the wavelength effectively scans the probing depth of the SAW.The amplitude of this motion is used as the sensitivity kernel in the proposed inverse method.

Fig. 3 .
Fig. 3. Dispersion curves obtained for samples with (a) stiff lower layer and (b) stiff upper layer.The corresponding inversions of the dispersion curve in (a) and (b) are shown in (c) and (d), respectively.

Fig. 4 .
Fig. 4. Images of a breast tissue phantom with a tumor-like inclusion.(a) B-mode ultrasound of a region with the inclusion (marked with dashed circle).SAW phase reconstructions of the surface over the inclusion are shown in (b) for excitation at 96Hz (Media 3) (c) for 168Hz (Media 4), and (d) for 228Hz (Media 5).The dashed circles in (b-d) correspond to the approximate lateral position of the inclusion.SAWs are perturbed by the inclusion only at low frequencies, when the SAW penetration depth is great enough to reach the inclusion.(e) Bmode ultrasound of the same phantom in a control region without the inclusion.SAW phase reconstructions of the surface in this region are shown in (f) for 96Hz (Media 6), (g) for 168Hz (Media 7), and (h) for 228Hz (Media 8).

Fig. 5 .
Fig. 5. Images of a raw chicken thigh with a bone through the middle.Ultrasound images are shown in (a), with the bone indicated by the dashed circle.SAW phase reconstructions of SAWs over the region containing the bone are shown in (b) for 84Hz (Media 9) (c) for 168Hz (Media 10), and (d) for 432Hz (Media 11).High frequency SAWs do not penetrate to the depth of the bone and appear more circular.