Improved patient-specific hyperthermia planning based on parametrized electromagnetic and thermal models for the SIGMA-30 applicator

Abstract Objective To create an improved planning method for pediatric regional hyperthermia (RHT) using the SIGMA-30 applicator (SIGMA-30). Materials and Methods An electromagnetic model of SIGMA-30 was generated for use with the finite-difference time-domain (FDTD) method. Applying special MATLAB-based algorithms, voxel models of a pediatric patient with pelvic rhabdomyosarcoma were created from Computed-Tomography (CT) contours for use with the FDTD method and the finite-difference (FD) method capable of using either temperature-independent or temperature-dependent perfusion models for solving the Bioheat Transfer Equation (BHTE). Patient models were parametrized regarding, first, the positioning in the applicator, second, the absorbed power range and, third, different perfusion models, resulting in the so-called Parametrized Treatment Models (PTMs). A novel dedicated optimization procedure was developed based on quantitative comparison of numerical calculations against temperature and power measurements from two RHT therapies. Results Using measured data, a realistic absorbed power range in the patient model was estimated. Within this range, several FDTD and BHTE runs were performed and, applying the aforementioned optimization scheme, the best PTMs and perfusion models were identified for each therapy via a retrospective comparison with measurements in 14 temperature sensor positions: 5 in the tumor, 8 in rectum and one in bladder. Conclusion A novel dedicated optimization procedure for identification of suitable patient-specific electromagnetic and thermal models, which can be used for improved patient planning, was developed and evaluated by comparison with treatment-derived measurements using SIGMA-30. The optimization procedure can be extended to other hyperthermia applicators and to other patient types, including adults.


Introduction
Along with increasing significance of regional hyperthermia (RHT) therapy for sarcomas and pediatric tumors [1][2][3][4], there is a growing need for patient-specific therapy planning procedures including detailed models of pediatric RHT applicators like the SIGMA-30 applicator (SIGMA-30). SIGMA-30 (formerly referred to as the Mini Annular Phased Array Applicator -MAPA) is used as part of BSD-2000 Hyperthermia System (Pyrexar Medical Inc., www.pyrexar. com) [5]. It was originally designed for hyperthermia of limbs and was tested with phantoms and amputee limbs [5][6][7].
SIGMA-30 was also used in RHT treatments of pediatric patients [3]. However, no systematic comparisons with measurements have been performed for pediatric patients until now.
During RHT therapies, temperatures must be monitored to control heat delivery to the target and to prevent overheating and damages to surrounding healthy tissues [8][9][10]. For temperature monitoring using temperature sensors, the gold standard is the invasive intratumoral thermometry, but alternative methods such as endocavitary thermometry are also applicable to tumors in the pelvic area [11,12]. On the other hand, development of patient-specific therapy planning has also been the subject of research for many years [13][14][15][16][17][18][19][20]. In the future, accurate temperature predictions could, first, make the need for intratumoral thermometry obsolete and, second, provide a solution for areas less accessible or even inaccessible by endocavitary thermometry like abdominal or extremity areas, which requires a high accuracy of planning methods and models. However, due to the high complexity of the human body and its thermoregulatory system under hyperthermia (e.g., temperaturedependent perfusion changes) special attention is required for creation of perfusion models.
The highly complex perfusion behavior in the human body under hyperthermia has been the subject of numerous investigations both experimental and numerical [19,[21][22][23][24][25][26]. Unfortunately, there is still a lack of not only robust on-line measurement methods under hyperthermia but also of accurate perfusion models which could account for all influence factors, including non-thermal ones, e.g., stress. Therefore, the exact perfusion values in a pediatric patient under hyperthermia are unknown.
Taking these facts into account, our approach, presented in this paper, is to search for the optimal perfusion models via a comparison of numerical calculations with retrospective clinical data, particularly with thermometry and power measurements from two therapies of a pediatric pelvic rhabdomyosarcoma treated with SIGMA-30. To this end, we developed a dedicated procedure based on the optimization of the so-called 'Parametrized Treatment Models' (PTMs), i.e., patient-specific electromagnetic and thermal models with different parameter sets. This optimization procedure can be extended to other patient types, including adults, and to different hyperthermia applicators.

Patient therapy setup
Partially resected pelvic rhabdomyosarcoma in a three-yearold pediatric patient (height 106 cm and weight 17 kg) was treated with hyperthermic chemotherapy under general anesthesia using SIGMA-30 (operating frequency: 123 MHz.) Temperature sensors (Bowman high-impedance thermistors; Pyrexar Medical Inc.) were placed into endocavitary catheters in rectum and bladder, and into an invasive tumor catheter. Temperatures were mapped along the tumor and rectum catheters. Measurements were taken every 1 cm along the catheter. With a tumor diameter of 5.5 cm, this provided a total of 5 measuring points in the tumor. The insertion depth of the rectal catheter was 7 cm, resulting in a total of 8 measuring points along the rectal catheter, including the position '0' at the anus. In contrast, due to the cuff fixation, the temperature sensor in the bladder catheter was fixated at one measuring point, 1 cm above the bladder neck. Consequently, temperatures in a total of N cath ¼14 measuring points in the pelvic region were recorded. We refer to these 14 catheter values as a (measured) Temperature Set. Moreover, several vital functions were monitored during anesthesia, including heart rate, blood pressure, oxygen saturation and sublingual temperature.
A total of eight RHT therapies were performed, but only the data of the 3rd and 4th therapy are discussed in this paper. During the 1st and the 2nd therapy, for preventing postoperative infections, a thick gauze bandage was used for the fixation of the tumor catheter, which locally heavily impeded the heat transfer through the skin into the cooling water bolus. For consideration of the heat accumulation underneath the gauze bandage a much more complicated thermal model is needed than that available for calculations performed in this paper. Therefore, the 1st and the 2nd RHT therapy are not discussed in this paper. The gauze bandage was removed before the 3rd therapy. Furthermore, due to therapy design, the tumor catheter was removed after the 4th therapy. For these reasons, the presentation of data is restricted only to the 3rd and the 4th therapy. As the 4th therapy exhibited higher power (and thus higher temperatures) than the 3rd therapy (see Results section), we refer in this paper to the 4th therapy as 'Higher-Power'-Therapy (HP_THER), and to the 3rd therapy as 'Lower-Power'-Therapy (LP_THER).
The patient is registered in the HT-Registry Study, which was approved by the Charit e's Ethics Committee (approval no.: EA2/114/14) with written informed consent obtained from the patient's legal guardian.
Overview of the measurement-based optimization procedure Figure 1 displays the general workflow of optimization procedure for improved patient-specific RHT planning. First, an electromagnetic model of SIGMA-30 was created for use with the finite-difference time-domain (FDTD) method [27][28][29][30][31][32][33][34]. A Perfectly Matched Layer (PML) formulation for dealing with open boundaries was used [30]. Second, applying our MATLAB-based (MATLAB R2016b) in-house segmentation algorithm [35], Computed-Tomography (CT) -based voxel patient models were generated for use with our in-house FDTD method [31][32][33][34] and with our in-house finite-difference (FD) method [36], the latter being capable of using either temperature-independent or temperature-dependent ('dynamic') perfusion models for solving the Bioheat Transfer Equation (BHTE), introduced by Pennes in 1948 [37]. Next, patient models were parametrized related to the positioning in the applicator ('FDTD-related parametrization') and several FDTD runs were performed, yielding relative SAR distributions. For each FDTD run, a certain 'channel setting' (i.e., sets of power amplitudes and phases in amplifier channels connected to the applicator antennas) was applied for simulating the real channel setting used for a therapy. Generally, the determination or choice of a particular channel setting is not part of our optimization procedure, but it acts as an 'external' input for both the therapy and the calculations (see the box 'channel setting' in Figure 1). As a next point, a realistic range of absorbed power in the patient model was derived from measurements, and due to these power values absolute SAR distributions were generated, necessary as inputs for BHTE calculations. On the other hand, BHTErelated parametrizations were performed, resulting in several perfusion models of a patient as other inputs for the BHTE. Combinations of all parametrizations (therapy, positioning, power, perfusion models) are referred to as PTMs. The total number of PTMs (which is identical to the total number of BHTE runs) is indicated in Figure 1 as 'n'. For each PTM a BHTE run was performed with an absolute temperature distribution as output, containing as a part of the solution also N cath ¼14 temperature values at the points along the catheters in locations which were defined by measurements (Temperature Set). In the next step, the best PTMs were found via a root-mean-square-deviation (RMSD)-based minimization using a comparison between the BHTE-calculated and the measured Temperature Sets (from therapies HP_THER and/or LP_THER). Finally, the best perfusion models were identified via dedicated evaluation approaches. The core parts of the optimization routine that are bounded by dashed lines in Figure 1 are described in more detail in Figure 2 in the paragraph Two-Step Optimization, further below. The computer requirements and time to perform the optimization procedure are discussed in the Supplementary material, in the paragraph 'Supplement Part 5'.
In the following, the main steps of the optimization procedure are presented in more detail.

Applicator and patient models
For performing FDTD calculations the two necessary initial steps are, first, creation of a SIGMA-30 model, and, second, of a CT-based patient model (Figure 1). SIGMA-30 is the smallest applicator of the SIGMA series of applicators (SIGMA-60, SIGMA-Eye, SIGMA-40, and SIGMA-30, Pyrexar Medical Inc.) used in the Department of Radiation Oncology at Charit e Campus Virchow Klinikum. It consists of eight flat copper-strip dipole antennas, connected via tuning/matching networks in pairs to four power amplifier channels. Electromagnetic models of applicators SIGMA-60 and SIGMA-Eye, used in our department, were previously created based on the FDTD method. These applicator models were validated using different methods including lamp matrix comparisons and temperature rise measurements in various homogeneous and inhomogeneous phantoms, as well as Magnetic-Resonance (MR)-thermography measurements [31,[38][39][40][41][42]. As a further development of these models, an FDTD model of SIGMA-30 was generated for comparisons including the pediatric case presented in this paper. Before these comparisons the SIGMA-30 model had been tested using a cylindrical agar phantom developed for this purpose (see 'Supplement Part 1').
For the creation of the patient model planning CT data of the patient, which included the postoperatively inserted tumor catheter, were acquired. Various structures in the CT data sets (tumor, organ-at-risk, etc.) were manually contoured using the Varian Treatment Planning System (Eclipse Version 11.0). A novel segmentation procedure was developed based on inhouse created MATLAB algorithms [35]. The segmentation output results in an optimal combination of clinically/manually contoured regions and HU-threshold-segmented portions. Figure 1. General workflow of the novel optimization procedure for improved regional hyperthermia (RHT) planning using optimization of the so-called 'Parametrized Treatment Models' (PTMs) via a root-mean-square-deviation (RMSD)-based comparison between calculated and measured temperatures along tumor, rectal and bladder catheters. An oval form of boxes refers to measurement parts of the optimization scheme. The core parts of the optimization routine that are bounded by dashed lines are described in more detail in Figure 2.
Further details concerning the model creation see 'Supplement Part 2' and 'Supplement Part 5'.

FDTD-related parametrization
In the optimization scheme in Figure 1, immediately before the step 'FDTD runs', a step called 'FDTD-related parametrization' is performed. In this step the user may define/ change numerical parameters that have impact on the relative SAR-distribution in the patient model. The total number of relative SAR distributions resulting from this kind of parametrization is indicated in Figure 1 as 'n_SAR'.
Numerous FDTD-relevant effects may be investigated like modifications in the model geometry and resolution, changes in material and electromagnetic properties, but also the use of different antenna or feeding models, or other effects like smoothing of SAR or fluctuations of channel settings (particular amplitude and phase values generated by amplifiers that are connected to the antenna feedpoints). In our paper, however, the focus is on the investigation of the effects of different perfusion models, resulting from the BHTE-related parametrization (see further below), and therefore we restrict the FDTD-related parametrization to the investigations of the influence of the adjustment of the patient model in the applicator model.
Concerning the longitudinal position the patient was positioned so that the central applicator transversal symmetry plane was situated in the central tumor plane, 7 cm caudally from the navel. Possible changes of longitudinal position are not investigated in this paper.
Regarding the transversal position, for avoiding any superficial or subcutaneous hotspots in the anesthetized patient, the default clinical patient positioning for both therapies was 'central' (as defined below) and a 'central channel setting' was set, too (equal power amplitudes and no phase delays in power channels). Apart from health hazard issues the use of the central channel setting was also reasonable because the tumor region was situated almost centrally. For investigating how the final solution may depend on patient positioning, we introduced the following parametrization related to transversal patient positioning: 'Central Patient Positioning' (CE_POS): the center of the patient model was aligned in the center of the applicator so that in the transversal symmetry plane the horizontal distances between patient surface and right/left applicator boundaries were equal (%5 cm) and the vertical distances between patient surface and top/bottom applicator walls were also equal (%10 cm). The shape of the pediatric patient in the transversal applicator symmetry plane was nearly elliptical with a major axis of 19 cm and a minor axis of 9 cm. 'Shifted Patient Positioning' (SH_POS): the patient model was shifted by a distance 'd' in the dorsal direction. It is often observed in pediatric patients during therapy in SIGMA-30 that the weight of the patient sinks into the special pediatric-size hammock-like couch, which is made of pliant fabric. Possible effects of such a shift were studied by using this positioning. In accordance with the small size of the pediatric patient, we investigated a dorsal shift of d ¼ 10 mm.
In conclusion, this particular FDTD-related parametrization results in only two different relative SAR distributions, i.e., n_SAR ¼ 2.

FDTD runs
The purpose of this step was SAR calculation via solving Maxwell's equations by using the FDTD method. For this, FDTD requires a CT-based voxel model of a patient and a model of SIGMA-30 as inputs ( Figure 1). Depending on patient positioning this resulted in two geometry models (CE_POS/SH_POS) which were downsampled to uniform 2.5 mm-grids representing electrical material properties. These grids were used for the staggered FDTD (Yee) lattice with 5 mm distance between E-field components thus providing a clear assignment of the electrical material properties to every E-field component. For each of so generated two FDTD models corresponding to patient positions CE_POS and SH_POS, respectively, a relative SAR distribution was calculated, simulating the central channel setting used for both therapies. For further details concerning FDTD runs see 'Supplement Part 2' and 'Supplement Part 5'.

Power range
As indicated in Figure 1, the input for every BHTE run was an absolute SAR distribution that can be obtained from the aforementioned relative SAR distribution via multiplication ('scaling') of its values in every point of the model grid with a (scalar) power scaling factor f ¼ P/P rel . P rel is the absorbed power resulting from integration of the relative SAR in the patient model. P is a desired absolute power value. Note that a given relative SAR distribution can be scaled with different scaling factors, f i ¼P i /P rel , corresponding to different power values P i , which can belong to a certain power range like that defined in Equation (1), below. The volume of integration was given by the FDTD model and was equal for this particular pediatric patient to $6.85 liters and the corresponding model mass was $7.13 kg (remember that the total mass of the child was 17 kg).
A crucial question is how to choose a particular value of P and thus how to scale the relative SAR for obtaining a realistic absolute SAR as input for a BHTE run. The exact value of absorbed power in the patient is unknown, and also it may vary in time during therapy. However, a certain physically consistent, time-averaged Power Range, (P_range), P_min P P_max, can be assumed for every patient model in use. P_range can be derived from power measurements (see Results section, paragraph Clinical measurements and derivation of absorbed power in the patient model). This step of the optimization workflow is indicated on the left side of Figure 1 as an oval box 'Power measurements'. P_range is a continuous interval, but we approximate it with a set of discrete power values P i (the total number of power values: n_power), as follows: P range % fP 1 , :, P i , . . . , P n power g i ¼ 1 . . . n power (1) with DP is an equidistant spacing between neighboring power values. The necessary condition is n_power >1. Thus, in the simplest case, n_power ¼ 2: P 1 ¼P_min and P 2 ¼P_max. For our calculations we used five power values, n_power ¼ 5, [see Results section, Equations (11a) and Equation (11b)]. Estimation and discretization of P_range according to Equation (1) is indicated as a box below the oval box on the left side of Figure 1 (Estimation of Power Range … ).

BHTE-related parametrization: perfusion models
As shown in Figure 1, before BHTE calculations, a parametrization was applied related to thermal properties, especially different perfusion values. We refer to this as parametrization related to different perfusion models. Perfusion models were region-based, i.e., inside a region/organ the perfusion value did not vary from point to point.
All investigated perfusion models of a patient were grouped into three 'Perfusion Behavior Types' (PBTs): 'homogeneous' (H) Perfusion Behavior Type (H_PBT): all organs/regions had the same perfusion value. (This is the simplest, unrealistic perfusion behavior, used here for fast and coarse estimation of the perfusion in a simplified patient model, as described in the next paragraph Two-Step Optimization). 'non-dynamic' (ND) Perfusion Behavior Type (ND_PBT): for every organ/region a temperature-independent perfusion value was defined. 'dynamic' (D) Perfusion Behavior Type (D_PBT): for every region a temperature-dependent ('dynamic') perfusion model was defined.
In the D-type perfusion models the perfusion was allowed to rise twofold or threefold (depending on region) from basal perfusion values, described by the so-called 'ramp functions' similar to Equation (3) in [36] (there used only for muscle and fat), dependent on the calculated temperature, T, as follows: where R ¼ 0.8 and S ¼ 3 for muscle and fat,and where R ¼ 0.4 and S ¼ 2 for the tumor, rectum, bladder and bones. Note that in our paper the basal values used in Equation (2) arein difference to Equation (3) in [36] not fixed, but are part of an optimization problem.
In the historical context, functions similar to Equation (2) were used, for instance, in [43] and [44]. Further details of the approaches in [43] and [44], as well as the rationale for using particular tissue properties (mainly based on [44]), are discussed in detail in [36], Section 'B. Thermal Tissue Properties'.
The total number of perfusion models of a patient (indicated as 'n_perf' in Figure 1), arising from the BHTE-related parametrization, can be formulated as follows: In Equation (3) the corresponding consecutive perfusion are the total numbers of H-or NDor D-type perfusion models of a patient.
In general, for each of the three PBTs the total number of perfusion models of a patient can be formulated as a combination of region perfusion models in the following way: n_reg is the total number of regions in a perfusion model of a patient.

Two-step optimization
According to Equation (4), J_TYPE(i_reg), i.e., the total number of perfusion values/models created per region can be extremely high or even theoretically infinite. Thus, the general aim is to limit the total number of perfusion models to a discrete number, but still yielding reasonable results within P_range. As demonstrated in Figure 2, in our optimization procedure this discretization of perfusion values is performed in two steps: Coarse Step (Pre-optimization) and Fine Step (Optimization). The coarse pre-optimization uses H_PBT, the Fine Optimization uses ND_PBT or/and D_PBT, respectively. The aim is to coarsely estimate a whole-body perfusion level (indicated in Figure 2 as 'm_best_H') that optimally corresponds to the applied power levels that were defined in Equation (1). For evaluation, we use the so-called 'Highest-Incidence-Index' approach (see the last paragraph in the section Materials and Methods). In the Fine-Optimization step, region perfusion values are varied, which is indicated in Figure 2 as 'm_best_H ±Dm'. For definition of the particular region perfusion values/models see the section Results.

PTMs
Due to the FDTD-and BHTE-related parametrizations and discretization of P_range several combinations of, 1, therapy, 2, positioning, 3, power level, and, 4, perfusion were created, referred to as PTMs. The total number of PTMs for each of the three PBTs can be determined as follows: n TYPE ¼ n therapies n SAR n power n perf TYPE (5) , n_therapies: total number of investigated therapies; n_SAR: total number of relative SAR distributions due to the FDTD-related parametrization; n_power: total number of power values representing P_range for one therapy (see Equation (1)); and, n_perf_TYPE as explained in Equation (3).
In our particular case in the section Results the values are: n_therapies ¼ 2, n_SAR ¼ 2, n_power ¼ 5, but, generally, other values are also possible.
In analogy to Equation (3) the total number of PTMs can be given as follows: where by n_H, n_ND, n_D are defined via Equation (5).

BHTE runs
For each of n PTMs a BHTE run was performed using our FD scheme, presented in the Appendix of [36]. The general BHTE formulation is not repeated here, it is presented, for instance, in Equations (1) and (2) in [36]. For detailed information about thermal and material parameters applied for solving the BHTE see 'Supplement Part 3'. The transient temperature behavior in every voxel of the BHTE model is not discussed in this paper, but it is similar to the 'steady imaging SAR' ('S-SAR') mode curve shown in Figure 4 in [36] with a temperature increase phase and a steady-state phase with a (maximum) constant temperature.
Only steady-state temperatures are used for comparisons of Temperature Sets in Equation (7) further below (see also again Figures 1 and 2). These temperature comparisons are described in the next paragraphs. For further details to BHTE runs see 'Supplement Part 3' and 'Supplement Part 5'.

RMSD-based optimization
The general purpose of this step and the following 'Evaluation' step ( Figure 1) is to find those PTMs and perfusion models which best match the measured data.
To this end all n PTMs (and some sets of PTMs, see below) were grouped together, due to the particular therapy (HP_THER/LP_THER), positioning (CE_POS/SH_POS), and power level (P 1 … P 5 ) (and additionally also due to the entire power range, see below). This resulted in a certain total number of cases, n_cases, (in our optimization scheme, the total number of n_cases is for each PBT identical). Then, for each case, i_case, a minimum RMSD value was identified via comparison of a total number of n_perf_TYPE perfusion models, as follows: where h i : measured temperatures shortly before power-off in i ¼ 1 … N points, and, T i : calculated steady-state temperatures In the simplest case the total number of squared differences, N, is equal to the total number of catheter points, N cath (here 14), but for some additional cases, if several PTMs are grouped together, it can also be equal to multiples of N cath (see below).
In this paper, those PTMs which exhibit minimum RMSD according to Equation (7) using N¼N cath are referred to as 'suitable PTMs'. Their total number is identical for each PBT and can be specified, using Equations (5) and (6), by the following expression: In the particular case presented in the section Results, n_suitable_ptm ¼ 225 ¼ 20.
Furthermore, we refer to a particular PTM that exhibits a minimum RMSD in Equation (7) among some specified group of suitable PTMs as the 'optimum' PTM.
Our dedicated PTM comparison procedure in terms of RMSD-based minimization in Equation (7) was composed of two parts.
In Part 1, the RMSD-based minimization was performed on 'single' PTMs (i.e., without grouping them together). Thus, the total number of squared differences, N, in Equation (7) was equal to N cath (N cath ¼14 in our particular case), resulting in the total number of results identical to the total number of suitable PTMs in Equation (8) (20 results in our particular case).

Part 2: 'ALL5-RMSD-Minimization'
In Part 2, those PTMs which were related to the same power range [Equations (1b) and (1c)] were grouped together into sets of PTMs, each set containing a total number of n_power PTMs (here n_power ¼ 5, hence 'ALL5'). Each set of PTMs was 'simultaneously' evaluated in a single RMSD expression, thus N ¼ n_powerN cath in Equation (7), resulting in our particular case in 70 (instead of 14 in Part 1) calculated temperature values T i in Equation (7). Finally, due to this RMSD-based minimization scheme including above Parts 1 and Part 2, the total number of cases used in Equation (7) can be calculated using the following expression (identical for each PBT): n cases ¼ n suitable ptm þ n therapies n SAR (9) whereby the value 'n_suitable_ptm' is specified in Equation (8). In our particular case presented in the section Results: We refer to those perfusion models that are associated with the cases in Equation (9) and optimized via Equation (7) as 'suitable perfusion models'. Their total number can be as high as n_cases. Usually, however, it is much smaller because a particular suitable perfusion model can be associated with more than one case in Equation (9).

Evaluation: identification of optimum perfusion models
We present below three Perfusion Evaluation Approaches for identifying optimum perfusion models of a patient among all suitable perfusion models associated with n_cases in Equation (9).

Perfusion evaluation approach 1: (SinglePTM approach)
This evaluation approach ('Approach 1') considers only the results of SINGLE-PTM-RMSD-Minimization. Among a defined group of suitable PTMs, the PTM with the minimum RMSD is determined (optimum PTM) and the perfusion model associated with this PTM is selected as the optimum perfusion model for Approach 1. This approach, however, should rather not be used for the prospective planning (i.e., prediction of the next therapy data) because there is rather a low probability that the aforementioned optimum PTMwith a particular power level and other particular parameters associated with this PTMwill be the best one for the next therapy, too. In contrast, this approach yields the best retrospective therapy information.

Perfusion evaluation approach 2: (ALL5-Approach)
This evaluation approach ('Approach 2') considers only the results of ALL5-RMSD-Minimization. Among a defined group of ALL5 cases, the case with the minimum RMSD is determined and the perfusion model associated with this ALL5 case is selected as the optimum perfusion model for Approach 2.
This approach is more general than Approach 1 (it includes a whole power range instead of one particular power value) and can be applied for prospective planning for the next therapy of the same patient when it is assured that the next therapy will be carried out similarly like the case associated with the optimum perfusion model (similar patient positioning, similar power range, the same channel setting, etc.).
This approach yields also good retrospective information about the therapy or therapies which were already performed, but, due to its averaging nature, not as exact as Approach 1.

Perfusion evaluation approach 3: ('Highest-Incidence-Index' approach)
This evaluation approach ('Approach 3') considers both the results of SINGLE-PTM-RMSD-Minimization and of ALL5-RMSD-Minimization. Approach 3 selects that perfusion model which was the most frequently occurring among a defined group of suitable PTMs and a defined group of ALL5 cases, whereby the descriptor is the so-called 'Incidence Index', P N, defined for a particular suitable perfusion model of a patient as follows: where n_single: total number of suitable PTMs containing a particular suitable perfusion model of a patient, n_range: total number of ALL5 cases containing this suitable perfusion model. (This means, with n_power ¼ 5, a fivefold weighting of those suitable perfusion models that are associated with ALL5-RMSD-Minimization). Thus, P N is equal to the total number of squared differences in Equation (7) associated with a particular suitable perfusion model of a patient (examples see Results section).
Approach 3 is the most general one because the identified perfusion model of the patient is associated with the highest number of different cases, e.g., originating from two therapies, from two patient positioning cases and from two different power ranges, like in our case in the section Results. Thus, it can be applied for the prospective planning for the next therapy of the same patient when, e.g., the power range and/or the patient positioning are not yet exactly known. This approach can also be used for different/new patients. However, in this case a reasonable preselection of parameters and generation of 'therapy-clusters' is necessary, which represent similar therapies/ cases. In contrast, due to the high number of involved different cases this approach yields still valuable but rather coarsely averaged retrospective information about the therapy or therapies which were already performed.
In conclusion, it is recommended to perform the retrospective analysis including all three perfusion evaluation approaches and to apply Approach 2 and/or Approach 3 for prospective planning. Moreover, only Approach 3 was used for the H_PBT-based coarse pre-optimization ( Figure 2). Further details of the particular use of these approaches can be found in Discussion.

Results
Clinical measurements and derivation of absorbed power in the patient model Note that in Equations (11a) and (11b) the equidistant spacing between neighboring power values, DP, defined in Equation (1), is DP ¼ 5 W for both HP_THER and LP_THER.
The derivation of particular values in Equations (11a) and (11b) using forward and reflected power measurements and FDTD calculations is described in detail in 'Supplement Part 4'.

Results of the PTM-based optimization procedure
The identification of suitable (and optimum) PTMs and perfusion models of the patient was performed applying the Two- Step Optimization scheme as depicted in Figure 2. The main optimization steps are presented below as Steps 1-6.
Step 1: initial coarse BHTE-related parametrization using H_PBT In this initial step, H_PBT was assumed to be valid for the whole patient [n_reg ¼ 1 in Equation (4)] with uniform thermal properties equal to those of muscle (2nd line from the bottom in Table 2). Seven perfusion levels were investigated (3)] with typical perfusion values of m ¼ 2.5, 5, 7.5, 10, 15, 20, and 30 ml/100g/ min), resulting in a total number of n_H ¼ 140 PTMs [n_H ¼ 2257 in Equation (5)] to be compared. The range of above seven basal perfusion values includes typical values used in the literature for muscle and fat tissue. Again, a good overview is given in [36], but other hyperthermia groups use values in this range, too. For instance, in [17,] the values for fat/muscle/tumor perfusion of fat are given as W b ¼ 1.1/3.6/1.8 kg/m 3 /s, which, with W b ¼qÁq b Á m, is equivalent to m % 7/20/10 ml/100g/min, respectively.
The optimization results for H_PBT are listed in the lines 5-10 of Tables 3 and 4, respectively (24 cases for HP_THER and for LP_THER, respectively). For every case, using the aforementioned short notation, the respective number of Table 1. Clinically measured power and temperature values for two investigated therapies of a pediatric pelvic rhabdomyosarcoma: Higher-Power-Therapy (HP_THER) in the middle column, and Lower-Power-Therapy (LP_THER) in the right column. In the last line the power range P_range is given, derived via a combined experimental-numerical power investigation, involving forward and reflected powers at the amplifier outputs, cable losses between amplifier outputs and antenna feedpoints, as well as and stray power losses. (7) for the best case in Table 3  each suitable perfusion model exhibiting minimum RMSD according to Equation (7) is also reported. The lines 5-9 refer to SINGLE-PTM-RMSD-Minimization (10 suitable PTMs per therapy, 5 for CE_POS and 5 for SH_POS), for each of the five power values separately. The 10th line refers to ALL5-RMSD-Minimization (2 ALL5 results per therapy, one for CE_POS and one for SH_POS).

Figure 4. A root-mean-squared-deviation (RMSD)-based minimization in Equation
Step 3: evaluation: identifying the optimum H-type perfusion model The purpose of this step was to identify a particular perfusion model (indicated as 'm_best_H' in Figure 2) as a starting point for the next step, i.e., the fine optimization. As this purpose was merely a prospective one (and no particular cases were to be preferred/excluded) only the most general perfusion evaluation approach, Approach 3 ('Highest-Incidence-Index' Approach) was applied. According to Tables 3 and 4 Tables 3  and 4 that are associated with H-4, we end up with the Incidence Index, P N, as high as 532 (1814 þ 470), and for H-5 we obtain only P N ¼ 28. Table 5 summarizes, in the right column, the results of Approach 3 displaying the corresponding Incidence Index P N. The two H-type perfusion models can be found in lines 2-3. Furthermore, all suitable perfusion models (i.e., those which are displayed in Tables 3 and 4) are reported in alphabetic order and all region perfusion values/models are listed.
Step 4: fine BHTE-related parametrization: creating NDtype and D-type perfusion models In this step, with knowledge of the optimum whole-body perfusion m ¼ m_best_H ¼ 10, as found above, perfusion values of different organs/regions were specified, which is indicated in Figure 2 as 'm_best_H ±Dm'. Apart from the perfusion, no other tissue-specific thermal properties were parametrized. Their (fixed) values as used for ND_PBT and D_PBT are indicated in Table 2.
The creation of ND-type and D-type perfusion models of the patient was performed for six regions, n_reg ¼ 6 in Equation (4). In the source code implementation six nested loops were written with a consecutive perfusion index or Table 4. Optimization results only for the 'Lower-Power'-Therapy (LP_THER) in terms of minimum RMSD values, displayed in an analogous way to  Table 5.   Table 3 are generally better than the respective ones in Table 4. Region perfusion values for every perfusion model reported in this table can be found in Table 5. Table 2. Tissue perfusion values in perfusion models exhibiting the nondynamic (ND) and dynamic (D) Perfusion Behavior Type (PBT) in the optimization procedure described in Figure 1, 2, and 3.  For ND_PBT the highest number for perfusion values, five, was allowed for the two tissues with inserted catheters, in the tumor and in rectum. For these regions two perfusion values were chosen beneath m ¼ 10, and two above. Following ND-type perfusion values were set for each of the six regions (a symbol 'j' is used as a delimiter): The total number of ND-type perfusion patient models of the patient n_perf_ND, can be calculated as a combination of region perfusion models defined in Equation (12), yielding, according to Equation (4), n_perf_ND ¼ 553332 ¼ 1350.
Almost analogically, the loop structure for D_PBT was realized. However, the use of D_PBT, with varying/dynamic perfusion values in each region, allowed a better 'covering' of the perfusion range, so we used only 4 tumor and rectum perfusion models instead of 5 in N_PBT. In turn, the number of fat and muscle perfusion models was increased in comparison to N_PBT, from 3 to 4. For simplicity, we describe dynamic region perfusion models in Equation (2) using a short notation with arrows like '5 ! 10' for a double perfusion increase, or '5!15' for the threefold increase. The following D-type region perfusion models were set for each of the six regions: The total number of D-type perfusion models of the patient, n_perf_D, can be calculated as a combination of region perfusion models defined in Equation (13), yielding, according to Equation (4), n_perf_D ¼ 443442 ¼ 1536.
Step 5: PTM creation and fine optimization using ND_PBT and D_PBT In this step, with the perfusion models identified in Equation (12), a total number of n_ND ¼ 27000 ND-type PTMs was created and 27000 BHTE runs were performed (Equation (5), Figures 2 and 3). Similarly, due to Equation (13), n_D ¼ 30720 PTMs and BHTE runs were applied for D_PBT, respectively. Next, all calculated Temperature Sets were compared with two measured Temperature Sets (for HP_THER and LP_THER) in terms of minimum RMSD (Equation (7)), using SINGLE-PTM-RMSD-Minimization and ALL5-RMSD-Minimization, resulting, in accordance with Equation (9), in an identical number of comparison cases, n_cases ¼ 24, for each of the In the last row, results of 'Approach 3' are reported by displaying the corresponding 'Incidence Index', RN, defined in Equation (9). The H-type, ND-type, and Dtype perfusion models with the highest RN values, respectively, are marked in bold.
two PBTs, respectively (20 single PTM cases and 4 ALL5 cases). In analogy to H_PBT optimization results of all 224 cases are listed in terms of minimum RMSD values (212 results for HP_THER in Table 3 and 212 results for LP_THER in Table 4, respectively). ND_PBT cases are displayed in the middle segment (lines 12-17) of Tables 3 and Table 4, respectively, whereas D_PBT cases are listed in the bottom segment of these tables (the last 6 bottom lines).
Generally, RMSD values for HP_THER (Table 3) were lower than those for LP_THER (Table 4) and CE_POS yielded overall better results than SH_POS. In addition, ND-type results were overall slightly worse than D-type results (but much better than H-type results).
In Figure 4 we illustrate the RMSD-based minimization via Equation (7) for the comparison between i_perf ¼ 1 … n_perf_D ¼ 1536 calculated Temperature Sets (indicated as 'MODEL-NUMBER' on the abscissa) and one measured Temperature Set for HP_THER. The comparison is made for the case that is associated with the optimum PTM (among all PTMs) parameter combination '(D_PBT K HP_THER K CE_POS K P ¼ 120 W)', which is marked in bold in Table 3 in the next paragraph. A minimum RMSD value (0.3882 C) was found for 'i_perf ¼ 1001', i.e., for a PTM that involves the perfusion model D-1001. High RMSD values (up to almost 5 C) were obtained for low perfusion model indices that are associated with obviously too low tumor and rectum perfusion values.
Step 6: evaluation of optimum ND-type and D-type perfusion models Perfusion evaluation approach 1 (SINGLE-PTM approach).
Among all 20 ND-type and 20 D-type suitable PTMs, the optimum PTM [with the lowest RMSD in Tables 3 and 4 (0.3882 C, marked bold in Table 3)] was a D-type PTM and was characterized by the following parameter combination: '(D-1001 K HP_THER K CE_POS K P ¼ 120 W)'. An optimum ND-type PTM (RMSD ¼ 0.5882 C, thus worse than D-1001, also marked bold in Table 3) was obtained for the following combination: '(ND-1011 K HP_THER K CE_POS K P ¼ 120 W)'. Thus, as a result of Approach 1 the dynamic perfusion model D-1001 and the non-dynamic perfusion model ND-1011 were found as optimum perfusion models. For particular perfusion values/models see Table 5. Both Approach-1 results are shown in Figure 5 in the next paragraph.
Perfusion evaluation approach 2 (ALL5 approach). The optimum ALL5 value (among all ALL5 values) was associated with a dynamic perfusion model D-994 (HP_THER K CE_POS, RMSD ¼ 0.5417 C, marked in bold in Table 3). For ND_PBT an optimum ALL5 result was obtained again for HP_THER and was characterized by the perfusion model ND-992 and CE_POS (RMSD ¼ 0.5882 C, marked in bold in Table 3). Thus, as a result of Approach 2 the dynamic perfusion model D-994 and the non-dynamic perfusion model ND-992 were found as optimum perfusion models. For particular perfusion values/models see Table 5.
In conclusion, all optimum perfusion models, found via Approach 1 and Approach 2, were obtained for HP_THER and CE_POS. Furthermore, the RMSD values in Table 3 are generally lower than those in Table 4, thus HP_THER was better approximated by the optimization procedure than LP_THER.
Perfusion evaluation approach 3 ('Highest-Incidence-Index' approach). As demonstrated in Table 5, the optimum perfusion models found via Approach 3 were D-906 for D_PBT ( P N ¼ 154, marked in bold) and ND-955 for ND_PBT ( P N ¼ 140, marked in bold, too), respectively. The dynamic perfusion model D-906 is present in Table 3 (HP_THER) for three power values for SH_POS and for one ALL5 value, for SH_POS, too. Furthermore, this perfusion model can be found in Table 4 (LP_THER) also for three power values for CE_POS. The non-dynamic perfusion model ND-955 is present in Table 3 (HP_THER) for two power values for SH_POS and for one ALL5 value for SH_POS, too. Moreover, in Table 4 (LP_THER), this perfusion model can be found for one power value for CE_POS and two power values for SH_POS.
Perfusion evaluation approaches 1 and 2 for LP_THER. If skipping Table 3 and comparing only LP_THER (Table 4), Approach 1 yields as an optimum perfusion model D-913 (LP_THER K CE_POS K 113 W), with the RMSD value as high as 0.61 C. This is much higher than 0.3882 C which was the best case in Table 3. The model D-913 is also the optimum perfusion model for the entire P_range (Approach 2) yielding 0.6358 C for RMSD. Thus, D-913 proves to be a useful perfusion model for LP_THER alone. For ND_PBT, moreover, Approach 1 yields as an optimum perfusion model ND-943 (LP_THER K CE_POS K 108 W) with RMSD ¼ 0.6308 C. Both Approach-1 results are shown in Figure 6 in the next paragraph.

Graphical comparison along the catheters for both therapies
In Figure 5 we graphically compare temperatures which were measured along the catheters in the steady-state phase of HP_THER shortly before 'power-off' (black dashed curve) with three calculated Temperature Sets which were associated with the optimum PTMs in Table 3, each one for each of the three PBTs, respectively (H_PBT, ND_PBT, D_PBT). A similar comparison is performed in Figure 6 for the optimum PTMs in Table 4 for LP_THER. For both therapies measured temperatures in the tumor were slightly higher than in rectum. Temperatures for HP_THER were on average higher than those for LP_THER by around 0.5-1 C (see also Table 1). The highest temperature value for HP_THER was measured in the bladder (black cross icon). For LP_THER the highest temperature value was in the center of the tumor. In the tumor, a decrease in temperature toward the lateral tumor margin (point 4) was observed, which may be explained by cooling effects originating from water bolus or by higher perfusion values for the tumor margin than for the tumor center. In the rectum, the temperatures were relatively constant, except for a decrease toward the caudal end (point 0). This can be explained by the influence of extracorporeal temperature conditions (cooling effects) at the anus.
The H-type PTMs exhibited for both therapies too high temperatures in the tumor and too low in rectum. Interestingly, improvement in inhomogeneous case was not due to lowering rectum perfusion in comparison to that of the tumor, but rather by lowering the fat perfusion and at the same time by equally increasing both rectum and tumor perfusion values (compare, in Table 5, the region perfusion values between the models H-4 and ND-1011 or between the models H-4 and D-1001).
In conclusion, although a substantial minimization of RMSD values was achieved (Figure 4), no one of the optimum PTMs exactly approximated the characteristic tumor temperature behavior that exhibited a temperature increase in the central tumor part, especially in the case of LP_THER ( Figure 6). In contrast, a good approximation of temperature behavior in rectum was achieved. Probably, for future comparisons, a more complicated subtumoral model is needed, for instance with a less perfused tumor center and a higher perfused margin.

Discussion
This paper describes a novel optimization procedure for standardized comparisons in patient models using clinical measurements (temperature profiles along the catheters shortly before power-off and power measurements). The procedure is presented using data of two hyperthermia therapies of the same pediatric patient treated in SIGMA-30.   Figure 6. Calculated results refer to those PTMs, which exhibited the minimum RMSD in Table 4, i.e. those which are the best for LP_THER only, for each perfusion behavior type (PBT), respectively. These PTMs are indicated as follows: 'H-4-centr-93 W' (red, circle symbols) (homogeneous PBT, CE_POS, p ¼ 93 W). 'ND-943-centr-108' W (blue, square symbols) (non-dynamic PBT, CE_POS, p ¼ 108 W), and 'D-913-centr-108W' (purple, triangle symbols) (dynamic PBT, CE_POS, p ¼ 108 W). Region perfusion values are displayed in Table 5.
Forward and reflected power measurements and FDTD calculations were used for identification of the absorbed power range for each therapy. For comparisons, involving both the FDTD and the BHTE, various combinations of, 1, therapies, 2, patient positioning in the applicator, 3, power levels/ranges, and, 4, perfusion models were created, referred to as Parametrized Treatment Models (PTMs). Three different perfusion behavior types (PBTs) were investigated, homogeneous (H), non-dynamic (ND) and dynamic (D). Perfusion models were region-based, i.e., inside a region/organ the perfusion value did not vary from point to point. Three approaches for evaluation of the best perfusion models for both retrospective analysis and prospective planning were presented.
In the PTM-related retrospective analysis, the PTM with the lowest RMSD (0.3882 C) was a D-type PTM which was characterized by the following parameter combination: '(D-1001 K HP_THER K CE_POS K P ¼ 120 W)'. However, it cannot be generally expected that the dynamic perfusion model D-1001 that was identified to be the optimum one for the retrospective analysis (Perfusion Evaluation Approach 1) will be the best one for the prospective planning, too. This is because it is rather a low probability that the aforementioned optimum PTMwith the particular power level of 120 W and other particular parameters associated with this PTMwill be the best one for the next therapy, too. This can rather be expected for the dynamic perfusion model D-994, the optimum one among all ALL5 results (Approach 2) which was associated with the parameter combination '(HP_THER K CE_POS, RMSD ¼ 0.5417 C)' when it is assured that the next therapy is carried out similarly like the case associated with the perfusion model D-994 [similar patient positioning CE_POS, similar power range to that specified in Equation (11a) and the same channel setting]. With other words, a perfusion model like perfusion model D-994, the optimum model for Approach 2, which optimally approximates 5 power values [values in Equation (11a)], is more likely to approximate the next power value (in the next therapy) than a model that approximates a single power value (for instance a power of 120 W that is 'required' for the perfusion model D-1001, the optimum model for Approach 1).
The Perfusion Evaluation Approach 3 identifies the best perfusion model via the Highest-Incidence Index, P N, defined in Equation (10). As explained in the section Materials and Methods, this approach is the most general one, and it can be applied for the prospective planning for the next therapy of the same patient and/or for different/ new patients. However, in the latter case it requires a reasonable preselection of parameters and generation of 'therapyclusters' representing similar therapies/cases. In contrast, due to the high number of involved different cases this approach yields still valuable but rather coarsely averaged retrospective information about the therapy or therapies which were already performed.
Note that all optimum perfusion models identified via Approach 1 and Approach 2 relate to HP_THER and CE_POS. Thus, firstly, the higher-power therapy was easier to approximate than the lower-power one, and, secondly, the real positioning of the patient in the applicator was indeed more similar to the central one than to the shifted one. Moreover, all optimum perfusion models identified via Perfusion Evaluation Approaches 1-3 exhibited the dynamic perfusion behavior type (D-1001 for Approach 1, D-994 for Approach 2, and D-906 for Approach 3). Thus, the dynamic perfusion behavior specified in Equation (2) seems to be a very useful model for hyperthermia treatment planning.
In general, presentation of the results in terms of all cases, like it is done in Tables 3 and 4, is strongly recommended. In this way, evaluation of parameter combinations which are subsets of n_cases, can be performed, for instance only those cases which are associated with one particular therapy, one particular positioning, etc., can be selected. These constellations of cases should be selected depending on the purpose. If, for instance, a third therapy of the same patient should be predicted from the retrospective data of two previous ones, then the perfusion evaluation should include data of both therapies (like HP_THER and LP_THER in Approach 3 resulting in the identification of the perfusion model D-906). If a completely new patient is planned, only data of similar patients should be included in the database.
Such prospective planning purposes can result in a further generalization of RMSD-based optimization, i.e., a further increase of the value of 'N' in Equation (7) after a reasonable selection of cases for the database, depending on the concrete purpose. Following this methodology, the next step beyond the RMSD-based minimization scheme presented in this paper could be an enlargement of N to 140 or even to N ¼ 280. In an 'optimization sense' this would be equivalent to searching for a single perfusion model that is related 'simultaneously' to 20 power values, two therapies and two patient positions (instead of searching for many optimum perfusion models case by case as presented in this paper). These 'N-enlargement cases' are not shown in this paper in detail, but, interestingly, they result in an identification of the same perfusion model D-906 that was also found using Approach 3, confirming the general nature of this perfusion model. Alternatively/additively, several cases originating from different therapies and different power ranges can be included in the evaluation scheme. Due to the general nature of Approach 3, it can be applied to other RMSD-based minimization schemes originating from any combinations of particular values of N in Equation (7) and the evaluation can consider any number of cases, applying weighting procedures similar to those in Equation (10). For proving if thesetheoretically plausibleperfusion evaluation approaches are useful in clinical practice, further investigations are needed in the near future.
The optimization procedure is presented for a pediatric patient in SIGMA-30, but it can be easily extended to other patient types, including adults, and it can be applied to other hyperthermia applicators. The necessary condition for obtaining realistic perfusion values is, however, estimation of a reasonable range of absorbed power in a given patient for a given applicator.
The future goal is to be able to predict any therapy using prospective analysis based on such patient-specific PTMs and by improving the accuracy of planning methods and modelsto reduce the need of thermometry. This requires a higher number of parameters than those investigated in this paper (several channel settings, more perfusion models, various patient positions, etc.). Moreover, data derived from local T-rise/-drop measurements in patients can also be included into future PTM comparisons, requiring, however, additional comparison formulas to that in Equation (7). For that, further retrospective investigations are needed to build a database, including a sufficiently large pool of different therapies based on optimization procedure principles developed in this paper, resulting in a cumulative learning effect from one therapy to the next. Moreover, collecting PTM data of different patient therapies to the aforementioned database would allow prediction of temperatures in similar patients in the future. To this end, this paper is the first step in this perspective.

Conclusions
A novel procedure for identification of patient-specific electromagnetic and thermal models with different parameter settings and parameters, which we refer to as Parametrized Treatment Models, was developed and evaluated by comparison with measurements obtained from two clinical hyperthermia therapies of a pediatric patient treated in SIGMA-30. The determining factor for accurate calculations of patientspecific therapy planning was the proper choice of the perfusion models according to the specified absolute values of the absorbed power in the patient model. For identification of optimum perfusion models, which can be used for improved patient planning, dedicated evaluation approaches were proposed. The optimization procedure can be extended to other patient types, including adults, and to different hyperthermia applicators. The future goal is to be able to predict temperature during any therapy using prospective analysis based on such patient-specific PTMs andby improving the accuracy of planning methods and modelsto reduce the need of thermometry.

Disclosure statement
I hereby declare that the author and coauthors of this manuscript have no conflicts of interest in the authorship or publication of this contribution.