Thermal modelling using discrete vasculature for thermal therapy: A review

Abstract

Reliable temperature information during clinical hyperthermia and thermal ablation is essential for adequate treatment control, but conventional temperature measurements do not provide 3D temperature information. Treatment planning is a very useful tool to improve treatment quality, and substantial progress has been made over the last decade. Thermal modelling is a very important and challenging aspect of hyperthermia treatment planning. Various thermal models have been developed for this purpose, with varying complexity. Since blood perfusion is such an important factor in thermal redistribution of energy in in vivo tissue, thermal simulations are most accurately performed by modelling discrete vasculature. This review describes the progress in thermal modelling with discrete vasculature for the purpose of hyperthermia treatment planning and thermal ablation. There has been significant progress in thermal modelling with discrete vasculature. Recent developments have made real-time simulations possible, which can provide feedback during treatment for improved therapy. Future clinical application of thermal modelling with discrete vasculature in hyperthermia treatment planning is expected to further improve treatment quality.

• Ablation
• discrete vasculature
• hyperthermia
• hyperthermia treatment planning
• optimisation
• RF heating

Introduction

Reliable temperature information during clinical hyperthermia and thermal ablation is essential for adequate treatment control, both in the tumour and at potential hot spot locations in normal tissue. Clinical thermometry is usually performed invasively, with thermometry probes inserted in closed end catheters. The number of invasive measurements are limited to avoid pain and infections. When heating pelvic tumours, intraluminal or intracavitairy thermometry probes are often placed in existing lumina/cavities close to the target (e.g. vagina, rectum or urethra). These standard sparsely sampled temperature measurements do not provide sufficient information to assure an adequate thermal dose as the temperature distribution is generally non-uniform. This is due to heterogeneity of dielectric and thermal properties of human tissue, which causes a non-uniform power distribution and temperature rise. Blood flow accounts for up to 90% of heat removal [1] and can increase up to a factor of 10 in response to hyperthermia, especially in normal tissue [2]. Non-uniform heat removal by blood flow further aggravates the temperature heterogeneity.

Non-invasive 3D thermometry techniques are under development. Most are still insufficiently accurate or practical for clinical application during hyperthermia treatments. Microwave radiometry has shown promise for future application to non-invasive monitoring of volume average temperatures to several cm depth [3,4]. For true 3D thermal imaging deep in the body, MRI techniques seem to have the best potential for clinical use [5–7]. MR-based temperature measurements are successfully applied for thermal ablation [8]. MR thermometry combined with hyperthermia has been used in monitoring abdominal and pelvic tumour sites [9,10] and has demonstrated good correlation with conventional measurements, which makes it feasible for validating the quality of heating. MR temperature measurements in extremity soft tissue sarcomas have also shown good correlation with conventional measurements, with a mean difference ≤1 °C and have shown promise for real-time feedback to treatment control [11]. Although results are very promising, MR thermometry cannot yet replace conventional invasive thermometry during hyperthermia, since the much smaller temperature range compared to thermal ablation requires a higher accuracy of the temperature measurements. Furthermore, some heating systems may contain large metal structures, which hamper combining heating with MR thermometry.

Hyperthermia treatment planning is another very useful tool to improve treatment quality, and significant progress has been made over the last decade. Prospective treatment planning has shown its value for improving the heating effectiveness during locoregional hyperthermia [12,13]. Because of the significant heat removal by blood flow, thermal modelling is a very important and challenging aspect of treatment planning. Several thermal models with varying complexity have been developed for this purpose.

Various bioheat transport models have been summarised by Bhowmik et al. [14]. The continuum model proposed by Pennes [15] is the most straightforward thermal model and remains the most frequently used. To model heat removal by blood flow, the Pennes bioheat equation applies a heat sink proportional to a non-directional volumetric average perfusion rate and local temperature rise. Baish et al. have shown that this model can adequately predict the thermal effect of flow in thermally poorly equilibrated parallel tubes in a dynamic tissue equivalent phantom [16]. The isotropic distribution and the poor equilibration of these tubes produce an isotropic heat sink similar to that of the Pennes equation. In the human body, however, blood vessels branch into increasingly smaller vessels which ultimately reach thermal equilibrium with the surrounding tissue. A drawback of the Pennes model is that it accounts neither for heat exchange between tissue and realistic vasculature, nor for the direction of blood flow. Several adaptations of this simple continuum model have been proposed to account more accurately for blood flow [17–19]. An enhanced effective tissue conductivity has been developed, in which the effective conductivity is described by a tensor, thereby including information about the direction of the local blood flow [18,20]. However, this model applies only for smaller vessels (<∼0.5 mm) and does not account for the presence of large blood vessels.

Thermal simulations are most accurately performed by modelling thermally significant vasculature discretely. In general, vessels are considered thermally significant when the thermal equilibration length (i.e. the length over which inlet blood vessel temperature equilibrates with the surrounding tissue temperature) is comparable to or larger than the actual vessel length [18,21]. The presence of a large thermally significant vessel or a large counter-current artery–vein pair entering a heated volume can have a significant cooling effect [22,23]. The first models including discrete vasculature studied the basic impact on heat transfer and were limited to straight vessels, but later also modelled more realistic branching and curving vessel architecture.

This review describes the progress in thermal modelling with discrete vasculature for the purpose of treatment planning for hyperthermia, starting with more basic vascular models suitable for research purposes and thermal ablation planning, then progressing to more sophisticated vascular models suitable for clinical hyperthermia treatment planning. Computation times, data acquisition and other challenges are discussed. Table I summarises the categories of discrete vessel models and the corresponding applications.

Table I. Overview of the categories of discrete vessel models and the corresponding applications. The pyramid at the left indicates that the number of models decreases with increasing complexity.

Thermal models with discrete vasculature

Basic discrete vessel models

Some of the discrete vessel models have been designed not for application in hyperthermia treatment planning, but to obtain basic insight into the heat transfer between large blood vessels and tissue. Therefore, these models are relatively simple with straight vessels represented as a tube with a specified diameter [22–25]. These basic models are also very useful to calculate temperature distributions induced by thermal ablation, where it is necessary to account for a single large blood vessel passing through a target region. Many researchers investigated the basics of heat transfer using basic discrete vessel models and a selection are discussed below.

The model developed by Lagendijk calculates the temperature distribution around a large blood vessel in the heated region during hyperthermia [23]. It was demonstrated that the cooling effect of inflowing blood on the surrounding tissue temperature can cause underdosing of parts of the tumour.

Mitchell and Myers developed an analytical model to analyse counter-current heat exchange phenomena [26]. It was observed that significant counter-current heat exchange occurs only if the conductance between arteries and veins is high. Additionally, the counter-current effect becomes more significant with decreasing blood flow rate.

Crezee and Lagendijk studied the impact of single and counter-current paired large vessels on the temperature uniformity during hyperthermia [22]. The vessel wall temperature determines the minimum tumour temperature, and it was shown that even when the blood temperature is low, the vessel wall temperature may reach therapeutic levels in specific cases, especially when the perfusion of surrounding tissue is high. Thus, a high tissue blood flow reduces the temperature gradient near large vessels, which implies that it is more difficult to achieve a uniform dose in poorly perfused than in well perfused tumours. They also showed that including a relatively large margin of normal tissue in the heated volume will improve the temperature homogeneity in the tumour region.

Rawnsley et al. included three to four counter-current paired blood vessels both in the Pennes bioheat equation and in the effective thermal conductivity equation [24]. Comparison with temperature profiles measured in canine thighs heated with scanned focused ultrasound showed that adding this small number of vessels yielded a significant improvement in predicted temperature profiles near large vessels.

Chen and Xu studied tissue temperature oscillations due to vascular thermoregulation in response to surface heating in an isolated pig kidney [27], both experimentally and using a numerical model. The model consisted of a counter-current artery–vein pair in a tissue unit. A good agreement was found between the simulations and experimental results. It was observed that the heating rate, the time delay and rate of the change in local flow are the main parameters affecting temperature oscillations.

Valvano et al. developed a 3D small artery model to describe the heat transfer contribution of perfusion in the canine cortex [25]. The finite difference model was based on the physical anatomy of the kidney cortex vasculature, with interlobular arteries running approximately parallel to each other. Arteries and veins were evenly distributed in the model with an overall density of 72 arteries and 72 veins per cm² of tissue. The steady-state and sinusoidal temperature response of self-heated thermistors were correctly predicted by the model.

Heat transfer coefficients used to calculate heat transfer to large blood vessels depend on, for example, entrance effects, tissue thermal resistance and vessel shape, and may vary circumferentially around the vessel during transient increases in temperature. Therefore, algorithms were developed to solve temperature distributions in tissue without using constant heat transfer coefficients that are valid only under specific theoretical conditions. An example is the model described by Moros et al. [28], in which heat transfer between vessels and tissue is modelled implicitly with one field equation. This is also less computationally intensive compared to a coupled set of two equations. They illustrated the model by simulation of a thermally significant counter-current vessel pair of 0.6 mm diameter, embedded in a tissue volume and feeding a tumour heated in the hyperthermia temperature range. These simulations showed that the temperature of the blood entering the heated tumour volume remained below therapeutic levels.

A model using a similar solution to avoid using heat transfer coefficients is described by Kolios et al. [29,30]. A cylindrical model with angular symmetry was solved with a variable grid size to concentrate the density of grid points near the large vessel. This cylindrical model provides relevant insight into the large temperature gradients that exist near large vessels in heated tissue [29]. Additionally, it was shown that the cooling effect of large vessels is largely determined by perfusion in the surrounding tissue, whether it was modelled by the Pennes equation or by effective conductivity [29,30].

Discrete vessel modelling with networks of straight vessels

In addition to the basic discrete vessel models, more advanced 3D models with branching vessel networks provide more insight in heat transfer processes in tissue heated by hyperthermia. Brinck and Werner developed a 3D vascular model in which convective heat transfer between branching counter-current vessel pairs and tissue was calculated [31,32]. The Nusselt number is the only assumption made in this model using a finite difference approach with an irregular grid, which allows a small voxel spacing around the vessels and in the space between vessels to improve accuracy. Heat transfer from the vasculature was calculated individually, combined with a continuum formulation using the local arterial blood temperature. The model was applied to the cross-sectional area of a human limb with pre-arteriole and post-venule vessels of the muscle layer incorporated in counter-current pairs consisting of eight generations. Simulation results showed that the temperature gradient around the artery is larger than the temperature gradient around the vein. This can be explained by the fact that the heat flow between arteries and tissue is larger than the heat flow between veins and tissue. The mean tissue temperature perpendicular to the vessels was thus closer to the temperature of the venous blood, while the mean tissue temperature in between the artery and the vein is represented by the mean arterial and venous blood temperatures.

Zhu et al. investigated the thermal interaction between counter-current supply artery–vein (SAV) pairs from 300 to 1000 μm and surrounding tissue for two vascular branching patterns under normal and hyperaemic conditions [33]. In the ‘pure branching’ pattern a vessel branches into two equal vessels in each successive generation, with no capillary bleed-off. In the ‘pure perfusion’ pattern a SAV vessel tapers as the smaller vessels branch from the main axial vessels. Analytical solutions showed that the temperature variations along these vessel pairs strongly depend on the branching pattern and the local perfusion rate. These models could be extended to investigate the thermal interaction under hyperthermic conditions.

Huang et al. [34] modelled a simple counter-current generic vessel network to demonstrate the potential usefulness of fully conjugated vessel network models. In this model, vessels are straight line segments running parallel to one of the Cartesian axes, which implies that vessels bend at 90 ° angles. Seven generations were modelled and the diameters decreased by a constant ratio at branching points. A finite difference scheme using an equidistant grid was applied. They demonstrated that the model could be used to predict the significance of heat removal by arteries and veins and to investigate thermal significance of lower generation vessels. Additionally, a more extensive 3D network of 8178 vessels with diameters ranging from 10 to 1000 μm was used to address the issue of thermally significant blood vessels in a counter-current network [35]. The results showed that the thermal significance of vessels should be determined on a case-by-case basis.

Blanchard et al. proposed an efficient finite element–finite difference method for modelling non-orthogonal vessel networks [36]. This method combines finite elements with a 1D finite difference method for the interior of the vessels. The method was applied to a four-level counter-current vessel network and comparison with a pure finite element method demonstrated that the hybrid technique combines the accuracy and flexibility of pure finite element techniques with computational speed of 1D finite difference computations.

More sophisticated discrete vessel models have been developed for hyperthermia treatment planning purposes. Lagendijk et al. [37] developed a 3D finite difference model to calculate temperature distributions in vascularised tissues during hyperthermia. The model combines discrete vasculature with a grid of rectangular tissue voxels where the effect of microvasculature is described either by Pennes’ bioheat equation or with the effective conductivity model. The discrete vessels are embedded at the centre of tissue voxels and divided into elements with a length equal to the voxel size in the finite difference modelling. Thus vessels always run through a straight row of voxels and bend at angles of 90 °. Blood vessels are inserted such that the wall temperature is easily computed in a finite difference scheme by dividing the tissue voxel around the vessel into four boundary elements. Each boundary element has a vessel diameter-dependent irregular shape such that the distance between the vessel wall and the voxel border is equidistant to the distance between the voxel centre and border in all neighbouring tissue voxels (Figure 1). The four vessel wall temperatures are assumed equal. The discretisation depicted in Figure 1 allows a relatively easy update of the heat balance for the blood element and the vessel boundary voxels.

Figure 1. (A) A blood vessel embedded in a tissue voxel with surrounding tissue voxels 1–8. The tissue voxel around the vessel has been divided into four boundary elements. The dotted lines indicate the lines along which the conductive heat flow in the tissue is calculated in the model. (B) shows the geometry around a vessel for two different vessel radii [37]. This picture has been reproduced from Lagendijk JJ et al. (A three-dimensional description of heating patterns in vascularised tissues during hyperthermic treatment. Phys Med Biol 1984;29:495–507. http://dx.doi.org/10.1088/0031-9155/29/5/002). © Institute of Physics and Engineering in Medicine. Published on behalf of IPEM by IOP Publishing Ltd. Reproduced by permission of IOP Publishing. All rights reserved.

The model has been experimentally verified in several phantom experiments, and an overall accuracy better than 0.3 °C was found. However, there are some limitations hampering modelling of detailed realistic vasculature in the human body. The maximum diameter of the modelled vessels is equal to the voxel spacing and the centres of two vessels should be at least three voxels apart. Another drawback of this model is the use of straight vessels. Extension of the model for more general cases of bending and branching vessels with arbitrary cross-sectional shapes turned out to be very complicated.

Discrete vessel modelling with semi-curved vessel networks

A more flexible algorithm in which the model geometry is subdivided into a vessel space and a tissue space has been developed by Mooibroek and Lagendijk [38] to obtain a more realistic, semi-curved representation of vessel networks for use in hyperthermia treatment planning. In this finite difference model, vessel segments are represented by connected strings of vessel nodes. The nodes with their centres closest to the vessel axis are considered vessel nodes. The description of a 3D vessel segment is used as a building block, which allows modelling of realistic complex vessel networks. A vessel-specific second discretisation step in time is performed to describe convective heat transfer within the vessel space. This makes it possible to incorporate vessels with different flows and diameters.

The implemented discretisation routine for the vessel segments yields an overestimation of the actual path length and thus of the flow and heat transport along a bending vessel trajectory. A correction factor has been implemented to minimise this effect. Additionally, correction factors were incorporated to correctly represent the cylindrically shaped vessels on a square grid. Simulations with a single vessel were compared to results obtained with a benchmark cylindrical model at flow velocities of 0.05 and 10 cm/s. The largest deviation for the low flow situation was 0.2 °C. Increasing flow to 10 cm/s reduced the maximal deviation to 0.12 °C. Furthermore, predicted thermal equilibration lengths were compared to theoretical values, and a reasonable agreement was observed. Comparison with phantom measurements showed an accuracy within the range of the experimental error. This shows that cylindrically shaped vessels can be modelled accurately using a square grid.

This model is a clear improvement compared to the preceding model by Lagendijk et al., since important features such as branching, bending, obliqueness, widening and tapering are preserved. However, because of the cubic subdivision of blood vessels, a one-to-one description of a true 3D vessel network cannot be realised with this model. The vessel discretisation yields stair-casing errors and the accuracy of the vessel modelling depends on the voxel size. This is a drawback for clinical applications in treatment planning, where modelling of more detailed 3D vessel networks is required for accurate thermal simulations.

Discrete vessel modelling with real 3D vasculature networks

To overcome the above-mentioned problems, a more sophisticated thermal model that allows thermal modelling with real detailed 3D discrete vasculature networks (DIVA) has been developed by Kotte et al. [39,40]. Similar to Mooibroek and Lagendijk [38], the model geometry is subdivided into a vessel space and a tissue space, but in the DIVA model the vasculature is described by 3D curves with a specified diameter. This approach allows taking into account all relevant blood vessels independent of the voxel size. Besides the clear advantage of tissue voxel resolution independence, modelling blood vessels as 3D curves makes the model compatible with MR/CT angiography vessel reconstruction software, which is essential for routine use in treatment planning.

The exact heat flow between a vessel segment and its surrounding tissue is difficult to calculate for a realistic situation with a heterogeneous vessel network, power and temperature distribution. To solve this problem, a method to estimate the heat flow using a simplified situation has been developed [39]. For a vessel segment embedded in a tissue cylinder, an analytical expression for the heat flow can be derived when cylindrically symmetric boundary conditions and a thermally fully developed flow are assumed. The latter is justified since the entrance length is considerably shorter than the equilibration length [22].

Vessel segments are discretised with a user-defined sample density (m⁻¹). A discrete vessel sample is labelled a ‘bucket’ (Figure 2A). To handle the tissue–vessel interaction two sets of voxels are defined for each bucket: the estimation set and the exchange set. The estimation set voxels are used for calculation of the tissue–vessel interaction. The exchange set voxels are then used to effectuate the heat exchange in tissue. Figure 2(B) illustrates the definition of the exchange and estimation sets for a vessel cross-section projected onto a tissue grid. For each voxel in the estimation set, an estimate for heat flow rate density is calculated. The average value is multiplied with the bucket surface to yield the heat flow rate into the bucket. The heat exchange with tissue is distributed evenly over the exchange set voxels. Thermal equilibration of the blood from the final generation of discrete vessels is completed by attaching user-defined sub-volumes to vessel end points and applying a local heat sink [40].

Figure 2. Projection of a vessel segment (A) and a vessel cross-section (B) onto a tissue grid, together with the associated estimation and exchange sets. (A) shows three discrete vessel samples, which are labelled as ‘buckets’ [39]. Heat flow towards the vessel is withdrawn from the tissue in the exchange set voxels. The estimation set voxels are used to calculate the thermal interaction between the vessel and surrounding tissue. This picture has been reproduced from Kotte et al. (A description of discrete vessel segments in thermal modelling of tissues. Phys Med Biol 1996;41:865–84. http://dx.doi.org/10.1088/0031-9155/41/5/004). © Institute of Physics and Engineering in Medicine. Published on behalf of IPEM by IOP Publishing Ltd. Reproduced by permission of IOP Publishing. All rights reserved.

Validation of the DIVA thermal model

During the development and validation of the DIVA model, special attention was given to the overall heat balance [40]. The correctness of the model has been validated by comparison with analytical solutions [39]. Simulations of a single vessel embedded in a coaxial tissue cylinder showed deviations <1% between simulated and theoretical thermal equilibration lengths. Simulated and theoretical axial blood temperature profiles also corresponded within 1%. The position and orientation relative to the tissue grid were varied, as well as vessel radius, vessel sample density and power deposition in tissue. Excellent agreement with theory was observed for all cases. The correctness of modelling vessel segment connections was tested by subdividing a single segment into three connected segments. Deviations between simulated and theoretical axial blood temperature profile were again <1%.

Next, the DIVA thermal model was validated experimentally using isolated perfused animal tissue [41]. An isolated bovine tongue was heated with three hot water tubes, and temperature profiles were measured at perfusion levels of 0, 6 and 24 mL (100 g)⁻¹min⁻¹. Cryo-microtome slices were used to reconstruct the experimental set-up, including the discrete vasculature down to 0.5 mm diameter. Comparison of measured and simulated temperature profiles showed excellent agreement, with deviations in the order of the experimental error.

Applications using the DIVA thermal model

The DIVA thermal model has been applied to investigate the thermal impact of detailed 3D discrete vasculature for a wide variety of applications. Craciunescu et al. derived discrete vasculature and perfusion maps using MR imaging for a patient with a high grade sarcoma [42]. Simulations were performed for hyperthermia with an annular phased array applicator. Temperatures were calculated using Pennes’ bioheat equation, perfusion maps and large vasculature extracted from MR angiography combined with perfusion maps and artificially generated vasculature using the algorithm by Van Leeuwen et al. [43]. The different simulation results were compared to MR thermometry and the simulation using tracked vasculature combined with perfusion maps showed the best correspondence.

Van den Berg et al. reconstructed the discrete vasculature in the pelvis from a dynamic contrast-enhanced computerised tomography (CT) scan of a prostate patient. Additionally, a detailed model of the prostate vasculature was reconstructed from a post-mortem prostate. Simulations with DIVA were compared to results with the conventional Pennes model for locoregional heating of the pelvis [44]. Differences in predicted temperatures were up to 1–2 °C, which clearly demonstrates the need for discrete vessel modelling in treatment planning for locoregional hyperthermia. Furthermore, this study showed that pre-heating of the prostate vessels differs from vessel to vessel, as illustrated in Figure 3.

Figure 3. The vessel network in a prostate, together with the simulated 40.5 °C iso-temperature surface for a homogenous SAR level of 20 W kg⁻¹. Pre-heating of the prostate vessels differs from vessel to vessel [44]. This picture has been reproduced from Van den Berg et al. (Towards patient specific thermal modelling of the prostate. Phys Med Biol 2006;51:809–25. http://dx.doi.org/10.1088/0031-9155/51/4/004). © Institute of Physics and Engineering in Medicine. Published on behalf of IPEM by IOP Publishing Ltd. Reproduced by permission of IOP Publishing. All rights reserved.

The DIVA model has also been used in applications other than hyperthermia treatment planning, for example to address radiofrequency safety issues. Van Leeuwen et al. applied discrete vasculature modelling to calculate the temperature change in the brain due to exposure to a mobile phone [45] and Flyckt et al. investigated the temperature rise in the human eye and orbit exposed to a mobile phone [46]. Van Lier et al. evaluated the temperature rise in the head after exposure to a 300 MHz radiofrequency field induced by a 7-T MRI coil [47]. The DIVA thermal model is also useful to evaluate the effect of hypothermia therapies. Van Leeuwen et al. performed thermal simulations using DIVA to investigate application of a cooling bonnet for hypothermia as a neuroprotection therapy after hypoxia-ischaemia in newborn infants [48].

Fast simulations of the steady-state temperature using DIVA

A drawback of the implementation of the DIVA model as described above is the finite difference approach, applying many small time steps to calculate the evolving temperature distribution until a user-defined end time. For treatment planning purposes, the steady-state temperature distribution is most relevant since it reflects the thermal dose in the tumour and the location and magnitude of hot spots in normal tissue. However, calculation of a steady-state temperature distribution can take up to many hours, depending on the level of detail of the vasculature included. This makes it impractical for simulations in routine clinical applications.

To overcome this problem, a new calculation method for fast simulation of the steady-state temperature distribution has been developed recently by Kok et al. [49]. The bioheat transfer and the interaction between vasculature and tissue were represented by a linear system Ax = b. The coefficients in the matrix A describe the thermal processes. The blood and tissue temperatures are represented by the vector x and the right-hand side vector b holds the boundary conditions and the power deposited by the heating system. The linear system was efficiently solved using an iterative method. To reduce computation times solving of Ax = b was implemented running in parallel on a graphical processing unit (GPU).

Simulations were performed for pelvic heating of a human anatomy with realistic 3D vasculature obtained from angiography [49]. Both the arterial and venous networks consisted of 174 segments and 93 end points with a diameter of 1.2 mm. Results showed that the steady-state solver was significantly faster than the original implementation of DIVA using time-stepping. Calculation of the steady-state temperature distribution was reduced from several hours to only ≤1.5 min [49], which makes it suitable for real-time treatment planning during hyperthermia treatment.

Treatment optimisation using discrete vasculature-based thermal modelling

During clinical hyperthermia and thermal ablation, a homogeneous tumour temperature is pursued. However, achieving this is a challenge due to heterogeneous perfusion. An additional difficulty for (locoregional) hyperthermia treatments are hot spots in normal tissue, which limit overall power deposition. Optimisation strategies that minimise hot spots by intelligent phase-amplitude steering can help improve the achieved temperature distribution. The vast majority of temperature-based optimisation methods described in the literature apply the Pennes bio-heat equation to calculate temperature distributions during the optimisation process. Taking discrete vasculature into account in these optimisation methods significantly increases the computational complexity.

Some optimisation methods have been proposed that account for a single large blood vessel passing through a target region, a situation relevant for both hyperthermia and thermal ablation. Huang et al. developed a fast and robust temperature-based adaptive power scheme to optimise tumour temperature, taking into account the thermal influence of a large thermally significant blood vessel in the tumour region [50]. Shih et al. investigated the cooling effects of thermally significant blood vessels on the extent of a thermal lesion after ultrasound thermal therapy [51]. Several different heating schemes were simulated for blood vessels with a diameter of 200 μm and 2 mm. Their results showed that short duration high intensity heating yielded the best target coverage for the case of a blood vessel with 200 μm diameter. For the case of a much larger 2-mm diameter blood vessel however, complete necrosis could not be achieved.

Optimisation becomes even more computationally intensive when including vessel networks, rather than single vessels. Huang et al. demonstrated the importance of taking the thermal influence of multiple blood vessels into account during optimisation of power and temperature distributions for hyperthermia treatments [52]. A counter-current network of straight vessels, consisting of seven generations was modelled. Optimisation was performed by continuously adjusting power deposition in the target region, thereby aiming for a homogenous tumour temperature. Results showed that high resolution power steering is needed to realise a homogenous temperature throughout the target region.

Clinical application of treatment planning can help improve the temperature distribution [12,13] and real-time application of temperature-based optimisation with detailed discrete vasculature can be expected to further improve treatment quality. However, optimisation requires many time-consuming recalculations of the temperature distribution, which is a limitation for real- time application. To overcome this problem Kok et al. recently devised a very efficient temperature-based optimisation method based on the DIVA thermal model [49]. During optimisation, the steady-state temperature distribution is calculated efficiently using superposition of pre-computed temperature distributions. This principle was first introduced by Das et al. [53] for the Pennes bioheat equation and extended to the DIVA model to include real 3D vasculature. The power, tissue temperatures and mixing cup blood temperatures are written as a vector-matrix-vector multiplication. The temperature matrix elements can be solved using the efficient steady-state solver developed by Kok et al. described earlier in this review [49]. Calculation of the temperature matrices can be performed in a pre-processing step, which makes the optimisation applicable for real-time feedback control applications.

Kok et al. compared results for optimisation based on Pennes bioheat equation with results obtained after optimisation based on DIVA for pelvic heating of a human anatomy with realistic 3D vasculature obtained from angiography [49]. Large differences were found both in optimised phase-amplitude settings as well as in predicted optimal tumour temperatures. This shows the necessity for temperature-based optimisation based on DIVA to improve tumour temperatures.

Vessel generation software

Patient-specific discrete vessel input for thermal simulations can be obtained using CT or MR angiography. However, not all thermally significant vasculature might be visible on the angiogram. To overcome this problem, vessel generation software has been developed to create vessel networks or add thermally significant vasculature to incomplete networks derived from an angiogram.

Details of the branching structure of vessel trees can be described accurately using mathematical relations. Several algorithms have been developed to realistically model the growth of blood vessels in 2D [54–57]. A 3D vessel network construction algorithm, modelling the diffusion of angiogenic factors in tissue has been developed by Baish [58]. This 3D algorithm is an extension of Gottlieb’s 2D algorithm [55] and the generated vessel trees closely match real vasculature in number and diameter of vessels as well as the degree of asymmetry in branching. However, a drawback is the alternating growth process, which hampers automatic generation of vasculature in real anatomies.

Van Leeuwen et al. developed a sophisticated 3D algorithm to construct vessel networks that may consist of both arteries and veins [43]. The algorithm also provides the possibility of generating counter-current networks. To create new vessel segments a potential function is defined in the tissue to control the paths by which points are connected to existing vessels. The presence of bone and cavities or implants can be taken into account during the vessel construction process to allow generation of vasculature superimposed into real anatomy. Furthermore, the desired perfusion distribution can be defined and the algorithm can be used to expand vessel networks obtained from angiography. Evaluation showed that physiologically realistic networks can be created with this algorithm [43]. Figure 4 shows an example of expanded vasculature in the pelvic region using this algorithm.

Figure 4. (A) Large arteries (red) and veins (blue) reconstructed from a CT angiogram of a human pelvis, together with the bony structures. (B) Expanded vasculature using the vessel generation algorithm devised by Van Leeuwen et al. [43].

A similar algorithm for vascular network construction has been devised by Prishvin et al. [59]. In this model the artificially generated vasculature is discretised into a semi-curved vessel network and combined with an efficient finite difference-based thermal model using a modified Pennes equation and a blood velocity vector field to model blood flow in smaller sized vessels [59].

Summary and future perspective

Substantial progress has been made towards clinical application of hyperthermia treatment planning over the last decade. Most applications so far use specific absorption rate (SAR) modelling, but this does not account for the significant thermal redistribution effects of perfusion and cooling by large vasculature, which hamper accurate prediction of thermal dose and treatment-limiting hot spot locations. For the future, thermal modelling that includes realistic heterogeneous tissue perfusion will be increasingly important for more accurate treatment planning that can further improve treatment quality.

A large number of basic models have been developed to provide fundamental insight into heat transfer between large blood vessels and tissue. These basic models are also very useful to calculate and optimise temperature distributions induced by thermal ablation, where it is necessary to account for a single large blood vessel passing through a target region. Furthermore, a combination of a continuum model and a small number of large blood vessels can be very useful to predict the most significant cooling effects by the vasculature for local and locoregional hyperthermia applications.

For more accurate planning of locoregional heating techniques more advanced discrete vessel models are required. To date, there has been significant progress in thermal modelling with discrete vasculature and recent developments allowing real-time thermal simulations during treatment using the DIVA thermal model with real 3D discrete vasculature. Furthermore, a very efficient temperature-based optimisation technique has been devised for the DIVA model, enabling real-time applications. To date the DIVA thermal model has been applied in a variety of applications and is expected to increase into widespread use for hyperthermia treatment planning. A simplified version of DIVA using boundary conditions to account for the effect of large blood vessels rather than modelling real 3D curved vasculature is already commercially available in SEMCAD X (SPEAG, Zurich, Switzerland). Although this implementation can be further improved, this is a first step towards the use of advanced discrete vessel models in commercial software. Actual clinical use of such models will lead to model refinement and is likely to motivate increased use of imaging to obtain more detailed vasculature data. When such detailed vasculature information is not (yet) available, vessel generation software can be used to provide realistic vasculature as input for thermal simulations. Smaller vasculature not modelled discretely can be described using a continuum model.

Accurate knowledge about blood flow and perfusion is important for reliable thermal simulations. Blood flow under normo-thermic conditions can be measured using dynamic contrast enhanced CT/MR imaging or arterial spin labelling with MRI, but these seriously underestimate perfusion during hyperthermia. On-line implementation of these imaging methods to assess hyperthermia-induced changes in perfusion remains difficult [60]. Therefore, literature values are used in treatment planning, though these are approximate since perfusion varies spatially and temporally, as well as with temperature and physical condition of the patient. Note that results of a recent study comparing the reliability of SAR-based and temperature-based optimisation techniques for cases with uncertain perfusion showed that temperature-based optimisation is superior to SAR-based optimisation even when uncertainty in both the absolute perfusion level as well as in the intra-tissue variation is taken into account [61].

Thus, despite the present uncertainty in perfusion, there is mounting evidence that ultimately thermal modelling is essential for clinical hyperthermia treatment planning. Prospective treatment planning using temperature-based optimisation has already proven its potential to improve the heating effectiveness during locoregional hyperthermia [12]. In future applications, adaptive treatment planning based on thermal modelling with discrete vasculature is expected to further improve tumour temperatures.

Conclusions

There has been substantial progress in thermal modelling with discrete vasculature. Many basic models have been developed to provide insight into the cooling effects of vasculature and temperature gradients around large blood vessels. Modelling of straight and semi-curved vessel networks has led to improved characterisation of heat transfer between vasculature and tissue. Fully realistic vascular DIVA-class models now allow thermal modelling with real detailed 3D discrete vessel networks for application in treatment planning, and recent developments have made real-time modelling during treatment possible. Sophisticated vessel generation software is available to expand thermally significant vasculature not visible on angiograms. Furthermore, a very efficient temperature-based optimisation technique has been devised for the DIVA model, enabling real-time applications. Future clinical application of DIVA-class thermal models in hyperthermia treatment planning is expected to further improve the quality of hyperthermia treatments.

Declaration of interest

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