Shear Wave Propagation in Soft Tissue and Ultrasound Vibrometry

1. 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.

In the past ten years, Oestreicher’s method was utilized to quantitatively measure both tissue elasticity and viscosity. Ultrasound vibrometry has been developed to noninvasively and quantitatively measure tissue shear moduli [10-16]. It induces shear waves using ultrasound radiation force [5, 6] and estimates the shear moduli using shear wave phase velocities at several frequencies by measuring the phase shifts of the propagating shear wave over a short distance using pulse echo ultrasound [10-16]. Applications of the ultrasound vibrometry were conducted for viscoelasticities of liver [16], bovine and porcine striated muscles [17, 18], blood vessels [12, 19-21], and hearts [22]. A recent in vivo liver study shows that the ultrasound vibrometry can be implemented on a clinical ultrasound scanner of using an array transducer [23].

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-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-21], and liver [8], when ultrasound vibrometry is used.

2. 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 c_s(ω) of the wave can be estimated by measuring the phase shift Δϕ over a distance Δz:

$$c_s(\omega )=\omega \Delta z/\Delta \varphi$$E1

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:

$$\frac{d^{2}S}{d{z}^2}+k^{2}S=0$$E2

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

$$k=k_r-i{k}_i$$E3

The solution of (2) is a standard solution of a homogeneous wave equation:

$$S=\widehat{x}{S}_0{e}^{-ikz}$$E4

where S₀ is the displacement at z = 0, $\widehat{x}$ is an unit vector in x direction. The plane wave is independent in y direction. The real time time-harmonic shear wave is:

$$S(\omega ,z,t)=\widehat{x}\mathrm{Re}\left\{S{e}^{i\omega t}\right\}=\widehat{x}{S}_0{e}^{-k_{i}z}\cos (\omega t-k_{r}z)$$E5

Although attenuation coefficient α = –k_i 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 $$d(\omega t-k_{r}z)/dt=0$$:

$$c_s(\omega )=\frac{dz}{dt}=\omega /k_r$$E6

The complex wave number k of the plane shear wave is a function of the frequency and the complex modulus of the medium [9]:

$$k=\sqrt{\rho {\omega }^2/\mu }$$E7

where ρ is the density of the tissue and the complex modulus that connects stress σ and strain ε:

$$\mu =\sigma /\epsilon ={\mu }_1+i\omega {\mu }_2$$E8

which describes the relationship between stress and strain in the Voigt tissue. The Voigt model consists of an elastic spring μ₁ and a viscous damper μ₂ connected in parallel, which represents the same strain in each component as shown in Figure 1.

Figure 1.

Voigt model consists of an elastic spring μ₁ and a viscous damper μ₂ connected in parallel.

The relation between stress σ and strain ε of the Maxwell tissue is:

$$\sigma ={\mu }_1\epsilon +{\mu }_2\frac{d\epsilon }{dt}$$E9

For a harmonic motion, (9) becomes:

$$\sigma =({\mu }_1+i\omega {\mu }_2)\epsilon$$E10

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):

$$c_s(\omega )=\sqrt{\frac{2({\mu }_1^2+{\omega }^2{\mu }_2^2)}{\rho ({\mu }_1+\sqrt{{\mu }_1^2+{\omega }^2{\mu }_2^2}}}$$E11

The elasticity μ₁ and viscosity μ₂ 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 c_s(ω) increases at the rate of square root of the frequency and there is no the upper limit for c_s(ω). 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 μ₁ and μ₂. In tissue characterization applications, μ₁ is often in the order of a few thousands and μ₂ is often less than 10. Thus, when the wave frequency is very low (less than a few Hz),

$${\mu }_1\approx \rho c_s^2(\omega )\forall \text{ very low }\omega .$$E12

When the frequency is very high (higher than a few tens of kHz),

$${\mu }_2\approx \rho c_s^2(\omega )/2\omega \forall \text{ very high }\omega .$$E13

Figure 2.

Plot of phase velocity of shear wave having μ₁=3 kpa and μ₂=1 pa.s in Voigt tissue

A broad frequency range is needed to accurately estimate both μ₁ and μ₂. (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.

Figure 3.

Maxwell model consists of a viscous damper η and an elastic spring E connected in series.

The relation between stress σ and strain ε of the Maxwell tissue is:

$$\frac{1}{E}\frac{d\sigma }{dt}+\frac{\sigma }{\eta }=\frac{d\epsilon }{dt}$$E14

For a harmonic motion, (14) becomes:

$$\frac{\sigma }{\epsilon }=\frac{i\omega \eta E}{E+i\omega \eta }=\mu =\frac{{\omega }^2{\eta }^{2}E}{E^2+{\omega }^2{\eta }^2}+i\frac{\omega \eta E^2}{E^2+{\omega }^2{\eta }^2}$$E15

which is the complex shear modulus of the Maxwell model. Unlike the Voigt model, real and imaginary components of (15) are functions of the frequency. When the frequency is fixed, the complex modulus is a function of $\eta$ and E. Substituting (15) into (7), the shear wave speed in Maxwell medium can be found from (6):

$$c_s(\omega )=\sqrt{\frac{2E}{\rho (1+\sqrt{1+E^2/{\omega }^2{\eta }^2}}}$$E16

Equation (16) can be also obtained by replacing μ₁ and μ₂ 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 c_s(ω) 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 η. $$E\approx \rho C_s^2(\omega )$$for a very large ω, $$\eta \approx \rho C_s^2(\omega )/2\omega$$ for a very small ω, c_s(ω) is zero for ω=0, and c_s(ω) approaches $$\sqrt{E/\rho }$$ when ω is very high.

Figure 4.

Plot of phase velocity of shear wave having E = 7.5 kpa and η = 6 pa.s in Voigt tissue

Figure 5.

Zener model adds an elastic spring E₁ to the Maxwell model (η, E₂) in parallel.

The Zener model adds an additional elastic spring, having the elasticity of E_1, to the Maxwell model (η, E₂) 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:

$$\mu =E_1+\frac{i\omega \eta E_2}{E_2+i\omega \eta }=E_1+\frac{{\omega }^2{\eta }^2{E}_2}{E_2^2+{\omega }^2{\eta }^2}+i\frac{\omega \eta E_2^2}{E_2^2+{\omega }^2{\eta }^2}$$E17

Substituting (17) into (7), the shear wave speed in Zener medium can be found from (6):

$$c_s(\omega )=\sqrt{\frac{2({\omega }^2{\eta }^2{(E_1+E_2)}^2+E_1^2{E}_2^2)}{\rho ({\omega }^2{\eta }^2(E_1+E_2)+E_1{E}_2^2+\sqrt{({\omega }^2{\eta }^2{(E_1+E_2)}^2+E_1^2{E}_2^2)({\omega }^2{\eta }^2+E_2^2)}}}$$E18

Equation (18) shows that $$E_1+E_2\approx \rho C_s^2(\omega )$$ for a very large ω, $$E_1\approx \rho C_s^2(\omega )$$ for a very small ω, η is proportional to the slop of the speed curve, and c_s(ω) approaches$$\sqrt{\left(E_1+E_2\right)/\rho }$$ when ω is very high. A numerical example of phase velocity of Zener tissue is shown in Figure 6.

Figure 6.

Plot of phase velocity of shear wave having E₁ = 4.5 kpa, η = 1.5 pa.s, and E₂ =7.5 ka in Zener tissue

3. 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-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,

$$d(t)=D\sin ({\omega }_{s}t+{\varphi }_s)$$E19

where ω_s=2πf_s 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]:

$$r(t,k)=\left|g(t,k)\right|\cos \left({\omega }_{0}t+{\varphi }_0+\beta \sin ({\omega }_s(t+kT)+{\varphi }_s)\right)$$E20

where T is the period of the push pulses shown in Figure 9 and the modulation index is:

$$\beta =2D{\omega }_0\cos (\theta )/c$$E21

where c is the sound propagation speed in the tissue, ω₀ is the angular modulation frequency of detection tone bursts, g(t,k) is the complex envelope of r(t,k), ϕ₀ 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]:

$$g\left(t,k\right)=I\left(t,k\right)+jQ(t,k)$$E22

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]:

$$s_A(t,k)={\tan }^{-1}(Q/I)-\text{mean of }{\tan }^{-1}(Q/I)=\beta \sin \left({\omega }_s(t+kT)+{\varphi }_s\right)$$E23

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.

Figure 7.

Block diagram of quadrature demodulator

Another motion detection method [14] uses a complex vector that is a multiplication between two successive complex envelops [29]

$$X\left(t,k\right)+\text{ i}Y\left(t,k\right)=g\left(t,k+1\right)g^{*}\left(t,k\right)$$E24

Thus, the motion velocity in slow time can be obtained,

$$s_B(t,k)={\tan }^{-1}\left(\frac{Y(t,k)}{X(t,k)}\right)=-2\beta \sin ({\omega }_{s}T/2)\cos \left({\omega }_s(t+kT)+{\varphi }_s+{\omega }_{s}T/2\right)$$E25

which is proportional to the velocity of the tissue harmonic motion for ω_sT/2 << 1. Thus, sin(ω_sT/2) ≈ ω_sT/2 and the velocity amplitude is $${\omega }_{s}T\beta$$, which is also $$D{\omega }_s=\beta {\omega }_{s}c/2{\omega }_0\cos (\theta )$$ because of the derivative relation between (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:

$$\Delta \mathrm{\varphi }={\mathrm{ }\mathrm{\varphi }}_{\mathrm{s}}\left(\mathrm{z}\right)-{\mathrm{ }\mathrm{\varphi }}_{\mathrm{s}}(\mathrm{z}+\mathrm{\Delta }\mathrm{z})$$E26

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 = -}k_rz. Given a tissue location (axial location) in fast time, all constant phases are removed by (26) except the phase shift $\Delta \mathrm{\varphi }$ in the lateral location.

Ultrasound vibrometry is developed to induce the shear wave described by (19) and detect the phase shift $\Delta \mathrm{\varphi }$ 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 f_v and harmonics nf_v. 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 f_v of the push pulses is in the order of 100 Hz, and pulse repetition frequency f_PRF of the detection pulses is in the order of 2 kHz.

Figure 8.

Array transducer for transmitting ultrasound radiation force and detecting shear wave propagation

Figure 9.

Interleaved push pulses for ultrasound radiation force and detection pulses

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].

4. 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.

Figure 10.

Schematic diagram of simulated tissue region (domain) and

All domains had the same material property of the Voigt tissue and all boundaries were set free to avoid reflections. The material parameters were: density of 1055 kg/m³, Poisson’s ratio of 0.499, and Voigt rheological model of the viscoelasticity model. The Voigt model was converted and represented in the form of Prony series. The store modulus and loss modulus were calculated using frequency response analysis for demonstrating the conversion of the Prony series. The complex shear modulus of the Voigt model is the same as (8):

$$\mathrm{\mu }\left(\omega \right)=\mathrm{ }{\mu }_1+i\omega {\mathrm{\mu }}_2$$E27

where elasticity modulus µ₁ and viscosity modulus µ₂ 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.

	Shear Wave Velocity (m/s)				Viscoelasitcity Estimation
100Hz	200Hz	300Hz	400Hz	µ1(kPa)	µ2(Pa*s)
Reference value	1.5574	1.9372	2.3470	2.7362	2	2
Measurement 1	1.46238	1.91972	2.37439	2.66872	1.69	1.90
Measurement 2	1.50648	1.94791	2.42833	2.86637	1.63	2.10
Measurement 3	1.52748	1.92955	2.46275	2.81089	1.74	2.10
Average	1.49479	1.9216	2.44359	2.77296	1.69±0.056	2.03±0.11
Std	0.02828	0.02455	0.05672	0.08517

Table 1.

Estimated Viscoelasticity of Voigt tissue having µ₁ = 2 kPa and µ₂ = 2 Pa*s

The shear wave velocities in red represent the theoretical values of wave speeds in Voigt tissue. The estimates of the speeds and viscoelasticity moduli of three simulations are shown by three sets of the measurement. Their average values are close to the theoretical values as shown in Figure 11, except the elasticity µ₁. Note that the differences between the average velocities and the reference velocities are less than 9% but the estimate error of µ₁ 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.

Figure 11.

Estimated shear phase velocities and set reference values

5. 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.

Figure 12.

Block diagram of the experiment system

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.

Figure 13.

Experiment setup with SonixRP.

Figure 14.

Ultrasound Research Interface (URI) of SonixRP

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.

Figure 15.

Timing sequence of the experiment system

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.

Figure 16.

Ultrasound RF echo (a) through the entire liver and (b) around the focus point 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.

Figure 17.

Displacements of the vibration and its frequency spectrum

The extracted displacement signal s_B(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.

Figure 19.

Curve fittings of three models with the estimates of the phase velocities of the liver tissue

Voigt Model, μ1, μ2, fitting error	4.10 kPa	1.51 Pa·s	0.019 m/s
Maxwell Model, E, η, fitting error	7.18 kPa	4.27 Pa·s	0.143 m/s
Zener Model, E1, E2, η, fitting error	4.07 kPa, 45.9 kPa	1.47 Pa·s	0.020 m/s

Table 2.

Estimated viscoelasticity moduli of three models

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². 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.

Figure 20.

Shear wave speeds in superficial and deep liver tissues

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.

Figure 21.

Model fittings for shear wave in the deep tissues (a), and superficial tissue (b)

Depth=0.4 mm	100 to 400 Hz	100 to 800 Hz	Changes
Voigt, μ1, μ2, fitting error	2.48, 2.00, 0.152	3.71, 1.46, 0.204	50%, 27%
Maxwell, E, η, fitting error	10.7, 2.50, 0.043	11.7, 2.36, 0.048	10%, 6%
Zener, E1, E2, η, fitting error	0.578, 9.033, 2.85, 0.028	1.34, 9.843, 2.56, 0.0569	132%, 9%, 10%
Depth=4.9 mm
Voigt, μ1, μ2, fitting error	2.74, 1.35, 0.108	3.59, 0.791, 0.151	31%, 41%
Maxwell, E, η, fitting error	5.68, 2.82, 0.016	5.90, 2.70, 0.021	4%, 4%
Zener, E1, E2, η, fitting error	1.49, 4.20, 2.44, 0.015	1.70, 4.25, 2.19, 0.018	14%, 1%, 10%

Table 3.

Estimates of Shear Moduli (elasticity in kPa, viscosity in Pa s)

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. The Zener model provides the best curve fitting with the minimum fitting error.

		3 sec.	6 sec.	9 sec.	12 sec.
Voigt Model	μ1	9.23	9.67	11.2	13.0
μ2	1.60	1.72	2.54	3.01
Error, (m/s)	0.103	0.114	0.114	0.121
Maxwell Model	E	18.3	19.6	32.3	39.2
η	3.60	3.73	3.93	4.48
Error, (m/s)	0.117	0.173	0.172	0.231
Zener Model	E1	7.68	8.40	9.60	11.9
E2	15.0	18.0	35.0	63.3
η	1.90	1.91	2.81	3.10
Error, (m/s)	0.029	0.034	0.0344	0.102

Table 4.

Shear Moduli of thermally damaged livers

6. 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.

7. 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.

Acknowledgement

This research was supported in part through a grant from National Institute of Health (NIH) of USA with a grant number of EB002167 and a grant from Natural Science Foundation of China (NSFC) with a grant number of 61031003.