2.2. Computer simulation

The modeling can be summarized as follows:

1. Each shot produces a pressure field p. The calculated pressure field depends on the physical and the tissue assumptions. This point will be addressed in section 2.3.

2. The energy of the beam is converted into heat. The tissue temperature evolution over time is obtained by solving Pennes’ bioheat transfer equation (BHTE) [7] which describes thermal processes in human perfused tissues. The BHTE has the following 3D+t form:

   $${\rho }_t{C}_t\frac{\partial T}{\partial t}=k_t{\nabla }^{2}T+V{\rho }_b{C}_b(T_a-T)+Q$$

   (1)

   where T is the temperature and ρ_t, C_t, k_t, V, ρ_b, and C_b are respectively the density, the specific heat and thermal conductivity of tissues, the perfusion rate per unit volume of tissues, the density and the specific heat of blood. This equation describes the variation of energy per unit of volume. In (1), it can be seen that this variation can be divided into three terms (from left to right): the thermal diffusion, the effect of perfusion and Q the acoustic power provided by the exposures. Q is directly deduced from the simulated pressure field (Q ∞ p²). A finite difference method was used to solve this equation [8]. The temperature of the cooling water serves as an iso-thermal boundary condition. As ebullition was not considered in the present model and for avoiding unrealistic temperatures; a threshold ensures that the temperature could never exceed 100 ° C. We take also the tissue necrosis into account: the perfusion rate V is set to zero when tissues are necrosed.

3. The concept of thermal dose [9], which is quite standard in tumor hyperthermia literature, was used to assess thermal damages. For each location and history of temperature, the thermal dose is defined as an equivalent time the temperature of reference of 43 ° C should be applied to obtain the same thermal damage. It can be computed directly from the temperature map:

   $$D_{43°\text{C}}(x,y,z,t)=\sum _{i=1}^t{R}^{(43-T(x,y,z,t))}\mathrm{\Delta }t$$

   (2)

   with R = 0.5 if T > 43 ° C and R = 0.25 otherwise.

4. Liver tissues are considered as irreversibly thermally damaged when D_43°C > 340 minutes [10].

2.3. Calculation of the pressure field

The determination of the pressure field is one major feature of the overall process. Numerical methods for estimating the pressure field can take into account the nonlinear propagation of high intensity waves [11, 12, 13] or remain purely linear. In our study, with a non-focused plane transducer, the propagation can be considered linear [14]. From this assumption, the pressure field can be exactly computed for the entire volume on the basis of the Rayleigh integral with O’Neil’s hypotheses [5]. The goal is thus to compute the acoustic pressure at a location point M; the transducer of surface S is sampled into elements of surface ΔS which size is negligible in comparison with the wavelength λ; M and one element of surface ΔS are connected by a straight segment ΔSM of length l sampled into elements Δl. The pressure at a location M is given by the discrete version of the Rayleigh integral:

with p₀ the pressure at the transducer surface (Pa), k the wave number (2π/λ), f the transducer frequency (Hz) and α_i, the tissue attenuation coefficient at location i along the line segment ΔSM (Np · cm⁻¹ · MHz⁻¹).

This twofold integration (over the surface of the transducer and the segment ΔSM) is very computation intensive. Some simplifications are usually made in order to accelerate the simulation [5]: the tissues are assumed to be homogeneous around the applicator and its physical properties (the tissue attenuation coefficient) constant in time. They have the three following benefits:

• The propagation of a wave from one element of surface ΔS to M is carried through only 2 media: the cooling water of the transducer (α₁, l₁) and the tissues (α₂, l₂)- Because α₁, is negligible and p₀ is estimated at the surface of the applicator, the equation (3) can therefore be simplified as:

  $$p(M)=\left|\sum _{S}j\frac{p_0}{\lambda }\mathrm{\Delta }S\frac{{exp}^{-\mathit{jkl}}}{l}{exp}^{-f{\alpha }_2{l}_2}\right|$$

  (4)

• The simulated air-backed transducer providing only forward waves and the homogeneous medium of propagation allow computing pressure in a quarter of a space in front of the propagation axis only for symmetry reasons.

• Regardless to the applicator orientation, the generated pressure field remains constant in the transducer coordinates system. So, the pressure field generated after the rotation of the applicator can simply be obtained by an angular shift of the first estimated pressure field matrix. The necrosis simulation can be summarized as: 1) estimation of the pressure field for one exposure; 2) temporal integration of this pressure field within the BHTE (1) after rotation by the planned angles.

Adaptative sampling of the transducer and volume of calculation allows a relatively fast pressure field computation time (8 minutes on a conventional PC: AMD Athlon 2200+ proc. with 1MB RAM).

2.4. Dynamic tissue attenuation coefficient variation

The attenuation coefficient associated to a given tissue is no more constant over time. This coefficient changes as temperature changes and as tissue necroses. The attenuation coefficient curves (ie coefficient variation vs. temperature) was measured for various tissues [15]. For liver, Connor et al. [13] made a 6th order polynomial approximation of the experimental data (see curve in Fig. 2).

Open in a separate window

Figure 2

Continuous and stepwise linearized tissue attenuation coefficient vs. temperature.

Taking this attenuation coefficient variation into account modifies the simulation framework described in section 2.2: 1) the pressure field is first computed and then incorporated into the BHTE; 2) according to the computed temperatures, the tissue attenuation coefficients at each location of the volume are updated; 3) a new pressure field is then estimated and so on.

The consequence of this time-varying behavior is that the hypothesis of homogeneity and time invariance of the tissue properties is no more true. The three previous described simplifications (analytical integration simplification of (4), symmetry and pressure field invariance) cannot be used anymore. The brute force algorithm based on (3) induces a tremendous computing time for the overall necrosis estimation. However, several simplifications can be proposed in order to reasonably reduce this computing time without loosing accuracy on the estimate of the necrosed volume:

• Within the BHTE, the temperature variation is relatively slow, so the tissue attenuation coefficient variation is slow too. Two different time steps can be considered: a fast one (Δt = 0.02s) for resolving the BHTE (1) and a longer one (2s) for updating the attenuation coefficients and the pressure field.

• The tissue attenuation coefficient curves is stepwise approximated (Fig 2). This allows to cluster the volume elements into larger regions with constant tissue coefficients. The segment ΔSM crosses now regions with a constant attenuation coefficient α_i over a distance l_i. Equation 3 can then be rewritten as:

  $$p(M)=\left|\sum _{S}j\frac{p_0}{\lambda }\mathrm{\Delta }S\frac{{exp}^{-\mathit{jkl}}}{l}{exp}^{-f{\sum }_{i=1}^N{\alpha }_i{l}_i}\right|$$

  (5)

  Where N denotes the number of crossed regions along the segment ΔSM.

• Segment coherence was also considered: adjacent ΔSM segments share almost the same attenuation properties. For nine adjacent connected transducer surface element ΔS, we compute the ${exp}^{-f{\sum }_{i=1}^N{\alpha }_i{l}_i}$ term of equation 5 for only the ΔSM segment which originates from the central surface element and apply this term to the segments which originate from the eight adjacent surface elements.

These simplifications, with the volume and transducer adaptive sampling mentioned in section 2.3, allow reducing the computation time to around 16h30. Some optimizations are still under study in order to keep diminishing the computation time.

This work is part of the national SUTI project supported by an ANR Grant (ANR05RNTS01106). The authors would like to thank J.-L. Coatrieux for his advices in conducting this research project and his contributions in writing this paper.