Shear Wave Propagation in Soft Tissue and Ultrasound Vibrometry

Studies have found that shear moduli, having the dynamic range of several orders of magnitude for various biological tissues [1], are highly correlated with the pathological statues of human tissue such as livers [2, 3]. The shear moduli can be investigated by measuring the attenuation and velocity of the shear wave propagation in a tissue region. Many efforts have been made to measure shear wave propagations induced by different types of force, which include the motion force of human organs, external applied force [4], and ultrasound radiation force [5].


Introduction
Studies have found that shear moduli, having the dynamic range of several orders of magnitude for various biological tissues [1], are highly correlated with the pathological statues of human tissue such as livers [2,3]. The shear moduli can be investigated by measuring the attenuation and velocity of the shear wave propagation in a tissue region. Many efforts have been made to measure shear wave propagations induced by different types of force, which include the motion force of human organs, external applied force [4], and ultrasound radiation force [5].
In past 15 years, ultrasound radiation force has been successfully used to induce tissue motion for imaging tissue elasticity. Vibroacoustography (VA) uses bifocal beams to remotely induce vibration in a tissue region and detect the vibration using a hydrophone [5]. The vibration center is sequentially moved in the tissue region to form a two-dimensional image. Acoustic Radiation Force Imaging (ARFI) uses focused ultrasound to apply localized radiation force to small volumes of tissue for short durations and the resulting tissue displacements are mapped using ultrasonic correlation based methods [6]. Supersonic shear image remotely vibrates tissue and sequentially moves vibration center along the beam axis to create intense shear plan wave that is imaged at a high frame rate (5000 frames per second) [7]. These image methods provide measurements of tissue elasticity, but not the viscosity.
Because of the dispersive property of biological tissue, the induced tissue displacement and the shear wave propagation are frequency dependent. Tissue shear property can be modeled by several models including Kelvin-Voigt (Voigt) model, Maxwell model, and Zener model [8]. Voigt model effectively describes the creep behavior of tissue, Maxwell model effectively describes the relaxation process, and the Zener model effectively describes both creep and relaxation but it requires one extra parameter. Voigt model is often used by many researchers because of its simplicity and the effectiveness of modeling soft tissue. Voigt model consists of a purely viscous damper and a purely elastic spring connected in parallel. For Voigt tissue, the tissue motion at a very low frequency largely depends on the elasticity, while the motion at a very high frequency largely depends on the viscosity [8]. In general, the tissue motion depends on both elasticity and viscosity, and estimates of elasticity by ignoring viscosity are biased or erroneous.
Back to the year of 1951, Dr. Oestreicher published his work to solve the wave equation for the Voigt soft tissue with harmonic motions [9]. With assumptions of isotropic tissue and plane wave, he derived equations that relate the shear wave attenuation and speed to the elasticity and viscosity of soft tissue. However, Oestreicher's method was not realized for applications until the half century later.
One of potential applications of the ultrasound vibrometry is to characterize shear moduli of livers. The shear moduli of liver are highly correlated with liver pathology status [24,25]. Recently, the shear viscoelasticity of liver tissue has been investigated by several research groups [23,[26][27][28]. The most of these studies applied ultrasound radiation force in liver tissue regions, measured the phase velocities of shear wave in a limited frequency range, and inversely solved the Voigt model with an assumption that liver local tissue is isotropic without considering boundary conditions. Because of the boundary conditions, shear wave propagations are impacted by the limited physical dimensions of tissue. Studies shows that considerations of boundary conditions should be taken for characterizing tissue that have limited physical dimensions such as heart [22], blood vessels [19][20][21], and liver [8], when ultrasound vibrometry is used.

Shear wave propagation in soft tissue and shear viscoelasticity
The shear wave propagation in soft tissue is a complicated process. When the tissue is isotropic and modeled by the Voigt model, the phase velocity and attenuation of the shear wave propagation in the tissue are associated with tissue viscoelasticity. Oesteicher documented the detailed derivations of the solution of the sound wave equation for Voigt tissue [9]. We extended the solution to other models [8] for the applications of ultrasound vibrometry [8]. In this section, we provide the simplified descriptions of the shear wave propagation in tissue modeled by Voigt model, Maxwell model, and Zener model.
Assuming that a harmonic motion produces the shear wave that propagates in a tissue region, the phase velocity cs(ω) of the wave can be estimated by measuring the phase shift Δϕ over a distance Δz: The phase velocity is associated with the tissue property, which can be found by solving the wave equation with a tissue viscoelasticity model. For a small local region, the wave is approximated as a uniform plane wave, which has a simple form in isotropic medium: where S is the phasor notation of the displacement of the time-harmonic field of the shear wave, z is the wave propagation distance which is perpendicular to the direction of the displacement of the shear wave, and the complex wave number is The solution of (2) is a standard solution of a homogeneous wave equation: where S0 is the displacement at z = 0, is an unit vector in x direction. The plane wave is independent in y direction. The real time time-harmonic shear wave is: Although attenuation coefficient α = -ki carries information of the complex modulus of tissue, the phase measurement is often more reliable because it is relatively independent to transducers and measurement systems. The phase velocity is the speed of the wave propagating at a constant phase, which is a solution of ( ) / 0 The complex wave number k of the plane shear wave is a function of the frequency and the complex modulus of the medium [9]: where ρ is the density of the tissue and the complex modulus that connects stress σ and strain ε: which describes the relationship between stress and strain in the Voigt tissue. The Voigt model consists of an elastic spring μ1 and a viscous damper μ2 connected in parallel, which represents the same strain in each component as shown in Figure 1. The relation between stress σ and strain ε of the Maxwell tissue is: For a harmonic motion, (9) becomes: which is the same as (8). Substituting (8) into (7) and finding the real part of the wave number, the phase velocity of the shear wave in Voigt tissue can be obtained from (6): The elasticity μ1 and viscosity μ2 are two constants and independent to the frequency.
A numerical example of phase velocity of Voigt tissue is shown in Figure 2. Equation (11) shows that cs(ω) increases at the rate of square root of the frequency and there is no the upper limit for cs(ω). As shown in the Figure 2, the phase velocity is determined by both elasticity and viscosity. Ignoring the viscosity introduces errors and biases for elasticity estimates. However, examining the velocities at the extreme frequencies is useful for understanding the model and obtaining initial values for numerical solutions of μ1 and μ2.
In tissue characterization applications, μ1 is often in the order of a few thousands and μ2 is often less than 10. Thus, when the wave frequency is very low (less than a few Hz), When the frequency is very high (higher than a few tens of kHz), A broad frequency range is needed to accurately estimate both μ1 and μ2. (12) and (13) are only useful for estimating initial values for the numerical solutions of (11) with measured velocities, and they should not be used for final estimates.
Equation (7) can be used for other models for the plane shear wave having a single frequency. The Maxwell model consists of a viscous damper η and an elastic spring E connected in series, which represents the same stress in each component, as shown in Figure 3. The relation between stress σ and strain ε of the Maxwell tissue is: For a harmonic motion, (14) becomes: fixed, the complex modulus is a function of and E. Substituting (15) into (7), the shear wave speed in Maxwell medium can be found from (6): Equation (16) can be also obtained by replacing μ1 and μ2 of (8) with the real and imaginary terms of (15).
A numerical example of phase velocity of Maxwell tissue is shown in Figure 4. Note that cs(ω) gradually increases to a limit that is proportional to the square root of the elasticity. As shown in the Figure 4, the phase velocity is determined by both elasticity and viscosity. However, examining the velocities at the extreme frequencies is useful for understanding the model and obtaining initial values for numerical solutions of E and η.
for a very small ω, cs(ω) is zero for ω=0, and cs(ω) approaches / E  when ω is very high.  The Zener model adds an additional elastic spring, having the elasticity of E1, to the Maxwell model (η, E2) in parallel. The Zener model combines the features of the Voigt model and the Maxwell models and describes both creep and relaxation. Based on the Maxwell model, the complex shear modulus of the Zener model can be readily obtained: Substituting (17) into (7), the shear wave speed in Zener medium can be found from (6): Equation (18) shows that for a very small ω, η is proportional to the slop of the speed curve, and cs(ω) approaches   when ω is very high. A numerical example of phase velocity of Zener tissue is shown in Figure 6.

Ultrasound vibrometry
Ultrasound vibrometry has been developed to induce shear wave in a tissue region, measure phase velocity of the shear wave, and calculate the tissue viscoelasticity based on (11), or (16), or (18). The basics of the ultrasound vibrometry are described in details in references [11][12][13][14][15][16][17]32]. Ultrasound vibrometry induces tissue vibrations and shear waves using ultrasound radiation force and detects the phase velocity of the shear wave propagation using pulse-echo ultrasound.
From the solution of the wave equation, equation (5) can be represented by a harmonic motion at a location, where s=2fs is the vibration angular frequency, the vibration displacement amplitude D and phase s depend on the radiation force and tissue property. (19) is another representation of (5). Applying detection pulses to the motion that causes the travel time changes of detection pulses and phase shift changes of the return echoes, the received echo becomes [11]: where T is the period of the push pulses shown in Figure 9 and the modulation index is: where c is the sound propagation speed in the tissue, 0 is the angular modulation frequency of detection tone bursts, g(t,k) is the complex envelope of r(t,k), 0 is a transmitting phase constant and  is an angle between the ultrasound beam and the tissue vibration direction.
Received echo r(t,k) is a two-dimensional signal. When one detection pulse is transmitted, its echo from the different depth of tissue is received as t changes. In medical ultrasound field, variable t is called fast time. When multiple detection pulses are transmitted, the multiple echo sequences are received as k changes. Variable k is called as slow time. r(t,k) in fast-time t is called as fast-time signal to represent the echo signal in beam axial direction or the depth location in the tissue. Its variation in slow-time k is called slow-time signal to represent the signals from one echo to another echo. If there is no tissue motion, r(t,k) will be the same for different k values. The tissue motion information is carried by modulation index β and phase s. A quadrature demodulator is used to obtain β and phase s.
As shown in Figure 7, a quadrature demodulator is applied to extract the motion information from r(t,k). The complex envelop consists of the in-phase and quadrature term [29]: Operating on the in-phase and quadrature components I and Q with input r(t,k), we obtain the tissue motion in slow time [11]: A phase constant can be added to the local oscillator of the demodulator [11] to avoid zeros in I. The signal extracted by (23) is proportional to the displacement of a harmonic motion induced by the push pulses.
Shear Wave Propagation in Soft Tissue and Ultrasound Vibrometry 9 Another motion detection method [14] uses a complex vector that is a multiplication between two successive complex envelops [29]  Thus, the motion velocity in slow time can be obtained,  (19) and (25)..
The slow-time signal s(t,k) represents the tissue motion at a particular location, its amplitudes and phases change over distances are described by (5). The measurements of amplitudes and phases at two locations are used to calculate attenuation and phase velocity. As shown in (1), the phase velocity is related to the frequency and inversely related to the phase difference Δϕ over a short distance Δd. Thus, estimating the phase differences is the key step of the ultrasound vibrometry. The phase difference can be obtained by comparing phases ϕs of the slow-time signals s(t,k) at two locations z and z+Δz: There are several methods to estimate the phases of slow-time signals: Fourier transform, correlation method, and Kalman filter [14]. The estimated phase of the slow-time signal at a location include some phase constants due to the tissue location t and different pulse k, and phase ϕs = -krz. Given a tissue location (axial location) in fast time, all constant phases are removed by (26) except the phase shift ϕ in the lateral location.
Ultrasound vibrometry is developed to induce the shear wave described by (19) and detect the phase shift ϕ described by (26) for characterizing the tissue shear property using (1) and (11), (14), and (16). Ultrasound virbometry uses interleaved periodic pulses to induce shear wave and detects the phase velocity of shear wave propagation using pulse-echo ultrasound. Figure 8 shows an application setup of the ultrasound vibrometry. An ultrasound transducer transmits push beams to a tissue region to induce vibrations and shear waves. The push beams are periodic pulses that have a fundamental frequency fv and harmonics nfv. During the off period of the push pulses, the detection pulses are transmitted and echoes are received by the transducer at lateral locations that are away from the center of the radiation force applied, as shown in Figure 9. In some of our applications, fundamental frequency fv of the push pulses is in the order of 100 Hz, and pulse repetition frequency fPRF of the detection pulses is in the order of 2 kHz.  There are different variations of the excitation pulses beside the on-off binary pulses: continuous waves [11], non-uniform binary pulses [15], and composed pulses or Orthogonal Frequency Ultrasound Vibrometry (OFUV) pulses [30,31]. The OFUV pulses can be designed to enhance higher harmonics to compensate the high attenuations of high harmonics. The OFUV pulses have multiple binary pulses in one period of the fundamental period [30,31]. Other variations of the ultrasound vibrometry include consideration of background motion and boundary conditions that require more complicated models of tissue motions [13] and wave propagations [22].

Finite element simulation of shear wave propagation
Simulations using Finite Element Method (FEM) were conducted to understand the shear wave propagation in tissue. The simulation tool is COMSOL 4.2. The simulated tissue region is a two-dimensional axisymmetric finite element model of a viscoelastic solid with a dimension of 100 mm × 100 mm, as shown in Figure 10. The size of domain Ω1 is 100 mm × 80 mm. The domain is divided to 25,371 mesh elements and the average distance between adjacent nodes is 0.95 mm. The schematic diagram shown in Figure 10 includes simulation domains (Ω1, Ω2, Ω3) and boundaries (B1,B2). A line source (with a length of 60 mm) in the left of the solid represents as an excitation source of the shear wave. where elasticity modulus µ1 and viscosity modulus µ2 were set to be 2 kPa and 2 Pa*s, respectively, in this simulation.
Transient analysis was used and the time step for solver was one eightieth of the time period of the shear wave. Uniform plane shear wave was produced by oscillating the line source with ten cycles of harmonic vibrations in the frequency range from 100 Hz to 400 Hz with a maximal displacement in the order of tens of micrometers. The displacements of the shear wave were recorded for post-processing at 8 locations, 1 mm apart, along a straight line that is normal to the line source. The phases of the wave were estimated by the Kalman filter and the average phase shifts were estimated using a linear fitting method [14]. The estimates of shear wave velocity and viscoelasticity are shown in Table 1.  Figure 11, except the elasticity µ1. Note that the differences between the average velocities and the reference velocities are less than 9% but the estimate error of µ1 is 15.5%. It is due to the fact that viscoelasticity moduli are proportional to the square of the phase velocity. Any small estimation errors of phase introduce large biases in the estimates of viscoelasticity, which is an intrinsic weakness of the ultrasound vibrometry, demonstrated by this example.

Experiment system and results
Experiments were conducted for evaluating ultrasound vibrometry. The diagram of an experiment system is shown in Figure 12. This system mainly consists of a transmitter to produce the ultrasound radiation force and a receiver unit using a SonixRP system. Two arbitrary signal generators were utilized to generate the system timing and excitation waveform. The waveform was amplified by a power amplifier having a gain of 50 dB to drive an excitation transducer for inducing vibrations in a tissue region. The SonixRP system was applied to detect the vibration using pulse-echo mode with a linear array probe. The SonixRP is a diagnostic ultrasound system packaged with an Ultrasound Research Interface (URI). It has some special research tools which allow users to perform flexible tasks such as low-level ultrasound beam sequencing and control. The center frequency of the excitation transducer was 1 MHz. The center frequency of the linear array probe was 5 MHz and the sampling frequency of SonixRP was 40 MHz. The excitation transducer and detection transducer were fixed on multi-degree adjustable brackets and were controlled by three-axis motion stages. The picture of experiment system setup is shown in Figure 13. The left lobe of a SD rat liver was embedded in gel phantom and placed in water tank. Before experiment, the SonixRP URI was run first to preview the internal structure of the liver. In the interface shown in Figure 14, the B-mode image and RF signal of a selected scan line were displayed together to help users selecting test points inside the liver tissue. The positions of the excitation transducer and the detection probe were adjusted to focus on two locations in the liver at the same vertical depth.  Computer programs based on the software development kit (SDK) of SonixRP were developed for detecting the vibrations and shear wave propagation. The programs defined a specific detection sequencing and timing that repeatedly transmit pulses to a single scan line and repeatedly receive the echoes with a PRF of 2 KHz. The timing of the excitation and detection pulses is shown in Figure 15. The pulse repetition frequency of the excitation pulses was 100 Hz. An example of the typical fast-time RF ultrasound signal acquired by the SonixRP is shown in Figure 16. Figure 16a shows the echo through the entire liver tissue region, while Figure  16b shows the echo around the focus point (75 mm in depth) in the liver tissue. The vibration of shear wave at a location was extracted from I and Q channels using the I/Q estimation algorithm described by equation (23). Figure 17a shows the vibration displacement and Figure 17b shows the spectral amplitude of the vibration. The extracted displacement signal sB(t,k) was processed by the Kalman filter [14] that simultaneously estimates phases of the fundamental frequency and all harmonics. Figure 18 shows estimated vibration phase shifts of the first four harmonics over a distance up to 4 mm. Linear regression was conducted to calculate the shear wave propagation speed for each frequency. Figure 18. Estimates of phase shifts over distances using vibration displacements and Kalman filter Figure 19 shows the phase velocities at different harmonics and the fitting curves of three models: Voigt, Maxwell, and Zener models. The fitting values are shown in Table 2. As shown by the figure and table, the Voigt model and Zener model fit the measurements of the phase velocity of the liver tissue better than the Maxwell model for this liver.  The second experiment was conducted to demonstrate the impact of boundary conditions. Because boundary conditions play very important roles in wave propagation, in vitro experiments were also conducted to investigate shear moduli of the superficial tissue of livers (0.4 mm below the capsule) and the deep tissue of livers (4.9 mm below the capsule). The excitation pulses were tone bursts having a center frequency of 3.37 MHz and a width of 200 μs for the binary excitation pulses and 100 μs or less for the OFUV excitation pulses. The pulse repetition frequency of the excitation pulses was 100 Hz. The broadband detection pulses had a center frequency of 7.5 MHz and pulse repetition frequency (PRF) of 4 kHz. Liver phantoms using fresh swine livers were carefully prepared so that the interface between the gelatin and the liver was flat. The thicknesses of liver samples were more than 2 cm and the areas were about 4×4 cm 2 . The phantom was immersed in a water tank.
The shear wave speeds were measured from 100 Hz to 800 Hz over a distance up to 5 mm away from the center of the radiation force application. Figure 20 shows the estimates of the shear wave speeds. Each error bar was the standard deviation of 30 estimates from five data sets of repeated measurements and six distances (1 to 4 mm, 1 to 5 mm, etc). The estimates from 100 Hz to 400 Hz were almost identical for the binary excitation pulses and the OFUV excitation pulses. Because the estimate errors using binary excitation pulses were too high for the frequency beyond 400 Hz, the estimates at 4.9 mm were based on the OFUV method. Figure 20 represents the trend of our experiment results that the shear wave speed in the superficial liver tissue is generally higher than that in the deep tissue. The results should be carefully examined. One of the possibilities is that we think it is caused by the liver capsule as we have verified it with Finite Element (FE) simulations, and another possibility is that the shear wave speeds of the gelatin are between 3 to 4 m/s from 100 to 800 Hz, higher than that in the liver tissue. The estimates of shear wave speeds at deep tissue of 4.9 mm and superfical tissue of 0.4 mm were used to numerically solve for the shear moduli of the three models. The curves generated by the models were compared with the measurements. As shown in Figures 21a and 21b, we find that the Voigt model may not always suitable for modeling liver shear viscoelasticity, at least for in-vitro applications with increased frequencies of shear waves in some of our studies.
On the other hand, we find that the Zener models matches the measurements very well with very small fitting errors as shown in the Figure 21 and Table 3. Table 3 shows the estimated shear moduli of different models with two different frequency ranges at two different depths in liver tissues based on our experiment data. Each modulus is an average of 30 estimates from 5 data sets and 6 distances. All elasticity has the unit of kPa and all viscosity has the unit of Pa·s. The fitting errors (m/s) are the deifferences between the measurements and calculated shear wave speeds using the models. The changes represent the variations of the estiamtes from one frequency range to another. The statistics are not conclusive because of the small number of samples. But this study indicates the variations of estimates and importance of the selection of tissue viscoelasticity models.   The third experiment was conducted to demonstrate the effectiveness of the ultrasound vibrometry to characterize the injury of liver tissue. Table 4 shows that the measured shear moduli of the livers thermally damaged by a microwave oven using different amount of cooking time (3, 6, 9, and 12 seconds). All estimates were from the superficial tissue region. It shows that the shear wave speeds estimated in the superficial tissue region are effective for indicating the damage levels of the livers. The errors are the standard deviations of the differences between the measurements and calculated speeds of the models.

Discussion
Shear moduli have very high dynamic ranges and are highly correlated with the pathological statues of human tissue. The solutions of the wave equation with constitutional models of tissue viscoelasticity show that the shear moduli of tissue can be estimated by measuring the phase velocity and attenuation of shear wave propagation in the tissue. However, it is a challenge to effectively and remotely generate vibrations and shear waves in a tissue region. It is also a challenge to measure shear wave because shear wave attenuates very fast as the propagation distance increases.
In the past fifteen years, the use of pulsed and focused ultrasound beams has been demonstrated as an effective method to remotely induce localized vibrations and shear waves in a tissue region. Several useful technologies have been developed for characterizing tissue viscoelasticity: Vibroacoutography, ARFI, Supersonic imaging, and ultrasound vibrometry, etc.
The ultrasound vibrometry is only technique that quantitatively estimates both tissue elasticity and viscosity. We found that the estimates of tissue elasticity by ignoring the viscosity are erroneous. Shear phase velocity are frequency dependent because the dispersive property of the biological tissue. Therefore, regardless of the usefulness of the viscosity, accurate estimates of tissue elasticity require the inclusion of the viscosity in the tissue models, as indicated by the solutions of the wave equation with three viscoelasticity models.
The ultrasound vibrometry transmits periodic push pulses to induce vibrations and shear waves in a tissue region, and detects the shear wave propagation using the pulse-echo ultrasound. The push pulses and detection pulses are interleaved so that one array transducer can be used for the applications of both pulses. The application of the array transducer allows the detection over a distance so that the phase velocities of several harmonics can be measured for calculating shear moduli.
Accurate estimates of shear moduli require an extended frequency range over an extended distance. The current technology is only effective for a few hundred Hertz in the frequency and a few mm in the distance away from the center of the radiation force applied. Shear wave having a high frequency attenuates very quickly as distance increases. Other vibration methods such as OFUV may be worth to explore.
We found that the shear wave speeds of livers are location dependent or dispersive in locations. Our experiment results indicate that the shear moduli estimated from a superficial tissue region and from a deep tissue region can be significantly different. Boundary conditions play a very important role in shear wave propagation and its phase velocity. The solution of the wave equation with boundary conditions should be considered for a tissue region that has a limited physical size. Some studies in this area have been done for myocardium and blood vessel walls.
The measurements of the ultrasound vibrometry are based on the assumption that tissue under the test is isotropic, which is not true for most tissues. Nevertheless, the measurements may be useful in clinical practices, which need to be evaluated in vivo experiments and clinical studies. On the other hand, the solutions of the wave equation with anisotropic tissue are needed.
Limited by the extensive contents in this chapter, we do not discuss the application of the Kalman filter in this work. The Kalman filter has great potential to include more complicated tissue models and motion models that are not fully explored yet, at least are not publically reported yet. On the other hand, Fourier transform and correlation method are also effective tools to calculate phases of the slow-time signals, if the motion model is simply sinusoidal.
Our experiments demonstrate that the ultrasound vibrometry can be readily implemented by using commercial medical ultrasound scanners with minimum alterations. Our experiment results also demonstrate that the ultrasound vibrometry is effective to characterize the stiffness and injury levels of livers.
We find that the Zener model fit the shear wave speeds of the livers better than the Voigt model and Maxwell model in almost all cases that include different frequency ranges, different locations, and different tissue conditions. Our study also indicates that the Voigt model is sensitive to the change of the observation frequency. Measurements at higher frequencies should be included when the Voigt model is used. In this case, the OFUV is useful to enhance the higher frequency components of the shear waves. The Zener model and Maxwell model appear to be less impacted by the frequency changes with our experiment data.

Conclusion
Tissue pathological statues are related to tissue shear moduli, which can be estimated by measuring the phase velocity of shear wave propagation in a tissue region. Ultrasound vibrometry is an effective tool to quantitatively measure tissue elasticity and viscosity. Ultrasound vibrometry induces vibrations in a tissue region using pulsed and focused ultrasound radiation force and detects the shear wave propagation using pulse-echo ultrasound. Experiment results demonstrate the effectiveness of the ultrasound vibrometry for characterizing tissue stiffness and liver damages.