Planning, optimisation and evaluation of hyperthermia treatments

Abstract

Background: Hyperthermia treatment planning using dedicated simulations of power and temperature distributions is very useful to assist in hyperthermia applications. This paper describes an advanced treatment planning software package for a wide variety of applications.

Methods: The in-house developed C++ software package Plan2Heat runs on a Linux operating system. Modules are available to perform electric field and temperature calculations for many heating techniques. The package also contains optimisation routines, post-treatment evaluation tools and a sophisticated thermal model enabling to account for 3D vasculature based on an angiogram or generated artificially using a vessel generation algorithm. The use of the software is illustrated by a simulation of a locoregional hyperthermia treatment for a pancreatic cancer patient and a spherical tumour model heated by interstitial hyperthermia, with detailed 3D vasculature included.

Results: The module-based set-up makes the software flexible and easy to use. The first example demonstrates that treatment planning can help to focus the heating to the tumour. After optimisation, the simulated absorbed power in the tumour increased with 50%. The second example demonstrates the impact of accurately modelling discrete vasculature. Blood at body core temperature entering the heated volume causes relatively cold tracks in the heated volume, where the temperature remains below 40 °C.

Conclusions: A flexible software package for hyperthermia treatment planning has been developed, which can be very useful in many hyperthermia applications. The object-oriented structure of the source code allows relatively easy extension of the software package with additional tools when necessary for future applications.

Keywords:

• Hyperthermia
• hyperthermia treatment planning
• optimisation

Introduction

Clinical hyperthermia, i.e. heating of tumour tissue to 40–43 °C for 1 h to enhance the cytotoxic effect of radiotherapy and/or chemotherapy, is applied to a variety of tumour locations. Examples are recurrent breast cancer, cervix uteri carcinomas, bladder cancer, melanomas and sarcomas [1–5]. Hyperthermia can be applied using several radiofrequency, microwave, water-filtered infra-red-A, ultrasound and capacitive heating techniques, depending on the tumour location. For example, locoregional hyperthermia using phased array systems with one or more rings of antennas surrounding the patient is applied to deep-seated tumours [6–8], whereas superficial hyperthermia is applied to tumours up to 4 cm depth from the skin [9–13]. Intraluminal or intracavitary heating can be applied when the tumour is accessible through an existing lumen or cavity [14,15]. Interstitial heating using an implant of small antennas can be combined effectively with brachytherapy (i.e. interstitial radiotherapy) [16].

Several studies have demonstrated a thermal–dose effect relationship for hyperthermia [17–19], which makes it important to pursue a sufficiently high and homogeneous temperature distribution. The strong heterogeneity of dielectric and thermal properties of human tissues complicates achieving a homogeneous temperature distribution. Especially blood flow causes large heterogeneities, because blood at a temperature close to the body core temperature entering a heated region will produce cold tracks [20,21].

To ensure a good treatment quality, it is important to have insight in the 3D power and temperature distribution and to have treatment planning tools to optimise the treatment. Some research and commercial software packages are available that can be used for simulation of heating patterns for user-defined antennas. Duke University has developed software to simulate SAR distributions for electromagnetic hyperthermia and applied it to predict the SAR for heating with the BSD Sigma60 locoregional system [22]. The University of California, San Francisco has developed dedicated treatment planning software for ultrasound hyperthermia [23]. Commercial packages developed for more general physics/electromagnetic simulations, such as COMSOL (www.comsol.com), Ansys High Frequency Structural Simulator (HFSS; www.ansys.com) and CST STUDIO SUITE (www.cst.com) can also be used for electromagnetic and basic thermal simulations in hyperthermia treatment planning. However, additional software is required for more specific hyperthermia treatment planning purposes. SEMCADX/Sim4Life (SPAEG, Zurich, Switzerland) is a commercially available package that does include additional tools for hyperthermia treatment planning, such as numerical optimisation techniques to optimise system settings for adequate heating. Sigma HyperPlan is a hyperthermia treatment planning system developed by the Konrad Zuse Institute. It is commercially available, though its application is limited to modelling locoregional hyperthermia systems of the BSD Sigma family [24].

Heat removal by blood flow is a very important mechanism influencing the achieved temperature distribution during hyperthermia. Therefore, advanced thermal simulations accounting for the thermal impact of large blood vessels should be included for patient-specific treatment planning [21,25]. The available software packages typically provide thermal simulations where perfusion is modelled as a heat-sink and where boundary conditions can be applied to account for the effect of large blood vessels rather than thermal simulations that include modelling of real 3D curved vasculature. This paper describes Plan2Heat, a dedicated in-house developed software package for hyperthermia treatment planning, supporting advanced thermal modelling using discrete vasculature. The software is voxel-based and applicable to several applications of hyperthermia, such as superficial, locoregional, interstitial and capacitive heating. Its use is illustrated by two examples, one for locoregional hyperthermia with treatment optimisation and one for interstitial hyperthermia including modelling of 3D vasculature.

Methods

The hyperthermia treatment planning software package Plan2Heat is voxel-based and has been developed in C++ and runs on a Linux operating system with an NVIDIA graphics card for efficient parallel implementation of electric field (E-Field) and temperature calculations. Exceptions are the complementary software tools jVolumetool and jTracktool, which have been written in java. JVolumetool is a software application providing a framework for using multiple data sets from different imaging modalities during the delineation phase of radiotherapy treatment planning and has been developed at the University Medical Center Utrecht [26]. jTracktool is comparable to jVolumetool, and is applied to delineate track paths and store the track coordinates.

The data loaded and produced during hyperthermia treatment planning consist of 3D volumes and geometrical objects. All 3D volumes in the treatment planning process are stored as hierarchical data format (HDF; https://www.hdfgroup.org/). The HDF format is very useful for this purpose, as it satisfies the following requirements for flexibility:

• Libraries: A library is available to access the data from C++.

• Geometrical information: It is possible to include geometrical information about a volume in the file.

• Data portability: The data is accessible without conversion when moved to a machine with different architecture. This implies that users do not need to worry about big-endian or little-endian byte ordering.

• Self-describing: Apart from the data itself, additional information can be stored to allow programmes reading the data to understand the structure.

• Inherent data type support: A number of built-in data types are supported, such as byte (char), short, int, float and double.

Structures that need to be modelled or manipulated, such as large blood vessels or catheters, are represented as geometrical objects. These geometrical objects are defined in world coordinates in ASCII files using the Generic Object Format (GOF) [27]. The GOF file is used to stream objects to and from a file or network connection. The geometrical objects establish a spatial relation between a volume and the geometrical objects in that volume. The objects are not discretized beforehand, which makes them independent of the resolution of the calculation grid. In GOF several geometrical structures are supported (e.g. block, sphere, cylinder, cone, tube, line, etc). Geometrical objects are composite objects, which act as a container for their components [28]. This means that they can be composed into tree-structures to represent part-whole hierarchies. The advantage of this composite pattern is that it allows to manipulate individual objects and compositions of objects uniformly. Operations, e.g. rotations and translations, can be applied to the composite object, which applies it to all its components. Especially when handling complex structures, such as vessel networks, a composite pattern is very flexible.

The software package Plan2Heat is a collection of software tools, rather than a monolithic software package, which makes it very flexible in expansion and modification. For each step in the treatment planning process a separate software tool is available, which is called manually from the Linux terminal. Most software tools are called with a parameter file as argument to specify the user input (e.g. the name of the input file and output file, etc.). Each planning process the user performs can be easily structured and documented using noweb literate programming (https://www.cs.tufts.edu/∼nr/noweb/). Additionally, the open-source software package ParaView (http://www.paraview.org) can be used to visualise 3D volumes. A preferred data format for ParaView is the visualization toolkit (VTK) format. For visualisation using ParaView HDF files can be converted to the VTK format using hdf2vtk. Similarly, structured vasculature can be converted for 3D visualisation using vessel2vtk. The process of treatment planning is discussed in the next sections and visualised in Figure 1.

Figure 1. General outline of the treatment planning process.

Preparation of the treatment planning

Create geometry for E-field and absorbed power calculations

To perform treatment planning, first a geometry for E-field and absorbed power calculations should be created. This process is illustrated in Figure 2. For simple phantom simulations, the phantom is defined in world coordinates using GOF. DrawInVolume then rasterises the structure in a resolution-independent manner in a volume defined by the user by specifying the number of voxels, origin and extent.

Figure 2. The process of creating a geometry for E-field and absorbed power calculations for phantom or patient simulations.

For clinical applications of treatment planning, a geometry of the patient anatomy is constructed by segmentation of a CT or MRI data set, preferably in treatment position. Delineations by the physician are stored in a DICOM structure set (e.g. the tumour definition). This imaging and delineation process is the same as for radiotherapy treatments, which allows a flexible integration in the clinical workflow. Next, the DICOM data set is converted to HDF by importing the data set into jVolumetool and saving it as HDF. The delineations are also imported into jVolumetool and stored as GOF. SuppressBackground can be used to remove the background of the scan. This tool identifies the largest object in the data set, such that the table and other objects that are small in comparison to the patient dimensions are automatically removed when a thin air layer is present between the patient and these objects (e.g. by using a towel). Segmentation of different tissues and organs can be done using the tool segment, which applies Hounsfield Unit based segmentation to distinguish muscle, fat, bone, lung and air. Index values are assigned to all segmented tissues, which refer to a MediaTable, i.e. a look-up table that specifies the dielectric and thermal tissue properties for all segmented tissues. As an alternative, or in addition to this automatic segmentation, jVolumetool can be used for manual delineation of specific tissues and organs. Possible artefacts (e.g. by metal objects or catheters) in the data set can also be removed by manual delineation in jVolumetool. These delineations are again stored as GOF and the voxels corresponding to the regions defined by the contours of all manual delineations are assigned a user defined index value using contour and inserted in the patient volume using insertHDF. To replace index values in a specific region when desired, a region growing algorithm, floodFill, is available, which can be applied either in 2D or 3D. AutoCrop can be applied to remove redundant air outside the geometry. For memory or calculation time purposes, the segmented anatomy can be downscaled using the tool downScale. Three different strategies can be selected for downscaling:

1. winner take all, which assigns the most frequently occurring tissue type of all the high resolution voxels corresponding to the new low resolution voxel.

2. averaging, which assigns the average tissue properties of the high resolution voxels corresponding to the new low resolution voxel.

3. anisotropic volumetric averaging, which applies a volumetric averaging in the three orthogonal directions.

After the user has created the complete phantom or patient geometry the antenna(s) need to be defined, again using GOF. For example, when the AMC-4 or AMC-8 system is modelled, rectangular waveguides are defined. Other antenna types, e.g. dipoles or patch antennas, can be modelled analogously. For interstitial heating using local current field electrodes, a dedicated electrode definition is available in GOF [27]. If convenient for positioning of the antenna(s), the origin of the CT/MRI data set can be redefined using setOrigin.

Create geometry for temperature calculations

For thermal simulations, the geometry used for the E-field and power calculations can be cropped around the patient or phantom using autoCrop, if necessary. Boundary conditions can be specified by replacing the indices of the tissue voxels at the boundaries by an index referring to a tissue type in the MediaTable with special properties, typically a fixed temperature or adapted thermal conductivity. This replacement of boundary voxels is done using extendHdf.

Define vessel structures

For advanced thermal simulations, discrete vasculature can be defined. First the vessel tracks need to be segmented using jTracktool. With this tool root points can be defined where the arteries or veins enter the heated volume. By clicking points the corresponding paths and branches of the vessel tracks can be indicated.

In a next step, these vessel tracks are combined into a vessel forest, which is also a GOF definition. This is done using the tool track2vesselForest. Using the specified diameter of each root segment and a user-defined estimation of the flow for each root segment, the diameters and flow are determined automatically for all vessel segments, using Murray’s law. Separate networks for arteries and veins can be reconstructed.

To generate realistic artificial vasculature or to complete incomplete vessel networks, dedicated vessel generation software is available: networkConstruct. The algorithm also provides the possibility to generate counter current networks and the presence of bone and cavities or implants can be accounted for. By adequate selection of construction parameters, physiologically realistic networks can be created with this algorithm [29]. 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. Furthermore, the desired perfusion distribution can be defined. To analyse a created vessel network in terms of branches and vessel lengths etc. of various generations, the tool trackStat is available. This tool assists the user to determine how realistic the vessel network is compared with the known characteristics from the literature [29].

The thermal equilibration in the vessel networks is usually not complete since small vasculature (typically <0.5 mm diameter in case of locoregional heating), where the equilibration takes place, is very numerous and will not be modelled discretely for practical reasons. Thermal equilibration of the blood from the final generation of discrete vessels should therefore be completed by the user by attaching sub-volumes to vessel end points to which a local heat sink is applied [30]. This can be done using the tool addPollGeom. The geometry of these sub-volumes is defined as a GOF structure, e.g. a sphere, a cube or an organ delineation. The sub-volume is automatically attached to all vessel end points in the vessel forest. It is also possible to model a bleed-off along the vessel tree by assigning the sub-volume to junctions where flow is lost to tissue [30]. For arteries, these sub-volumes are used as local sink sets, over which the blood from a terminating artery is distributed, representing the local perfusion. For veins, these sub-volumes are used as local sample sets, and the inflowing blood temperature is the average temperature in the sample set. The user should take the heat balance into account when defining these sink and sample sets: for a correct heat balance the average tissue temperature in the artery sink volumes should be equal to the average tissue temperature in the vein sample volumes. This requirement is automatically met when modelling vasculature in a specific volume or organ and using the complete volume/organ geometry as a sink and sample set.

Define ferromagnetic seed sets

The software can also model heating by ferromagnetic seeds, i.e. cylindrically shaped implants of a ferromagnetic alloy, which are placed inside the tumour [31]. An externally applied oscillating magnetic field induces an electric current within the seeds, which causes a temperature rise. The seeds have a self-regulating control mechanism, determined by the Curie temperature of the ferromagnetic alloy [32]. Because of the conductive nature of this hyperthermia technique, no SAR calculations are required. For treatment planning, a single seed can be defined as a GOF structure specifying the coordinates of the seed track, the diameter, seed density and specific heat, the internal conductivity of the seed and the base power production of the seed in W m⁻³. Additionally, a number of layers between the seed core and the surrounding tissue can be defined to model coatings of the seeds. The inner radius and thermal conductivity of each layer can be specified. The seeds in the implant are then combined into a seed set.

Treatment planning: E-field and power calculations

For radiofrequency or microwave heating, the electromagnetic field distribution is calculated using the Finite Difference Time Domain (FDTD) method [33], implemented to run in parallel on the graphics card (FDTDGPU). The FDTD method solves Maxwell’s equations for a specified operating frequency. A perfectly matched layer is applied to avoid reflections of the electromagnetic fields [34]. In case of multiple antennas, the E-fields can be calculated for each antenna separately with unit amplitude and zero phase. This allows calculation of the power or specific absorption rate (SAR) distribution for any arbitrary phase-amplitude settings by superposition, using the superposition tool batchFadder.

For capacitive heating, the quasi-static approximation of Maxwell’s equations is solved to determine the potential distribution. From this potential distribution the power distribution (P) is calculated using: (1) where σ (S m⁻¹) is the electrical conductivity. Solving the quasi-static approximation of Maxwell’s equations was implemented in quasar. An initial potential or current distribution can be specified. In case of interstitial heating with local current field electrodes, the diameter of the electrodes is typically quite small compared to the spacing of the simulation grid. A correction method using an analytical solution was implemented in quasar to correct the calculated power distribution for the size of the electrode [35]. This feature allows adequate modelling of irregular implants; electrodes do not need to be aligned with the rectangular grid.

Treatment planning: temperature calculations

Thermal simulations implemented in heatran are based on Pennes’ bio heat equation [36]. (2) with c the specific heat capacity (J kg⁻¹ °C⁻¹) and ρ the tissue density (kg m⁻³). The term ∇ ċ (k_tis ∇ T) represents the heat conduction in tissue, with k_tis (W m⁻¹ °C⁻¹) the thermal conductivity. The second term on the right hand side models the perfusion, with c_b the specific heat capacity of blood, W_b (kg m⁻³s⁻¹) the volumetric perfusion rate and T_art the local arterial or body core temperature. During hyperthermia, the perfusion increases with increasing temperature. To account for this increase usually enhanced perfusion values as expected at steady-state hyperthermic temperatures are applied for calculation of steady-state temperatures. Another option is to use a temperature-dependent perfusion model [37]. This yields larger calculation times, but would be particularly useful when transient temperatures are investigated. The power density deposited by the heating system is represented by P (W m⁻³).

When DIscrete Vasculature (DIVA) is included in the thermal simulations, the reconstructed vessel networks are modelled as 3D curves, separate from the tissue grid [38,39]. The advantage of this approach is that it allows for realistic modelling of vasculature, independent of the resolution of the tissue grid. Also closely spaced counter-current vessel pairs can be modelled adequately with this implementation [40]. Vessel segments are discretized into “buckets” with a user-defined sample density (m⁻¹). To handle the tissue–vessel interaction, two sets of voxels are automatically defined for each individual vessel sample: the estimation set and the exchange set. The estimation set voxels are used for calculation of the tissue–vessel interaction. For each voxel in the estimation set, an estimate for heat flow rate density is calculated. The average value is multiplied with the vessel sample surface to yield the heat flow rate. The exchange set voxels are used to effectuate the heat exchange in the neighbouring tissue. The exchange set voxels have their centre inside the bucket or, for small vessels, contain part of the heart line of the bucket and are adjacent to at least one voxel of surrounding tissue. After calculating the temperature distribution voxels belonging to an exchange set can have pseudo temperatures that do not represent the local tissue temperature. Therefore, in the final temperature distribution the values of these voxels are replaced by the wall temperature of the corresponding blood vessels, calculated from the mixing cup blood temperature and the tissue temperature of its surroundings.

When sink/sample sets are included in the simulation the perfusion in these sets should be matched to the blood flow in the terminal branches to which they are connected. This can be done for example using a flow defined approach, in which the blood flow and the volume of the sink set determine the perfusion. A perfusion defined approach allows to use a perfusion map, e.g. derived from MR measurements, to determine the perfusion for each local sink set and subsequently the blood flow in all branches.

Thermal calculations can be performed by a finite difference scheme using small time steps until a user-defined end-time is reached. However, for hyperthermia applications often the steady state temperature distribution is most relevant and can be calculated fast using heatranSST, which solves a linear system of equations by an iterative method [39], either with or without discrete vasculature taken into account.

For interstitial heating using ferromagnetic seeds, the model also applies a separation of the tissue grid and the modelled seeds [31]. The implementation is similar to discrete vasculature. To calculate the interaction of the seeds with the surrounding tissue again estimation sets and exchange sets are defined. Seeds are discretized with a sample density and for each estimation set voxel an estimate of the heat flow is calculated. This heat flow is then divided equally over the exchange set.

Treatment planning: treatment optimisation

When phased array systems are used, the principle of superposition can be applied to combine the electromagnetic fields of each individual antenna. Treatment planning can then be used to predict the best antenna settings for adequate tumour heating by optimising the phases and amplitudes. This phase-amplitude optimisation can be based on SAR or temperature [41–45].

SAR optimisation in our treatment planning software is performed efficiently by solving an eigenvector problem for a linear transformation which is characterised by a positive definite Hermitian matrix for tumour and normal tissue as described by Bardati et al. [41]. So-called influence matrices are generated for target tissue and normal tissue separately by the tool influence. For each tissue type, a different weight factor can be chosen to prevent heating of critical tissues. As a goal function the ratio between the total SAR deposited in the tumour and normal tissue is maximised. Constraints to antennas can be applied that define the minimum or maximum percentage of the total power to be delivered by an antenna. A standard steepest descent method is applied by the optimiser optInf [46].

Temperature-based optimisation requires frequent calculation of the temperature distribution during evaluation of the goal function and constrains, which is computationally intensive. For efficient computations, the method proposed by Das et al. has been implemented [43]. This method uses temperature super-positioning to calculate the temperature distribution for any arbitrary set of phases and amplitudes by a vector-matrix-vector multiplication. The vector represents the amplitudes and phases and the matrix is a complex Hermitian matrix, from which the elements are pre-calculated by decomposition of the bio heat equation. To this end first the absorbed power is decomposed using createPmatrix, then the temperature matrix elements are solved at steady state using heatranComplex for each individual element, or heatranComplexSST. Inclusion of discrete vasculature is supported in heatranComplexSST for advanced optimisation based on DIVA [39]. Constrains are applied to normal tissue temperatures. The temperature in critical structures can be suppressed more actively by using separate constrains to those tissue types. If necessary, even individual voxels can be defined and assigned an individual constraint to suppress hot spot regions. As for the SAR-based optimisation, constraints to antennas can be applied that define the minimum or maximum percentage of the total power to be delivered by an antenna. Different objective functions have been implemented in optimize:

1. Maximise T₉₀ in the tumour,

2. Minimise difference between goal temperature T_G (typically 43 °C) and actual temperature: min ∑_{i ∈ tumour} (max (T_G - T_i, 0))²

3. Maximise tumour temperature, i.e. the sum of the temperature over all tumour voxels:max ∑_{i ∈ tumour} T_i

4. Maximise ratio between the sum over all tumour and normal tissue temperatures:

Other objective functions can be added easily if desired. The actual optimisation is performed using the CFSQP package (C routines for Feasible Sequential Quadratic Programming) [47]. To minimise the risk of finding a local optimum, a number of runs is performed with random initial amplitudes and phases and the solution with the best value for the objective function is then selected as the optimum.

Computation costs are associated with the number of constraints, but in the optimised temperature distribution only a relatively small fraction of the normal tissue voxels reaches the constraint temperature. Therefore, an under-sampling strategy has been developed and implemented to speed up calculations [37]. Rather than including constraints to all normal tissue voxels the grid is under-sampled, such that one voxel is included in the constraint set for every block of i × j × k voxels, with i, j and k user-defined integers. After completion of the optimisation the constraint temperature will be exceeded in some tissue voxels. These voxels are added to the constraint set and the optimisation is restarted. This process is repeated until all normal tissue voxels satisfy the user-defined constrains. Another option to speed up the optimisation is by using element grouping [43,48]. Element grouping uses average temperatures of selected voxel sets instead of single voxels. Elements which achieve their maximum heating potential for approximately the same phase-amplitude setting are grouped. To form groups, eigenvalues and eigenvectors of the temperature matrices are used [43].

Treatment planning: post-treatment evaluation

Some tools are available for retrospective analysis of hyperthermia treatments. hdfHist creates a histogram or cumulative histogram of a distribution and can be used to determine the statistics, e.g. minimum, average and maximum values, overall and per tissue type. This can be used for example to create a temperature volume histogram to analyse the quality of tumour heating. Indexed temperatures T10, T50 and T90, i.e. the temperature at least achieved in 10, 50 and 90% of the volume are calculated, but the user can specify any other percentage of interest to be calculated. The output is written to a standard text file, which can be imported into any programme suitable for analysis (e.g. Excel, gnuplot or Matlab). hdfHist can also be applied to normalise a power distribution, such that a user-defined amount of power is absorbed in specified tissues.

Values of a calculated distribution along specified profiles can be extracted and stored in a standard text file using hdfprofile. Either tracks along a straight line, defined by two coordinates, can be extracted or tracks along paths delineated using jTracktool. This is especially useful to compare simulated SAR or temperatures along catheter tracks with measurements during treatments.

Validation of treatment planning

Validation of treatment planning software is very important to guarantee reliable results of simulations and therefore extensive validation studies have been performed for this software package [35,49–52]. The quasi-static implementation of Maxwell’s equations in quasar has been validated by comparing simulations and analytical solutions for absorbed power in tissue [35]. The implementation of the current source electrodes has been validated by comparing different discretizations of an electrode configuration, as well as by comparing numerically estimated electrode impedances with analytically calculated impedances [35].

Calculations of electromagnetic fields using FDTD have been validated by comparing SAR measurements and simulations for heating homogeneous and inhomogeneous tissue-equivalent phantoms with the 70 MHz locoregional AMC-4 system [51,52]. The differences between measurements and simulations were generally in the order of the measurement accuracy. i.e. 10%. Next, in vivo validation of SAR calculations using FDTD has been performed in 7 cervical cancer patients heated with the AMC-4 system, showing good correlations between measured and simulated SAR profiles along thermometry trajectories [50].

The implementation of vasculature in the DIVA thermal model has been validated by comparing theoretical and simulated equilibration lengths and blood temperature profiles for a vessel inside in a coaxial tissue cylinder [38,53]. Next, experimental validation of DIVA was performed in an isolated perfused bovine tongue, heated with hot water tubes [49]. Comparison of measured and simulated temperature profiles for different controlled perfusion levels showed very good agreement, with deviations of the order of the experimental error [49].

These validation studies show that the software package can accurately perform SAR and temperature predictions when dielectric and thermal tissue properties are known. Further validation is required to determine the clinical value of treatment planning given the present uncertainties in dielectric tissue properties and perfusion. In view of this clinical validation, we evaluated 78 treatment sessions and showed that there is a correlation between measured and simulated changes in average SAR after adapting antenna settings [54]. This indicates that treatment planning can be useful to assist the operator in determining alternative antenna settings to suppress hot spots and increase tumour temperatures based on qualitative information. This clinical validation is part of ongoing research at our department.

Examples

The use of the software will be illustrated using two examples, one for locoregional hyperthermia combined with phase-amplitude optimisation and one for interstitial hyperthermia combined with temperature simulations using DIVA.

Locoregional heating

In this example, a pancreatic cancer patient is treated with locoregional hyperthermia delivered by the 70 MHz AMC-8 waveguide system that can be used as a single-ring or a double-ring system [6]. For this patient a single ring of four antennas was used. Based on a CT scan in treatment position, a segmented anatomy was created on a 2.5 × 2.5 × 2.5 mm³ resolution, using the process described above (Figure 2). The patient had a large seroma close the tumour, which was delineated manually to model convection inside the fluid during temperature calculations by using an enhanced thermal conductivity [55]. The aorta was also delineated manually and kept to a constant temperature of 37.5 °C. Cooling by the water boluses was taken into account by a fixed temperature of 12 °C. Tissue properties were taken from literature and are listed in Table 1. The segmented patient anatomy combined with the antenna model is shown in Figure 3. Electric fields for the individual antennas were calculated using FDTDGPU and the SAR and temperature distribution according to standard treatment settings was calculated. Next, SAR-based and temperature-based optimisation were applied to determine the improvement in heating quality that can be realised using treatment planning.

Figure 3. Transversal and sagittal slices of the CT scan with delineations and the segmented patient anatomy combined with the antenna model for the 70 MHz AMC-8 system, using a single ring of four antennas.

Table 1. Values of the dielectric and thermal properties for different tissue types at 70 MHz, used in the simulations; conductivity (σ [S m⁻¹]), relative permittivity (ɛ_r [–]), density (ρ [kg m⁻³]), specific heat capacity (c [J kg⁻¹°C⁻¹]), thermal conductivity (k [W m⁻¹°C⁻¹]) and perfusion (W_b [kg m⁻³s⁻¹]).

	σ (S m^−1)	ɛr (–)	ρ (kg m^−3)	c (J kg^−1 °C^−1)	k (W m^−1 °C^−1)	Wb (kg m^−3 s^−1)
Air	0	1	1.29	10 000^a	0.024	0
Aorta	0.83	75	1025	3740	0.58	0
Bone	0.05	10	1595	1420	0.65	0.12
Fat	0.06	10	888	2387	0.22	1.1
Lung	0.2	20	750	3000	0.25	3.6
Muscle	0.75	75	1050	3639	0.56	3.6
Seroma	0.75	75	1050	3639	10.0	0
Tumour	0.74	65	1050	3639	0.56	1.8

^a The value of c used for air was ten times too high in order to allow larger time steps in thermal computations. This has a negligible effect on the steady state temperature.

Interstitial heating

This example considers a spherical tumour model (Figure 4) with a diameter of 4 cm and twelve implanted polyethylene catheters containing dual electrode 27 MHz capacitively coupled electrodes [56]. The electrodes in the outer ring were 8 mm long; the central electrodes were 14 mm long and the spacing between the positive and negative electrode was 5 mm; the modelling resolution was 1 mm. The tumour was embedded in a 6 cm wide block of muscle tissue. Outside this block adiabatic boundary conditions were applied for power calculations to model an extending volume. The thermal conductivity outside this block was set to 1/10 of the value for muscle tissue to simulate a fixed temperature boundary condition at a larger distance [57]. Applied tissue properties were literature-based and are shown in Table 2. The power distribution was calculated using quasar, solving the quasi-static form of Maxwell’s equations. Realistic 3D arterial and venous networks were constructed by first manually creating counter-current root segments connected to two child segments in the positive quadrant of the sphere. Next, detailed vasculature was added using networkConstruct with 75 arterial and 75 venous end points of 0.1 mm diameter, taking into account the presence of the catheters. The arterial and venous root diameters were 0.436 mm and 0.438 mm, respectively. Similar networks were defined for the other quadrants of the sphere. The inflow temperature was assumed to be 37 °C. The tumour volume was used as sink/sample set for the vessel end points, ensuring a correct heat balance. Temperature distributions were compared for calculations based on the Pennes bio heat equation (Equation (2)(2) ) and the DIVA thermal model, with a blood flow such that the perfusion in the tumour is the same for both simulations. During interstitial heating blood flow increases substantially and an enhanced tumour perfusion of W_b = 7.5 kg m⁻³s⁻¹ was applied [58].

Figure 4. 3D representation of the modelled tumour volume with a 12 dual-electrode implant, together with the initial discrete vasculature (left) and the complete artificially generated detailed vasculature (right).

Table 2. Values of the dielectric and thermal properties for different tissue types at 27 MHz, used in the simulations; conductivity (σ [S m⁻¹]), relative permittivity (ɛ_r [–]), density (ρ [kg m⁻³]), specific heat capacity (c [J kg⁻¹°C⁻¹]), thermal conductivity (k [W m⁻¹ °C⁻¹]) and perfusion (W_b [kg m⁻³s⁻¹]).

	σ (S m^−1)	ɛr (–)	ρ (kg m^−3)	c (J kg^−1 °C^−1)	k (W m^−1 °C^−1)	Wb (kg m^−3 s^−1)
Muscle	0.71	102	1050	3639	0.56	3.6
Polyethylene	1.2·10^−6	2	920	2300	0.335	–
Tumour	0.6	100	1030	4000	0.6	7.5

Results

Locoregional heating

Based on clinical experience standard phase-amplitude settings for treatment in supine position were selected, yielding phase setting top:bottom:left:right= 0°:40°:20°:20° and power ratios top:bottom:left:right =1.0:0.67:0.8:1.0. The simulated temperature is shown in Figure 5. The power deposition was scaled such that the maximum temperature in the patient did not exceed 45 °C. This resulted in a total power of 258 W absorbed in the patient, of which 1.5 W in the 48.7 cc tumour volume and a T90, T50 and T10 of 40.0 °C, 40.9 °C and 42.7 °C, respectively.

Figure 5. Simulated temperature distributions for heating of a pancreatic cancer patient with the AMC-8 locoregional hyperthermia system using standard clinical phase-amplitude settings, SAR-based optimised phase-amplitude settings and temperature-based optimised phase-amplitude settings.

The heating quality improved after SAR-based optimisation to maximise the SAR ratio between the tumour and the rest of the patient volume. The power in the tumour increased while the total absorbed power decreased. Optimized phase settings were top: bottom: left: right =0°:110°:96°:127° and power ratios were top: bottom: left: right =1.0: 0.05: 0.11: 0.22. The total absorbed power was 185 W, of which 2 W was absorbed in the tumour, yielding a simulated T90, T50 and T10 of 40.4 °C, 42.5 °C and 43.5 °C, respectively. Temperature distributions are shown in Figure 5 and again, the total power was scaled such that the maximum temperature did not exceed 45 °C.

When temperature-based optimisation was performed a goal temperature of 43 °C and normal tissue constraints of 45 °C were defined. The T90, T50 and T10 were 41.3 °C, 42.6 °C and 43.6 °C, respectively (Figure 5). This was achieved with phase settings top:bottom:left:right =0°:31°:124°:144° and power ratios top:bottom:left:right =0.47:0.56:0.42:1.0. The total absorbed power in the patient volume was 527 W, of which 2.34 W reached the tumour; an increase of more than 50% compared to the standard settings. Although the T50 and T10 are fairly comparable to the result after SAR optimisation, the simulated T90 is almost 1 °C higher. This is very relevant since clinical outcome is correlated to the T90 achieved [17,59]. This example demonstrates that treatment planning tools can be useful to focus heating to the tumour.

Adequate tumour coverage is very important in clinical hyperthermia and might be improved by focussing alternatingly onto different parts of the tumour. In this example, the tumour region was sub-divided into three non-overlapping zones and three separate temperature-based optimisations were performed to determine antenna settings that yield maximal heating in each of these zones. To this end, the tumour was split into three different target regions with equal length (2 cm). To assign a separate target index to these regions floodFill was applied, first in 2D to define the target boundaries and then in 3D to replace the index of the target region. When applying the three different optimal antenna settings alternatingly with a short duty cycle (e.g. 30 s per setting) the T90, T50 and T10 became 41.4 °C, 42.8 °C and 43.9 °C. For this example, the gain of using alternating antenna settings is limited because heating with a single setting is already very effective with a T50 of 42.6 °C, but in general such an advanced optimisation can help to improve tumour coverage and reduce treatment limiting hot spots.

Interstitial heating

The characteristics of vessel networks as constructed in Figure 4 can be specified by the branching ratios for the number of branches (RB) and the branch lengths (RL) for each order i, defined as: with N the number of vessels and L the vessel length. Calculations using trackStat showed that for the generations counting at least three branches the mean RB was 3.9 (range 3.8–3.9) and 3.3 (range 2.9–3.8) for the arteries and veins respectively. The mean value of RL was 1.4 (range 1.3–1.5) and 1.4 (range 1.3–1.6) for the arteries and veins respectively. In the literature values of RB between 2.3 and 4.1 are reported and the physiologically feasible range for RL is 1.3–1.7 [60]. This implies that the generated artificial vessel networks are physiologically realistic [29].

After calculation of the power density distribution due to the electrode implant by quasar, the temperature distribution was calculated using Pennes’ bio heat equation, without taking the vasculature into account. The power distribution was scaled such that the maximum tissue temperature did not exceed 45 °C. Indexed temperatures in the tumour volume were T90 = 39.6 °C, T50 = 41.5 °C and T10 = 43.4 °C. The minimum temperature was 38.4 °C. Next, the same power deposition was used to calculate the temperature distribution using the DIVA model accounting for the discrete vasculature. Figure 6 shows opaque isosurfaces of the temperature distribution with and without discrete vasculature, visualising the volume with a temperature with 40 °C or higher. It can be seen that the blood at body core temperature entering the heated volume produces cold tracks, where the temperature remains below 40 °C. This also has some impact on the indexed tumour temperatures: T90 = 40.0 °C, T50 = 41.3 °C and T10 = 43.1 °C. The minimum and maximum temperature in the tumour were 37.5 °C and 44.5 °C.

Figure 6. Opaque iso-temperature surface of the simulated temperature distribution for the 12 dual electrode implant in Figure 4, resulting from the Pennes’ bio heat model and the DIVA model. The picture for the DIVA model also shows the temperature distribution in the vasculature, responsible for the cold tracks where the temperature remains below 40 °C.

Discussion

In this work Plan2Heat, a flexible hyperthermia treatment planning software package, was described, which can be used for a variety of heating techniques. The geometry of the applicator model can be defined by the user by means of GOF structures. This allows to model radiative as well as capacitive heating systems for locoregional, superficial, intraluminal or interstitial hyperthermia. The software also provides the possibility to model ferromagnetic seeds for controlled interstitial heating [31]. This implementation could be extended to modelling of magnetic nanoparticles, a technique that is currently more actively investigated [61–63]. The solver for the quasi-static Maxwell’s equations can be used to simulate capacitive heating methods [64], but is also capable to model the temperature development and distribution during electroporation of tissue; a treatment that is presently investigated by several research groups [65–67].

The applications for which this package can be used vary in scale and the user should select the modelling resolution based on the specific application, such that relevant phenomena as local SAR and temperature peaks can be predicted accurately. For example, as demonstrated in earlier studies, prediction of hot spots during locoregional hyperthermia requires sufficient anatomical detail [48,68]. In the example in this paper a resolution of 2.5 mm was chosen, ensuring both relatively good anatomical detail and suppression of noise after Hounsfield Unit based segmentation. Applications on a smaller scale, such as interstitial heating or electroporation, would require a higher resolution. Currently, simulations are performed on fixed resolution meshes. When necessary, relatively thin or small structures compared with the voxel size can be modelled accurately by implementing additional correction strategies similar to the correction for the small diameter of current source electrodes modelled in quasar [35,69–71]. Implementation of odd-shaped or irregular structures can require more sophisticated techniques. Dedicated solutions for modelling of e.g. helical antennas with FDTD have been published, for example using a weighted equivalent source approximation [72]. Another solution approximates a circular helix with a rectangular helix of wires, which does not require modification of the FDTD source code [73].

The modularity of the software is a great advantage for maintenance and development. On the other hand, it might in principle result in higher demands on user experience. Therefore, for routine simulations that are part of our clinical applications we have structured the consecutive steps and required parameter files in standard noweb files mentioned in the methods section. These noweb files can be executed with minimal user interaction, also limiting the chance of operator errors. This way we have the advantages of modularity without requiring extensive user experience for standard applications.

Thermal modelling is a very important and challenging aspect of hyperthermia treatment planning. There has been significant progress in thermal modelling with discrete vasculature and several models exist that account for straight or semi-curved vessels[74–77], but the possibility to include 3D discrete vasculature networks is a unique feature of the Plan2Heat software package. First steps towards the use of advanced discrete vessel models are taken by companies distributing commercial software, but currently boundary conditions are applied to account for the effect of large blood vessels rather than modelling real 3D curved vasculature. The use of 3D vasculature models in personalised treatment planning requires further research, particularly regarding the required level of detail. When limited vasculature data are available for individual patients, realistic artificially generated vasculature or an atlas-based vessel description can be used, which might have an impact on the reliability of temperature predictions. Future research should reveal down to which detail patient-specific vasculature data are required for personalised treatment planning. For example, a patient-specific description of the largest vessels combined with a more or less generic atlas-based vasculature representation or artificial vasculature could be sufficient for reliable temperature predictions outside the tumour region. Tumours on the other hand usually have very heterogeneous and chaotic vessel networks and the reliability of temperature predictions in the tumour might be limited when using such generalised or artificial vasculature. The example presented in Figure 4 confirmed that the artificially generated vasculature is physiologically realistic by comparing the branching parameters with the ranges reported in the literature. This will be sufficient for realistic modelling of general heat transport by the vasculature in hyperthermia simulations, but a more thorough comparison might be desired for specific applications or organs where the artificial vasculature should represent the real vasculature in detail.

Ongoing progress in research requires continuous updates, adaptations and extensions of treatment planning software. The design using separate software tools, rather than a monolithic software architecture, and the object-oriented structure of the source code of Plan2Heat allows relatively easy extension with additional tools when necessary. Currently, research is performed at our department to create an MR-based patient-specific 3D map of dielectric tissue properties for E-field and power calculations to replace the Hounsfield Unit based segmentation as used in our standard workflow. Hounsfield Unit based segmentation is adequate to distinguish between muscle-like tissue, fatty tissue, bone, lung tissue and air, i.e. tissues with large dielectric contrast, but there is a strong inter and intra-patient variation in dielectric tissue properties, which can amount up to 50% [78]. A patient-specific 3D map of dielectric properties would reduce uncertainties in absorbed power calculations and thus improve the reliability of treatment planning. Future research should also focus on patient-specific perfusion values to further improve simulation accuracy for temperature predictions. Another improvement currently investigated is to add a convective model for the bladder [79]. This is especially relevant when the bladder is the target of heating or located close to the target. The bladder is usually modelled as solid tissue, but during hyperthermia convection will occur in the liquid bladder filling, which has a substantial impact on the temperature distribution [79,80]. Accurate modelling of convection inside the bladder would improve the reliability of treatment planning for bladder cancer patients. These (and other) new features can be included relatively easily in the treatment planning platform because of the object-oriented structure of the software.

Applications of the hyperthermia treatment planning software are numerous, both clinical as well as preclinical. Treatment optimisation to improve clinical outcome is an important application, but because of the uncertainties in tissue properties treatment limiting hot spots might still occur after pre-treatment optimisation. Some research groups explore the possibility to combine treatment planning with MRI feedback [81–83], but this requires MR-compatibility of the heating system. Therefore also on-line adaptive use of treatment planning is currently investigated [84–86]. The feasibility has been demonstrated by a study, which showed that there is a correlation between measured and predicted changes in SAR after phase steering [54]. First clinical applications demonstrated that target coverage improves by using treatment planning during the hyperthermia treatment. Ongoing research to derive patient-specific dielectric tissue properties mentioned above, will improve the quality of pre-treatment optimisation, thereby yielding a better starting point for the treatment. These developments will probably further improve the value of adaptive treatment planning.

Another important application of treatment planning is to compare simulated heating patterns for different hyperthermia systems, either for development of a new prototype or to compare different heating techniques for a specific patient or tumour site [87–94]. For this application, the qualitative reliability of pre-treatment planning can be used to compare heating characteristics of different hyperthermia systems to select the optimal system for tumour heating. Applicator selection for individual patients and prototype development based on treatment planning can thus be very helpful to improve the treatment quality.

Hyperthermia treatment planning presently does not take into account the synergy with radiotherapy and this radiosensitizing effect can be expressed as a local increase in (tumour) dose [95,96]. Ongoing research focuses on integration of radiotherapy and hyperthermia treatment planning to realise thermoradiotherapy treatment planning [97]. In the future, this framework could be extended to model also more complex multi-modality treatments including other combinations with radiotherapy and hyperthermia, such as chemotherapy, halogenated pyrimidines or poly(ADP-ribose) polymerase-1 (PARP1) inhibitors [98,99].

Although in this paper applications were focussed on hyperthermia, the software package is more widely applicable for electromagnetic and thermal simulations, for example in research on high field MRI, RF ablation, mobile phones or RF safety in general [100–107]. With the widespread use of MR investigations and the increasing field strength of MR scanners, safety aspects become very important and the risk of MR-induced hot-spots can be evaluated using simulations [106]. Furthermore, applications of electromagnetic fields increase rapidly and devices should satisfy certain standards for radiofrequency EM radiation. The actual standards are formulated for a maximum tolerable SAR. Advanced computer simulations can help to obtain more insight in these safety issues [100,104,105].

Conclusions

A flexible software package for hyperthermia treatment planning has been developed, which can be very useful in many hyperthermia applications. The object-oriented structure of the source code and the design using separate software tools, rather than a monolithic software architecture allows relatively easy extension of the software package with additional tools when necessary for future applications.

Disclosure statement

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

Additional information

Funding

This work was financially supported by the Dutch Cancer Society (grant UvA 2012-5393).