1.

Introduction

The behavior of low-frequency surface waves in biological tissues has received only limited attention in the literature. Yet tracking low-frequency surface acoustic wave propagation in tissues provides a noninvasive means of investigating the mechanical behavior of the tissue. Such an approach may yield diagnostic methods based on local differences in mechanical behavior of tissues that may be indicative of certain pathologies. Indeed, this is the basis of manual palpation.

Potts, Chrisman, and Buras¹ investigated the propagation and attenuation of low-frequency waves in human skin over a frequency range from near zero to 1 kHz. Their results indicated that at frequencies below a few hundred hertz, the waves propagated primarily along the skin surface. At higher frequencies, the waves tend to propagate as bulk waves. They inferred from these results that at low frequencies, the waves provide information from the surface layer (i.e., the stratum corneum), and at higher frequencies, information about the dermis is provided. Karzakov and Klochkov² investigated the mechanical properties of skin tissues in the human arm over a range of 10 to 100 Hz. One of their key results was that at these low frequencies, the speed of propagation was relatively slow, ranging between 21.0 and 1.6 m/s, respectively. They also noted that there is a very large interindividual variation in the speed of propagation at a given frequency. These results are consistent with those of Potts, Chrisman, and Buras.¹

Low-frequency surface wave propagation has been measured in other tissues besides skin. Ganesan, Man, and Lai-Fook³ evaluated surface waves in lung tissue as a function of transpulmonary pressure as a means of evaluating lung edema. At frequencies <10 Hz, they observed propagation velocities between 2.87 and 1.8 m/s, depending on the transpulmonary pressure. Kirkpatrick, Hinds, and Duncan⁴ evaluated the propagation and decay of 1-Hz acoustic waves in bovine nuchal elastin as a means of noninvasively determining the viscoelastic properties of this tissue for applications as a tissue engineering scaffold.

The purpose of this study was to evaluate low-frequency surface wave propagation and decay in porcine skin ex vivo and to directly relate the wave characteristics to the viscoelastic behavior of the skin. The goal of this study was to provide both a theoretical and empirical framework for future studies on the use of low-frequency acousto-optical elastography for potential clinical and biomechanical applications.

1.1.

Theoretical Development

1.1.1.

Physical model

Skin, while grossly anisotropic, can be considered to be locally isotropic. Thus, for this model, we assume a locally isotropic half-space subjected to a harmonic external normal force S_z=S₀ exp(iωt) over an area of radius a (Fig. 1). The time-dependent mechanical behavior of the solid material that fills the half-space can be described adequately by the constitutive relation

Eq. (1)

$${\mathit{\sigma }}_{\mathit{i}\mathit{j}}(\mathit{t}){\displaystyle {\int }_0^{\mathit{t}}\mathit{\lambda }}(\mathit{t}-\mathit{\tau })\mathit{\delta }\mathit{i}\mathit{j}\frac{\mathit{d}{\mathit{\epsilon }}_{\mathit{k}\mathit{k}}}{\mathit{d}\mathit{\tau }}\mathit{d}\mathit{\tau }+{\displaystyle {\int }_0^{\mathit{t}}2\mathit{\mu }(\mathit{t}-\mathit{\tau })\frac{\mathit{d}{\mathit{\epsilon }}_{\mathit{i}\mathit{j}}}{\mathit{d}\mathit{\tau }}\mathit{d}\mathit{\tau }}\text{.}$$

where σ_ij is stress, λ and μ are the two independent Lame´ constants, δ_ij is the Kronecker delta (δ_ij=1 if i=j; δ_ij=0 if i≠j) , ε_kk=ε₁₁+ε₂₂+ε₃₃ is strain, and τ is the time variable of integration. Note that the two viscoelastic functions may not have the same time dependence. Alternatively, engineering functions may be considered instead of the tensorial ones:

Eq. (2)

$$\begin{array}{c}G=\mu ,\\ E=\frac{G(3\lambda +2G)}{(\lambda +G)},\\ v=\frac{\lambda }{2(\lambda +G)},\\ B=2G\frac{1+v}{3(1-2v)},\end{array}$$

where G is the shear modulus, E is the elastic modulus, ν is the Poisson’s ratio, and B is the bulk modulus.

Figure 1

Isotropic, viscoelastic half-space subjected to a sinusoidal loading over an area of radius a.

Following Zhang et al.;⁵ and letting μ=μ₁+(∂/∂t)μ₂ and λ=λ₁+(∂/∂t)λ₂ where subscript 1 denotes the elastic component and subscript 2 denotes the viscous component, then for the harmonic excitation at angular frequency ω, we have

Eq. (3)

$$\begin{array}{l}\lambda ={\lambda }_1+i\omega {\lambda }_2,\\ \mu ={\mu }_1+i\omega {\mu }_2,\end{array}$$

Thus, from Eqs. (1), (2), and (3), the material is viscoelastic in nature, and the harmonic excitation will result in a nonzero imaginary part of E. Recalling that the half-space is subjected to a harmonically varying force S_z, the stress on the material varies sinsuoidally in time as

Eq. (4)

$$\mathit{\sigma }(\mathit{t})={\mathit{\sigma }}_{\mathit{o}}\exp (\mathit{i}\mathit{\omega }\mathit{t}).$$

The strain will also be sinusoidal, but will lag behind the stress by a phase, or loss, angle δ:

Eq. (5)

$$\mathit{\epsilon }(\mathit{t})={\mathit{\epsilon }}_{\mathit{o}}\exp (\mathit{i}\mathit{\omega }\mathit{t}-\mathit{\delta }).$$

The phase angle is related to the time lag between the sinusoids by

Eq. (6)

$$\mathit{\delta }=\frac{2\mathit{\pi }\Delta \mathit{t}}{\mathit{T}}$$

where T is the period of the waveforms. As a result of this phase lag, the complex dynamic modulus E^* is represented as a complex number:

Eq. (7)

$$\mathit{E}*=\frac{\mathit{\sigma }*}{{\mathit{\epsilon }}_0^{*}}=\mathrm{Re}(\mathit{E})+\mathrm{lm}(\mathit{E})\mathit{i}.$$

The first term on the right-hand side (RHS) of Eq. (7) is the in-phase component and is known as the storage modulus, while the second term on the RHS is the out-of-phase component and is commonly referred to as the loss modulus. The phase angle is a dimensionless measure of the viscoelastic damping of the material and represents the amount of energy lost (e.g., dissipated as heat) during each loading-unloading cycle. The ratio Im(E)/Re(E) is commonly referred to as the loss tangent or loss factor and is related to δ by the equation

Eq. (8)

$$\frac{\mathrm{lm}(\mathit{E})}{\mathrm{Re}(\mathit{E})}=\mathit{\omega }{\mathit{\tau }}_{\mathit{r}}=\mathrm{tan }\mathit{\delta },$$

where τ_r is the relaxation time of the material and is the ratio of the viscous-to-elastic components of the overall behavior of the material.

In the present scenario, the existence of δ should manifest as: 1. a decay in the amplitude of the propagating surface wave along $\overrightarrow{x}$ and 2. a significant lag between the imposed acoustic stress wave and the resulting speckle shift. The two manifestations are not unrelated to one another. To see this, it is of interest to relate tan δ with the characteristic decay distance, or relaxation distance x_δ, which is the propagation distance in which the magnitude of the propagating surface wave (which in the current case is measured by the magnitude of the laser speckle shift), decays to 1/e of its original value, A₀. The relaxation distance is given by

Eq. (9)

$$\mathit{A}={\mathit{A}}_0\exp (-\mathit{x}/{\mathit{x}}_{\mathit{\delta }}).$$

Recalling Stokes’ relaxation time formulation for viscosity, we can redefine tan δ in terms of density ρ₀ and dynamic viscosity η as

Eq. (10)

$$\mathrm{tan }\mathit{\delta }=\frac{\mathit{\omega }\mathit{\eta }}{{\mathit{\rho }}_0{\mathit{V}}_0^2},$$

where V₀ is velocity. Without any loss in generality, we can redefine the relaxation time τ_r as

Eq. (11)

$${\mathit{\tau }}_{\mathit{r}}=\frac{\mathit{\eta }}{{\mathit{\rho }}_0{\mathit{V}}_0^2}$$

with the denominator of Eq. (11) reflecting the well-known relationship between Young’s modulus E density, and wave velocity for an isotropic elastic body E=ρV². Thus, the viscous relaxation time and tan δ can simply be related to each other through Eqs. (8) and (11) as

Eq. (12)

$$\frac{\mathrm{tan }\mathit{\delta }}{\mathit{\omega }}=\frac{\mathit{\eta }}{{\mathit{\rho }}_0{\mathit{V}}_0^2}={\mathit{\tau }}_{\mathit{r}}.$$

τ_r is simply the time it takes for a wave to travel a distance x_δ. Tan δ and x_δ are then related to each other through τ_r. This point is discussed further later on.

1.1.2.

Surface acoustic (Rayleigh) waves

Consider again the situation of Fig. 1. The harmonic excitation launches a wave polarized in the xz (sagittal) plane with surface normal along $$\overrightarrow{z}$$ and propagating in the x direction. If a plane-strain situation is assumed, then the displacement and velocity components are in x and z directions and there is no coupling to the transverse waves with displacement along y [i.e., shear horizontal (SH) mode], perpendicular to the sagittal plane, so there is no mode conversion or reflection. Assume that the half-space is a semi-infinite solid with a free surface, such that the boundary conditions are that the tangential and normal stresses are zero on the surface at z=0.

Because the imposed displacements and therefore strains are small, we adopt the Lagrangian definition of stress and can express the general form of the displacements and the stress components, respectively, as

Eq. (13)

$$\begin{array}{c}{u}_x=\frac{\partial \varphi }{\partial x}-\frac{\partial \psi }{\partial z},\\ u_z=\frac{\partial \varphi }{\partial z}-\frac{\partial \psi }{\partial r},\\ {\sigma }_{zz}=\lambda \left(\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial z^2}\right)+2\mu \left(\frac{{\partial }^2\varphi }{\partial x^2}-\frac{{\partial }^2\psi }{\partial x\partial z}\right),\\ {\sigma }_{xz}=\mu \left(\frac{{\partial }^2\varphi }{\partial x^2}+2\frac{{\partial }^2\varphi }{\partial x\partial z}-\frac{{\partial }^2\psi }{\partial x^2}\right),\end{array}$$

where ϕ and ψ are potentials for the longitudinal z and transverse x directions, respectively, and λ and μ are the first and second Lame´ constants. The origin of Eq. (13) can be seen by defining a scalar and vector potential for the plane-strain situation described earlier:

Eq. (14)

$$\overrightarrow{\mathit{u}}=\overrightarrow{\nabla \mathit{\varphi }}+\overrightarrow{\nabla }\times \overrightarrow{\mathit{\psi }}$$

The potentials satisfy the wave equations

Eq. (15)

$$\begin{array}{l}{\nabla }^2\varphi =\frac{1}{c_L^2}\frac{{\partial }^2\varphi }{\partial t^2},\\ {\nabla }^2\overrightarrow{\psi }=\frac{1}{c_L^2}\frac{{\partial }^2\psi }{\partial t^2},\end{array}$$

for the longitudinal and shear waves, respectively, with the corresponding wave speeds, c_L and c_T. ⁶

Since we have assumed a plane-strain case, the only nonzero component of $$\overrightarrow{\psi }$$ is in the y direction. Recalling the harmonic time dependence, then the plane wave solutions for ϕ and ψ satisfying Eq. (15) are given by

Eq. (16)

$$\begin{array}{c}\varphi =A_1\exp \left[j(\omega t-\beta x-{\gamma }_{L}z)\right],\\ \psi =-j{A}_2\exp \left[j(\omega t-\beta x-{\gamma }_{T}z\right],\end{array}$$

where A₁ and A₂ are arbitrary constants and the wave numbers are defined by

Eq. (17)

$$\begin{array}{l}{\gamma }_L^2={\beta }^2-k_L^2,\\ {\gamma }_T^2={\beta }^2-k_T^2,\end{array}$$

Here k_L and k_T are the bulk wave numbers and β=k_R is the wave number of the harmonic wave traveling along the x axis.⁷

Noting that in the plane-strain case, the normal stresses vanish at the surface and thus σ_zz|_z=0=σ_xz|_z=0=0. Setting the determinant of the coefficients A₁ and A₂ equal to zero, then the characteristic Rayleigh equation may be obtained:⁷

Eq. (18)

$$({\mathit{\beta }}^2-{\mathit{\gamma }}_{\mathit{T}}^2)+4{\mathit{\beta }}^2{\mathit{\gamma }}_{\mathit{L}}{\mathit{\gamma }}_{\mathit{T}}=0.$$

Rewriting Eq. (18) in polynomial form with the definitions⁷

Eq. (19)

$$\begin{array}{l}\eta \equiv \frac{k_T}{\beta }=\frac{c_R}{c_T},\\ \xi =\frac{k_L}{k_T}=\frac{c_T}{c_L},\end{array}$$

we arrive at

Eq. (20)

$${\mathit{\eta }}^6-8{\mathit{\eta }}^4+8(3-2{\mathit{\xi }}^2){\mathit{\eta }}^2-16(1-{\mathit{\xi }}^2)=0.$$

Depending on Poisson’s ratio ν of the material as defined in Eq. (2), several roots to Eq. (20) arise. One of them, the real root with the smallest magnitude, corresponds to the existence of a Rayleigh surface wave. An approximation to this root is given using Bergmann’s formula⁷

Eq. (21)

$${\mathit{\eta }}_{\mathit{r}}=\frac{0.87+1.12\mathit{v}}{1+\mathit{v}}.$$

Thus, over the allowable ranges of 0<ν<0.5, the Rayleigh wave velocity varies from 0.87c_T⩽c_R⩽0.96c_T.

Using the definitions of Eq. (17), the characteristic equation, Eq. (18) may be rewritten in terms of wave speeds as:⁸

Eq. (22)

$$\mathit{r}-4\mathit{q}\mathit{s}=0,$$

where

Eq. (23)

$$\begin{array}{c}q={\left[1-{\left(\frac{c_R}{c_L}\right)}^2\right]}^{1/2},\\ s={\left[1-{\left(\frac{c_R}{c_T}\right)}^2\right]}^{1/2},\\ r=2-{\left(\frac{c_R}{c_T}\right)}^2.\end{array}$$

In this fashion, solutions for the particle displacements maybe obtained in an efficient form:⁸

Eq. (24)

$$\begin{array}{l}{u}_z=\left[A\exp (-qz)-2sq\exp (-sz)\cos \beta (x-c_{R}t)\right],\\ u_x=Aq\left[r\exp (-qz)-2\exp (-sz)\right]\sin \beta (x-c_{R}t),\end{array}$$

where A is the amplitude, β=2π/λ_R is the Rayleigh wave number, and λ_R is the Rayleigh wavelength. The Rayleigh wavelength is given by the standard equation

Eq. (25)

$${\mathit{\lambda }}_{\mathit{R}}=\frac{{\mathit{c}}_{\mathit{R}}}{\mathit{f}},\mathrm{where}\mathit{f}=\frac{\mathit{\omega }}{2\mathit{\pi }}.$$

Inspection of Eq. (24) reveals that there is an overall elliptical particle motion in the presence of Rayleigh waves. It is this surface particle motion that we wish to track by observing the translating laser speckle. Also, inspection of the equation reveals that u_z and u_x decay exponentially with z, and that since there is no frequency term in Eq. (20), the root C_R is nondispersive in an elastic body.

1.1.3.

Estimating speckle shift

The details of the numerical algorithm used to estimate the speckle pattern shift as a function of acoustic stress imposed by the propagating Rayleigh waves have been given elsewhere,⁴ so only the key features of the algorithm are presented herein.

The crux of the acousto-optical method is to determine the shift ζ of the speckle pattern as a function of record number (time) and acoustic stress. This shift was calculated using a maximum likelihood approach developed in earlier works.^{9 10} We adopted a frozen speckle model, much like Taylor’s frozen turbulence hypothesis, and assume that over a time on the order of a couple sequential exposures of the camera, the structure of the speckle pattern is fixed. The only change with time is its lateral motion. Thus, the speckle motion can be modeled as

Eq. (26)

$${\mathit{g}}_{\mathit{j}+1}({\mathit{x}}_{\mathit{i}})={\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}}-\mathit{\zeta }),$$

where the subscript i denotes the pixel (spatial dimension) and the subscript j represents the record (temporal dimension). We assume that the shift ζ is small compared to a pixel, so that Eq. (26) can be approximated as

Eq. (27)

$${\mathit{g}}_{\mathit{j}+1}({\mathit{x}}_{\mathit{i}})\approx {\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})-\mathit{\zeta }{\mathit{g}}_{\mathit{j}}^{\text{'}}({\mathit{x}}_{\mathit{i}}).$$

This is simply the first two terms of the Taylor series expansion for g. To introduce a degree of symmetry into the problem, we inspect the two speckle records on either side of the record of interest,

Eq. (28)

$$\left[{\mathit{g}}_{\mathit{j}+1}({\mathit{x}}_{\mathit{i}}),{\mathit{g}}_{\mathit{j}-1}({\mathit{x}}_{\mathit{i}})\right].$$

We then determine the ζ that minimizes the error,

Eq. (29)

$$e_j^2={\displaystyle \sum }_{i=1}^N{\left[g_{j+1}(x_i+\zeta )-g_{j-1}(x_1-\zeta )\right]}^2,$$

where the summation is over all pixels in the two records. We are thus seeking the ζ that brings these two records into registration. If we make use of the approximation in Eq. (27), then differentiating with respect to ζ and rearranging yields the formula

Eq. (30)

$${\widehat{\mathit{\zeta }}}_{\mathit{j}}=\frac{{\mathit{\Sigma }}_{\mathit{i}=1}^{\mathit{N}}\left[{\mathit{g}}_{\mathit{j}+1}({\mathit{x}}_{\mathit{i}})-{\mathit{g}}_{\mathit{j}-1}({\mathit{x}}_{\mathit{i}})\right]\left[{\mathit{g}}_{\mathit{j}+1}^{\prime }({\mathit{x}}_{\mathit{i}})-{\mathit{g}}_{\mathit{j}-1}^{\prime }({\mathit{x}}_{\mathit{i}})\right]}{{\mathit{\Sigma }}_{\mathit{i}=1}^{\mathit{N}}{\left[{\mathit{g}}_{\mathit{j}+1}^{\prime }({\mathit{x}}_{\mathit{i}})-{\mathit{g}}_{\mathit{j}-1}^{\prime }({\mathit{x}}_{\mathit{i}})\right]}^2}.$$

The term in the first square bracket in the numerator is proportional to the first central difference approximation to the derivative;

Eq. (31)

$$\frac{\partial {\mathit{g}}_{\mathit{j}}({\mathit{x}}_{\mathit{i}})}{\partial {\mathit{t}}_{\mathit{j}}}\approx \frac{{\mathit{g}}_{\mathit{j}+1}({\mathit{x}}_{\mathit{i}})-{\mathit{g}}_{\mathit{j}-1}({\mathit{x}}_{\mathit{i}})}{2}.$$

The spatial derivatives may be approximated similarly:

$$\frac{\partial {\mathit{g}}_{\mathit{j}+1}}{\partial {\mathit{x}}_{\mathit{i}}}\approx \frac{{\mathit{g}}_{\mathit{j}+1}({\mathit{x}}_{\mathit{i}+1})-{\mathit{g}}_{\mathit{j}+1}({\mathit{x}}_{\mathit{i}-1})}{2},$$

Eq. (32)

$$\frac{\partial {\mathit{g}}_{\mathit{j}-1}}{\partial {\mathit{x}}_{\mathit{i}}}\approx \frac{{\mathit{g}}_{\mathit{j}-1}({\mathit{x}}_{\mathit{i}+1})-{\mathit{g}}_{\mathit{j}-1}({\mathit{x}}_{\mathit{i}-1})}{2}.$$

Note that the speckle shift parameter ζ is the time rate at which the speckle pattern shifts. The units are pixels/record. This estimator was derived in terms of the shift that minimizes the mean-square error (MMSE) between adjacent records. It is also the maximum likelihood (ML) estimate of the shift.^{9 10}

2.

Methods and Materials

The fundamental method employed in this study was to use an acoustic speaker to supply a low-frequency, sinusoidally varying acoustic stress to the skin. Simultaneously, specific points on the skin were illuminated with coherent laser light and the motions of the backscattered speckle patterns as a function of the driving acoustic stress wave were recorded. The magnitude and the time-resolved pattern of the speckle shift were calculated using the maximum likelihood estimator described before.

All of the experiments were performed using samples of fresh porcine flesh obtained from a local abattoir. The tissue samples consisted of skin, adipose tissue, and underlying muscle tissue. Typical dimensions of the tissue sample were 20×10×5 cm in the x, y, and z dimensions, respectively. Tissues were kept refrigerated until approximately $$\frac{1}{2}$$ h before testing, at which time they were removed from refrigeration, surrounded on all sides, except for the skin surface, with damp paper towels, and allowed to equilibrate to room temperature. To apply the acoustic stress wave, a small speaker (87 W, 7.0 cm diam) was placed in direct contact with the skin at one end of the sample. The speaker was sinusoidally driven at 1 Hz by a sine wave originating from a function generator. In addition, the speaker was assumed to have remained in solid contact with the skin during excitation. Optical illumination was supplied by a 5-mW green HeNe laser (λ=543 nm) coupled through single-mode fibers. This wavelength was chosen because it restricts our interrogation depth to depths on the order of 100 μm or so, thus the measurements are primarily surface measurements.

The backscattered speckle pattern was imaged with a linear array charge-coupled device (CCD) camera that consisted of 5000 pixels with a 7-μm pitch. To reduce the amount of data acquired, only the middle 512 pixels were used. A 1:1 telecentric lens was mounted on the camera. The camera was triggered at 50 Hz. Each 1-D exposure of the camera was sequentially arranged into a 2-D array so that exposure number, or time, was along the vertical axis, and pixel number, or space, was along the horizontal axis. We have previously referred to such arrays as stacked-speckle histories.^{9 10 11}

2.1.

Rayleigh Wave Velocity and Wavelength

A standard two-beam configuration was employed to determine the velocity and wavelength of the 1-Hz surface acoustic waves (Fig. 2). The speaker was placed at one end of the skin samples and the HeNe laser was used to illuminate two spots on the skin of approximately 1 mm diam, so that the laser spots and the speaker formed a straight line along the long axis of the skin. The distance of the spot closer to the speaker was chosen rather arbitrarily at 2.5 cm. This distance was chosen with the awareness that at distances too close to the speaker, the skin is overdriven and the speckles show substantial decorrelation during the acoustic loading-unloading cycle, and at too great of distances, the magnitude of the speckle shift is restrictively small. The distance between the two laser spots varied between approximately 1.4 and 2.0 cm. The purpose of this was to ensure that regardless of the distance between the two spots, the calculated velocity and wavelength did not vary. The linear array CCD camera was aligned with the array parallel to the line formed by the speaker and the illuminated spots. A slight misfocus of 1 cm was purposely induced to ensure that the speckles translated in space as opposed to simply boiling, as would be the case for a perfectly focused system.¹²

Figure 2

Acousto-optical arrangement for determining the Rayleigh wavelength and velocity.

The Rayleigh wave velocity C_R was determined by the equation

Eq. (33)

$${\mathit{C}}_{\mathit{R}}=\frac{\mathit{\omega }\mathit{D}}{{\mathit{\varphi }}_1-{\mathit{\varphi }}_2},$$

where ω is the angular frequency of the driving waveform, D is the distance between the two laser spots, and ϕ₁−ϕ₂ is the difference between the phases of the sinusoids that describe the motion in the stacked speckle histories generated from each of the illuminated spots [Eqs. (26) through (32)]. Once the velocity C_R was determined, the Rayleigh wavelength λ_R could then be estimated from Eq. (25).

Finally, once both C_R and λ_R were determined, the relative particle motions in $$\overrightarrow{x}$$ and $$\overrightarrow{z}$$ were estimated using Eq. (24).

This procedure was conducted on three different locations each on five separate skin samples.

3.

Results

3.1.

Rayleigh Wave Velocity and Wavelength

A typical dual-stacked speckle history is shown in Fig. 3. This figure shows the sinusoidal motion of the backscattered speckle patterns from the two illuminated spots on the porcine skin sample. The spots were separated by a distance D=2.0 cm, and because of this separation, the sinusoids describing the motions are very slightly out of phase with each other (Fig. 4). This phase difference is the denominator of Eq. (33).

Figure 3

Dual stacked-speckle history showing the translation of the speckle patterns at two distinct spots on the skin separated by 2.0 cm due to the sinusoidal acoustic stress wave.

Figure 4

Cumulative speckle shift as a function of time for each of the speckle histories shown in Fig. 3. The two waves are very slightly out of phase with each other. For the dashed line, the peak shift of −2.25 pixels in the image plane corresponds to a peak shift of 0.16 μm of the skin surface.

Based on Eq. (33), the mean velocity C_R of the 1-Hz surface acoustic waves was estimated to be 2.9±1.8 m/s. No difference (student’s t -test, p>0.05) in C_R was found, regardless of the separation distance between the two spots, lending a degree of confidence to the measurements. Figure 5 shows the mean and standard deviation of the velocity measurements at each of the four separations investigated. Using the value C_R=2.9±1.8 m/s, the Rayleigh wavelength λ_R was then calculated using Eq. (25). Based on this equation, λ_R=2.9±1.8 m for a 1-Hz surface wave. This long wavelength indicates that our measurements are all in the acoustic near-field.

Figure 5

Velocity and standard deviation of the surface wave at four different separation distances between the two laser spots. There was no significant difference between the different separation distances in the measured velocity. Velocities are in m/s. The overall mean velocity C_R was 2.9±1.8 m/s.

With knowledge of C_R and λ_R, the relative particle displacements in u_x and u_z as a function of depth can be estimated via Eq. (24), allowing ν=0.45 (a reasonable number for Poisson’s ratio of biological tissue) and C_L=1500 m/s. The results of this theoretical prediction are plotted in Fig. 6. The primary result here is that both components of particle displacement have decay constants ∼1λ_R, indicating that the wave is in fact a surface wave. This depth, however, is substantially greater than the dimensions of skin, indicating that subsurface features may influence the behavior of the wave as observed on the surface. However, it should be appreciated that the actual displacements, even on the surface, are very small, typically submicron, and the attenuation of the wave is rapid. For example, from Fig. 5, which plots the cumulative speckle shift in the image plane for two illuminated spots located 2.5 and 4.5 cm, respectively, from the speaker, we can approximate the motion on the skin (object plane). For the spot located closest to the speaker (2.5 cm away, dashed line in Fig. 5), a peak shift of ∼2.25 pixels was observed for this specific experiment. For a camera pitch of 7 μm, this corresponds to a shift of 15.8 μm in the image plane. Since a 1:1 telecentric lens was used with a misfocus distance of 1.0 cm, this motion relates to a peak-to-peak shift in the object plane (i.e., on the skin) of only 0.16 μm. This peak shift was over $$\frac{1}{2}$$ of an acoustic cycle. Since the sample rate was 50 Hz, the 0.16-μm shift occurred over 25 samples, or 0.5 s, yielding a shift of ∼6.4 nm per record.

Figure 6

Relative amplitudes of the x and z components of displacement of the skin Rayleigh wave as a function of depth based on the theoretical Eq. (24). The Poisson’s ratio was assumed to be 0.45 and C_L was assumed to equal 1500 m/s.

3.2.

Decay in Amplitude of Surface Acoustic Wave

The relative amplitude of the surface acoustic wave, as measured by the magnitude of the speckle motion, decayed exponentially as a function of distance from the acoustic source. The results are plotted in Fig. 7. The line shows the fit to Eq. (34), and the characteristic decay distance x_δ≅2.9 cm. The equation describing the relative decay of the amplitude was

Eq. (35)

$$\mathit{y}=\exp [(-\mathit{x}/2.88)],$$

where y is the peak-to-peak magnitude of the speckle shift and the decay constant is in centimeters.

Figure 7

Relative decay in speckle shift as a function of distance from the speaker. The decay followed a negative exponential with x_δ≅2.9 cm, as indicated by the solid diamond.

4.

Discussion

The slow velocity C_R of the surface waves and corresponding long wavelengths measured in this study are consistent with the surface wave velocities measured in previous studies by other authors using very different methods.^{1 2 3} Thus, it is likely that the actual velocity of 1-Hz surface acoustic waves in skin is on the order of 3.0 m/s, and that the Rayleigh wavelength is of the same order. Vexler, Polyansky, and Gorodetsky¹³ examined the propagation velocity of shear waves in human skin at an acoustic frequency of 5.6 KHz and found that the velocities ranged between approximately 20 and 80 m/s, depending on the direction of wave propagation in relation to the anisotropic behavior of the skin. They used these velocities as indirect measurements of the viscoelasticity of the skin for determining the physiological status of skin.

As another confirmation of the value of C_R, Eq. (12) can be used to solve for τ_r. Using Eq. (12), and allowing tan δ=0.14, the viscous relaxation time τ_r=0.02 s. Physically, this can be interpreted that it takes 0.02 s for the surface acoustic wave to travel the distance x_d=2.9 cm. Using these values to solving for velocity C_R=1.45 m/s, which is lower than the directly measured velocity of 2.9±1.8 m/s, but still within the standard deviation of the measurements. A possible explanation for this discrepancy may be that the acoustic field is probably not quite flat, but is probably very slightly expanding. Recall that the diameter of the speaker was only 3 cm less than the width of the tissue samples. Thus, the measured decay in the relative amplitude of the acoustic stress wave with distance from the speaker is the result of some combination of wavefront expansion and the viscous nature of the tissue. If the wavefront was truly flat, the measured and calculated values of C_R would likely be even closer.

The decay in the amplitude of the acoustic wave with distance from the speaker is indicative of the viscous nature of the tissue. Although not investigated here, the characteristic decay constant should change with a change in frequency, as indicated by the description of the physical model from before. Furthermore, it is very likely that this decay constant will change with tissue conditions as well, including location of the skin sample on the body, hydration levels, age and physiological (i.e., healthy or diseased) state. These are certainly issues that need to be investigated.

The variation in tan δ values between samples is consistent with the variation seen between individuals in Karzakov and Klochkov’s² study on the low-frequency mechanical properties of skin tissues in the human arm. Potts, Chrisman, and Buras¹ also noted large variations between individuals. This variation simply reflects the natural variation in tissues, and the fact that the mechanical properties of skin are very dependent on location on the body and physiological parameters. It has been shown, however, that the values of tan δ of skin can be varied in a consistent fashion through the local injection of cross-linking agents (decreases tan δ) or enzymatic agents to partially digest the tissue (increases tan δ) .¹⁴

The long wavelength of the Rayleigh wave (∼3 m) brings into question the assertions that the skin sample can be modeled as a semi-infinite medium, and that the method evaluates the viscoelastic properties of the skin only, and not the entire tissue slab, including fat and muscle tissue. However, careful consideration of the results and of the optical sensing of the acoustic wave supports both assertions. First, the rapid attenuation of the surface wave as demonstrated by the short decay constant in the x and z directions by Fig. 6, along with the very small surface motions observed (see Results in Sec. 3), leads to the observation that the surface wave has attenuated dramatically before reaching a boundary, and reflections from boundaries are not a significant factor. Indeed, the smooth sinusoidal motion of the speckle patterns (Figs. 3 and 4) supports this observation. Reflections of the acoustic wave back into the measurement area would lead to period speckle decorrelation and nonsinusoidal shifts in the pattern of the speckle shift. Neither of these phenomena was observed. Thus, there is no reason to model the tissue slab as a waveguide, and it is better modeled as a semi-infinite medium, as was done here. Second, making the reasonable assumption that porcine skin has similar optical properties to human skin,¹⁵ then the approximate roundtrip probing depth of the 543-nm HeNe laser was on the order of 150 μm. Thus, the measurements of the viscoelastic behavior of the tissue were limited to about that depth, and not the entire thickness of the tissue slab.

Clearly, the variation in the measured mechanical properties of skin makes it somewhat impractical to rely solely on absolute values of the mechanical properties for diagnostic purposes. However, it seems reasonable that relative hydration levels of the skin could be determined in this fashion. On the other hand, variations in the local tissue mechanical properties may provide a means to interrogate the tissue. For example, if a large area of skin, including a suspected pathological site, were to be illuminated and a sequence of speckle images were taken as the skin is acoustically stressed, then neighborhood variations in the speckle shift across the images could be determined by extending the speckle shift algorithm to three dimensions (the third dimension being time).¹⁶ In this manner, the tissue surrounding the suspect tissue would serve as an internal control. Stiffer tissue areas would exhibit less speckle motion, while less stiff areas would exhibit greater motion under the same acoustic stress, indicating local variations in mechanical properties.