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 Buras1 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 Klochkov2 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.1
Low-frequency surface wave propagation has been measured in other tissues besides skin. Ganesan, Man, and Lai-Fook3 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 Duncan4 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.
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 Sz=S0 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σ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=ε11+ε22+ε33 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: G is the shear modulus, E is the elastic modulus, ν is the Poisson’s ratio, and B is the bulk modulus.
Following Zhang et al.;5 and letting μ=μ1+(∂/∂t)μ2 and λ=λ1+(∂/∂t)λ2 where subscript 1 denotes the elastic component and subscript 2 denotes the viscous component, then for the harmonic excitation at angular frequency ω, we have1), (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 Sz, the stress on the material varies sinsuoidally in time as 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: 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 τ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 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, A0. The relaxation distance is given byρ0 and dynamic viscosity η as V0 is velocity. Without any loss in generality, we can redefine the relaxation time τr as 11) reflecting the well-known relationship between Young’s modulus E density, and wave velocity for an isotropic elastic body E=ρV2. Thus, the viscous relaxation time and tan δ can simply be related to each other through Eqs. (8) and (11) as τ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.
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 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, asz 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: cL and cT. 6
Since we have assumed a plane-strain case, the only nonzero component of is in the y direction. Recalling the harmonic time dependence, then the plane wave solutions for ϕ and ψ satisfying Eq. (15) are given byA1 and A2 are arbitrary constants and the wave numbers are defined by kL and kT are the bulk wave numbers and β=kR is the wave number of the harmonic wave traveling along the x axis.7
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 A1 and A2 equal to zero, then the characteristic Rayleigh equation may be obtained:718) in polynomial form with the definitions7 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 formula7 0<ν<0.5, the Rayleigh wave velocity varies from 0.87cT⩽cR⩽0.96cT. 8 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 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 uz and ux decay exponentially with z, and that since there is no frequency term in Eq. (20), the root CR is nondispersive in an elastic body.
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,4 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 asi 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 g. To introduce a degree of symmetry into the problem, we inspect the two speckle records on either side of the record of interest, 27), then differentiating with respect to ζ and rearranging yields the formula 9 10
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 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
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.12
The Rayleigh wave velocity CR was determined by the equationD is the distance between the two laser spots, and ϕ1−ϕ2 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 CR was determined, the Rayleigh wavelength λR could then be estimated from Eq. (25).
Finally, once both CR and λR were determined, the relative particle motions in and were estimated using Eq. (24).
This procedure was conducted on three different locations each on five separate skin samples.
Decay in Amplitude of Surface Acoustic Wave
The physical model described before predicts an exponential decay in the magnitude of the acoustic wave with distance from the acoustic source [Eq. (9)]. To determine the characteristic decay distance, or relaxation distance xδ, which is the propagation distance in which the magnitude of the propagating surface wave as measured by the magnitude of the speckle shift decays to 1/e of its original value A0, an experimental arrangement similar to that described earlier was employed. However, only a single illumination beam was used, and the location of this illuminated spot was moved sequentially further away from the speaker. The magnitude of the speckle shift as a function of distance from the speaker was determined at each point distant from the speaker. As before, the speaker was driven sinusoidally at 1 Hz. Measurements were made at distances ranging from approximately 3.3 to 6.4 cm from the edge of the speaker. Beyond this latter distance, the shift in the speckle histories became prohibitively small (due to attenuation of the surface wave) to a point where the shift could no longer be analyzed with any degree of confidence.
Stacked speckle histories were generated for each point, and the magnitude of the speckle shift [Eqs. (26) through (32)] was taken to be proportional to the magnitude of the acoustic wave at each distance from the speaker. The magnitude of the speckle shift was plotted against distance, and the data were fit in a least-squares sense to an equation of the formx is the distance from the edge of the speaker along and the characteristic decay distance xδ was solved as a free variable.
This procedure was repeated on 3 separate samples of skin.
Tan δ of Porcine Skin
While the earlier velocity measurements relied on evaluating the phase difference between oscillating speckle patterns from two illuminated spots, the measurement of tan δ relies on determining the relative phase difference between the driving acoustic sinusoid and the sinusoidal shift in the speckle pattern from a single illuminated spot. In this case, the illuminated spot was located 5 cm from the edge of the speaker.
The speaker was again driven with a 1-Hz sinusoid from a function generator and the backscattered speckle pattern was imaged with the linear array CCD camera. As before, stacked speckle histories were generated and the speckle shift was plotted along with the driving waveform as a function of time (or record number).
Based on the measured Rayleigh wave velocity CR, a correction was made for time of flight of the surface acoustic wave and the phase difference δ between the driving acoustic wave. The resulting speckle shift was calculated. The phase of each waveform was determined as the arctangent of the ratio between the imaginary and real parts of the waves. The tangent of the phase difference was taken to be the desired variable, tan δ. The phase delay due to the electronics (i.e., the time delay between the time the signal is sent to the speaker and the coil actually moves) was considered to be negligibly small and was ignored in this analysis.
This procedure was conducted on three different locations each on five separate skin samples.
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).
Based on Eq. (33), the mean velocity CR 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 CR 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 CR=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.
With knowledge of CR and λR, the relative particle displacements in ux and uz 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 CL=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 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.
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 wasy is the peak-to-peak magnitude of the speckle shift and the decay constant is in centimeters.
Tan δ of Porcine Skin
Based on the phase differences between the driving acoustic sinusoid and the sinusoid describing the time-resolved speckle shift, the mean value of tan δ was 0.14±0.07, or a phase angle of approximately 8 deg. Table 1 shows the actual values for each skin sample and their standard deviations. A phase angle of 8 deg physically means that the strain in the skin sample lagged behind the imposed acoustic stress by 8 deg when the skin was subjected to a dynamic loading regime at 1 Hz.
|Measured tan δ values for five skin samples. Each value is the mean of three measurements.|
|Skin sample||Tan δ (mean±s.d.)|
The slow velocity CR 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 Gorodetsky13 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 CR, 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 xd=2.9 cm. Using these values to solving for velocity CR=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 CR 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’s2 study on the low-frequency mechanical properties of skin tissues in the human arm. Potts, Chrisman, and Buras1 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 δ) .14
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,15 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).16 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.
This grant was funded in part by NSF grant numbers BES-0196172 and BES-0201841, and NIH grant number 5R24EB000224-04. The authors wish to thank Deborah Baker and Jon Adams for their assistance in making some of the tan δ measurements.