Modelling of endoluminal and interstitial ultrasound hyperthermia and thermal ablation: Applications for device design, feedback control and treatment planning

Abstract

Endoluminal and catheter-based ultrasound applicators are currently under development and are in clinical use for minimally invasive hyperthermia and thermal ablation of various tissue targets. Computational models play a critical role in device design and optimisation, assessment of therapeutic feasibility and safety, devising treatment monitoring and feedback control strategies, and performing patient-specific treatment planning with this technology. The critical aspects of theoretical modelling, applied specifically to endoluminal and interstitial ultrasound thermotherapy, are reviewed. Principles and practical techniques for modeling acoustic energy deposition, bioheat transfer, thermal tissue damage, and dynamic changes in the physical and physiological state of tissue are reviewed. The integration of these models and applications of simulation techniques in identification of device design parameters, development of real time feedback-control platforms, assessing the quality and safety of treatment delivery strategies, and optimisation of inverse treatment plans are presented.

• Endoluminal ultrasound
• hyperthermia
• interstitial ultrasound
• theoretical model
• thermal ablation
• treatment planning

Introduction

Endocavity, intraluminal, and interstitial ultrasound devices can deliver hyperthermia and thermal ablation from direct placement within the body [1,2]. The ultrasound energy source (transducer array) of these catheter-based applicators is positioned directly within or adjacent to a deep target volume, via placement within a body cavity, lumen, or by direct insertion. These approaches are preferable for sites where energy localisation from external devices is difficult or where localisation of all power and energy propagation within the target tissue is critical. In contrast to large fixed focus or phased arrays typical of external ultrasound technology, the dimensions of interstitial (1–10 mm diameter) and endocavity (20–30 mm diameter) ultrasonic transducer assemblies are constrained by limits imposed by anatomy, endoscopic and intraluminal access, and clinical utility. Typical operational frequencies between 3–12 MHz range translate to smaller wavelengths and higher attenuation/absorption, and result in energy penetration limited to ~3 cm, with the exception of endorectal high intensity focused ultrasound. Transducer assemblies typically consist of single or multiple element arrays of planar, curvilinear, or tubular transducer elements operated asynchronously with collimated beam output, although miniature phased arrays with electronic beam forming have also been investigated. Temperature regulated water-cooling and coupling of the applicator/tissue interface is commonly applied to protect intervening tissues and enhance penetration characteristics. By proper design of the applicator and procedures of use, these small devices demonstrate enhanced penetration with precise and predictable spatial control of power deposition. Independent real-time power control to transducer elements, frequency modulation, and/or rotation or translation of the transducer assembly can be used for dynamic control over the heating patterns, including magnetic resonance temperature imaging (MRTI) feedback control. Thus, in implementation, many of these ultrasound devices can selectively direct or conform the heating volume to a specified target area while protecting or avoiding other tissues. Some examples of designs and clinical implementation of this broad technology include: transurethral prostate ablation with MR guidance [3–6], endoluminal arrays for esophageal [7,8], cervix [9], and biliary tract tumour treatment [10,11], intravascular and cardiac ablation [12,13], and percutaneous/interstitial for ablation of soft tissues including prostate, liver, and brain [14–21].

A myriad of both complex and simplified acoustic and thermal models have been applied to the design of these type of devices, development of treatment strategies and control schemes, as well as applications to treatment planning. Theoretical models provide an invaluable tool during the design and optimisation of thermal therapy devices, providing an inexpensive and convenient method to isolate promising designs for fabrication and experimental evaluation. Models are also applied to aid the design and evaluation of feedback control strategies for thermal therapy procedures. Model-based feedback control can be integrated into transient models and used to explore strategies for modulating energy delivery during a procedure [3,22–25]. Treatment planning tools use predictions of the induced temperature and thermal dose or damage distributions to determine a set of treatment parameters that optimise the delivery of thermal therapy for individual patients [26]. In clinical scenarios where volumetric thermometry is not available, models may also be used as a 3D thermal dosimetry tool to augment the limited data available from invasive thermometry [20]. Varying requirements on the accuracy and computational efficiency of models for these applications have led to the development of models optimised for specific applications.

The objective of the following work is to provide a brief review of aspects of acoustic and biothermal modelling specific to endoluminal and catheter-based ultrasound. Models of ultrasound propagation and energy deposition in tissue, and associated relevant physical parameters, are described along with applicable approximations. Numerical techniques for modelling energy deposition and bioheat transfer, heat-induced changes to tissue properties, and calculating resultant temperature distributions are described, including methods for validating these models in tissue phantoms, ex vivo tissue, in vivo animal models and human clinical studies. Finally, this review concludes with a summary of current applications and research directions.

Modelling acoustic power deposition in tissue from endoluminal and interstitial devices

The time-averaged energy deposited in a material by an acoustic field is given by Equation 1: where α [Np/m/MHz] is the acoustic absorption coefficient, I [W/m²] is the acoustic intensity, p [Pa] is complex pressure profile, ρ [kg/m³] is density, and c [m/s] is sound velocity. Acoustic pressure calculations are incorporated in models and simulations to calculate the intensity of energy deposited in tissue. Computation of these pressure fields often involves numerical evaluation of the Rayleigh-Sommerfeld integral [27]. The surface of a transducer is assumed to comprise point-source radiators surrounded by a rigid baffle. The complex acoustic pressure at a point in space is then approximated as the summation of acoustic pressures along the entire surface of the transducer, given by Equation 2. The most straightforward strategy is to numerically evaluate Equation 2 by segmenting a transducer face into multiple point sources, spaced significantly less than a wavelength apart [28–30], and integrate over the entire surface for each field point. This is a computationally intensive approach for 3–12 MHz endoluminal and interstitial ultrasound devices that produce sub-millimetre wavelength pressure distributions with rapid spatial variations.

Several numerical techniques and approximations have been developed to expedite the computation of pressure profiles from practical ultrasound transducers. The rectangular radiator method approximates the pressure profile of a transducer as the summation of pressures radiated by several small rectangular sub-elements [31]. With this method, iterative calculations are needed similar to point source evaluations, but the number of calculations is substantially reduced because of the larger dimensions of the rectangular sources. This technique has been applied by several investigators for the design of interstitial and endoluminal ultrasound devices for thermal therapy [3,22,32–34]. To model near-field pressure distributions immediately adjacent to the transducer, requirements on the sub-element dimensions decrease and the number of sub-elements increases, resulting in computation times that vary quadratically with the number of sub-elements. Fast near-field methods have been developed by McGough et al. for circular, spherical and rectangular pistons [35,36]. These methods reduce two-dimensional surface integrals to a sum of one-dimensional integrals where computational burden varies linearly with number of elements. Hybrid angular spectrum (HAS) approaches are useful for fast calculation of pressure distributions in large 3D volumes [37]. These employ a combination of frequency and spatial domain methods to propagate a pressure wave through a defined medium, implemented efficiently following a fast Fourier transform (FFT). These approaches may be useful for modelling situations where propagation through complex interfaces is required.

To expedite pressure computations, some simplifying assumptions can be made by taking advantage of transducer geometry. Bessel functions can be used for circular or spherically focused transducers [38,39], and Fresnel approximations [40] can be used for rectangular apertures to compute an approximate Rayleigh integral using a quasi-analytic form. Several designs of interstitial or percutaneous applicators for ablation and hyperthermia utilise planar, tubular, or sectored tubular transducers arranged in a linear array. Assuming collimated output normal to the transducer across the entire transducer surface, the acoustic intensity in tissue can be approximated with analytical expressions incorporating exponential attenuation as shown in Equation 3a for a planar transducer [41,42], and inclusive of 1/r radial decay from the transducer surface for a tubular transducer as described in Equation 3b [43]. This approximation provides a computationally efficient method to simulate ultrasound power deposition from planar and tubular sources, and as demonstrated in Figure 1, this approximation agrees well with more detailed computations using the Rayleigh-Sommerfeld solution and experimental measurements. Multiple transducers on a single applicator are often driven asynchronously, and this power deposition from individual collimated elements can simply be summed [44–49]. Depending upon the transducer dimensions and acoustic output patterns, simplified expressions describing the shape of the intensity output normal across the surface of a transducer can be applied [42], including active sector angle directional tubular transducers [49].

Figure 1. A comparison between measured (hydrophone measurements) and simulated (Rayleigh-Sommerfeld and simplified models) acoustic intensity and pressure distributions is shown for (A) a sectored tubular transducer along the angular position, and a planar transducer along (B) the axial range and (C) longitudinal distance. (B) and (C) are reproduced from Lafon C et al. (Cylindrical thermal coagulation necrosis using an interstitial applicator with a plane ultrasonic transducer: in vitro and in vivo experiments versus computer simulations. Int J Hyperthermia 2000;16:508–22). The transducer tubes were sectioned longitudinally to have two active sectors (120°). Thermal simulations are shown for transurethral prostate hyperthermia/ablation by (D) planar and (C) tubular transducer computed using finite difference (FDTD) methods (recreated from Wootton JH et al. (Prostate thermal therapy with high intensity transurethral ultrasound: The impact of pelvic bone heating on treatment delivery. Int J Hyperthermia 2007;23:609–22)). Acoustic energy in panels (D) and (E) were calculated using the rectangular radiator method and the geometric approximation in Equation 2, respectively.

Many modelling studies account for ultrasound attenuation and ultrasound absorption by using a complex valued wavenumber while evaluating the Rayleigh-Sommerfeld integral. Tissue heterogeneity can also be incorporated in this formulation along with acoustic properties of the applicator housing, including catheters and coupling balloon materials. Acoustic refraction due to the curved oesophageal wall was modelled by Yin et al. [30] in their simulation study for endocavity cardiac ablation with an ultrasound phased array. Sinelnikov et al. devised an endovascular cardiac ablation device which consisted of parabolic reflectors for directing acoustic energy at the target tissue [50]. In a recent study investigating bone tumour ablation, Scott et al. [51] modelled the impact of including reflection and refraction phenomena at the bone-soft tissue interface during ablation with interstitial ultrasound applicators.

Typically, a linear relation is assumed between pressure and density variations due to wave propagation. Non-linear effects can be modelled using the KZK equation, however, these are generally not required for unfocused catheter-based ultrasound devices [52]. Acoustic cavitation and effects are typically neglected for the interstitial and endoluminal applicators operating at higher frequencies. During thermal ablation, elevated temperatures in excess of 90 °C may lead to tissue water vaporisation and boiling, which can shield ultrasound propagation causing shallow lesions with tight tissue coagulation [40]. Typically in practice and simulation, the applied power levels can be controlled to limit peak tissue temperatures below the thresholds for water vaporisation [48].

Thermal models

The acoustic energy deposition pattern due to a particular device is coupled with a bioheat transfer model to estimate the transient or steady-state temperature distributions. The Pennes bioheat transfer equation (BHTE) is the most widely used model, given by: where ρ is density, c is specific heat capacity, T is temperature, t is time, k is thermal conductivity, Q_US is the acoustic power deposition (given by Equation 1), m_bl is the blood mass perfusion rate, c_bl is the specific heat capacity of blood, T_bl is the temperature of blood, and Q_m is the metabolic heat production. The Q_m term is negligible for hyperthermia and thermal ablation applications where Q_US is several orders of magnitude larger. The Pennes model approximates heat transfer between small blood vessels and heated tissue with a volumetric blood perfusion term which acts as a heat sink [53]. Another widely used bioheat transfer model is the effective thermal conductivity model, which lumps together the effects of thermal conduction and blood perfusion into a single effective thermal conductivity term [54]. While neither the Pennes nor the effective thermal conductivity models account for heat transfer due to discrete vessels, they remain a reasonable approximation for modelling heat transfer in tissue and are widely used for evaluating performance of hyperthermia and ablation technology [3,21,24,41,55–61]. Detailed models that account for heat transfer due to thermally significant vessels and vessel trees have been developed [62,63].

Models of standard 30–60 min hyperthermia treatments often employ a steady-state solution of the bioheat transfer model (left hand side of Equation 4 is set to zero), since treatments are planned to hold a constant temperature profile and the time taken to reach the target temperature is small compared to duration of the treatment [26,44,64]. Transient models are, however, required if thermal dose is to be calculated. Further, steady-state models are not suitable for short duration applications such as thermal ablation, where the tissue temperature remains transient over the course of the heating period. Physical changes induced by heating (e.g. increases in tissue absorption due to coagulation, changes in blood perfusion) can be incorporated into transient models by modifying tissue parameters as a function of the time-temperature history, as discussed in later sections of this review. At elevated temperatures >90 °C where tissue water vaporisation may occur, the heat sink associated with the latent heat of water vaporisation can be incorporated to define an effective specific heat capacity [48,65].

Models of tissue thermal effects and damage

Thermal damage to tissue is a complex function of the time–temperature history during heating. Several models of thermal injury, effects, and destruction due to heating have been reported in the literature and have been reviewed by Dewhirst et al. [66], Dewey [67], and Pearce [68]. The Arrhenius thermal damage model and the thermal isoeffective dose model are commonly applied for describing effects of hyperthermia or thermal ablation in general, and commonly for ultrasound applicators. Thermal damage (Ω) can be calculated as a rate process (Equation 5), or related to the ratio of damaged to undamaged cells after heating for duration t (Equation 6): The parameters ΔE, activation energy, and A, frequency factor, are specific to tissue type, and need to be established experimentally. He et al. [69] have reviewed techniques for measuring the Arrhenius parameters, ΔE and A, using different indicators for assessing thermal damage.

Sapareto and Dewey [70] formulated a thermal isoeffective dose model, derived from the Arrhenius relationship as detailed in [69], where non-isothermal heating can be compared to isothermal heating at a reference temperature, which is typically 43 °C. The cumulative equivalent thermal dose of heating at 43 °C is given by Equation 7. For hyperthermia applications as an adjunct to radiation or chemotherapy, effective thermal doses greater than t₄₃ >  5–10 min are required [71–74]. For thermal ablation, the threshold for coagulative necrosis is within the range 100 < t₄₃ < 1000 min [67]. Numerous in vivo studies have shown a threshold of t₄₃ ≥ 240 to be in good agreement with the extent of coagulative necrosis after thermal ablation [67,75–77]. Typical values applied for safety thresholds to represent non-lethal exposure are T < 43 °C, t₄₃ < 5 min [66].

Prakash et al. conducted a theoretical study to compare the use of critical temperature, thermal isoeffective dose, and Arrhenius thermal damage thresholds for estimating the extents of necrosis after thermal ablation specific to catheter-based ultrasound applicators [48]. They found that for the 5–15 min ablations typical with catheter-based ultrasound applicators, the extent of the ablation zone as indicated by t₄₃ > 240 min and Ω > 4.6 boundary were all in good agreement. The T > 54 °C and T > 50 °C boundaries were in good agreement with the thermal dose and damage measures for short (~5 min) and long (~15 min) ablations, respectively. This study indicated that, provided treatment duration was taken into account, temperature-based thresholds may serve as a computationally cheaper surrogate for estimating thermal injury.

Thermal injury thresholds are also used to dynamically adjust physical properties of tissue in theoretical models of ablation. While experimental studies indicate that changes in tissue properties are not just a function of temperature, but the time–temperature history, models that include temperature-based adjustments are less computationally intensive and may be better suited to treatment planning applications.

Modelling physiological and physical effects of hyperthermia and thermal ablation

Hyperthermia and thermal ablation induce several physiological and physical effects on tissue that affect energy deposition and bioheat transfer. Comprehensive models of tissue heating with endoluminal and catheter-based ultrasound devices can account for these heating-induced changes in tissue parameters. Moderate levels of heating may modify tissue vasculature leading to an increase in blood perfusion within the hyperthermic temperature range. Sustained heating at elevated temperatures in excess of 50 °C leads to coagulation and microvascular stasis, and regions of tissue desiccation. At temperatures in excess of 90 °C, tissue desiccation and charring may be observed. While many of these changes are not relevant for modelling of hyperthermia, physical changes to tissue may dramatically alter energy deposition and heat transfer during thermal ablation. For modelling ablation with catheter-based ultrasound applicators, the most important physical parameters are the acoustic absorption and attenuation coefficients and blood perfusion rates.

Several studies have reported the impact of heating on acoustic absorption and attenuation coefficients of various tissues, in both in vivo and ex vivo tissue models [78–81]. Damianou et al. measured the acoustic attenuation and absorption coefficients of ex vivo canine kidney, liver, and muscle [79]. Acoustic attenuation increased at thermal dose values in excess of t₄₃ > 100–1000 min and T > 50 °C, reaching maximal values up to 2–3 times that of nominal values. They also noted a dependence of the increase in attenuation with the rate of heating. Similar increases in attenuation at ablative temperatures were measured in porcine liver and human prostate by other groups [78,80]. Changes in the attenuation coefficient are primarily attributed to structural changes due to protein coagulation and denaturation at elevated temperatures [78–81]. Since heating-induced attenuation changes below 50 °C are not significant, models of hyperthermia typically do not account for dynamic changes in tissue properties. While several studies have reported attenuation data in ex vivo tissues, there is limited data available for the temperature-dependent changes in attenuation for tissue in vivo.

Models of dynamic changes in perfusion as functions of temperature and thermal dose have been applied for simulating thermal ablation and hyperthermia. Temperature-based thresholds have been applied to model perfusion changes during hyperthermia based on experimental data [82,83]. The dynamic perfusion model of He et al. derived from experimental data and based upon an Arrhenius model, can be applied to estimate the degree of vascular stasis (DS) based on the transient temperature profile typical of ablation [84]. Simple maximum temperature and thermal dose thresholds of vascular stasis have also been employed [85].

Tyreus and Diederich developed a theoretical model of thermal ablation with interstitial ultrasound devices incorporating thermal dose-dependent changes in attenuation and microvascular perfusion [49]. This comprehensive model yielded close agreement between estimated extents of ablation zone boundaries and temperatures profiles as measured in vivo. The study by Prakash and Diederich employed five models of dynamic tissue changes of varying complexity to simulate feedback-controlled ablation with interstitial and transurethral ultrasound applicators [46]. Tissue models incorporating dynamic changes in blood perfusion and acoustic attenuation as a function of the time-temperature history predicted ablation zone volumes within 2.5% of each other. The proper modelling of acoustic absorption/attenuation is critical (see Figure 2); models using constant values for acoustic attenuation predicted ablation zone volumes up to 50% larger in prostate (low perfusion) and liver (high perfusion) tissue, respectively. Theoretical studies of prostate ablation with transurethral ultrasound applicators indicate that models including dynamic changes in attenuation predict longer treatment times and larger maximum temperatures, than those that do not include heating-induced increases in attenuation [32,33]. Mast et al. [40] and Tyreus and Diederich [49] showed close agreement between theoretical models and experimentally observed ablation zones in liver tissue when using models that accounted for heating-induced changes in attenuation and tissue water vaporisation.

Figure 2. (A) Temperature profile through central slice of a dual-sectored transurethral ultrasound as computed by a transient temperature model. Boundaries of the predicted ablation zone are shown when using a constant attenuation (dashed line) and thermal dose-dependent increase in attenuation. Models that assume a constant attenuation coefficient predict larger ablation zones and lower maximum temperatures. (B) MR image of ex vivo vertebra ablated with an interstitial ultrasound applicator and (C) temperature profile computed with theoretical modelling employing parameters similar to the experimental case. The theoretical model included dynamic changes in tissue attenuation and accounted for reflection and transmission of energy at the soft tissue–bone interface. The computed extent of region experiencing 20 °C rise (solid line) is in good agreement with experimental measurements (dashed line) obtained using MRTI. (B) and (C) are reproduced from Scott SJ, et al. (Interstitial ultrasound ablation of tumors within or adjacent to bone: Contributions of preferential heating at the bone surface. Proc SPIE, Energy-based treatment of tissue and assessment VII, 2013:85840Z–Z).

Numerical techniques

Numerical strategies are commonly employed for evaluating acoustic energy deposition and bioheat transfer models. For interstitial or endoluminal treatments, heating applicators are deployed within or adjacent to tumour targets. It is important to model tissue inhomogeneity arising from differences in properties of tumours, healthy tissue, organs, intervening structures such as lumen walls and bones. For certain organs such as liver, heart, lung and kidney it may be necessary to individually model large blood vessels, ducts and/or airways. Interstitial ultrasound applicators are acoustically coupled to the targets with flowing water. It is critical to include additional cooling in the models because of this water flow. Dynamic temperature-based changes in tissue properties and perfusion also play an important role in the models. Computation of 3D temperature distribution is a complex and computationally odious task. Several analytical and numerical techniques are available to accomplish these calculations. Some approaches are finite element methods (FEM), finite differences, and analytical methods based on Green's functions. We note that the application of finite differences and FEM techniques for evaluating the BHTE are not unique to modelling thermal therapy with catheter-based ultrasound applicators, and have been widely employed for modelling ablation and hyperthermia with other energy sources [86,87].

Finite difference methods have been used in several interstitial or endoluminal ultrasound modelling studies due to their ease of implementation [21–23,30,32,33,49,88–94]. With this technique, the partial derivatives of the BHTE (Equation 4) are discretised with finite differences. Explicit time-stepping methods are seldom employed due to the requirement for extremely small time steps to maintain stability. Rather, implicit schemes such as the Crank-Nicholson method and the alternating direction implicit (ADI) method are frequently employed to allow for arbitrarily large time steps, while maintaining stability. Finite difference (FDTD) methods for simulating heat transfer from endoluminal and interstitial sources have been implemented in Cartesian and cylindrical coordinate systems, with variable spatial discretisation. A well-known difficulty with finite difference methods is the efficient discretisation of complex geometries, such as those that may arise when constructing patient-specific models. Conformal corrections have been developed and employed for hyperthermia modelling to limit the staircasing errors that arise due to poor spatial discretisation of complex geometries with standard gridding techniques. More details on applying finite difference methods to thermal therapy simulations can be found in Chato et al. [95] and Neufeld et al. [96].

Finite element methods (FEM) are attractive tools in developing models of these catheter-based devices within heterogeneous or complex tissue environments, and for treatment planning applications. Tumour targets and anatomical structures can have complex irregularly shaped boundaries. In FEM implementation the simulation domain is divided into small constituent basis elements which often take the shape of triangles or tetrahedrons, allowing efficient high-resolution capture of anatomical shapes. FEM methods have been successfully employed to model bioheat transfer during hyperthermia and ablation with interstitial ultrasound [9,26,44,46–48]. In these studies, critical organs were segmented from radiological scans and then converted into FEM meshes using commercially available FEM software. The heat source for the BHTE was calculated using geometric approximations or rectangular radiator methods, and then interpolated on to the FEM grid. Most models of endoluminal and interstitial ultrasound thermal therapy have employed the Galerkin formulation of the FEM. Implicit transient solvers integrated within commercial FEM packages were utilised to compute 3D transient and steady-state temperature distributions. More details on the FEM technique as applied to bioheat transfer modelling can be found in Chato et al. [95] and Neufeld et al. [96].

While FEM and FDTD methods are most commonly used, other methods include FFT-based implementation and Green's function method. Dillenseger employed FFTs with incremental time steps while evaluating the Pennes equation for liver ablation with interstitial ultrasound [97]. The heterogeneous tissue and applicator structure, cooling interfaces and near-field heating typical of endoluminal and catheter-based ultrasound may make it difficult to accurately or efficiently model heating with these approaches.

Model validation

Comparison of experiment to theoretical predictions is essential for validating and defining the limits of the models, specifically as applied to a particular device. For catheter-based ultrasound devices, the primary experimental measurements of interest are the acoustic radiation pattern of transducers, transient temperature profiles induced in tissues or phantoms, and thermal necrosis and coagulation in tissue. Several studies have compared experimentally measured acoustic intensity profiles with numerically evaluated intensities computed by models of varying complexity (see Figure 1) [16,22,41,98,99]. Other studies have demonstrated good agreement between experimentally measured temperature profiles and coagulation patterns, with estimates from numerical models [41,49,57]. Further, modelling and ex vivo experiments were conducted by Prakash et al. to investigate clustered arrays of interstitial ultrasound applicators for liver ablation [48]. Dimensions of the observed ablation zones were statistically equivalent with thermal dose contours estimated by modelling (t₄₃ > 600 min). MRTI can be applied for 3D volumetric temperature measurements within phantoms or in vivo. For example, Burtnyk et al. compared 3D model-based temperature calculations with MR thermometry measurements during ablation of tissue mimicking gel phantoms with a transurethral ultrasound ablation device [33,91,99], demonstrating the models were able to predict temperatures within 1.6 °C across multiple MRTI slices (Figure 3). Scott et al. reported good agreement between 3D FEM calculations and MRTI data during heating of ex vivo bovine vertebrae surrounded with muscle [51]. Ultimately, thermometry available from in vivo clinical data can be used to validate theoretical models. As an example, temperature measurements obtained during a pilot study of interstitial ultrasound hyperthermia in prostate [20] were compared to 3D temperature profiles estimated using FEM-based patient-specific models [26,47]. For these treatments, the measured and simulated T₅₀ values in the hyperthermia target volume were in good agreement ranging between 40.1–43.9 °C and 40.3–44.9 °C, respectively, and transient changes in temperature correlated consistently between measurement and simulation (Figure 3).

Figure 3. Comparison of patient-specific 3D FEM thermal models with clinical thermometry data: (A) transient temperatures and (B) T₅₀. (C) Comparison between MRTI data and model-based temperature measurements is shown in representative axial slices during ablative heating in tissue mimicking gel phantoms. (C) is reproduced with permission from Burtnyk M, et al. (3D conformal MRI-controlled transurethral ultrasound prostate therapy: Validation of numerical simulations and demonstration 306 in tissue-mimicking gel phantoms. Phys Med Biol 2010;55:6817–39. DOI: 10.1088/0031-9155/55/22/014).

Application of models for device design, treatment planning and feedback control

Theoretical models are widely used to inform the design and performance of catheter-based ultrasound applicators and systems for controlled delivery of thermal therapy. With respect to endoluminal and interstitial applicators, theoretical models have been used to select parameters such as: power requirements, operating frequency, transducer geometry, number of transducers, applicator dimensions, coupling-fluid temperature and flow-rates, and catheter materials [9,21–23,45,58,85,90,94,100]. As part of this design process, the temperature distributions, thermal dose distributions, and estimated zones of thermal coagulation can be used for characterisation and optimisation. For systems employing multiple applicators, models have been used to compare the performance of applicator configurations [48,101]. The technique for calculating energy deposition and heat transfer should be selected based on the type of applicators/transducers, the level of desired accuracy, and computational resources available. For tubular and planar sources, the geometric approximations described offer a good balance between accuracy and computational efficiency. However, for lightly focused transducers or phased arrays, numerical techniques that numerically approximate the Rayleigh-Sommerfeld integral are required.

Treatment planning platforms employ theoretical models of energy deposition and heat transfer calculated on patient-specific anatomies to quantitatively assess and optimise treatment delivery strategies. Forward planning tools allow the user to adjust treatment parameters and observe the outcome on the treatment, and are used to iteratively compare the performance of treatment parameters such as applicator positioning and applied power levels. Inverse planning tools couple theoretical models with numerical optimisation techniques to determine a set of treatment parameters that optimise a defined treatment objective. Since inverse planning involves repeated model evaluation using different sets of treatment parameters, computationally efficient models are desirable to facilitate evaluation of treatment plans within clinically relevant timeframes. A common approximation is to assess treatment quality based on temperature profiles computed using steady-state solutions of the BHTE rather than thermal dose profiles.

Typical objective functions specify minimal/maximal acceptable values for the temperature inside the target, while specifying constraints on the maximal acceptable temperatures in critical structures, as well as other tissues. Chen et al. [26] developed a patient-specific treatment planning platform for catheter-based ultrasound hyperthermia of pelvic targets adjunct to HDR brachytherapy (Figure 4) where transducer powers were determined based on prescribed temperature distributions. The objective function was to maximise the volume of the target region above the therapeutic temperature, while restricting the maximal tissue temperature and the temperature in normal (i.e. non-targeted) tissue: The constrained optimisation problem was cast as an unconstrained optimisation problem through the introduction of penalty parameters, and was solved using a quasi-Newton approach (gradient-based optimisation). When using gradient-based techniques, which converge to local optima, a common approach is to perform several iterations of the optimisation procedure using various starting points in order to find the best local optimum [26]. Derivative-free optimisation techniques have been employed for optimisation of thermal therapy procedures; however, such techniques typically require a greater number of function evaluations to converge, without providing any guarantees of optimality [102]. To increase the computational efficiency of their treatment planning procedure, Chen et al. employed a coarse spatial discretisation with the BHTE during optimisation. The optimal solution was then evaluated with a fine spatial discretisation to verify its suitability. Similar approaches that vary the discretisation of the simulation space during optimisation have been applied to determine optimal applicator placement during RF ablation [103].

Figure 4. Patient-specific thermal models employed for optimisation and planning of hyperthermia treatments for cervical (top row) and prostate tumours where heating applicators with multiple transducer elements were deployed by endoluminal (cervix) and interstitial (prostate) placement. Case geometries and FEM meshes were constructed from CT scans which were segmented to delineate critical organs, tumour targets and applicator positions. Reproduced from Chen X et al. (Optimisation-based thermal treatment planning for catheter-based ultrasound hyperthermia. Int J Hyperthermia 2010;26:39–55).

Salgaonkar et al. [47] developed a temperature superposition technique for rapidly computing the temperature profile induced by catheter-based ultrasound applicators. This method used precomputed solutions of the temperature profile induced by an individual ultrasound transducer to estimate the temperature profile induced by an array of applicators employing multiple transducers. The superposition technique yielded a 4–7× computational speed-up factor, while predicting temperatures within 0.4 °C of those computed by FEM models, thereby enabling rapid computation of optimal treatment plans that can be readily integrated into the clinical workflow. With this technique it was possible to generate optimised treatment plans within 10–15 min, which is practical for routine clinical planning. When choosing a method for evaluating the objective function for inverse planning, the associated computation time must be taken into consideration. The amount of time afforded for treatment planning, as well as the number of objective function evaluations required for convergence to an optimal solution, place constraints on the choice of the computational model. Further improvements in computational speed may be achieved with hardware acceleration and the use of parallel architectures [96].

The choice of tissue parameters employed in a simulation has a large impact on the outcome. Since tissue physical properties for individual patients are seldom available a priori, models used for treatment planning should not be expected to provide accurate predictions of the temperature profiles induced in tissue. Rather, they are valuable for comparative evaluation of several candidate treatment configurations. Wootton et al. assessed the impact of uncertainties in tissue parameters by comparing treatment plans optimised with various sets of tissue parameters [44]. There is a need to develop stochastic simulation techniques that provide statistical measures of treatment outcome.

Catheter-based ultrasound systems incorporate closed-loop feedback to dynamically control the energy deposition profile through power and frequency modulation, and rotation/translation of the transducer assembly. Theoretical models provide a powerful tool for determining optimal control parameters and comparative evaluation of control algorithms. Burtnyk et al. [33,91,92,99] developed a patient-specific FDTD modelling platform to simulate feedback-controlled transurethral prostate ablation with an applicator consisting of a linear array of planar ultrasound transducers. An algorithm to adjust rotation rate and applied power levels based on MRTI feedback was designed, and controller parameters determined using the theoretical model [33,56]. Simulations on patient-specific anatomies were used to assess the safety profile of the transurethral ablation system and optimise strategies for treating prostate glands of varying sizes [91]. These results were used to guide delivery of transurethral prostate ablation in a clinical trial [4].

Patient-specific models have also been applied as a computational dosimetry tool for estimating temperature and thermal dose profiles delivered during hyperthermia. In scenarios where MRTI is not available, catheter-based temperature sensors provide only limited information about the volumetric temperature profile. A recent study of hyperthermia with interstitial ultrasound devices compared temperature profiles computed by theoretical models with those measured by multi-junction thermocouples, and used this tool to estimate the physical parameters of the targeted tissue [20]. These patient-specific parameters were then applied to theoretical models to estimate the 3D temperature and thermal dose profiles to characterise the quality of hyperthermia treatments. Similar techniques have been employed with other hyperthermia energy modalities for correlating clinical response with model-derived estimates of temperature profiles.

Conclusion

Theoretical models of thermal therapy with endoluminal and interstitial ultrasound applicators are widely used to assess the feasibility of thermal therapy with devices in specific targets, optimise the design of devices, evaluate feedback control strategies, and for patient-specific treatment planning. This review detailed techniques for modelling the ultrasound energy deposition, bioheat transfer, and thermal damage profiles after thermal therapy with catheter-based ultrasound devices. Comprehensive models account for the dynamic changes in tissue properties induced during heating. For time-sensitive applications such as patient-specific treatment planning, computationally efficient techniques are desirable to enable rapid computation for integration into the clinical workflow.

Declaration of interest

We gratefully acknowledge support from US National Institutes of Health grants R01CA122276, R44CA112852, and P01CA159992. The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.