The Westervelt equation with a causal propagation operator coupled to the bioheat equation.

The Westervelt wave equation is frequently used to describe non-linear propagation of finite amplitude sound. If one assumes that the medium can be treated as a thermoviscous fluid, a loss mechanism can be incorporated. In this as in previous work the authors replaced the typical loss mechanism incorporated in the Westervelt equation with a causal Time Domain Propagation Factor (TDPF) which incorporates the full dispersive effects (both frequency dependent phase velocity and attenuation) in the numerical solution while remaining in the time-domain. In the present work we investigate heat deposition due to finite amplitude propagation through a dispersive medium (e.g., human tissue). To this end, the Westervelt equation with and without the TDPF is coupled to the Pennes bioheat equation and the coupled equations are solved using the method of finite differences to determine the resulting heat deposition. We show that non-linear effects are large and that proper treatment of dispersion results in significant changes as compared to modeling the medium as a thermoviscous fluid.


1.
Introduction. Modeling acoustic propagation in tissue presents a number of mathematical, physical, and numerical issues. To begin with, the propagation of sound in such a dispersive medium must be accurately described, and since the problem is inherently non-linear, finite-amplitude effects must be taken into account. Further, in many applications it will be of interest to determine the heat deposition, i.e., the process of heat accumulation or buildup, due to absorption of acoustic energy in the biological medium. This deposition traditionally is modeled using the bioheat equation [9] or some related heat diffusion equation. While this problem is intrinsically interesting, it has implications for ultrasound imaging as well.
In previous work [8,11] we have studied the propagation of sound in a dispersive medium by solving the Westervelt equation in the time domain, treating the dispersion by modeling the medium as a thermoviscous fluid, or, alternatively, by employing the time domain propagation factor (TDPF) [14,15,16,17]. It was found that employing the causal TDPF yields results that differ significantly from employing the Westervelt equation with the traditional loss mechanism, which is due to the fact that the latter fails to take into account the full dispersive characteristics of the medium. In this work we couple the Westervelt equation (with the viscous loss term as well as with the TDPF) to the Pennes bioheat (diffusion type) equation [9] in order to study the heating of the medium due to acoustic energy.

GUY V. NORTON AND ROBERT D. PURRINGTON
Solovchuk and Sheu [13] recently employed the Westervelt equation in modelling the heating of blood and tissue due to absorption of highly focused acoustic energy in ultrasound ablation of liver tumors. No attempt was made to incorporate dispersive effects and the treatment of the frequency dependence of the absorption is very different from the present work. Earlier modelling of non-linear acoustic propagation by Jing and Cleveland [5] utilized the KZK equation and did not examine heating effects.
Human tissue is inherently dispersive and when acoustic energy is propagated through tissue, causality requires that the propagating field be attenuated, and that the attenuation will be frequency-dependent. If the source is time limited, i.e. an acoustic pulse, its frequency content will be appropriately broad, and to accurately treat such a situation the numerical model will have to be able to take into account this frequency-dependent attenuation and phase velocity. Further, if the spectrum of the acoustic source is in the ultrasound region and of sufficiently high intensity then finite amplitude effects can lead to the production of non-linearly generated harmonics. Traditionally, the Westervelt equation has been used to model acoustic propagation in inhomogeneous media such as tissue, treating the medium as a thermoviscous fluid. As mentioned previously, Norton and Purrington [8,11] have replaced the traditional (thermoviscous) loss mechanism in the Westervelt equation with a causal TDPF that takes into account the full dispersive nature of the medium, (i.e. both frequency dependent velocity and attenuation).
In non-linear propagation, acoustic energy will be introduced into higher harmonics, a fact that has been exploited in recent years to improve medical imaging. The technique commonly known as harmonic imaging [18] has the benefit that when imaging at the higher harmonics one does not need to exclude or time-gate out the energy from the source. A further aspect of finite-amplitude (non-linear) ultrasound propagation in a dispersive medium involves heat deposition. While heating will always take place if the medium is lossy, i.e. dispersive, in the case of non-linear propagation, heat deposition is modified because acoustic energy is directed into the higher harmonics.
As indicated above, in this work the deposition of heat will be described by a diffusion type equation coupled to the Westervelt equation. In earlier work, Hallaj and Cleveland [4] have coupled the bioheat equation to the Westervelt equation containing a thermoviscous loss term. But if dispersion is not treated properly it is quite clear that the description of heat deposition will be inadequate. In this paper we will compare the deposition of heat as developed using the traditional Westervelt equation with that obtained using the causal TDPF. Since the acoustic field acts as the source term for the bioheat equation, differences in the acoustic field will lead to measurable differences in the deposition of heat. Specifically we show that the deposition of heat based on the acoustic field developed by the causal TDPF can yield results that differ significantly from employing the Westervelt equation with the traditional loss mechanism.
We begin by briefly reviewing the Westervelt equation containing the traditional loss term and then replace this term with the causal TDPF. Next we describe the bioheat equation and address how it is coupled to Westervelt equation. This is followed by a description of the finite-difference formalism and implementation of both equations. Next a numerical experiment is carried out in one dimension, using both forms of the Westervelt equation, each coupled to the bioheat equation.
Finally the differences and similarities observed with regard to the deposition of heat associated with the two approaches are discussed.

2.
Westervelt equation for non-dispersive and dispersive media. The development below follows from Purrington and Norton [11]. The propagation of finite amplitude sound in a lossless medium can be described by the Westervelt equation [6]: where p is the acoustic pressure, ρ 0 and c 0 are the ambient density and sound speed, β = (1 + B/2A) is the non-linearity coefficient for the fluid, and B/A is the non-linearity parameter. The third term in the equation describes non-linear distortion of the wave due to finite-amplitude effects. If the medium is assumed to be a thermoviscous fluid, 1 takes the following form: The added (third derivative) term is a loss term, which is due to the thermal conduction and the viscosity of the fluid. Here δ is the diffusivity of sound; see 13 below. In a thermoviscous fluid subject to a harmonic signal, the classical absorption coefficient α 0 is related to δ and the angular frequency ω = 2πf by where α 0 is a constant under classical Stokes-Kirchhoff theory. Note the 1/ω 2 dependence. The final form of 2, hereafter referred to as the "Viscous equation," takes the following form with the attenuation explicitly shown, If the attenuation in the medium follows a frequency power law, as in Sabo [14], then 4 generalizes to the following form [8,11,14,15,16,17]: where Γ (t) represents the causal TDPF and the * denotes a convolution. We will call 5 the "TDPF equation." The causal TDPF introduces both the causal attenuation and variable phase velocity and is defined as the inverse Fourier transform of the complex propagation factor γ (ω), The complex propagation factor is written as, where α(ω) is the frequency dependent attenuation, β is the frequency dependent wavenumber (hence dispersive component), and τ is the retarded time τ = t − r/c 0 . The attenuation coefficient and wavenumber are causally related through the Kramers-Krönig relations [7]. The connection is expressed in the following fashion: where H represents the Hilbert transform, thus leading to the following expression for the TDPF: here 1 + (τ ) represents the step function defined as [2] 3. Bioheat equation. The medium of interest in this experiment will be mammalian tissue. The heating of tissue due to the attenuation of finite amplitude acoustic waves has often been described by the bioheat equation, originally formulated by Pennes [9]. The bioheat equation for tissue at temperature T , takes the form where ρ t , C t , and κ are the density, specific heat, and thermal conductivity of the tissue, respectively. This is just the classic heat conduction equation with two terms added: Q 1 is the perfusion term, representing heat transfer as a result of blood flow through capillaries, and Q 2 incorporates the contributions from metabolism and the absorption of acoustic energy. In Q 2 , which is the term that couples to the acoustic propagation equation, we neglect the contributions from metabolism. The thermal source term Q 1 is given by Newton's law of cooling/heating, namely, where T b and C b are the ambient temperature and specific heat of the blood, and W b is the volumetric perfusion rate (cooling by blood flow). The resulting equation is, assuming κ does not vary spatially, is Pierce [10] has shown that the absorption term Q 2 , representing the attenuation of acoustic energy by the tissue, is expressible as where c t is the ambient speed of sound in the tissue and the attenuation is given by α = ω 2 δ/c 3 t , as in 3. The diffusivity of sound, δ is shown by Pierce [10] to have the form: where µ is the viscosity, C p is the specific heat capacity of blood at constant pressure, and γ is the usual ratio of specific heats. The final form of the Bioheat equation is then: 4. Finite difference formulation. A one-dimensional numerical solution of both the Viscous and TDPF equations were performed using a finite-difference time domain (FDTD) representation that was explicitly fourth order in both time and space. The second partial derivative of the pressure with respect to time was extended to fourth order, but because the usual fourth-order finite difference representation of the second partial derivative would lead to an unconditionally unstable scheme, a technique given by Cohen [3] based on the "modified equation approach" to obtain fourth-order accuracy in time was used. This technique, while improving the accuracy in time, preserves the simplicity of the second-order accurate time-step scheme. This was done in order that all derivatives have the same order of accuracy as the boundary conditions, which are also fourth order. The reader is directed to Ref. [3] for further information and explanation. The temporal derivatives were calculated as follows: The first derivative term is associated with the non-linearity and the third derivative term is the loss term. The spatial differences were computed using a fourth-order accurate centered differencing scheme For the Bioheat equation the temporal and spatial derivatives were calculated as follows: To solve the Bioheat equation and thus to determine the heating of the tissue, the first time derivative of the acoustic pressure in 15 are required at each grid point and for four previous time steps 15. The complementary operator method (COM) was utilized in the development of the absorbing boundary conditions (ABC). The COM method, which was introduced in electromagnetics, where it was shown to yield excellent results [12], is a differential equation-based ABC method, and differs from the other common approach of terminating the grid with the use of an absorbing material [1]. The details of the derivation of the COM operators are beyond the scope of this paper and the interested reader should refer to Ref. [12].
5. Numerical experiment. The numerical experiment is carried out in one dimension. The input signal is a 1MHz sinusoidal burst of six cycles modulated by a Gaussian envelope in time. The numerical grid consists of 1500 points with ∆X = 0.1mm. The time increment is ∆t = 10.681 ns. The source is located at X = 3.5 cm. The medium simulated is mammalian liver tissue in which the sound speed in the non-dispersive case is 1600 m s −1 and the density of the tissue is taken to be 1100 kg m −3 . The physical constants required by the Bioheat equation are listed in Table 1, [19]. The attenuation is assumed to have the following power law frequency dependence: α (ω) = α 0 ω y = 2.0976 × 10 −7 Np m −1 rad −1 s −1 ω 1.13 .
(21) Figure 1 displays the attenuation in dB m −1 versus frequency in Hz. The causal phase velocity that is associated with this attenuation is shown in Fig. 2. The truncated source, which is a sinusoid of constant frequency but finite time duration, has a significant frequency content, as shown in Fig. 3. The source has a band width of approximately 1 MHz with a center frequency of 1 MHz. The use of the TDPF in 5 introduces the frequency-dependent attenuation and phase velocity. When solving the TDPF equation a maximum of 1000 points where used for the convolution, and in the Viscous equation, the attenuation appropriate to the peak frequency (1MHz) is used, that is, 87.5 dB m −1 or approximately 10.1 Np m 1 (see Fig. 1).

Figure 1. Attenuation versus frequency.
In order to compare the effect that the two loss terms have on the propagating signal and hence the heating of the tissue, the effects of non-linearity are removed either by setting the initial pressure p 0 = 1 Pa or by excluding the non-linear term altogether. In order to determine whether the implementation of the Viscous and TDPF equations were properly performed, both the attenuation and phase velocity were extracted from the numerical signal recorded at two locations and then compared to the analytic values. The reader is directed to Ref. [8] for the details. Figure 4 compares the analytic attenuation, 21 to the attenuation extracted from the Viscous and TDPF equations. The extracted attenuation from the TDPF equation is essentially indistinguishable from the analytic attenuation indicating that the time domain propagation factor is introducing the correct amount of loss at each frequency. The attenuation as a function of frequency based on the Viscous equation crosses the analytic curve at the center frequency of 1 MHz as expected. At frequencies below the peak frequency, the viscous term produces less attenuation than the causal operator, with the situation reversed above the peak frequency.  There is a frequency squared dependence on attenuation as one would expect (see 3 ), and indeed for any medium the viscous equation imposes a frequency squared dependence on the attenuation.
The phase velocity associated with the Viscous equation is frequency independent, and in the present case has a value of 1600 m s −1 , in contrast to the frequencydependent phase velocity (Fig. 2)     order to determine computationally whether the difference in the treatment of attenuation makes an important difference in the observed heating of the tissue, and to sort out the importance of non-linear effects, the following relevant cases were examined: 1) the Viscous equation, with and without the non-linear term, and 2) treatment of attenuation using the TDPF equation, again with and without nonlinearity. Figure 5 shows a simulation in which the initial pressure p 0 was chosen to be 15 MPa. The temperature response of the tissue obtained from using both the Viscous and TDPF equations coupled to the Bioheat equation, along with the linear and non-linear responses are shown. Note that the initial temperature deposition is insensitive to the presence of non-linear effects, and that for the linear case, both the Viscous equation and the TDPF equation show nearly identical results. This is clearly because in the absence of the non-linear term, no acoustic energy is introduced into the higher harmonics. And since the attenuation has the same value at the peak frequency (1 MHz) in all cases (see Fig. 4), the energy is dissipated in a nearly identical fashion. The gradual decay of the acoustic signal naturally results in a gradual decrease in the heat deposition. Comparing the Viscous and TDPF (non-linear) equations results, the different treatment of attenuation (Viscous and TDPF) translates into a large difference in the heating of the tissue. While the initial heat deposition at the source location (X = 0.035 m) is the same for both cases and is identical to the linear result, the TDPF equation results in more heating and peaking later (further down range) when compared to the Viscous equation results, and the heating in the two cases has a completely different range dependency compared to the linear results. If we consider the effect of the different treatment of attenuation, we note that for the Viscous equation, the peak deposition occurs at a distance of 0.0126 m from the source (R1 = 0.0476 m) while when employing the TDPF equation, the peak deposition is delayed, occurring at a distance of 0.0177 m (R2 = 0.0527 m). Furthermore, the heat deposition resulting from the TDPF equation is larger than that from the Viscous equation. As we have noted, in the case of finite-amplitude (non-linear) propagation, acoustic energy is introduced to higher harmonics. The attenuation of this more complicated distribution of acoustic energy has a dramatic effect on heat deposition. As an aid in understanding this, we determine the frequency spectrum of the acoustic signal at the two peak heat deposition locations (R1 (Fig. 6) and R2 (Fig. 7)) as a function of the treatment of attenuation (lossless, TDPF, and Viscous). Note that at the first location, R1 (Fig. 6), the amplitude associated with the higher harmonics decreases with increasing frequency in all cases, something that is not true of the lossless signal at location R2 (Fig. 7), where the amplitude is greatest at 3 MHz. It is evident that less energy is deposited in the higher harmonics using the Viscous equation than using the TDPF equation, and both are below the lossless case. In order to study the dependence of heat deposition on the initial pressure p 0 , the numerical experiment was also carried out for a lower pressure of 10 MPa. Figure 8 compares the resulting heat deposition versus range with the earlier case of an initial pressure of 15 MPa. As expected the overall heat deposition or temperature change is reduced. Although the shape of the temperature change vs. range curve using the Viscous equation is quite different, the location of the peak heat deposition is essentially unchanged ( 0.0124 m or X = 0.0474 m). Using the TDPF equation, however, while the heat deposition curve is still strongly peaked, the change is significant, with the peak occurring later (.0215 m or X = 0.0565 m). In other words, correct treatment of attenuation significantly changes the nature of heat deposition in the tissue. We now compare the peak heat deposition locations using the Viscous and TDPF equations to the approximate location of shock formation for the two source pressures to see if there is a correlation. In the non-linear case, the approximate location of shock formation   can be determined from the following expression [10]:  It follows that for an initial pressure p 0 = 15 MPa the location of shock formation occurs at X = 0.0437 m and for 10 MPa it is X = 0.048 m, that is, further from the source. At the higher peak pressure the onset of shock occurs before the peak heat deposition (for both the Viscous and TDPF equations), while at 10 MPa this depends on how the attenuation is modeled. For the TDPF equation the peak of heat deposition occurs further from the source than shock formation for both peak pressures, while using the Viscous equation the situation is more complicated, with maximum heat deposition occurring at essentially the same location, independent of source pressure [See Table 2]. Thus there is little correlation between shock formation and peak heat deformation for the Viscous equation. It is, however, worth noting that for the TDPF equation the distance between peak heat deformation and shock formation is nearly independent of the initial pressure p 0 . In either case, Viscous or TDPF, the effects of non-linearity are dramatic, both in the magnitude of the effect and in the functional dependence of the heating on range when compared to the linear case. This is clearly due to energy being fed into the higher harmonics of the original signal. And the energy in these higher frequencies is attenuated more strongly, resulting in additional heating. Clearly there is a complicated relationship between the effects of dispersion (varying phase velocity and attenuation) and those of non-linearity, on heat deposition. 6. Concluding remarks. We have previously shown that implementing the causal propagation operator in a modified Westervelt equation gives a more faithful description of acoustic propagation in dispersive media. Therefore, we can safely assume that the calculation of the heating effects in tissue will also be given more accurately using this approach, whether non-linear effects are included or not. The fundamental question is whether non-linear effects will swamp any difference due to the more accurate treatment of dispersion. This obviously will depend on the amplitude of the pressure wave, but we expect that because the non-linearity introduces harmonics of the original signal, that dispersion increases with frequency, and that attenuation necessarily accompanies dispersion, that the components with higher phase velocity will be attenuated more than those with lower phase velocity. Since it is attenuation that produces the heating of the tissue, we expect that heating effects will peak earlier in the non-linear case, and that as the pressure increases, maximum heating will occur earlier. This is what is observed.