Advanced patient-specific hyperthermia treatment planning

Abstract

Hyperthermia treatment planning (HTP) is valuable to optimize tumor heating during thermal therapy delivery. Yet, clinical hyperthermia treatment plans lack quantitative accuracy due to uncertainties in tissue properties and modeling, and report tumor absorbed power and temperature distributions which cannot be linked directly to treatment outcome. Over the last decade, considerable progress has been made to address these inaccuracies and therefore improve the reliability of hyperthermia treatment planning. Patient-specific electrical tissue conductivity derived from MR measurements has been introduced to accurately model the power deposition in the patient. Thermodynamic fluid modeling has been developed to account for the convective heat transport in fluids such as urine in the bladder. Moreover, discrete vasculature trees have been included in thermal models to account for the impact of thermally significant large blood vessels. Computationally efficient optimization strategies based on SAR and temperature distributions have been established to calculate the phase-amplitude settings that provide the best tumor thermal dose while avoiding hot spots in normal tissue. Finally, biological modeling has been developed to quantify the hyperthermic radiosensitization effect in terms of equivalent radiation dose of the combined radiotherapy and hyperthermia treatment. In this paper, we review the present status of these developments and illustrate the most relevant advanced elements within a single treatment planning example of a cervical cancer patient. The resulting advanced HTP workflow paves the way for a clinically feasible and more reliable patient-specific hyperthermia treatment planning.

Keywords:

• Hyperthermia treatment planning
• dielectric imaging
• EPT
• convection modeling
• thermal modeling
• discrete vasculature
• biological modeling

Introduction

Clinical hyperthermia aims at heating the tumor to a temperature of 39–43 °C in order to increase the therapeutic effect of radiotherapy and/or chemotherapy. Several clinical trials have demonstrated that complementing these two standard oncological treatments with hyperthermia improves the tumor response and patient survival for different tumor sites (e.g., cervix [1,2], breast [3], bladder [4], soft tissue sarcoma [5]). In general, a therapeutic thermal dose of 43 °C for 1 h is pursued in the tumor [6,7], in view of the demonstrated thermal dose-effect relationship [8–10]. This desired tumor thermal dose is often not achieved because of treatment limiting hot spots exceeding 45 °C in normal tissue. Achieving optimal tumor heating while avoiding hot spots in healthy tissues is a clinical challenge. Hyperthermia treatment planning (HTP) can be helpful to improve treatment quality, especially for locoregional hyperthermia as applied to deep-seated (pelvic) tumors, for which the use of HTP is recommended by the most recent quality assurance guidelines [11]. Locoregional hyperthermia is usually applied by heating systems consisting of phased-arrays of radiofrequency (RF) antennas, operating at frequencies between 60 and 130 MHz for deep-seated pelvic tumors [12,13], between 300 and 600 MHz for brain tumors [14] and at 434 MHz for head and neck tumors [15]. Thus, some extent of normal tissue heating is inevitable and several degrees of freedom (i.e., the amplitudes and phases of the individual antennas) are available to optimize the temperature distribution. In this respect, HTP can be very useful to predict optimal settings for effective tumor heating during thermal therapy delivery [16–18]. In this introduction, we review the current practice in HTP workflow and the recent research progress made to improve HTP accuracy.

In the typical clinical HTP workflow, an anatomical patient model is generated by segmenting a pretreatment CT or MR image [18–20]. For each segmented tissue, dielectric (or electrical) and thermal tissue properties are normally assigned from tabulated literature data. Then, electromagnetic and thermal simulations are performed to obtain the electric field, power absorption and temperature distribution in the patient’s body for selected phase-amplitude settings. Thermal modeling is generally based on Pennes’ bioheat equation [21], which models the heat transfer from blood to tissue as a heat sink term proportional to the local volumetric blood perfusion and the difference between local tissue temperature and the assumed blood temperature. Subsequently, automatic optimization techniques [16,22–25] predict the phase-amplitude settings for the individual antennas to realize optimal heat deposition in the tumor while preventing overheating of normal tissue. Notwithstanding various studies reported good qualitative agreement between HTP simulations and clinical observations during treatment [26–29], present hyperthermia treatment plans lack quantitative accuracy because of uncertainties in dielectric [30–32] and thermal [25] parameters and because of approximate heat transport modeling of blood flow and heat convection in fluids [33–35]. Furthermore, these hyperthermia treatment plans report only power and temperature distributions, but do not allow easy interpretation of the increased radiosensitization and therefore of the added clinical benefit to the primary treatment. Over the past years, substantial effort has been put into research to improve the reliability of hyperthermia treatment planning.

An important advancement is the measurement and inclusion of patient-specific electrical tissue properties. Indeed, uncertainties (up to 50%) in dielectric properties, in particular in the electrical conductivity, can cause average errors of 20% in the predicted power deposition [30], therefore leading to suboptimal prediction of tumor heating and hot spots outside the target region. Permittivity and conductivity maps were initially obtained from the water and fat content determined using standard MRI sequences [36,37]. A shortcoming of these methods is that they are calibrated on literature values and thus they neglect the impact of heterogeneities in ionic concentration, blood content and tumor pathophysiology in the calculation of the electrical properties. Subject-specific maps of the electrical conductivity were later reconstructed with another MR-based technique, called MR Electrical Properties Tomography (EPT), which relies on measurements of the RF magnetic transmit field, also called B₁⁺, generated within the body during an MR scan [38,39]. Balidemaj et al. [40] reported a substantial effect of including EPT conductivity values in the patient models for HTP, compared to using literature values: differences in the predicted temperature ranged from 0.6 °C to 1.5 °C in five patients. However, the accuracy of EPT at tissue interfaces and its robustness against noise are limited. Lately, Mandija et al. [41] have pioneered an EPT technique exploiting Deep Learning (DL-EPT). In DL-EPT, electrical properties of a body site at a specific field strength are reconstructed by a convolutional neural network trained on large datasets of simulated B₁⁺ distributions. Very recently, a novel DL-EPT method has been developed for reconstruction of in vivo 3D conductivity maps in the pelvic region [42]. Both the in silico validation and the in vivo comparison with a conventional EPT technique demonstrated that this DL-EPT method is able to reconstruct pelvis conductivity maps with high noise-robustness and accuracy at tissue interfaces.

Sophisticated modeling of the thermal impact of blood flow, urine and cerebrospinal fluid has been developed to calculate more accurate temperature distributions during hyperthermia for various tumor sites. Since the substantial heat exchange between individual large blood vessels and the surrounding tissues is ignored in the Pennes heatsink model, thermal models including discrete vessels have been proposed [43–45]. The discrete vasculature model (DIVA) by Kotte et al. [46,47] efficiently represents 3D patient-specific vessel trees within a 3D thermal model. Raaymakers et al. [48] experimentally validated the DIVA temperature prediction, demonstrating that temperature profiles in the vicinity of a large vessel in the heated area corresponded well to the measured temperature profiles, in contrast to the temperature profiles predicted by the simple Pennes model [33,34,49]. Sumser et al. [49] reported maximum temperature differences of ∼1 °C between DIVA and Pennes models for target volumes close to large vessels. In case unbound body fluids, e.g., urine and cerebrospinal fluid, are present in the heated region, convection will influence the temperature distribution. Schooneveldt et al. developed convective fluid models of the bladder [35,50] and cerebrospinal fluid [51] to accurately represent the fluid dynamics during hyperthermia. In the bladder, for example, more realistic and homogeneous temperatures were predicted by the complete fluid thermodynamic model. Compared to the default model, in which the bladder was modeled as a solid perfused (muscle) tissue, mean deviations between temperature predictions and measurements during treatment reduced from 4 °C to 0.6 °C [35].

Several phase-amplitude optimization techniques have been developed to efficiently calculate the effective phase-amplitude steering settings for modern phased-array systems. These techniques are normally referred to as SAR-based [16,23,52,53] or temperature-based [25,54] optimization strategies, depending on whether their optimization principle is applied to the SAR distribution or the temperature distribution. SAR-based optimization typically maximizes the quotient of the absorbed power in the tumor and a weighted energy norm outside the tumor [22]. Additional weight functions can be defined to account for potential hot spot locations, or regions which are known to have a significantly different perfusion. Nonetheless, temperature hot spots can still occur after SAR-based optimization. Temperature-based optimization allows for direct constraints to normal tissue temperatures, but fast thermal computations are essential since many temperature evaluations need to be performed during optimization. For this purpose Cheng et al. [55] proposed a concept called ‘virtual sources’, which uses a weighted combination of pre-calculated phase-amplitude configurations. Das et al. [56] formulated a method based on superposition of pre-calculated temperature distributions, based on Pennes’ bioheat equation. By extending the concept of temperature superposition to DIVA thermal modeling, Kok et al. [57] enabled efficient temperature-based optimization accounting for discrete vasculature.

Recently, biological modeling has been introduced to quantify the radiosensitization effect of hyperthermia. In this regard, Kok et al. [58] presented a method to predict the effect of hyperthermic radiosensitization in terms of equivalent enhanced radiation dose, i.e., the radiation dose yielding the same biological effect as the combined radiotherapy plus hyperthermia treatment. In this model, the equivalent dose calculation is based on an extended version of the linear-quadratic (LQ) model with temperature-dependent radiosensitivity parameters α and ß reflecting the radio-sensitizing effect of hyperthermia. Literature data on the radiobiological parameters of different tumor and normal tissues at hyperthermic temperatures are scarce, therefore α and ß are normally determined based on few studies reporting the α and β values for some specific tumor cell lines at specific temperatures [59,60]. Myerson et al. [61] and van Leeuwen et al. [62] went a step further: Myerson et al. [61] modeled the heat-induced radiosensitization as changes in the α parameter based on results from in vitro hyperthermia studies; van Leeuwen et al. [62] derived the α and β behavior from mathematical modeling based on biological experiments on cervical cancer cell lines. Furthermore, Franken et al. [60] and Datta and Bodis [63] concluded that the rise in tumor temperature during hyperthermia significantly changes the tumor α/β ratio, based on cell survival experiments and on clinical data. Studies in in vivo tumor models [64,65] also demonstrated that hyperthermia is tumor-selective as it does not cause significant radiosensitization in normal tissues, provided radiotherapy and hyperthermia are given sequentially. Hence, normal tissue α and ß values are assumed to be much less temperature-dependent. In general, results from biological modeling in prostate [58], cervical cancer [66], and recurrent pediatric sarcoma patients [67] demonstrated considerable dose escalations for the combined treatment, highlighting the importance to take the synergistic action of radiotherapy and hyperthermia into account during the treatment planning.

Thus, substantial improvements on various aspects of hyperthermia treatment planning have been realized, but these elements should be integrated into a single HTP workflow to fully exploit the benefits of this research. Hence, in this paper we demonstrate how to integrate the most relevant advanced developments in this field into a single workflow, thereby illustrating an example for advanced patient-specific hyperthermia treatment planning.

HTP workflow

We demonstrate the advanced workflow for hyperthermia treatment planning on a cervical cancer patient. The patient underwent a diagnostic MRI examination, after which she was diagnosed with stage IIIB adenosquamous carcinoma of the cervix. The patient received a combined radiotherapy (RT) + hyperthermia (HT) treatment. An external beam radiation dose of 46 Gy (23 × 2 Gy) was given, followed by a 24 Gy brachytherapy boost using pulse dose rate (PDR) brachytherapy. Five weekly hyperthermia treatments were delivered after an external beam RT session, by using the 70 MHz AMC-4/ALBA-4D locoregional hyperthermia system. According to the clinical protocol, a planning CT scan was made for both external beam radiotherapy and hyperthermia. The external beam dose distribution was calculated using Oncentra External Beam (Elekta AB, Stockholm, Sweden). Hyperthermia treatment planning was performed using Plan2Heat [68].

Hereafter, we outline first the clinical imaging acquisition protocols, then we describe the methodologies used for advanced hyperthermia treatment planning, which include DL-EPT patient-specific tissue conductivity values, modeling of large vasculature, an effective thermal conductivity accounting for convection within the urine in the bladder, phase-amplitude optimization and biological modeling. The advanced hyperthermia treatment planning workflow is also schematically illustrated in Figure 1.

Figure 1. Workflow describing schematically all the methodological steps performed for advanced HTP. In HTP preparation, the patient anatomical model is obtained from CT segmentation. The vasculature, reconstructed from an MR scan and registered to the CT, is also included in this model. The dielectric properties are assigned to each tissue type based on the median DL-EPT values. The thermal properties are assigned from literature values, but an effective thermal conductivity accounting for convective heat transfer in the bladder is assigned. Next, treatment planning is performed, in which the electric field distribution is calculated and phase-amplitude settings for effective tumor heating are optimized. Finally, biological modeling predicts the radiosensitizing effect of hyperthermia in terms of an equivalent enhanced radiation dose.

Imaging protocols

MR

The diagnostic MRI examination was carried out with a 3T clinical scanner (Ingenia, Philips Healthcare, Best, the Netherlands) with a 26-channel torso receive coil. A 3D spoiled gradient echo was acquired with flip angle (FA) = 10°, TR = 5.6 ms, echo train length = 3, SENSE accelerator factor = 1.5 for both phase encoding directions, FOV = 350 × 350 × 98 mm³, resolution = 1.22 × 1.22 × 1.3 mm³, scan time = 1:51 min, and processed according to a generalized two-point Dixon’s method [69]. The resulting water-only image was used for vascular segmentation.

Additionally, a 3D AFI sequence [70] and a 2 D multi-slice spin echo (SE) sequence were performed to reconstruct the amplitude of the B₁⁺ field, |B₁⁺|, and the transceive phase, ϕ^±, respectively, which are necessary for EPT reconstruction. AFI was performed with FA = 65°, TR₁ = 30 ms, TR₂ = 120 ms, TE = 2.5 ms, total scan time = 3:30 min. SE parameters were: imaging FA = 90°, TR = 1200 ms, TE = 6.2 ms. SE was performed twice with opposing readout gradient polarities to acquire an eddy current-free transceive phase map (total scan time = 6:57 min). A FOV = 400 × 280 × 120 mm³ and resolution = 5 × 5 × 7.5 mm³ were used for both AFI and SE.

The duration of the entire MR examination for this patient was ∼50 min. This includes 5–10 min of patient preparation time and 27 min of standard diagnostic MRI. The additional MR-EPT protocol took 10 min.

CT

The planning CT image for external beam RT treatment (RT-CT) was obtained with a Philips Gemini TF 16 PET/CT scanner (Philips, Best, The Netherlands). The RT-CT image was acquired with FOV = 600 × 600 x 363 mm³ and resolution = 1.17 × 1.17 × 3.0 mm³ with the patient in supine position. As part of the clinical protocol, a bowel contrast agent was used during RT-CT scan session in order to allow bowel delineation and sparing during the following radiotherapy treatment planning.

The planning CT image for hyperthermia treatment (HT-CT) was acquired mimicking the hyperthermia treatment setup, i.e., with the patient in supine position on a water bolus and mattresses and with the intravaginal plastic pelotte for temperature monitoring during treatment in situ. The HT-CT image was acquired with a GE LightSpeed RT16 CT scanner (GE Healthcare, USA) with FOV = 650 × 650 × 680 mm³ and resolution = 1.27 × 1.27 × 2.5 mm³.

The durations of the RT-CT and HT-CT examinations were 10–15 min each. For logistic reasons, the HT-CT scan was performed in a separate imaging session after the first hyperthermia treatment. Both the RT-CT and HT-CT could be obtained within a single session; in that case, the HT-CT image should be acquired first, to avoid Hounsfield unit-based tissue segmentation (see next section) being complicated by the presence of bowel contrast.

Preparation of the patient model for HTP

Tissue geometry

The HT-CT scan was semi-automatically segmented based on Hounsfield units with an in-house developed software package, Plan2Heat [68], resulting into fat, muscle, bone and air segmentations. Gross tumor volume (GTV), bladder and the intra-vaginal plastic pelotte were segmented manually (delineation time: ∼20 min).

Vascular structures with diameters ≥∼3 mm were included explicitly in the patient model. External/internal iliac arteries and veins, as well as small blood vessels (with diameter ≥∼3 mm), were manually delineated on the water-only Dixon image, by exploiting the high contrast of vessels with respect to surrounding tissues. Subsequently, these delineations were mapped onto the HT-CT image by deformable registration between the Dixon and CT images. The multi-modal EVolution registration algorithm [71] was used, which maximizes edge alignment between the images being registered. This algorithm requires only contrast between anatomical tissues, preferably present in both image modalities, and assumes deformable/elastic tissues or volumetric tissue changes. The multi-modal EVolution registration algorithm provides a high accuracy in mapping contours between CT and MR images, as demonstrated in a previous abdominal imaging study [72].

The contours of the registered vessels were slightly adjusted on the CT to match the vessel shape in the CT. Furthermore, the vessels registered on the CT only partially covered the patient’s pelvis volume due to the smaller FOV of the water-only Dixon image (∼100 mm in cranio-caudal direction) with respect to the FOV of the HT-CT scan. Thus, the aorta, the renal arteries, the left inferior mesenteric artery, and the common iliac arteries which were identifiable only on the HT-CT image were delineated manually. The common iliac vessels were then manually connected to the registered delineations of the subordinate iliac arteries and veins.

Manual vessel delineation on the MR water-only image takes about 1.5–2.5 h, whereas vessel delineation and adjustment on the HT-CT image may require 3–5 h. Exact times depend on the quality of the image contrast, the image FOV and resolution, the desired minimum diameter of the vascular structures to include in the patient model and the operator’s experience.

Electrical properties of tissues

Recently, Gavazzi et al. [42] have developed a DL-EPT approach for pelvis conductivity mapping where a convolutional neural network infers the conductivity of pelvic tissues from |B₁⁺| and ϕ^± maps measured with MRI, after being trained on in silico dielectric pelvic models with realistic conductivity values at 128 MHz. For implementation details, validation and comparison with a conventional EPT method, the reader is referred to reference [42]. The conductivity is reconstructed in ∼10 s on an NVIDIA Tesla P100 GPU with 16 GB RAM.

Electrical properties were assigned a single value for each tissue type in the segmented patient model. The conductivity of fat, muscle, bone, tumor and bladder corresponded to the median DL-EPT conductivity value calculated in manually delineated 3D regions of interest (ROIs) per tissue. Tissue ROI delineation took approximately 30 min. These conductivity values, measured at 128 MHz (the Larmor frequency of the 3T MR system), were decreased by 3% to obtain conductivity values at 70 MHz (the operating frequency of our hyperthermia system) [40]. The tissue permittivity was taken from literature for two main reasons: first, retrieving the permittivity with EPT at 128 MHz has proven challenging for conventional EPT [38,73], and has not been shown yet in the pelvis with DL; second, even large changes in permittivity (e.g., in muscle) affect very little the temperatures predicted by HTP [40]. The electrical properties of blood were also assigned from literature [74,75]. This is justified, since ex vivo electrical properties of blood can be expected to be rather similar to in vivo values. Electrical properties of tissues are reported in Table 1, along with thermal properties.

Table 1. Electrical and thermal properties of pelvic tissues at 70 MHz.

	Conductivity σ (Sm^–1)^a	Permittivity εr (rel. units)	Density ρ (kg m^–3)	Specific heat capacity ct (J kg^–1 K^–1)	Thermal conductivity kt (W m^–1 K^–1)	Perfusion Wb (kg m^–3 s^–1)
Air	0	1	1.29	10000^b	0.024	0
Bladder	1.22	50	1024	4178	5.6^c	0
Blood	0.83	84.1	1050	3617	0.52	–
Bone	0.12	16.4	1908	1313	0.32	0.12
Fat	0.04	13.4	911	2348	0.21	1.1
Muscle	0.91	70.8	1090	3421	0.49	3.6
Cervical tumor	1.02	70.8	1090	3421	0.49	1.8

Bold font indicates the parameters that were introduced in this study and differing from literature values (in normal font).

^a The values of σ obtained with DL-EPT at 128 MHz were decreased by 3% to determine the values at 70 MHz.

^b The value of c_t for air is ten times higher to permit larger time steps in thermal computations. This has a negligible effect on the steady state temperature.

^c The value of k_t for bladder was ten times higher than the literature value to account for convection phenomena in the bladder.

Thermal properties of tissues

Thermal properties of tissues were derived from the literature [74,76]. A constant perfusion model was used, with enhanced perfusion values to account for the effect of hyperthermic temperatures, based on experimental data [77]. For the bladder a dedicated modeling approach was adopted. Recent studies have shown that modeling fluid convection has a large impact on the final temperature distribution, yielding more accurate hyperthermia treatment planning [35,50,51,78]. However, modeling convective heat transfer is computationally time-consuming at this moment, which makes it impractical for routine clinical use. A practical and time-efficient alternative is to increase the effective thermal conductivity of the fluid. This strategy was shown to approximate the thermal effect of fluid dynamics fairly well [51,78]. We included this strategy to account for convective heat transfer in the urinary bladder, which we modeled to have zero perfusion and thermal conductivity 10 times higher than the literature value [35,51,78].

Hyperthermia treatment planning

Simulation of the electric field

The dielectric patient model was down-sampled to a resolution of 2.5 × 2.5 × 2.5 mm³ based on the winner-takes-all principle [79]. The patient model was then placed within a model of the ALBA-4D hyperthermia system (Medlogix, Rome, Italy) [80], consisting of 1 ring with 4 waveguides operating at 70 MHz (Figure 1). Finite Difference Time Domain (FDTD) simulations of the electric field (E-field) were performed with perfectly matched layer boundary conditions [81]. The E-field was computed separately for each antenna with unit amplitude and zero phase and takes 1–1.5 h per antenna on an NVIDIA GeForce GTX 760 GPU. The exact computation time depends on the computational volume, which is determined by the patient size. The electric field maps were used as input for the optimization of amplitude-phase antenna settings. Superposition can be used to calculate the SAR distributions for specified antenna settings.

Thermal modeling and temperature-based optimization

Thermal modeling was based on Pennes’ bioheat equation [21]: (1) $$c_t{\rho }_t\frac{\partial T}{\partial t}=\nabla ·\left(k_t\nabla T\right)-\omega c_b\left(T-T_{\mathit{art}}\right)+Q+\rho ·\mathit{SAR}$$(1)

In Equation (1)(1) $$c_t{\rho }_t\frac{\partial T}{\partial t}=\nabla ·\left(k_t\nabla T\right)-\omega c_b\left(T-T_{\mathit{art}}\right)+Q+\rho ·\mathit{SAR}$$(1) the temperature change over time in the tissue lies on the left‐hand side, with c_t (J kg⁻¹ K⁻¹) the tissue's specific heat per mass, ρ_t (kg m⁻³) its mass density, T (K) the temperature, and t (s) the time. On the right side of Equation (1) the first term represents the thermal conduction, with k_t (W m⁻¹ K⁻¹) the tissue thermal conductivity; the second term models the perfusion, with c_b (J kg⁻¹ K⁻¹) the specific heat per mass of blood, ω (kg m⁻³ s⁻¹) the volumetric perfusion rate, and T_art the arterial blood temperature; Q (W m⁻³) is the metabolic heat rate, which is negligible during hyperthermia; the last term is the power absorption induced by the heating system, which is proportional to the SAR.

The delineated blood vessels, representing arteries and veins visible on the water-only Dixon and HT-CT images, have a large diameter (≥∼3mm) and thus also a high flow velocity (≥15 cm s⁻¹(82)). This implies that these vessels do not reach thermal equilibrium with surrounding tissue when reaching and traversing the heated region [82]. Therefore, modeling a constant temperature boundary condition is justified and a 37 °C temperature was imposed to the vessel delineations during the thermal computation.

Based on this thermal model, phase-amplitude antenna settings were optimized to realize a favorable 3D steady state temperature distribution for optimal tumor heating. For efficient calculations during optimization, temperature super-positioning was used to calculate the temperature distribution for any arbitrary set of phases and amplitudes by a vector-matrix-vector multiplication [56,57]. The vector represents the antenna settings and the matrix elements are pre-calculated by decomposition of Equation (1)(1) $$c_t{\rho }_t\frac{\partial T}{\partial t}=\nabla ·\left(k_t\nabla T\right)-\omega c_b\left(T-T_{\mathit{art}}\right)+Q+\rho ·\mathit{SAR}$$(1) . Quadratic programing was used to maximize the cost function, which was defined as the tumor T90 (i.e., the temperature exceeded in 90% of the tumor volume), while constraining the maximum temperature for normal tissues to 45 °C. Further constraints were imposed on the antenna contributions to avoid clinically unrealistic power settings: in particular, each antenna was set to deliver at least 10% and at most 40% of the total power. Calculation of the 11 ((N²+N)/2 + 1; N number of antennas, with N = 4) temperature distributions needed for superposition takes about 30 s each, using an NVIDIA GeForce GTX 760 GPU. The actual optimization takes less than 5 min on an Intel Xeon^® E5-1650 v3 3.5 GHz running CentOS 6.8.

Biological modeling

Biological modeling was performed to obtain insight into the radiosensitizing effect induced by hyperthermia, assuming hyperthermia is applied for 1 h once a week, directly after external beam radiotherapy. This effect was expressed in terms of equivalent radiation dose EQD_RT, i.e., the radiation dose yielding the same effect as the combined RT + HT treatment. The EQD_RT calculation is based on an extended linear-quadratic (LQ) model, which includes temperature-dependent radiosensitivity parameters α(T) and β(T) [58,66]. Calculations were performed using the in-house developed software package, X-Term [83], which determines the EQD_RT per voxel from the radiotherapy dose distribution and the temperature distribution. Calculation of the EQD_RT is described more in detail in the following sections.

Preparation

Because the equivalent radiation dose is calculated on a voxel-by-voxel basis from the RT dose distribution and from the temperature distribution, these two maps need to be matched on the same reference frame. Therefore, the RT-CT and HT-CT scans were registered with an intensity-based deformable image registration algorithm embedded in Velocity Medical Solutions (Varian Medical Systems, Palo Alto, CA). First, the bony anatomy was rigidly registered between both scans, followed by a deformable registration as described by Crezee et al. [66]. The resulting deformation vector field was then applied to resample the RT dose distribution onto the HT-CT image. This process takes up to 0.5–1 h, including corrections and fine-tuning of the (structure-guided) deformable image registration.

Equivalent radiation dose calculation

The equivalent radiation dose for RT + HT treatments was calculated instantly, by using the following equation [62]: (2) $${\mathit{EQD}}_{RT}=D·\frac{\alpha \left(T\right)+G·\beta \left(T\right)·D}{\alpha \left(37\right)+d_{\mathit{ref}}·\beta \left(37\right)}+\frac{7.38·{10}^{13}·\left(T+273.15\right)·\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{\Delta S}{2}-\frac{\Delta H}{2·\left(T+273.15\right)}\right)}{\alpha \left(37\right)+d_{\mathit{ref}}·\beta \left(37\right)}$$(2) where α(T) and β(T) are the temperature-dependent radiosensitivity parameters, α(37) (Gy⁻¹) and β(37) (Gy⁻²) are the radiosensitivity parameters at 37 °C, G is the Lea-Catcheside protraction factor accounting for repair during irradiation, D (Gy) is the total dose, d_ref (Gy) is the fraction dose of the RT treatment (i.e., 2 Gy), ΔS (cal/°C/mol) is the entropy of inactivation and ΔH is the inactivation energy of the critical rate-limiting molecules causing cell death.

The temperature-dependent tumor LQ parameters α and β were applied only within the GTV region, whereas α and β were kept at their normothermic value in the surrounding normal tissue regions. Similar equivalent dose calculations could be applied to normal tissue, but here we calculated only the EQD_RT in the tumor region to evaluate the possible effect of the combined treatment. This was motivated by the fact that radiosensitization in normal tissue is much lower compared to the tumor region, not only when hyperthermia is applied within 1–4 h after radiotherapy [64,65], but also for time intervals shorter than 1 h [84].

Similar to [62], the temperature-dependent LQ parameters were obtained from extensive SiHa cervical cancer cell line experiments and increased exponentially with the temperature according to the following expressions$:$ (3) $$\alpha \left(T\right)=\alpha \left(37\right)·\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{T-37}{41-37}·ln\left(\frac{\alpha \left(41\right)}{\alpha \left(37\right)}\right)\right],$$(3) (4) $$\beta \left(T\right)=\beta \left(37\right)·\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{T-37}{41-37}·ln\left(\frac{\beta \left(41\right)}{\beta \left(37\right)}\right)\right]$$(4)

After inserting Equations (3)(3) $$\alpha \left(T\right)=\alpha \left(37\right)·\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{T-37}{41-37}·ln\left(\frac{\alpha \left(41\right)}{\alpha \left(37\right)}\right)\right],$$(3) and (4)(4) $$\beta \left(T\right)=\beta \left(37\right)·\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{T-37}{41-37}·ln\left(\frac{\beta \left(41\right)}{\beta \left(37\right)}\right)\right]$$(4) within Equation (2)(2) $${\mathit{EQD}}_{RT}=D·\frac{\alpha \left(T\right)+G·\beta \left(T\right)·D}{\alpha \left(37\right)+d_{\mathit{ref}}·\beta \left(37\right)}+\frac{7.38·{10}^{13}·\left(T+273.15\right)·\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{\Delta S}{2}-\frac{\Delta H}{2·\left(T+273.15\right)}\right)}{\alpha \left(37\right)+d_{\mathit{ref}}·\beta \left(37\right)}$$(2) and substituting β(37) = α(37)/[α(37)/β(37)], the equivalent radiation dose becomes: (5) $${\mathit{EQD}}_{RT}=D·\frac{\frac{\alpha \left(37\right)}{\beta \left(37\right)}·\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{T-37}{41-37}·ln\left(\frac{\alpha \left(41\right)}{\alpha \left(37\right)}\right)\right]+G·D·\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{T-37}{41-37}·ln\left(\frac{\beta \left(41\right)}{\beta \left(37\right)}\right)\right]}{\frac{\alpha \left(37\right)}{\beta \left(37\right)}+d_{\mathit{ref}}}+\frac{\frac{\alpha \left(37\right)}{\beta \left(37\right)}·\left[7.38·{10}^{13}·\left(T+273.15\right)·\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{\Delta S}{2}-\frac{\Delta H}{2·\left(T+273.15\right)}\right)\right]}{\alpha \left(37\right)·\left(\frac{\alpha \left(37\right)}{\beta \left(37\right)}+d_{\mathit{ref}}\right)}$$(5) where the values for each parameter in Equation (5)(5) $${\mathit{EQD}}_{RT}=D·\frac{\frac{\alpha \left(37\right)}{\beta \left(37\right)}·\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{T-37}{41-37}·ln\left(\frac{\alpha \left(41\right)}{\alpha \left(37\right)}\right)\right]+G·D·\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{T-37}{41-37}·ln\left(\frac{\beta \left(41\right)}{\beta \left(37\right)}\right)\right]}{\frac{\alpha \left(37\right)}{\beta \left(37\right)}+d_{\mathit{ref}}}+\frac{\frac{\alpha \left(37\right)}{\beta \left(37\right)}·\left[7.38·{10}^{13}·\left(T+273.15\right)·\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{\Delta S}{2}-\frac{\Delta H}{2·\left(T+273.15\right)}\right)\right]}{\alpha \left(37\right)·\left(\frac{\alpha \left(37\right)}{\beta \left(37\right)}+d_{\mathit{ref}}\right)}$$(5) are summarized in Table 2. The tumor α/β ratio at 37 °C was 17.9 Gy. A confidence interval for the equivalent dose was also calculated as explained by van Leeuwen et al. [62].

Table 2. Parameters and 95% confidence intervals (within brackets) used in the EQD_RT calculation expressed in Equation (5(5) $${\mathit{EQD}}_{RT}=D·\frac{\frac{\alpha \left(37\right)}{\beta \left(37\right)}·\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{T-37}{41-37}·ln\left(\frac{\alpha \left(41\right)}{\alpha \left(37\right)}\right)\right]+G·D·\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{T-37}{41-37}·ln\left(\frac{\beta \left(41\right)}{\beta \left(37\right)}\right)\right]}{\frac{\alpha \left(37\right)}{\beta \left(37\right)}+d_{\mathit{ref}}}+\frac{\frac{\alpha \left(37\right)}{\beta \left(37\right)}·\left[7.38·{10}^{13}·\left(T+273.15\right)·\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{\Delta S}{2}-\frac{\Delta H}{2·\left(T+273.15\right)}\right)\right]}{\alpha \left(37\right)·\left(\frac{\alpha \left(37\right)}{\beta \left(37\right)}+d_{\mathit{ref}}\right)}$$(5) ).

Parameter (measure unit)	Value
α(37) (Gy^–1)	0.386 (0.364 − 0.409)
α(37)/β(37) (Gy)	17.9 (14.6 − 22.5)
α(41)/α(37) (–)	1.73 (1.64 − 1.83)
β(41)/β(37) (–)	0.41 (0.28 − 0.59)
ΔS (cal/°C/mol)	392.08 (391.62 − 392.46)
ΔH (cal/mol)	147907.8 (NA)^a

These values were obtained from SiHa cervical cancer cell line experiments [62]. ^aΔH has no confidence interval since it is a fixed parameter.

HTP workflow results

HTP preparation results

Figure 2 presents the MR-based reconstruction of the patient-specific conductivity by using DL-EPT. The |B₁⁺| and ϕ^± maps given as input to the trained convolutional neural network are shown along with an MR magnitude image. These maps displayed characteristic patterns in the pelvis at 3T. The reconstructed DL-EPT conductivity map exhibited clear tissue conductivity contrasts, noise-robust estimation and sharp muscle/fat interfaces.

Figure 2. DL-EPT. Transverse and coronal cross-section of MR magnitude image from AFI sequence, |B₁⁺| map, transceive phase (ϕ^±) map, conductivity map obtained with DL-EPT (σ_DL). |B₁⁺| and ϕ^± maps are given as input to the trained network, which infers σ_DL. The tumor contour is shown in red. The tumor was delineated by a radiation oncologist on an ADC map (part of the diagnostic MRI protocol) [85,86]. The tumor delineation was then rigidly transferred to the frame of reference of the AFI image.

Figure 3 illustrates the conductivity histograms obtained within the 3D tissue ROIs. Narrow conductivity histograms were observed in muscle and fat, indicating a precise estimation. Median conductivity values for muscle and fat were 0.933 Sm⁻¹ and 0.044 Sm⁻¹, respectively. Moreover, muscle conductivity agreed with the in vivo mean value found in a previous study where 20 patients were analyzed with Helmholtz-based EPT [87], whereas it was 29% higher than the tabulated value by Gabriel et al. [75]. Fat conductivity was 19% higher than the literature value for non-infiltrated fat [75]. Figure 3 also depicts the conductivity histograms for bladder, bone and tumor. Bone presented a median conductivity of 0.123 Sm⁻¹, a value in between the average literature conductivity values of cortical and cancellous bone (0.067 Sm⁻¹ vs 0.180 Sm⁻¹, respectively) [75]. However, the bone conductivity histogram was skewed toward lower values, possibly reflecting the heterogeneity of bone tissues. The tumor and bladder median conductivities were 1.050 Sm⁻¹ and 1.255 Sm⁻¹, respectively. A similar conductivity value for cervical cancer was found by Balidemaj et al. (1.02 Sm⁻¹) [87]. The small side tails of tumor conductivity histogram were likely caused by inaccuracies due to motion artifacts, which entered the 3D tumor delineation. Analogously, the broader histogram for bladder content could be explained by breathing ghosts and bowel motion artifacts. Hence, the bladder median conductivity was 28% lower than previously reported urine conductivity [87,88], but still close to the range reported for 10 different patients (1.35 Sm⁻¹ to 2.60 Sm⁻¹⁾ [87]. Such variation in bladder electrical conductivity predominantly depends on urine salinity concentration. Note that the median conductivity values obtained with DL-EPT were obtained at 128 MHz (the Larmor frequency at 3T MR proton imaging). These values were converted to values at 70 MHz and then assigned to the corresponding tissue in the patient model (Table 1).

Figure 3. Histograms of DL-based conductivity in pelvic tissues. Conductivity histograms are reported for each tissue along with the median value (black dashed line). Conductivity histograms were based on conductivity values within manually delineated tissue ROIs, depicted on the MR magnitude image (bladder: yellow; bone: blue; fat: magenta; muscle: green; tumor: red).

Figure 4 depicts the vessels reconstructed from the water-only Dixon MR image and the HT-CT scan. Figure 4 clearly shows that iliac arteries and veins have high MR signal intensity with respect to the surrounding muscle in the water-only image, which was therefore exploited for vessel delineation. This benefit of using water-only Dixon scan for vascular imaging was also demonstrated in a previous study, highlighting the better image quality and ∼30% higher vessel-to-background contrast provided by water-only Dixon-based contrast-enhanced MR angiography than a conventional MR angiographic technique [89]. After the multi-modal deformable registration with EVolution registering the MR image to the HT-CT frame of reference, the iliac vessel contours were manually adjusted to accurately fit the vessel shape on the HT-CT image. This is represented in the transverse view of HT-CT image. The coronal views of both MR and CT scans in Figure 4 also reveal the different cranio-caudal dimensions in both scans. Because of this, vessels that were outside the MR image FOV, such as the aorta, the renal arteries and common iliac arteries, were manually delineated directly on the HT-CT image. A 3D rendering of the final reconstructed vessel tree is also reported in Figure 4.

Figure 4. Vessel definition on MR and CT scans. Vessel delineations are shown on the transverse and coronal cross-sections of water image generated with the DIXON method (left) and on the HT-CT scan (right). These delineations were manually adjusted after MR to CT registration (not shown). At the bottom, a 3D rendering of the final vessel tree is shown.

Hyperthermia treatment planning results

Figure 5 shows the SAR and temperature distribution after temperature-based optimization. The optimized phase settings were top/bottom/left/right = 0°/–23°/92°/72°, and the power ratios were top/bottom/left/right = 0.23/0.10/0.28/0.39. The total power absorbed by the patient was 314 W, of which ∼3.9 W were absorbed by the bladder (149 mL volume), whereas the tumor absorbed 1.77 W (82.6 mL volume). This is caused by the higher bladder electrical conductivity with respect to the tumor conductivity. As far as temperature is concerned, the predicted tumor T90, T50 and T10 (i.e., the temperature exceeded in 90%, 50% and 10% of the tumor volume) were 39.0 °C, 39.8 °C and 40.5 °C respectively. The temperature map shows that the average bladder temperature was 42.15 °C, which is 2.35 °C higher than the average tumor temperature, reflecting the negative impact of urine on the achieved tumor thermal dose. Higher tumor temperatures could not be achieved without overheating normal tissue, since hot spots were arising outside the target region, in particular at muscle/bone and/or muscle/fat interfaces (see, for example, close to the ischium in the sagittal view of Figure 5), where tissue temperatures reached the 45 °C constraint. Moreover, the cold tracks induced by the flowing blood in the vessels did not influence the tumor temperature because of their distant location. However, the external iliac artery and vein, as well as the internal iliac artery (visible in the transverse view), contributed to lowering the temperature in the bladder periphery. Furthermore, the external iliac vessels, running down toward the legs (coronal view of Figure 5), affected the temperature of potential hot spots occurring around the pubic bone.

Figure 5. Patient model, SAR and Temperature maps. Transverse, sagittal and coronal cross-sections of the segmented model (first row), showing all the tissues included in the patient model, the SAR distribution (second row) and the temperature distribution (third row) in the patient. The GTV contour is indicated in black.

Biological modeling results

Figure 6 displays the dose distribution related to external beam radiotherapy treatment alone, the temperature distribution and the equivalent radiation dose distribution for the combined treatment (with an assumed fraction dose of 2 Gy, as for RT). The mean RT dose in the GTV was 46.6 Gy. The equivalent radiation dose to the tumor increased on average by 5 Gy when hyperthermia was added, leading to a mean dose value of 51.4 Gy. This 10% increase in dose is attributed to the relatively mild temperature increase achieved in the tumor during HTP. The Dose Volume Histogram (DVH, Figure 7) in the GTV indicates that the D95% (i.e., the dose obtained in 95% of the tumor volume) was 45.45 Gy for radiotherapy alone, and 48.9 Gy for the combined RT + HT treatment, with an uncertainty range of 48.6 − 50.3 Gy reflecting the uncertainty in temperature-dependent LQ parameters. The D0.01 cc (the maximum dose received by at least 0.01 cm³ of the tumor volume) was 53.8 Gy and 60.8 Gy for RT alone and RT + HT combination, respectively. In addition, Figure 7 displays the DVH of the effective biological equivalent dose predicted when a routine HTP workflow was performed (RT + HT_lit, with literature-based electrical and thermal properties [74–76], bladder modeled as perfused muscle tissue, and no inclusion of blood vessels). This literature-based prediction led to a ∼1 Gy higher mean dose with respect to the advanced RT + HT prediction. The D95% and D0.01 cc for RT + HT_lit were 50.0 Gy (uncertainty range: 48.6 − 51.8 Gy) and 58.9 Gy, respectively. The impact of the advanced elements introduced in the proposed HTP workflow is thus also reflected in the predicted equivalent dose. The dose to normal tissue and organs at risk did not change, since radiosensitization in normal tissue was considered negligible [84].

Figure 6. Equivalent dose distributions in radiotherapy and radiotherapy plus hyperthermia. Transverse, sagittal and coronal cross-sections of the radiotherapy (RT) dose distribution (first row), the temperature distribution obtained from hyperthermia treatment planning (second row), and the equivalent dose distribution derived from combining radiotherapy and hyperthermia (third row). The equivalent radiation dose was calculated within the GTV, which is denoted with blue contours. NB: The high-dose region adjacent to the GTV reflects an additional lymph node boost (with mean dose of 54.59 Gy).

Figure 7. Dose Volume Histogram (DVH) reflecting the dose distribution in the GTV for radiotherapy alone (RT, dashed black line) and the combined radiotherapy and hyperthermia treatment as predicted with the advanced HTP workflow (RT + HT, solid black line). For comparison, the DVH predicted using a routine HTP workflow based on literature dielectric properties is also shown (RT + HT_lit, solid magenta line). The shaded area around the DVH for RT + HT and RT + HT_lit corresponds to confidence intervals resulting from uncertainty in LQ parameters.

Discussion

We presented a framework for advanced patient-specific hyperthermia treatment planning (HTP) in which we integrated, for the first time, the clinical HTP workflow with all currently available advanced features of HTP: (i) in vivo MR-based conductivity information, (ii) heat transfer modeling of the vasculature, (iii) modeling of heat convection in urine and (iv) modeling of the hyperthermia treatment effect in terms of equivalent radiation dose. As demonstrated in earlier research papers [34,40,48,49,58,66,83], each of these advanced modeling steps influences HTP predictions and ultimately contributes to reliable patient-specific treatment planning. This makes it important to include all developed advanced features into a single clinical workflow. The presented framework fulfilled this scope without excessive additional computational effort compared to the current state-of-art planning procedure. The E-field calculations of 4–6 h are unchanged and remain the most time-consuming computational part of the workflow. Also the computation time of other parts of the HTP process, e.g., temperature based optimization, did not change. For this advanced workflow, an additional 10 min of MRI scanning time is required for the MR-EPT protocol and the following DL-EPT reconstruction takes only a few seconds. Manual delineation of tissues to generate the dielectric model (i.e., for tissue segmentation on the HT-CT image or for DL-EPT median value calculation) did not exceed 30 min. Manual delineation of blood vessels is more time-consuming (on the order of hours) and varies predominantly with image quality, desired vessel diameter and operator’s experience. These manual delineations could be automated in the future, as discussed later. The workflow was illustrated for locoregional hyperthermia of a cervical cancer patient, but a similar pipeline including patient-specific conductivity reconstruction, vasculature modeling and equivalent dose calculation could also be applied for superficial hyperthermia treatment planning.

Including DL-EPT to obtain the patient-specific conductivity represents an important progress toward reliable hyperthermia treatment planning. We have recently shown that DL-based pelvis conductivity maps display high reconstruction accuracy at tissue interfaces, robustness against noise in measured B₁⁺ maps and acceptable sensitivity to conductivity variations [42]. In contrast to DL-EPT, the present state-of-art method for in vivo pelvis conductivity mapping [38,40,87], i.e., conventional Helmholtz-based EPT (H-EPT), suffers from high sensitivity to noise and so-called ‘boundary errors’ at tissue interfaces, which limit its practical use. Among the numerous EPT approaches [38], inverse methods such as CSI-EPT demonstrated reduced noise sensitivity and good accuracy, also at tissue boundaries, in pelvis conductivity simulations [90,91]; nonetheless, in vivo pelvis conductivity maps have not been reported yet. In principle, electrical properties could be obtained with any EPT reconstruction method, but the quality of conductivity maps demonstrated in vivo by DL-EPT is appealing for direct use as input to HTP, thus eventually eliminating the need for Hounsfield-based thresholding of CT scans to create the patient dielectric model. Nevertheless, direct input to HTP is not yet feasible because bowel peristalsis, breathing and pulsatile flow occurring in the pelvis during the MRI acquisition produce imaging artifacts which affect the |B₁⁺| and ϕ^± maps (Figure 2). At locations where these artifacts appear, DL-EPT retrieves an erroneous conductivity value. To date, no EPT method is capable to counterbalance these motion-induced artifacts. Therefore, we implemented a hybrid alternative where median DL-EPT conductivity values were calculated in tissue 3D ROIs and then assigned to a Hounsfield unit-based segmented model. A deterministic EPT approach as the one suggested by Bevacqua et al. [91], where the tissue segmentation serves as representation basis of the unknown electrical properties, may conceptually be simpler than our hybrid strategy if only median EPT values are desired. In this regard, a convolutional neural network could probably be trained to combine segmentation with learning of an average EPT value for each tissue, for the specific use in HTP. Note, however, that our DL-EPT method was originally developed as an EPT approach aiming at a wide range of clinical scopes [38], with HTP as one of the possible applications. Thus, our network was trained to directly reconstruct the whole 3D electrical conductivity distribution. To aim for full 3D conductivity maps as input for HTP, future work in DL-EPT should focus on compensating motion patterns. Furthermore, an in vivo pelvis permittivity map would also be necessary for fully MR-based HTP. EPT-based in vivo permittivity mapping has been proven more challenging than conductivity mapping in the frequency range of 64–128 MHz, corresponding to clinical MR field strengths [38,73]. However, preliminary results from a recent study [41] have shown that a good quality permittivity map could be reconstructed in the brain at 3T with a DL-EPT method. Further research is required to establish whether DL-EPT could reproduce good quality permittivity maps also in the pelvis.

A previous study [40] showed that the true in vivo electrical conductivity values, which are higher than literature values used in clinical HTP, have a profound impact on the power absorbed by the tumor and its achieved thermal dose with optimized antenna settings. This demonstrates the clinical relevance of using DL-EPT in vivo conductivity values in this workflow (Figure 5). In particular, the tumor SAR was ∼25% lower than the bladder SAR, which was caused by the bladder having a 20% higher electrical conductivity than the tumor. In the present example, the bladder absorbed twice as much power as the tumor, due to its relatively large volume and high conductivity. For a conservative fixed bladder electrical conductivity (1.80 Sm⁻¹), Balidemaj et al. [40] showed that increasing bladder volumes negatively affected the tumor SAR. They further reported that the resulting predicted tumor temperature T90 was diminished by 0.57–1.37 °C from the temperature predicted with the conventional HTP for bladder volumes between 33–497 mL. Considering that similar conditions apply for the patient examined here (augmented bladder electrical conductivity and 149 mL of volume), a negative impact on the tumor temperature of a similar magnitude is expected. Moreover, the in vivo electrical conductivity of the surrounding muscle tissue, which differed by −8% from that of tumor, negatively influenced the tumor SAR. The resulting decrease in tumor temperature can be estimated around 0.40 °C [40].

Another important advantage of using patient-specific conductivity information yielding trustworthy SAR predictions is that the location and magnitude of hot spots will be predicted more accurately, since temperature is correlated with SAR. The electrical conductivity of muscle, being higher than the average literature value, contributed to increase the dielectric contrast at muscle/fat and/or muscle/bone interfaces, therefore aggravating the severity of energy deposition at these locations. Figure 5 illustrates that high SAR values, resulting in potential hot spots, occurred at such tissue interfaces outside the heated region. These predicted arising hot spots, reaching the 45 °C constraint for normal tissue during phase-amplitude optimization, limited a further increase of the tumor temperature.

The effect of convective heat transfer in the urinary bladder was integrated to account for its impact on the tumor temperature distribution, owing to the bladder’s vicinity to the cervical tumor. For a clinically acceptable computation time, fluid heat convection in the urine was accounted for by a higher bladder thermal conductivity, which we set to a value 10 times higher than the literature value. Previous works on bladder HTP showed that using a ten-fold thermal conductivity value approximates the thermal effect of fluid heat convection that can be modeled by a computationally intensive thermophysical model of the bladder [35,78]. As a result of this improved modeling, the mean predicted bladder temperature increases with respect to current clinical HT plans where the bladder is typically modeled as solid perfused tissue. Modeling the bladder as solid tissue was shown to yield errors of 3.8 °C in the median predicted bladder temperature [35,40]. Furthermore, 1 °C temperature differences can be expected also at 1–1.5 cm distance from the bladder [35,50], therefore affecting adjacent tissues such as the cervical tumor and muscle. Considering the realistic bladder electrical conductivity and convection, hyperthermia delivery for cervical cancer patients could be expected to benefit from an empty bladder. A previous study modeling only the electrical in vivo conductivity of urine already suggested this [40]. However, during a 1-h treatment some bladder filling would be unavoidable.

Incorporating a patient-specific network of large blood vessels traversing the heated region in Pennes’ bioheat transfer model represents a further advance in HTP. The Pennes model is commonly used and proved to be valuable but neglects the direction of blood flow and the local heat exchange of vasculature, resulting in deviations in temperature estimations of ∼1–2 °C [34,57]. Explicit modeling of large blood vessels is thus required for accurate temperature predictions. Here we combined manual delineations of blood vessels ≥ ∼3 mm in the pelvic region based on the contrast of water-only Dixon scan with vessel delineations outside the pelvic region directly on the HT-CT image to enable modeling of the largest thermally significant arteries and veins. Further improvements would also include smaller thermally significant blood vessels using the more sophisticated DIVA thermal model [46,47]. This requires a patient-specific discrete vasculature model, knowledge on blood flowrates and directionality and a 3D perfusion map in input. Hence, for future routine applications, additional vessel visualization methods should be acquired, for example angiograms which are typically realized with contrast-enhanced MRI, and (automatic) reconstruction of vessel networks from these angiograms should be incorporated in the clinical workflow [49]. If necessary, available realistic vessel generation software could be used to update incomplete vessel branches automatically [92]. Moreover, flow imaging could be included in the workflow, for example within the MRI protocol, to retrieve functional information on blood flow parameters. Ideally, a quantitative 3D map of the perfusion in the patient would also be required for accurate temperature predictions. It was previously shown that a 25% uncertainty in perfusion parameters of muscle and fat affected the optimized tumor T90 by 0.2–0.4 °C [31]. Furthermore, van den Berg et al. [34] found that the temperature in the prostate was ∼1–2 °C lower when a 3D heterogeneous perfusion was modeled, compared to a homogeneous perfusion model. However, perfusion maps acquired with perfusion CT or MR imaging techniques cannot be obtained during hyperthermia and therefore the increase in perfusion levels at hyperthermic conditions, which can be up to a factor of 10 [93], is not accounted for. Alternatively, thermal decay measurements with standard clinical thermometry probes could be used to estimate the perfusion level [94], but these measurements can only be performed off treatment (i.e., while the power is off, e.g., right after the treatment). The disadvantage of thermal decay measurements is that they provide an indication of the perfusion only at a certain point in time and at the specific locations of these probes (generally within the bladder, the cervix and the rectum).

A temperature-based optimization strategy was applied to predict the antenna settings for optimal tumor heating, while limiting the incidence of treatment-limiting hot spots [54,56,57]. However, the proposed workflow can also be combined with SAR-based optimization, if desired [16,23,24,52,53]. SAR-based optimization techniques generally yield a lower computational cost compared to temperature-based optimization techniques, since the latter typically involve many temperature evaluations for which the Pennes bioheat equation should be solved. Nevertheless, an efficient temperature-based optimization method based on superposition of pre-calculated temperature distributions [56] was used in this workflow. This method allows for temperature-optimization within a few minutes. Here, a temperature-based optimization was preferred as it accounts for physiologically important cooling mechanisms such as perfusion and conduction, as well as bolus cooling. Previous HTP studies in the pelvic region demonstrated that temperature-based optimization is robust to perfusion uncertainty [25] and that simulated temperature changes after phase-amplitude steering correlate with measured temperature changes during HT [95]. Moreover, since tumor temperature and treatment-limiting normal tissue hot spots, which are related to treatment effect, are inherently accounted in the optimization cost function, a temperature-based optimization was also a straightforward choice in light of the subsequent biological modeling.

Biological modeling of the combined radiotherapy and hyperthermia treatment is a beneficial step toward a patient-tailored treatment plan for this treatment combination. We assessed the radiosensitization induced by hyperthermia in terms of equivalent dose by exploiting the extended LQ model [58]. The equivalent dose was calculated for the GTV only, because the radiosensitization in the tumor was shown to be higher than in normal tissue for hyperthermia delivery immediately after radiotherapy [84]. Figure 7 shows that the tumor temperature achieved in this patient escalated the equivalent radiation dose in the GTV by, on average, an extra 5 Gy. Besides, Figure 7 indicates that the advanced modeling steps influencing HTP predictions, as discussed above, also influence the predicted equivalent dose. The hyperthermic radiosensitization effect is reflected in the temperature-dependent α and β parameters. These parameters, obtained from extensive cell line experiments [62], relate to the inhibiting effect of hyperthermia on DNA damage repair and to direct cell kill, but do not account for other mechanisms such as reoxygenation [62,65]. Further (in vivo) research is needed to acquire knowledge on how these physiological effects influence the α and β values [84]. Moreover, LQ parameters for normal tissues should still be investigated, for example by extensive in vivo measurements at different radiation doses, temperature and time intervals [84], especially for organs at risk such as bladder, rectum and bowels. A previous planning study concluded that the hyperthermic radiosensitization in these organs at risk is negligible compared to the tumor region, also because the temperatures were lower [84]. However, a significantly higher bladder temperature is predicted by the advanced HTP. Although in clinical studies for pelvic tumors no enhanced bladder toxicity is observed [96], these higher temperatures might influence the equivalent bladder dose to some extent. In the future, reliable equivalent dose computations in both tumor and normal tissues and prediction of the potential cytotoxicity to organs at risk would become possible with more extensive knowledge on the behavior of LQ parameters under hyperthermic conditions. A next, though very challenging, step in patient-tailored treatment planning would be to directly optimize the equivalent dose distribution of the combined treatment, rather than optimizing the individual treatment modalities separately.

Even though there are still some uncertainties and challenges left in HTP, as discussed above, the advancements achieved with this sophisticated HTP framework will improve the reliability of pretreatment plans compared to the present clinical HTP procedure. Therefore, they are also likely to improve the treatment quality. A previous study demonstrated that, even for clinical pretreatment plans which are only qualitatively reliable, the treatment plan can be adapted on-line in response to patient’s hot spot-induced complaints without significantly diminishing the thermal dose to the tumor [97]. This treatment plan re-optimization is referred to as adaptive HTP. Adaptive HTP can be valuable for both hyperthermia treatments guided by standard thermometry probe feedback and locoregional MR-guided hyperthermia (e.g., using the BSD Sigma Eye MR or the Sigma 30 MR applicator). Instead of thermometry probe feedback, MR-guided hyperthermia uses more detailed MR thermometry information to re-optimize treatment settings. This is feasible for large fixated tumors, as shown for example for treatment of extremity sarcomas [98]. In case of pelvic hyperthermia, nevertheless, MR thermometry information is less reliable due to motion artifacts and standard probe thermometry feedback is still preferable [99]. Ongoing research aims to close the gap toward quantitative temperature monitoring in the pelvic region [100]. In general, it is expected that adaptive HTP will benefit from advanced pretreatment plans, since these plans more accurately predict the magnitude and locations of potential hot spots, which means that the pretreatment plan will be closer to the real optimum. Furthermore, the correlation between the simulated and measured temperature changes during hyperthermia as a result of phase-amplitude steering would be improved [95], which will render adaptive HTP strategies more effective.

Conclusion

A patient-tailored advanced hyperthermia treatment planning workflow was presented, integrating the following novel research elements in the current clinical practice of hyperthermia treatment planning: (i) patient-specific in vivo tissue conductivity reconstructed with DL-EPT from MR measurements; (ii) patient-specific modeling of the thermal influence of large blood vessels; (iii) modeling of the effects of heat convection in urine; (iv) patient-specific assessment of the radiosensitizing effect induced by hyperthermia in equivalent radiotherapy dose. The proposed workflow ensures more accurate predictions of the power deposition and temperature in the patient, as well as a more realistic prognosis of the expected tumor control (through the equivalent radiation dose of the combined treatment) without requiring excessive computational effort.

Acknowledgments

The authors would like to thank Thomas Willigenburg, MD (Radiotherapy Department, University Medical Center Utrecht) for assistance with blood vessel delineation.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The present research was financially supported by Dutch Cancer Society [KWF; project UVA 2014‐7197].