Optimisation-based thermal treatment planning for catheter-based ultrasound hyperthermia

Abstract

A patient-specific optimisation-based hyperthermia treatment planning program for catheter-based ultrasound technology was developed for a priori evaluation of proposed applicator implant strategies and determination of initial applied power settings. The interstitial and endocavity heating applicators, designed for delivering 3-D controllable hyperthermia within High Dose Rate (HDR) brachytherapy implants, consist of linear and sectored arrays of ultrasound transducers with variable power control in both length and angle. A 3D biothermal model, which incorporates relevant anatomical structures and implant geometries based upon HDR treatment planning, has been developed to simulate the temperature distributions induced by these ultrasound applicators within the catheter implants. A temperature-based constrained optimisation algorithm was devised and integrated within the finite-element thermal solver to determine the optimal applied power levels. A temperature-expressed objective function and constraints were employed to limit maximum temperature (T_max), maximise target coverage (T_target), and minimise thermal exposure to normal tissue and surrounding organs. The optimisation-based treatment planning was applied on representative examples of clinical HDR implants for endocavity treatment of cervix (n = 3) and interstitial treatment of prostate (n = 3). Applicator positioning and orientation, T_max, and T_target, were varied, and temperature volume and thermal dose volume histograms calculated for each plan. The optimisation approach provided optimal applied power levels (4–24 independent transducer sections) leading to conforming or tailored temperature distributions for all cases, as indicated with improved temperature index T₉₀ in the target volume and negligible temperature and thermal dose (t_43,max < 1 min) exposure in surrounding non-targeted tissues, such as bladder and rectum. The precision of the optimised power estimates was shown to be within <5% for a range of starting levels and were similarly convergent. The execution times of this optimisation (<16 min) and forward thermal treatment planning (<22 min) is sufficiently fast to be integrated into the clinical setting. This optimisation-based treatment planning platform for catheter-based ultrasound applicators is a useful tool to provide feedback for applicator selection (sector angle, number of transducer sections along length), positioning (angle or orientation), optimal initial power settings, and has potential to significantly improve the delivery of hyperthermia in conjunction with HDR brachytherapy.

• hyperthermia
• optimisation
• ultrasound
• finite element modelling
• treatment planning

Introduction

There is ample biological and clinical evidence to demonstrate the efficacy of hyperthermia as an adjunct to radiation therapy for the treatment of recurrent and/or locally advanced tumours 1–4. Hyperthermia can be delivered to target sites using external or minimally invasive applicators typically based upon radio-frequency, microwave, or ultrasound energy sources 3, 5. Nevertheless, this combination use of hyperthermia and radiation therapy in clinical practice can be limited by the inherent deficiency of heating devices in controlling the heat delivery, where inadequate temperature exposure in tumour and overheating in surrounding healthy tissue remain a problem 6.

Catheter-based interstitial and endocavity ultrasound applicators are under development as an approach to deliver 3D controllable hyperthermia from within HDR brachytherapy implants 7–11. The key design advantage of these applicators includes the configuration of linear arrays of axially distributed and circularly sectored ultrasound transducers, which may allow for improved energy penetration and 3D spatial control of heating patterns in both length and angle within the target volume. This controllability compares favourably to existing clinical RF and MW interstitial and intrauterine heating systems, which are limited in ability to tailor temperature profiles 5, 12. The multi-electrode RF current source is an exception, which has demonstrated power control along the length of an interstitial implant to improve resultant temperature distributions 13, 14. In order to take advantage of the controllability of catheter-based ultrasound devices in the clinical setting, selection of applicator orientation within the implant and determination of optimised power levels specific to each independently controlled transducer section could be applied to best maximise therapeutic heating and thermal dose coverage to the target region while minimising thermal exposure to surrounding non-targeted tissue.

Hyperthermia treatment planning platforms can be used to guide applicator placement, set-up configurations, and applied power levels prior to delivery of a treatment or for post-procedure analysis of temperature distributions 15. Some examples of clinical hyperthermia treatment planning include the superficial hyperthermia treatment of melanoma on lower leg and exophytic tumour16, intracavitary treatment of oesophageal cancer 17 and cervical carcinoma 18, reperfusion treatment of malignant limb cancer 19, interstitial treatment of brain and prostate 20, and in deep-regional clinical treatment of the patients with cervical, rectal, and prostate carcinoma 21.

Various treatment planning methods tend to pursue a reliable decision-making procedure to avoid the expertise-based intuitive guesswork. The first type of planning is based on forward simulation using finite-difference or finite-element method 21–29. The second type is optimisation, either power-based or temperature-based. The power-based optimisation pursues to optimise the power deposition in the tumour. For instance, the specific absorption rate (SAR) was optimised for phased-array radiofrequency hyperthermia 30–32. The temperature-based optimisation is to optimise the configuration of power input based on the temperature-expressed objective and constraints, i.e., to localise the temperature elevation within the tumour or meet a minimum threshold temperature within the tumour while with constrained temperatures in healthy tissue. For example, the optimisation of amplitude/phase settings of electromagnetic phased-arrays 17, 33, 34 or applicator position and power level settings 35.

In this paper, we developed a patient-specific temperature-based treatment planning platform for the delivery of catheter-based ultrasound hyperthermia which is composed of the following: (1) 3D geometric planning to accurately reconstruct the critical anatomical structures, clinical target volumes, and implant/applicator trajectories, which were acquired from sequential multi-slice 2D images obtained during HDR treatment planning; (2) an anatomy-based biothermal and acoustic model to simulate the dynamics of temperature distribution in tissue induced by ultrasound heating; (3) a temperature-based optimisation algorithm was devised and integrated within the finite element thermal solver to determine a priori the optimal applied power levels and the resulting 3D temperature distribution according to the temperature-expressed objective function and constraints on the reconstructed target and non-target tissues.

The optimisation-based treatment planning system was evaluated using six representative cases from clinical HDR implants for treatment of cervix (n = 3) and prostate (n = 3). The resultant hyperthermia treatment plans were characterised by reconstructed 3D temperature and thermal dose clouds overlaid on 3D anatomy, calculations of thermal dose volume histograms (TDVH), temperature volume histograms (TVH), T₉₀ in target, and a small level of thermal dose t₄₃ in the surrounding tissues. The impacts of different target and normal tissue constraints, variable blood perfusion, stability of convergence with different initial guesses, and computation time were examined.

Materials and methods

Catheter-based ultrasound applicators

Two types of catheter-based ultrasound devices are modelled in this work, following the generalised design schema depicted in Figure 1, with dimensions and operating parameters listed. Typical construct consists of multiple axially distributed and circularly sectored cylindrical piezoceramic ultrasound transducers, with fabrication and performance specifics described in detail elsewhere 7, 9, 11, 36, 37. Each of the transducers can be modified during the applicator fabrication process to emit an angularly directive energy pattern with an inactive zone, or with dual sectoring and each zone active and separately controllable. Interstitial applicators are placed within percutaneous implanted 13g brachytherapy catheters, while the endocavity applicator is placed within the tandem catheter inserted within the uterine cervix. The insertion position of these applicators within the catheters is variable, and can be locked in place once in final position. Cooling water is circulated between the catheter wall and the transducer array during power application. Each applicator configuration can activate up to four 1-cm long transducer sections, with pre-shaped angular heating fields between 90°–360° available and selected a priori. The endocavity applicator, a larger diameter configuration, can accommodate control of four sectors per tube (dual 180°) providing dynamic and independent control in angle and length. The angular directivity of each device is determined and fixed during fabrication (Figure 1). In clinical practice, a family of applicators with different angular heating patterns (e.g., 360°, 270°, and 180°) is available to be selected to best fit the implant geometry and heating requirements. The specific placement of the plastic implant catheters is dictated by the HDR brachytherapy surgical procedure. In practice, the applicators will be inserted and hyperthermia delivered after one of the fractions of HDR brachytherapy, during an interval when the patient is disconnected from the remote afterloader. For the examples shown in this work, the orientation and directionality of multiple applicators was selected for each planned case based upon intuition, and not optimised. Where applicable, directional applicators were selectively placed within peripheral catheters, directing heating energy to best cover the intended treatment zone and avoid non-targeted tissues such as rectum.

Figure 1. Simplified schematic diagram of a catheter-based ultrasound applicator composed of four axial transducers within an implant catheter (R is the radius of the catheter, Re is the radius of the transducer and Le is the length of the transducer). Each transducer can be circularly sectored to produce directional heating patterns; the sector angles applied herein include 360°, 180°, 270°, 2 × 180°, 120°, and 90° (arrows represent the direction of acoustic energy and heat delivery). The larger endocavity applicator can be dual-sectored. The parameters of the interstitial and endocavity applicators applied in this paper are listed.

Patient-specific 3D geometry

Sequential multi-slice CT or MR images of the HDR implant are typically obtained post-surgery and prior to radiation treatment planning for the brachytherapy patient. Following standard HDR treatment planning procedures, the outline of critical anatomical structures including prostate or cervix, bladder and rectum, treatment targets, and catheter trajectories were manually segmented by the physician. The coordinates of these points in 3D space were exported by the radiation planning software and used to reconstruct patient-specific 3D bio-thermal models within the thermal planning programme. The ultrasound applicators are positioned within selected implant catheters and reconstructed during the geometry planning. Figure 2 shows an example of the geometric planning and applicator placements performed for the implants of (a) uterine cervix and (b) prostate. For the purpose of finite element modelling, the structures were automatically meshed using first order adaptive tetrahedron elements to discretise the complicated structures. For solution stability and accuracy, operator-imposed restrictions were used to constrain the minimal element size to ∼0.1 mm within the catheter and adjacent tissue regions, and limited to a maximum element size at the surface of the control volume of 4 mm. These mesh dimensions and choice of first order 3D elements were found by numerical simulations to provide a stable temperature solution as compared to the use of higher order elements and finer mesh.

Figure 2. Examples of geometry planning and anatomy-based 3D thermal modeling of catheter-based ultrasound hyperthermia treatment using (a) an intrauterine applicator for cervix and (b) multiple interstitial applicators for prostate. CT scans (slice thickness 2.5 or 3 mm) were overlaid by segmented contours of rectum, bladder and target.

Biothermal and acoustic model

A 3D anatomy-based thermal model has been developed to simulate the temperature distribution dynamics and thermal dose coverage in tissue heated by an array of multi-element catheter-based interstitial or endocavity ultrasound applicators using the transient bio-heat transfer equation 38:where the transient temperature T(r, t) is a function of r, the three-dimensional position vector, and t, time. The first term on the right side represents the rate of heat conduction in tissue in which k is the thermal conductivity. The second term characterises heat loss due to blood perfusion in which ω₀ is the perfusion rate coefficient. The third term, q, represents the acoustic power deposition produced by an array of catheter-based multi-element interstitial or endocavity ultrasound applicators. The model parameters for generalised soft tissue and the applicators are listed in Table I, cited from 39–44. The initial and Dirichlet boundary conditions were indicated to be the basal temperature T_b = 37°C. The convective boundary condition occurs on the inside catheter surface, where the ultrasound applicators are internally cooled to reduce the local temperature adjacent to the heat source with H_c = 100–300 W cm⁻²°C⁻¹ 45, 46 and T_∞ = 30°C. When applicable for prostate treatment in this study, a urethral cooling catheter is modelled as the endocavity cooling catheter but without the ultrasound power elements.

Table I.  Acoustic and thermal parameters of tissue used for the treatment planning optimisation and simulations.

Parameters	Definition and unit	Value
ρt	Density of tissue (kg m^−3)	1050
ct	Specific heat of tissue (J kg^−1 °C^−1)	3639
k	Thermal conductivity of tissue (W m^−1 °C^−1)	0.56
ωb	Blood perfusion rate (kg m^−3 s^−1)	0.5, 2.5†, 3.5, 5.0
cb	Specific heat of blood (J kg^−1 °C^−1)	3639
μ	Attenuation coefficient, (Np m^−1 MHz^−1)	5.3 (43.9 catheter)
α	Absorption coefficient, (Np m^−1 MHz^−1)	3.5
Hc	Convective cooling coefficient (W m^−2°C^−1)	100, 300
D	Catheter diameter (mm)	2.4 (i), 6.0 (e)
T∞	Temperature of cooling catheter flow (°C)	30
f	Operating frequency (MHz)	6.5, 7.0
δ	Wall thickness of catheter or balloon (mm)	0.5(i), 0.01 (e)

†nominal value.

(i): interstitial applicator; (e): endocavity applicator.

The overall power deposition from the ultrasound applicators was calculated by superimposing the power distribution provided by each single applicator 9, 47, 48, as shown in Equation 2:where α is the absorption rate, μ is the attenuation rate, f is the operating frequency in MHz, is the direction vector of the ith catheter, P_e is the electrical power at the surface of the eth ultrasonic transducer (W cm⁻²), η is the acoustic efficiency of the conversion from electrical power to acoustic power output, θ_i is the sweep angle within the active sector angle of the ith catheter, θ_i,0 is the active sector angle of the ith catheter, N and M is the total number of catheters and transducer elements per catheter. U is the unit step function for setting the active acoustic sector angle, M is the transformation matrix for converting the local coordinate system to global coordinate system. For the test cases evaluated in this paper: the number of applicators was one for the cervix cases, and three or six for the prostate cases; the number of optimal power levels to be determined (active transducer elements) varied from four up to 24 for the cervix and prostate cases, respectively.

Temperature-based power optimisation

A temperature-based optimisation algorithm was devised and integrated to a static anatomy-based thermal modelling to iteratively estimate the optimal applied power levels according to a series of temperature-expressed objective function and constraints that were established on the well-delineated target and non-targeted structures.

For the interstitial or endocavity heating, the objective is to ensure that the temperature on the outer surface of the target volume is above a selected temperature: T_target > = T_targetmin (41°C, 42°C). A series of inequality constraints were considered as follows: (1) The global maximal temperature is controlled below T_max (T_max = 44°C, 46°C) so that the temperature level in target volume tends to be uniform, i.e., T < T_max; (2) The temperature on the surface of the surrounding normal tissue and organs are kept below the therapeutic level to prevent overheating, i.e., T_normal < T_normalmax (T_normalmax = 41.5°C).

To achieve the above objective and constraints, an optimisation algorithm was developed based on an auxiliary function 49, Equation 3, which is composed of objective function and a series of weighted constraints:where the objective function is given by:which is subject to the conditions:where Φ is the objective function on target surface, C₁ is overall maximum temperature constraint, and C₂ is the maximum temperature constraint on normal surface, where N_t and N_nt are the numbers of selected discrete points on the target surface and non-target organs, including rectum and bladder. In order to reduce computation cost, the total number of control points was limited to 50 for the cases evaluated in this work, although this number could be increased if the resultant temperature profiles or solutions are unacceptable. During planning the clinician or physicist would subjectively choose the control points surrounding the edges of the target (N_t) volume, and on the intervening or immediately adjacent surfaces of non-targeted tissues (N_nt) next to the target volume that are a concern, such as the posterior bladder or anterior rectum. The criterion of selection is subjective, and based upon intuition and induction of the approximate temperature profiles and locations of tissues of concern. β_is are the penalty parameters for penalty functions C_is representing the constraints. The penalty function increases as the square of the difference from each constraint and proportional to the weighting factor. The estimates of power levels P_e(i) were updated iteratively to search for the minimal point of the auxiliary function Φ*. The temperature-based optimisation was performed using the quasi-Newton-based iterative method 49.

Implementation

The bioacoustic-thermal model described above was implemented and solved using COMSOL Multiphysics v3.4 (COMSOL, Palo Alto, CA) finite-element software using the PARDISO direct solver 50, 51. Reconstruction of patient anatomy, evaluation of objective functions and constraints, and the quasi-Newton-based optimisation procedure 52 as described above were implemented in the MATLAB (Mathworks, Natick, MA) environment. As a means to reduce the execution time of this treatment planning, the power optimisation was performed using steady-state temperature calculations with a coarsely meshed model (<1e5 elements). The advantage of this strategy is that the computation time can be significantly reduced due to fewer computational elements, because only ordinary differential equations need to be solved for every iterative step, and only the temperature distribution in steady state needs to be considered in the optimisation algorithm. By using this approach, the more complex calculations of transient changes of temperature or thermal dose are not required to perform the power optimisation. The final temperature and thermal dose calculations were then subsequently performed using a full or high resolution (>1e6 elements) transient solution using the previously optimised power levels. A comparison of optimised power values obtained with fine and course models indicated insignificant differences (<5%), validating this approximation as a computational acceleration scheme. The cumulative thermal dose (t₄₃, equivalent minutes at 43°C) distributions were calculated from the transient temperature distributions using the Sapareto-Dewey formulation 53, 54. The planning software ran on a Dell Precision T5400n Workstation, dual quad core Xeon 2.5 GHz processors, 16 GB RAM, with Red Hat Linux operating system (Dell, LaVergne, TN).

The initial values of an optimisation, for example, P_GYN1-GTV,0 is a vector of power intensity levels of transducers of a single endocavity applicator. The elements of P_GYN1-GTV,0 were set as identical small values, i.e., P_GYN1-GTV,0 = [P0(1) = P₀(2) = ··· P₀(4)= const]. In this paper, the constant value, const, was set in the range of 0 < P₀ < 0.5 W cm⁻² to ensure convergence. This upper range of 0.5 W cm⁻² is suitable for generation of hyperthermia temperature distributions with these ultrasound applicators, whereas values approaching or greater than ∼5 W cm⁻² have been utilised for high-temperature ablation 37, 55. In all cases, sufficiently large penalty parameters (β = 106) combined with stringent termination criteria (ε = 10⁻², where ε is the difference between subsequent solutions of the objective function) were chosen by numerical experiments to facilitate the best convergence performance and avoid premature terminations. To evaluate whether the final power levels were readily estimated with acceptable precision, a series of outputs including the value of objective function, iteration number, gradients, and condition number were calculated.

Evaluations on clinical examples

This optimisation-based treatment planning was applied on six representative cases of clinical implants for HDR treatment of cervix and prostate. The hyperthermia targets to be examined included two uterine cervixes (clinical target volume–CTV), one gross tumour volume (GTV) within cervix, two prostates, and one dominant intraprostatic lesion (DIL) in prostate, which were labelled with the acronyms of GYN1, GYN2, GYN1_GTV, PROS1, PROS2 and PROS1_DIL respectively. In each of these cases, 3D geometry planning, anatomy-based thermal modelling and temperature-based optimisation was carried out, following the methods described in the previous sections. The selection of these cases cover a range of anticipated target dimensions, from large volumes in GYN1 to small more easily controlled volumes in GYN1_GTV and PROS1_DIL.

Treatment plans were calculated for each example case under different sets of constraints (T_max/T_min = 46/42°C, 46/41°C, 44/42°C and 44/41°C) and compared. The 3D temperature distributions, 41°C temperature cloud, and thermal dose cloud t₄₃ = 10 min are displayed to evaluate plans. These thresholds for temperature (41°C) and thermal dose clouds are based upon values proven by clinical hyperthermia studies to be efficacious 2, 4, 56, 57. For each optimised treatment plan the thermal dose volume histograms (TDVH) and temperature volume histograms (TVH) were calculated. TDVH is a plot of the thermal dose coverage volume (V) in the target as a function of thermal dose (t₄₃) accrued over the 60 min treatment interval 58. TVH is a plot of the isothermal volume (V) in target as a function of temperature (T) obtained at steady-state. At treatment steady-state, a critical index T₉₀, representing the temperature covering 90% of target volume was indicated.

For further evaluation, this optimisation was tested in the following aspects: (1) to test the impact of different initial guesses, we tested the optimisations on different initial settings of power levels; (2) since it is challenging to measure the precise blood perfusion of patient-specific anatomical structures, a constant and uniform perfusion was assumed, referring to the literature 59. Considering the uncertainty and heterogeneity of perfusion, the deviation of required powers for different perfusion values was calculated for PROS1 case.

Results

Optimised treatment plans for endocavity applicators in cervix

The dual-sectored endocavity ultrasound applicators were positioned within the tandem catheter centred within the target volume as previously described, with the initial applied power P_e to each transducer segment set to 0.25 Wcm⁻². A four-segment applicator was applied for GYN1 and GYN2, each with a large hyperthermia target volume defined from the radiation planning clinical target, which broadly encompasses the uterus and cervix. A two-segment applicator is used for GYN1-GTV, with a smaller hyperthermia target volume defined from the gross tumour volume within the distal cervix. The constraints were set at T_max = 46°C and T_targetmin = 42°C. T_{normal max} = 41.5°C. The optimisation of applied power levels, and forward transient simulations for 1 hour treatment duration were performed. The average and range of optimised power values to be applied for each case are shown in Table II. The resultant 3D temperature distributions, 41°C temperature volume, and thermal dose cloud t₄₃ = 10 min are shown in Figure 3. The TDVH and TVH for the same cases are shown in Figure 4 for four different constraints based on the combinations of the maximal temperature T_max of 44°C or 46°C and the target boundary temperature T_min of 42°C or 41°C. T₉₀ for each set is indicated by an arrow and solid lines. The maximum temperatures and thermal dose exposure to the rectum or bladder is summarised in Table III for each case.

Figure 3. Optimised hyperthermia treatment plans calculated for GYN1, GYN2 and GYN1-GTV cases for intrauterine treatment of cervical cancer: the 3D geometry and anatomical models, temperature distributions, isothermal clouds (T ≥ 41°C) and isothermal dose clouds (t₄₃ ≥ 10 min) are shown, using constraints of T_max = 46°C and T_targetmin = 42°C.

Figure 4. Thermal dose volume histograms (TDVH) and temperature volume histograms (TVH) for optimised treatment plans of intrauterine hyperthermia for cases GYN1, GYN2 and GYN1-GTV. The optimisations were carried out following the constraints as shown in the legend.

Table II.  Summary of final optimised applied power levels for each GYN and prostate treatment plan (T_max = 46°C, T_targetmin = 42°C). Average values and the range from minimum to maximum are reported.

	Applied power (W)
	PROS1	PROS2	PROS1_DIL	GYN1	GYN2	GYN1_GTV
Average ± SD	1.62 ± 0.62	1.45 ± 0.36	1.59 ± 0.48	2.45 ± 0.24	2.26 ± 0.36	1.77 ± 0.38
Range	0.61–3.07	0.77–2.05	1.02–2.36	2.22–2.76	1.72–2.78	1.51–2.33

Table III.  The maximum temperature and maximum thermal dose calculated within the surrounding rectum and bladder by the optimization-based treatment plans.

	Tmax (°C)/(t43,max (min))
	PROS1	PROS2	PROS1_DIL	GYN1	GYN2	GYN1_GTV
Rectum	41.0 (0.9)	41.5 (2.1)	41.0 (2)	37.4 (0)	38.8 (0.2)	37.8 (0)
Bladder	38.5 (0.2)	41.0 (2.3)	37.0 (0)	40.6 (0.3)	41.0 (2)	37.4 (0)

The thermal dose coverage (t₄₃ ≥ 10 min) and temperature profile (T ≥ 41°C) corresponded to the target volume fairly well in all cases (Figure 3). As expected, there were outer target regions in GYN1 beyond the therapeutic range of the ultrasound applicators (>2.5 cm radial) given a maximum temperature threshold of 46°C. In contrast, GYN2 and GYN1_GTV, it was shown that the optimised temperature profile was tailored in the depth and angle direction to conform more closely to the target. The TDVH curves for GYN2 (40.5 cc target volume) and GYN1-GTV (11cc target volume) demonstrated wider plateau than GYN1 (114 cc target), indicating better thermal dose coverage for the smaller targets (Figure 4). The value of T₉₀ is 41°C for GYN-GTV, 40.5°C for GYN2, and 39°C for the much larger GYN1. The T₉₀ could be significantly increased if the maximal temperature constraints could be released from 44°C to 46°C. The maximal width and length dimension of the 41°C temperature cloud for the studied CTVs are shown in Table IV for each case, indicating optimisation of power levels with these catheter-based ultrasound applicators can achieve good therapeutic coverage of intended target volumes.

Table IV.  Summary of the maximal width and length of the 41°C isothermal cloud and corresponding target volumes for the GYN and prostate cases.

	Width (cm)/Length(cm)
	GYN1	GYN2	GYN1_GTV	PROS1	PROS2	PROS1_DIL
41°C cloud	4.6/5.5	3.4/4.5	2.2/3.7	4.5/4.9	2.3/3.2	3.5/4.6
Target volume	4.9/6.0	3.3/4.5	2.0/3.5	4.8/5.0	2.2/3.0	3.3/4.5

Optimised treatment plans for interstitial applicators in prostate

PROS1 and PROS2 were implanted with six interstitial applicators (135°, 180°, 270°, 360° directivity), each with four transducers (n = 24), and directed within and toward the geometric centre of the whole prostate target volume. PROS1_DIL was implanted with three directional applicators in the posterior margin (n = 12, 180° directivity). A transurethral cooling-only applicator was placed within the urethra to simulate a urethral cooling balloon. The optimised treatment plans were calculated for PROS1, PROS2, and PROS1_DIL, with the resultant 3D temperature distributions, 41°C temperature volume, and thermal dose cloud t₄₃ = 10 min shown in Figure 5. The constraints were set at T_max = 46°C and T_min = 42°C. The TDVH and TVH are summarised in Figure 6 for constraints T_max of 44°C or 46°C and the target boundary temperature T_min of 42°C or 41°C.

Figure 5. Optimised hyperthermia treatment plans calculated for PROS1, PROS2 and PROS1-DIL cases for interstitial treatment of prostate target regions: the 3D geometry and anatomical models, temperature distributions, isothermal clouds (T ≥ 41°C) and isothermal dose clouds (t₄₃ ≥ 10 min) are shown, using constraints of T_max = 46°C and T_targetmin = 42°C.

Figure 6. Thermal dose volume histograms (TDVH) and temperature volume histograms (TVH) for optimised treatment plans for optimised interstitial hyperthermia treatment of PROS1, PROS2 and PROS1-DIL cases. The optimisations were carried out following the constraints, in which the T_max and T_targetmin are listed in the legend.

The optimal applied power values (Table II) produced thermal dose coverage (t₄₃ ≥ 10 min) and temperature distribution (T ≥ 41°C) that conformed well to the target volume (Figure 5) for these given implant configurations. The TDVH curves for PROS2 (40 cc) and PROS1-DIL (11 cc) demonstrate a wide plateau with fairly uniform and therapeutic heating throughout the target volumes. The TDVH curves for PROS1 (93.5 cc) target demonstrate therapeutic heating throughout a majority of the target volume. The value of T₉₀ is 43°C for PROS1-DIL, 40.5°C for PROS1, and 42.4°C for PROS2. Similar to the GYN cases, the effect of maximal temperature constraint is more dominant to affect the optimal results than the minimal temperature constraint. The maximum thermal dose t_43,max and maximum temperature T_max for each case are listed in Table IV, and indicate that the thermal exposure to rectum and bladder were negligible. The maximal width and length dimension of the 41°C temperature cloud for the studied CTVs are shown in Table IV for each case, demonstrating temperature coverage throughout the clinical target volumes for these specific catheter implants.

Convergence and stability of optimisation

The uniqueness and stability of the optimisation algorithm was tested by applying four different initial power settings for planning implants of GYN1-GTV and PROS1. The initial applied power intensities were selected as 0.1, 0.2, 0.3 or 0.5 W cm⁻². A boxplot of estimated power levels (P_GYN1-GTV), corresponding to GYN1-GTV, is shown in Figure 7a and demonstrates the deviation of the final estimates due to the difference in initial guesses, with upper and lower sides representing the first quartile (P_0.25) and third quartile (P_0.75) of the estimates. The power estimate of No.1 transducer has a large whisker and inter-quartile region, (2.6, 2.7) W cm⁻², signifying a relatively larger deviation. On the contrary, the estimate of No.2 transducer was more precise with a small inter-quartile region showing no significant difference induced by varying initial guesses. The log-scaled convergence curves for four different initial settings were plotted versus iteration numbers in Figure 7b.

Figure 7. Precision and convergence of optimisation for GYN hyperthermia with different initial guesses. (a) The boxplot of optimised power values of two dual-sectored (2 × 180°) ultrasonic transducers for GYN-GTV, (b) the log-scaled convergence plot for optimisations starting with four different initial settings.

Similarly, the range of optimised final power levels for PROS1 plans are shown in Figure 8a for the six interstitial applicators with 24 transducers total. The ranking of power estimates show that the 7th transducer was of a relatively lower precision while the 8th, 9th, 10th, 12th, 13th, 21st and 22nd transducer have a higher precision. The convergence curves for four different initial settings were shown in Figure 8b.

Figure 8. Precision and convergence of optimisation for prostate hyperthermia with different initial guesses. (a) The boxplot of optimised power values of six interstitial applicators with 24 ultrasonic transducers for PROS1, (b) the convergence plot for optimisations starting with four different initial settings.

To test the sensitivity of optimised power levels induced by a variation in the estimate of blood perfusion, the power levels were optimised at a high perfusion (5.0 kg m⁻³ s⁻¹) and low perfusion (0.5 kg m⁻³ s⁻¹) level and compared to values calculated for the nominal perfusion of 2.5 kg m⁻³ s^−1. The final estimates of power levels were regulated up as the perfusion was larger, and vice versa, as depicted in Figure 9. Of the 24 transducers, the range of power estimates was significant for transducer sections deeply located in the target tissue, such as 45% difference for the transducers of applicator 4 and 5, and the last three transducers of applicator 1. On the contrary, the range of estimates as a function of perfusion were not significant for those transducers located either near the constrained target surfaces or very close to each other such as the first transducer of applicator 1–3, all transducers of applicator 6, and the second transducer of applicator 3. Specifically, the first transducer of applicator 3 is partially out of the target, so the power input was optimised to zero regardless of perfusion variation.

Figure 9. (a) Range of optimised applied power levels for 24 transducer sections on six applicators as calculated for different blood perfusion values (0.5–5 kg m⁻³ s⁻¹) for treatment planning of PROS1; (b) Diagram demonstrating the relative positioning of the applicators, sequentially labelled a1through a6, and the corresponding transducers as indexed in (a). The transducer index of each applicator starts on the proximal end (towards the apex of the prostate) and increases toward the distal end (prostate base), with four transducers per applicator.

With the identical termination criteria, initial guesses, mesh rate, and penalty parameters throughout the cases, the execution time of the optimisation with a coarsely meshed static model and the forward simulation with a finely meshed transient model to simulate a 60 min treatment session are listed in Table V. The fastest computational times were for the GYN_GTV case, with only four sector power levels to be optimised and hence shorter iteration times. In general, computational times increase with the power levels to be optimised or with the number of applicators within the implant requiring additional fine meshing within the thermal model.

Table V.  Summary of the computation time required to generate each optimised treatment plan, distinguishing between time for power optimisation and forward treatment simulation.

	Time (s/case)
	GYN1	GYN2	GYN1_GTV	PROS1	PROS2	PROS1_DIL
Optimisation	180	135	58	902	950	890
Simulation	93	80	23	1280	1320	1002

The sensitivity of the optimisation to the uncertainty or changes in blood perfusion from values used in initial planning was investigated for each of the prostate and GYN test cases. The power inputs specific to each individual transducer or sector were determined by optimising with a basal or nominal perfusion, i.e., ω = 2.5 kg m⁻³ s⁻¹. The forward temperature calculations were then performed for a perfusion decrease (ω = 0.5 kg m⁻³ s⁻¹) and increase (ω = 5 kg m⁻³ s⁻¹) for comparison to temperature distributions at the nominal value. The maximum temperature T_max and T₉₀ for each case and perfusion change is summarised in Table VI. As anticipated, as perfusion increases from the planned value the T₉₀ and T_max are reduced; as perfusion decreases from the planned value the T₉₀ and T_max are larger.

Table VI.  The impact of changes in perfusion from nominal values, as used in the optimisation-based planning, on T_max and T₉₀ calculated within the target volume for the prostate and GYN test cases. Transducer power levels were optimised for a nominal perfusion (ω* = 2.5 kg m ^{− 3} s⁻¹), and applied for an increase (ω = 5 kg m⁻³ s⁻¹) and a decrease (ω = 0.5 kg m⁻³ s⁻¹) in perfusion values for comparative purposes.

	ω* = 2.5(kg m^−3 s^−1)		ω = 0.5(kg m^−3 s^−1)		ω = 5(kg m^−3 s^−1)
	Tmax	T90	Tmax	T90	Tmax	T90
PROS1	46.1	42.4	51	45	43.6	40.5
PROS2	46.2	43.0	51.2	44.5	43.9	41
PROS1-DIL	46.2	42.3	50.4	45.2	43.6	40.8
GYN1	46.2	39.0	49.4	39.3	44	37.8
GYN2	46.2	42.1	49.7	45	44	40.8
GYN1-GTV	46.4	41.4	49.7	43.6	44	40

Discussion

Ultrasound and MECS-RF catheter-based hyperthermia technology can provide the means to tailor or produce more conformal heating patterns when properly applied, typically in conjunction with HDR brachytherapy. The advantages of these devices and their ability for true 3D control are slightly tempered by the added complexity, especially with the requirement of power control and/or selection to a large number of independent sources within either single or multiple catheters. Typical power deposition or SAR-based optimisation schemes which have been developed for deep heating technologies 24 are not suitable for optimising power levels for catheter-based sources, mostly due to the steep radial fall-off of energy and small wavelengths of sources, which necessitate the need for temperature-based optimisation. Towards this end we have developed a hyperthermia treatment planning platform using temperature-based optimisation and constraints for enhancing hyperthermia delivered with interstitial and endocavity ultrasound applicators within an HDR brachytherapy implant.

Results of our study of this planning strategy demonstrate that 3D patient-specific thermal models can be generated from manual reconstruction of 2D image sets typically taken immediately after the surgical implant and prior to the planning of brachytherapy (Figure 2). Critical anatomy in the treatment zone, such as rectum and bladder, as well as the clinical or hyperthermia target volume can be defined within this model based upon clinical insight. Ultrasound applicator configurations, as intended to be inserted within the implant, can be positioned within desired catheters (along catheter length and angle direction if applicable). Temperature constraints can be imposed along the target boundary, maximum temperature thresholds set, and maximum temperature thresholds in boundary tissues prescribed. The temperature-based optimisation with these patient-specific 3D models has been shown to effectively determine power levels which can best tailor resulting temperature distributions and thermal dose coverage to the intended target for a given implant configuration, as demonstrated for multiple cases of both endocavity treatment of cervix (Figures 3–4) and interstitial treatment of prostate (Figures 5–6). By changing the T_max and T_target, better temperature coverage and greater T₉₀ of the target volume could be obtained (Figures 4, 6) for all cases while the maximal thermal dose coverage (t_43,max) in the rectum and bladder was kept at negligible levels as shown in Table II. The effect of maximal temperature constraint is more dominant to affect the optimal results than the minimal temperature constraint for large volume cases such as GYN1 and 2; in this scenario the planned target volume was defined as the full extent of the clinical radiation target volume which is much larger than a typical hyperthermia target volume or gross tumour volume, such as shown in GYN1-GTV.

The optimisation-based treatment planning provides a priori the best estimate of applied power levels for initial settings in delivering hyperthermia with catheter-based ultrasound in the clinical setting. In practice the blood perfusion, acoustic and thermal tissue properties, and real temperature distributions are difficult to predict accurately. Furthermore, the perfusion and tissue properties are often heterogeneous and dynamic during the hyperthermia, with blood perfusion the most dominant and least predictable parameter that often increases with heating. As shown in Table VI, the resultant temperature distributions from using optimised power level settings can change significantly as a function of perfusion changes. However, there is utility in providing an “optimal” setup for the complex interstitial or endocavity procedures following the best assumptions for the parameters. Additionally, optimised plans could be generated over an anticipated range of perfusion variation (Figure 9) and adjusted accordingly. The optimised starting power levels can be applied, and if necessary, the overall power to the applicator(s) can be altered based upon measured temperatures while maintaining the relative power weighting between transducer segments. Although not performed in this study, generating multiple optimised plans with different positioning of applicators within the implant catheters and/or various applicator configurations may be used to empirically determine an optimal applicator implant (position, applicator type) as well as starting power levels to produce tailored heating within each specific patient prior to treatment set-up. Similar analysis and investigations based upon forward simulations with patient-specific anatomy have been used, for example, to optimise patient-specific applicator positions for electromagnetic deep heating 60 or evaluation of microwave applicator designs 61. Although beyond the objectives of this study, future efforts to improve clinical integration and accuracy of this planning system could include the following: incorporation of optimisation techniques to determine the position of applicators and the angle or alignment of the angular sectors 35, implement optimisation algorithms which can better account for uncertainty or variability in tissue properties 62, 63, and inclusion of patient-specific discrete vessel networks and perfusion maps 64.

The sufficient precision and convergence of this optimisation was demonstrated with different initial guesses (Figures 7, 8). When the low initial guess of power input was selected, the temperature will start from basal level and iteratively approach to the optimal points. The robustness of the penalty function against violating the constraints is closely dependent on the penalty parameters, β_i, which determines the relative weight between the objective function and the constraint violations. In this paper, a penalty parameter (β > 10⁶) and terminated criteria (ε < 10⁻²) were numerically tested and chosen as a means to facilitate the best convergence performance and avoid premature terminations.

A maximum of eight transducer sectors (four 10 mm long tubes, dual 180° patterns) for the single endocavity applicator and up to 24 transducer segments within six interstitial applicators (each with four 10 mm tubes of various directionality) were modelled and optimised in this work for treatment of cervix and prostate, respectively. The applicator models were based upon practical catheter-based ultrasound configurations as under development within our group or applied previously in the clinic. The overall computation times required for both the optimisation and forward transient simulation was less than 5 min for the eight element (sectors) endocavity applicator, as determined on our specific computer workstation (Table V). For the prostate cases, the time for power optimisation of 24 transducers was less than 16 min, and time for the forward treatment simulation less than 22 min (Table V). Although the larger number of transducers and increased model complexity (number of nodes and degrees of freedom) required for the interstitial implants significantly increased the computational times, the times remain practical for determining a hyperthermia plan in between the HDR brachytherapy planning session and the associated HDR fraction (e.g., up to 3–4 fractions delivered over the course of two days). The implementation of a coarsely meshed static model for the optimisation component of this planning platform significantly improved computation times and did not significantly degrade the resultant final computed treatment history, as demonstrated for multiple treatment plans (Figures 3–6). If planning times become an issue, the steady-state optimisation can be applied alone for determination of optimal power levels based upon temperature constraints alone, without consideration of the thermal dose.

In summary, components of an optimisation-based treatment planning platform have been developed and shown feasible for planning applied power levels and predicting temperature distributions for multi-element catheter-based ultrasound applicators within an HDR implant. This optimisation-based treatment planning platform can be a useful tool to allow comparison of possible temperature distributions for various applicator selections (sector angle, number of transducer segments) and positioning (angle or orientation) under consideration for a specific patient, and calculate optimal estimates of initial power settings prior to the treatment setup and delivery; the implementation of this information has potential to reduce complexity of treatment delivery and control, and provide a means to deliver more conformal hyperthermia and greater thermal dose coverage.

Acknowledgement

This work was supported by NIH grant R01CA122276.

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