|Home | About | Journals | Submit | Contact Us | Français|
Shear wave propagation techniques have been introduced for measuring the viscoelastic material properties of tissue, but assessing the accuracy of these measurements is difficult for in vivo measurements in tissue. We propose using the Kramers-Kronig relationships to assess the consistency and quality of the measurements of shear wave attenuation and phase velocity. In ex vivo skeletal muscle we measured the wave attenuation at different frequencies, and then applied finite bandwidth Kramers-Kronig equations to predict the phase velocities. We compared these predictions with the measured phase velocities and assessed the mean square error (MSE) as a quality factor. An algorithm was derived for computing a quality factor using the Kramers-Kronig relationships.
Wave propagation methods have been used in the past two decades to measure the elasticity of soft tissue. Techniques that use shear, surface, and Lamb waves have been developed to noninvasively measure the viscoelastic properties of tissue (Yamakoshi et al., 1990; Muthupillai et al., 1995; Sarvazyan et al., 1998; Chen et al., 2004; Palmeri et al., 2008; Catheline et al., 1999; Kanai, 2005; Zhang and Greenleaf, 2007). The commonalities between all of these techniques are that some form of mechanical excitation is used to create low frequency waves, and the resulting propagating waves are measured using ultrasound or magnetic resonance (MR) methods. The motion or propagation velocity of the waves are used in conjunction with a model to extract viscoelastic material properties.
Particularly in measurements made in vivo, there is no way to assess the accuracy of the estimate of the viscoelastic tissue parameters. A quality factor for evaluating the accuracy of a set of measurements would be highly desirable, especially in clinical situations.
In the product FibroScan® produced by Echosens (Paris, France) two quality criteria are used to evaluate the success of a measurement (Kettaneh et al., 2007). The system's software assesses the quality of the vibration induced by the vibrator and the quality of the ultrasonic measurements of the motion. When a measurement or “shot” does not meet the given criteria, no value is reported. Another quality standard has been reported which involves repetition of the measurement ten times and reporting the median, in addition to requiring that the success rate of obtaining measurements is greater than 60%, and also the interquartile range of the successful measurements was required to be less than 30% of the median value reported (Lucidarme et al., 2009).
In a method proposed by Palmeri, et al., shear waves are induced using short tonebursts of ultrasound. The shear wave is tracked at positions laterally displaced from the ultrasound beam's axis. The time at which the peak of the shear wave occurs at different lateral positions is recorded using a method called Lateral Time-to-Peak (TTP) (Palmeri et al., 2008), and then a linear regression on the time of the maximum displacement of the shear wave versus measurement distance is performed to find the shear wave velocity and the shear modulus is calculated. This regression result is rejected if the goodness of fit metrics R2 > 0.8 and 95% confidence interval < 0.2 are not met.
Deffieux, et al., reported measuring shear waves of different frequencies and calculated the standard error of the slope of the linear fit for phase data over different lateral locations (Deffieux et al., 2009). When the slope changed significantly due to inclusion of points that did not lie along the same line the standard error increased and that point was excluded. This was performed for data at different frequencies, and at higher frequencies the cutoff distance was reduced because the shear wave had attenuated and the phase data became unreliable for velocity estimation.
The three quality metrics employed for the aforementioned techniques are empirically chosen and are not directly derived from any fundamental physical phenomenon. In most cases, rejection of a data set or measurement would occur if gross motion, noise, or other artifacts interfered with the measurement.
It is desirable to begin with the fundamental physical phenomena of wave propagation and derive a quality factor related to these relationships. Soft tissue is viscoelastic; therefore the modulus describing the tissue is modeled as a complex quantity. This infers that the modulus is complex and dispersive, or varies with frequency. If it is assumed that the measurements are causal and linear, then the phase velocity and attenuation of propagating waves can be related using Kramers-Kronig relationships (Waters et al., 2005). These relationships dictate that if the dispersive attenuation can be measured, then the phase velocity can be predicted and vice versa. This is useful because one independently measured parameter can be used to check another.
For the specific case of evaluating the viscoelastic parameters of soft tissue, we focus on the use of ultrasonic methods to induce and measure motion. In our method called Shearwave Dispersion Ultrasound Vibrometry (SDUV) (Chen et al., 2009), focused ultrasound produces a radiation force that causes the tissue to deform. Pulse-echo ultrasound combined with signal processing techniques are used to extract the induced motion. Using this motion information, phase information is extracted and used to estimate phase velocities at multiple frequencies. The dispersive wave velocity is fit to a model such as the Voigt model to find the shear elasticity and viscosity (Chen et al., 2004; Chen et al., 2009). This model has been found to fit previously measured data very well over the selected frequency range.
The purpose of this paper is to report a quality metric that is derived from the fundamental relationship between the complex parts of the shear wave propagation defined by the Kramers-Kronig relationship.
Szabo and Wu developed a model for wave propagation in viscoelastic media that is based on the absorption obeying a power law (Szabo and Wu, 2000). They detail how the Kramers-Kronig relationships can be utilized with measured attenuation to determine the phase velocity of the propagating waves. In the most general case, the propagating wave can be described by
where ω is the angular frequency, γ is a complex wave number, the M subscript describes the mode whether it is longitudinal, M = L, or shear, M = S, and r is propagation distance. For the purposes of this study we focus on the shear mode only, so M = S, and the subscript will not be carried any further. The complex wave number is
where α is the attenuation, and
The absorption function is modeled as a power law (Szabo and Wu, 2000)
The measured attenuation is fit to the model in (4). The original Kramers-Kronig relationships assume that data are available at all frequencies. However, measured data are restricted to finite bandwidths. To account for this, differential forms of the Kramers-Kronig equations have been developed (Szabo and Wu, 2000; Waters et al., 2003). If the exponent y in (4) is an even integer or a noninteger,
where c(ω0) is the measured phase velocity at a reference frequency, ω0. If y is an odd integer,
Using different reference frequencies will bias the predictions up or down. We will describe later how to select ω0. These relationships can be used with measured attenuation data to predict the phase velocity and compare the agreement between the data and the prediction.
The theory presented above has been used to compare longitudinal wave attenuation and phase velocity at frequencies above 1.0 MHz (Szabo and Wu, 2000; Waters et al., 2003). Szabo and Wu evaluated shear wave data in plastics at frequencies above 2.0 MHz. The frequency range that is useful for evaluation of shear waves in soft tissue is typically from about 20-1000 Hz. In SDUV, the Voigt model is used to describe the viscoelastic parameters of the tissue. The shear wave attenuation and velocity is modeled as (Catheline et al., 2004)
where ρ is the tissue density, μ1 is the shear elasticity and μ2 is the shear viscosity.
The attenuation was simulated for a frequency range of 50-600 Hz, a tissue density of 1000 kg/m3, shear elasticity of 2, 5, and 10 kPa, and shear viscosity of 0.5, 2, and 4 Pa·s. The frequency resolution of the simulations was 50 Hz. Figure 1 shows the simulated attenuation values. The open symbols represent the simulated values, and the solid lines are fits to (4) for each set of simulated points.
The exponents for the curve fits in most cases were nonintegers so (5) was used for prediction of the shear wave velocity. The shear wave velocity was simulated using (8) and plotted in Figure 2. The reference frequency in (5) was varied and the velocity prediction was compared with the simulated values and the mean square error (MSE) was calculated
where N is the number of frequencies evaluated, ĉs (ωn : ω0) is the predicted velocity at the nth frequency for reference frequency ω0, and s (ωn : ω0) is the measured value at the nth frequency for reference frequency ω0, and the units of the MSE are m2/s2. The reference frequency that gave the smallest MSE was deemed as the optimal frequency. The phase velocity prediction corresponding to this reference frequency was plotted for each combination of μ1 and μ2 in Figure 2. For most cases, the predictions matched the simulated values very well. The results of the curve fitting of simulated shear wave attenuation values, the reference frequency and the associated MSE are summarized in Table 1. The wave attenuation coefficients listed in Table 1 will give values of attenuation of Np/cm when used in (4).
A flow diagram for the algorithm introduced above is shown in Figure 3. We take radiofrequency (RF) data and perform motion detection. We use the displacement amplitude to estimate the shear wave attenuation at a set of frequencies. For the same set of frequencies, the displacement phase is used to estimate the shear wave velocity. The attenuation data are fit to (4) and the parameters are used to predict the shear wave velocity. The predicted and estimated shear wave velocities are compared, and the MSE is calculated. We also introduced an avenue to optimize data selection with an iterative loop to identify outliers. Lastly, we calculate the viscoelastic properties of the tissue using both the data and the predictions and report the results.
An experiment was performed to test how measured data could be used with the developed method. A section of porcine tissue was embedded in a gelatin phantom made with 300 Bloom gelatin (Sigma-Aldrich, St. Louis, MO) with a concentration of 10% by volume and glycerol (Sigma-Aldrich, St. Louis, MO) with a concentration of 10% by volume. The porcine tissue was chosen because of the clear and consistent orientation of the muscle fibers so that it was easy to determine the shear wave polarization with respect to the fiber direction. Ultrasound radiation force was produced using 3.0 MHz ultrasound tonebursts of length 200 μs repeated at a rate of 100 Hz. This produces a multifrequency radiation force with components at harmonics of this repetition frequency. Tissue motion was measured using pulse-echo ultrasound at 7.5 MHz (ECHO Ultrasound, Lewistown, PA). Both transducers were co-focused before the experiment. The excitation transducer was moved to change the propagation distance while the pulse-echo transducer remained fixed. The first excitation point was placed 2.0 mm from the measurement point to avoid any interference in the motion tracking from the toneburst that produces the radiation force. A block diagram of the experimental setup is shown in Figure 4. The echo data were processed with previously mentioned methods described by Hasegawa and Kanai (Hasegawa and Kanai, 2006) and Zheng, et al (Zheng et al., 2007). Values of μ1 and μ2 were calculated from measured shear wave propagation velocity in the muscle tissue. Measurements were made along and across the muscle fibers.
The wave velocity was estimated by measuring the phase at two different locations and computing the velocity using cs(ω) = ω·Δr/Δ, where Δr is the distance between measurement locations and Δ is the phase difference. For this experiment Δr = 4 mm. The excitation consists of repeated tonebursts which produces shear waves at the repetition rate of 100 Hz and harmonics of that repetition rate. Motion was tracked for frequencies 100-600 Hz in 100 Hz increments.
To find the attenuation, the normalized displacement amplitude data for each frequency was fit to an exponential function of the form y = a exp(−br) to find the attenuation coefficient b. However, before this attenuation estimation, a correction must be made for the geometric attenuation associated with the cylindrical wave produced by the radiation force. Chen, et al., modeled the amplitude of the displacement of the radiation force induced displacement in the axial or z-direction, uz(r), as (Chen et al., 2004)
where is the zeroth order Hankel function of the first kind, and . At distances where hr 0, the relationship in (10) is approximated by
Therefore, to correct for the geometric attenuation, the displacement values were multiplied by and then normalized for fitting to the decaying exponential function. The attenuation fits for measurements along the muscle fibers for frequencies 100-600 Hz are shown in Figure 5. The attenuation versus frequency is shown in Figure 6 along with the fit to (4) and the fit is αs (ω) = 0.00012ω1.18 + 0.0045 (Np/cm). The attenuation curve in Figure 6 was then used to predict the shear wave velocity along the muscle fibers. The set of velocity predictions versus reference frequencies are shown in Figure 7.
Figures 8, ,9,9, ,1010 show the attenuation estimation curve fits, the attenuation dispersion and fit to (4), and the shear wave velocity predictions for measurements made across the muscle fibers, respectively. The attenuation fit for propagation across the muscle fibers was αs (ω) = 0.0019ω0.85 − 0.44 (Np/cm). The MSE values for the shear wave velocity predictions for measurements made along and across the muscle fibers for different frequencies are given in Table 2. The reference frequencies with the lowest MSE were chosen for reporting. In this experiment, 100 Hz and 600 Hz were chosen as the reference frequencies for the measurements made along and across the fibers, respectively. Figure 11 shows the viscoelastic material property estimates for measurements of propagation along and across muscle fibers using the shear wave velocity in (8). We performed the fit to the raw velocity data values and compared them to the fit to the Voigt model to the estimated shear wave velocity predictions made using the Kramers-Kronig relationships. The values of μ1 and μ2 are very consistent between the two methods. The viscoelastic properties along the muscle fibers evaluated from the data and the predictions were μ1 = 12.65 kPa and μ2 = 2.91 Pa·s and μ1 = 12.84 kPa and μ2 = 2.73 Pa·s, respectively. The viscoelastic properties across the muscle fibers evaluated from the data and the predictions were μ1 = 5.32 kPa and μ2 = 1.05 Pa·s and μ1 = 5.12 kPa and μ2 = 0.98 Pa·s, respectively.
In the simulation results with the Voigt model it was found that using the wave attenuation to predict the wave velocity worked well for the ranges used for μ1 and μ2. The fits to the attenuation values from (7) are good. When (5) is used to fit to the calculated values from (8), the fits are poorer for combinations of high μ1 and μ2.
When experimental data were used, the measured normalized displacement amplitude was used to assess the wave attenuation at each frequency. Those attenuation values were then used to fit the model in (4) and the fit coefficients were then used to predict the wave velocity. Varying the reference frequency served to bias the wave velocity predictions. This is predicted by the first terms of (5) and (6). For the most part, the overall trends of the predicted wave velocity agreed with the data as evidenced by the general agreement among the MSE values for different reference frequencies. If the bias associated with using different reference frequencies was significant, the MSE values would vary more widely.
One of the more important points is that independent measurements of propagating shear waves displacement amplitude and phase provide data to check the quality of the data as a whole, and provide a guide as to which points should be potentially excluded from analysis.
These results provide evidence that the Kramers-Kronig relationships can be used to assess the quality of shear wave data in soft tissue. In this paper, the MSE served as a quality factor. The MSE did not vary drastically with using different reference frequencies and it was much lower in the measurements made across the muscle fibers.
A more sophisticated quality factor for the data could be formulated and could include multiple steps to insure the integrity of the final viscoelastic property estimates. We could incorporate some of the quality control steps outlined in the introduction for Fibroscan®, the Lateral TTP method, or the method described by Deffieux, et al. An advantage of this quality factor is that it would apply empirical measures as well as the fundamental Kramers-Kronig equations that relate the wave attenuation and velocity.
There are a few parameters that could influence the quality factor derived from the Kramers-Kronig relationships. Sufficient motion amplitude is necessary for adequate estimation of the attenuation at each frequency. The goodness-of-fit for the estimates of attenuation at each frequency could be used to assess the quality of the data. This may affect which frequencies are used to fit the model in (4) because at higher frequencies, the displacement will decrease. We may choose a certain bandwidth where motion is adequate and estimates are more likely to be accurate.
The reference frequency is also an important parameter. The reference frequency did not change the MSE significantly for the data presented in this study, but this parameter should not be treated as a trivial point. When the wave velocity values are estimated using different reference frequencies, the variance of the different estimates could be assessed and MSE values of the predictions compared to the data could be computed and used. Data outliers which increase the MSE significantly could be separated based on comparisons with the wave velocity and excluded from analysis.
We could fit the wave velocity predictions and the data for μ1 and μ2 and compare the results, and use the material property estimates if the values are within some tolerance. All these parameters could be aggregated in some way to formulate an overall quality factor that could qualify the confidence in the results which could then be evaluated by a user to accept or reject the measurement. Thresholds for the different parameters would be application specific and established with preliminary tests.
The ability to assess data quality of propagating waves is not specific to any particular excitation or measurement method. As long as the wave is propagating and can be generally modeled using (1), then this method could be employed to formulate a quality factor. Therefore, surface, longitudinal, Lamb, etc., waves could be used in this format. The method is also not limited to any rheological model used to find the viscoelastic material properties. The method or model used to evaluate the material properties should not be critical, but has been shown to work well with the Voigt model for the bandwidth and ranges of μ1 and μ2 used.
The method is also not restricted to using a power law frequency approximation for the attenuation. A different model could be formulated and the fundamental Kramers-Kronig relationships could be used to compare the two parts of the complex wave number.
Some biases may be present that are associated with measurement of wave propagation in a viscoelastic, anisotropic medium such as the skeletal muscle used in this study. As far as anisotropy is considered, we limited ourselves to a transversely isotropic material for analysis. The particular muscle used had very well delineated fiber structure so we were more confident in trying to isolate the attenuation to two cases, along and across the muscle fibers. The transducer that we are using to create the shear waves will also introduce complexities associated with diffraction from the radiation force source. This has been addressed previously (Bercoff et al., 2004b; Chen et al., 2004). Regarding the bias due to spatial source, we know that measurements made in the excitation region or very close (< 2 mm) at times can be biased. For this reason, we usually start our measurements at 1-2 mm from the center of the excitation beam. Lastly, there may be some error due to the rheological model used, in this case the Voigt model. The data seem to be well described by the model for the tissue examined in these experiments. Different models could be applied, but as the Voigt model is casual and linear, it should be suitable for the application of the Kramers-Kronig relationships.
The Kramers-Kronig relations do assume that information for the full bandwidth is known; however, in practice this will not be possible. A way to optimally choose the reference frequency used for predictions with finite bandwidth measurements may need to be devised. This has been addressed in studies performed for ultrasound attenuation and phase velocity measurements made in the megahertz range and may be adapted for use in shear wave propagation (Mobley et al., 2003; Waters et al., 2003; Waters et al., 2005). We will devote future theoretical and experimental studies to understanding the effects of the reference frequency has on the agreement between predicted and measured wave velocity values.
We presented results for a one-dimensional shear wave propagation experiment. This method could be applicable for other elastography applications in liver, muscle, and other types of tissue that utilize shear wave propagation (Sandrin et al., 2003; Palmeri et al., 2008; Gennisson et al., 2003). The only requirement is that shear wave attenuation and velocity data could be obtained over an adequate bandwidth for the Kramers-Kronig relationships to be valid and applicable. There should not be any fundamental barriers for extending this method to multiple dimensions such as supersonic shear imaging (Bercoff et al., 2004a), as long as dispersive measurements necessary for applying the Kramers-Kronig relationships can be made.
We have presented a quality factor based on the Kramers-Kronig relationships for low-frequency shear wave propagation in soft tissue. We used the shear wave attenuation to predict the shear wave velocity at frequencies in a finite bandwidth. In an ex vivo experiment with skeletal muscle, we found good agreement between the predicted shear wave velocity and the measured values. The mean square error was reported as a quality factor. The fundamental Kramers-Kronig relationships provide a unique basis for assessing the quality of shear wave data used to characterize the viscoelastic material properties of soft tissue.
The authors are grateful to Dr. Shigao Chen and Randall Kinnick for experimental support, Tom Kinter for computer support, and Jennifer Milliken for administrative support. This work was funded in part by grants EB002167 and EB002640 from the National Institutes of Health. Dr. James Greenleaf has financial interest in some of the methods described.