The role of mathematical modelling in thermal medicine

Introduction

The use of thermotherapies in treating diseases has its origin in antiquity. There are historical accounts in ancient cultures of the use of thermal therapies for treatment and healing. However, it was more of an art than a science. Modern scientific research unveiled that the deposited thermal energy may cause physiological changes due to molecular or cellular responses to thermotherapies. The human race has come a long way in using thermal effects for medical treatment by designing special purpose devices to control the delivery of thermal energy to targeted regions. Today, commonly used modalities of thermal medicine include cryotherapy, hypothermia, hyperthermia, and thermal ablation depending on the temperature range of operation. The delivery technology may involve laser, microwave, radio-frequency (RF), and focused ultrasound. Cutting edge delivery systems may also involve nanoparticle-mediated thermal therapy. As part of medical practice it is clear that the relationship between dosages and effects need to be established. Obviously it is an impossible task to determine that by experience. It is necessary to establish mathematical models based on underlying human and biology biophysics to predict the effects by deposited thermal energy via computer simulation.

Brief historical review

To predict tissue response due to thermotherapies there are at least two major mathematical components to be developed: heat transfer in tissues (commonly referred as bioheat transfer), and tissue damage (or cell damage in some literature) due to deposited thermal energy. Proposed in 1948, the Pennes bioheat equation was developed based on empirical observations [1]. It is now commonly used in computing temperature field in thermal therapies. The original form of the governing equation has three major assumptions: 1) the model properties such as mass density, heat capacity and thermal conductivity of tissue as well as metabolic heat parameter are constant, 2) the rate of heat transfer from blood to tissue is proportional to the temperature difference between tissue and blood in arteries, and 3) volumetric blood perfusion rate and arterial blood temperature are constant throughout the tissue. The advantages of the Pennes bioheat equation are the ease of implementation and that the perfusion rate can be calibrated to fit the experimental data using rather simple tests. However, tissue heterogeneity was not considered in the original publication and the perfusion term in the equation is a scalar, thus non-directional. The major criticism though is that the Pennes equation is unable to take into account thermal effect of large vessels in the tissue domain of interest.

Over the years, attempts were made to improve Pennes’ model and develop more comprehensive bioheat transfer models that account for tissue heterogeneity and effects of large blood vessels. In 1974, Wulff [2] published a modified version of Pennes’ equation by considering directional perfusion using so-called blood flux. In the same year, Klinger [3] further improved Pennes’ model by considering both spatial and temporal variations in convection term and heat source. The directional perfusion is also considered. However, the blood flow is treated as incompressible.

To determine the significance of blood flow in different sized blood vessels, Chen and Holmes [4] published a paper in the 1980s and separated perfusion into three modes: 1) the equilibrium of blood and tissue temperatures near large vessels, 2) the heat convection within tissue, and 3) thermal fluctuation. In this model the interface between blood and surrounding tissue is clearly defined. Effective thermal conductivity was derived by considering both solid tissue and blood. Also, the directionality of blood flow and simple vasculature are considered. However, these three perfusion modes, which could occur at different locations, were treated to coexist in the equation. Jiji et al. in 1984 [5] considered artery, vein and tissue separately and obtained a three-energy equation model. Their assumptions are: 1) small arteries and veins are parallel and the flow direction is countercurrent and 2) isotropic blood perfusion between vessels is assumed. This model was derived based on the applications of skeletal muscle tissue. Moreover, it was only a particular case of vascular configuration being chosen.

Xuan and Roetzel in 1997 [6] used a volume averaging approach with assumptions of porous media and local thermal nonequilibrium between blood and tissue. The difference of local temperatures of blood and tissues is considered. However, local heat transfer coefficient resulting from interface integral is difficult to compute unless more assumptions are introduced, which may or may not be physical. Roemer and Dutton in 1998 [7] introduced the volume averaged tissue energy equation. They used generic control volume which also resulted in interface integrals that are hard to evaluate unless simplification or approximation are introduced. Moreover, many model parameters need to be determined.

Nakayama and Kuwahara in 2008 [8] also used a volume averaging approach. The assumptions included: 1) tissues are filled with interstitial blood flow as porous media, 2) local thermal nonequilibrium exists between blood and tissue, and 3) interface integrals need to be simplified. The advantage of this model is that the difference between local temperatures of blood and tissue was considered. The mass and momentum equations are provided to complete the equation system. However, the information (such as diameters, locations of embedded vessels, and flow rate) for each averaging voxel volume is needed to determine the coefficients, as a result of deriving interface integrals, which are known to be hard to evaluate. Moreover, the resulting blood perfusion terms lack physical interpretation. Shrivastava and Vaughan in 2009 [9] used a volume averaged energy equation with two assumptions: 1) blood stays in vasculature, i.e. there is clear interface between blood and tissue, and 2) interface integrals were simplified. Both tissue and blood are considered. Similarly, information (such as diameters, locations of embedded vessels and flow rate) for each averaging voxel volume is needed to determine the coefficients, resulting from interface integrals as well. To date, a general bioheat model capable of considering the coupling effect of heat and transfer in tissues with embedded realistic blood vessels is still in need. For discussions of thermal damage modelling, which is the second component to model the outcome of thermotherapies, we refer to the articles by Pearce [10] and Whitney et al. [11] in this issue.

Summary of this special issue

The aim of this special issue is to present the latest developments in mathematical modelling with respect to the applications in thermal medicine. The articles included in this issue are summarised as follows.

Pearce [10] has reviewed theoretical models that have been used to predict thermal damage. He concludes that calculation of transient temperature fields alone is certainly not sufficient. Thus, useful models ought to predict tissue responses to adequately represent experiment analysis or to evaluate the effectiveness of new and/or existing modalities since heating time and the particular process under study are also dominantly important considerations.

Whitney et al. [11] have presented a study that aims to determine the effects of heating methods on the Arrhenius parameters and spatial viability response prediction. They demonstrate that heating methods have a greater effect on the Arrhenius parameters. This study recommends the Arrhenius parameters will be determined from thermally and biologically representative systems in order to obtain accurate spatial predictions of viability loss.

Prakash et al. [12] review the use of mathematical modelling for endoluminal and interstitial ultrasound hyperthermia and thermal ablation for device design, feedback control and treatment planning. Their main conclusion is that rapid computation techniques should be employed for integrating mathematical modelling in clinical settings.

Chiang et al. [13] have presented a paper on computational models for microwave thermal ablation. This paper includes a detailed review of the theory, and a numerical method for modelling thermal ablation with microwave ablation. These authors show how modelling assists in the design of devices and techniques to improve ablation procedures.

Reddy et al. [14] investigated the applicability of the thermal dose model to predict thermal damage of hepatocellular carcinoma cells subjected to 45–60°C for 0–32 minutes. Their results suggest that thermal dose model is not applicable for predicting thermal damage in this temperature range, in vitro. This report is of importance when considering the use of thermal dose model to predict the size of the necrosis zone following thermal ablation therapy.

Fahrenholtz et al. [15] have used the generalised polynomial chaos method to evaluate parameter uncertainties in modelling interstitial laser thermal therapy. The calculated temperatures are compared to real-time measurements of MR thermometry in vivo. This study shows that conservative estimates of the thermal field could be used for treatment planning, with uncertainties in optical properties and approximation of the Pennes bioheat equation.

Kok et al. [16] describe the progress made in the inclusion of discrete blood perfusion in mathematical modelling of hyperthermia. This paper presents a software package that includes the effect of detailed 3D discrete vessel networks on heat dissipation during hyperthermia. The authors conclude that this model could be used for treatment planning and optimisation of hyperthermia.

Paulides et al. [17] have reviewed simulation tools and techniques that were developed for clinical hyperthermia. The authors evaluate the current status computer models and their adequacy for clinical applications. They conclude that hyperthermia treatment planning has shown promise in improving treatment outcomes, and new non-invasive measurements techniques are expected to further improve accuracy for personalising thermal therapy, in the future.

This special issue covers the major advancement in the use of mathematical modelling in thermal medicine. The models allow us to study the effect of thermal dose on the response to thermal therapies. The advancement in computer simulations has enabled major improvements in the design of devices for thermal therapy. The inclusion of radiological imaging and real-time simulation and measurements is expected to enable the research community to tailor thermal therapies to individual patients.

Declaration of interest

The authors report no conflicts of interest. The authors alone are responsible for the content and writing of this article.