Iterative time-reversal for multi-frequency hyperthermia

TR focusing is achieved by placing a so-called virtual source at the desired focal spot and by recording the impulse response at the antenna locations (figure 1(a)). The recorded signals are then back-propagated simultaneously by all the antennas, generating an interference pattern that exhibits full phase coherence at the original source location (figure 1(b)). Together with the desired peak, however, secondary lobes appear due to both the limited aperture of the array and the anatomy of the patient. In the frequency domain, the TR operator corresponds to a complex conjugation (^*):

$$\begin{eqnarray}&&{p}_{c}(f)={\left[{p}_{c}^{R}(f)\right]}^{* },\end{eqnarray} \tag{ 1 }$$

where ${p}_{c}^{R}$ is the complex phasor recorded by the array channel c at frequency f during source excitation, while p_c is the complex steering parameter for the same channel to focus at the source location. Under this perspective, a discrete set of simultaneous operating frequencies can be selected for treatment, rather than wideband pulses. In practice, an applicator can be operated to switch between single frequencies in finely interleaved time slots, if the duration of each slot is small compared to the biological response (Trefná et al 2017). This simplifies the complexity of the hardware needed for feeding the HT applicator (Trefná et al 2012). The rest of the analysis will therefore be carried out assuming that n_f simultaneous operating frequencies can be independently controlled in phase and amplitude for each of the n_c channels of the array applicator.

Zoom In Zoom Out Reset image size

TR focusing requires in principle only one EM simulation per virtual source, which makes it attractive in terms of computational resources. Unfortunately, this single simulation is not sufficient to assess the quality of the resulting treatment planning. The HTP quality indicators must in fact be evaluated over the focused SAR distribution (figure 1(b)). Moreover, since the immediate TR solution for HTP might not be optimal or clinically viable, further adjustments might be needed, which would require additional simulations. In view of this, it is more convenient to run a separate simulation for each antenna/channel c in the array, as is the case for EV or PS. The TR solution at any point can then be determined from this set of E-field distributions via a generalized and source-free TR method described later in section 2.2. Denoting with ^t the transpose operation (without conjugate), the complex E-field distributions and steering parameters relative to all channels at a certain frequency can be joined into column vectors for compactness:

$$\begin{eqnarray}&&{{\boldsymbol{p}}}_{f}={[{p}_{1,f}\,{p}_{2,f}\,...\,{p}_{{n}_{c},f}]}^{t}\quad \quad {{\boldsymbol{E}}}_{f}={[{E}_{1,f}\,{E}_{2,f}\,...\,{E}_{{n}_{c},f}]}^{t}.\end{eqnarray} \tag{ 2 }$$

A number of challenges arise when implementing TR as HTP. First of all, the choice of the source location, or focal spot. Intuitively, the source should be placed at the center of the tumor, but this is often hard to define, especially for irregular shapes. The center of mass of the tumor can be picked as an initial guess, but this choice often leads to inhomogeneous or insufficient heating of deep or peripheral tumor areas. Cold-spots in the target volume should be avoided, as they reduce the treatment's efficacy. Better coverage can be achieved, for example, by translating the focal spot progressively towards the internal, deeper part of the target volume (Takook et al 2017). Second, HSs arising in healthy tissues must be suppressed, to allow a higher deposition in the target. In order to shift the focus away from these locations, the natural approach under the TR perspective would be to place secondary sources, record their impulse response, and subtract these secondary solutions from the primary one. Taken together, these considerations suggest that an initial TR solution with source at the tumor center can be iteratively improved by shifting the set of complex steering parameters away from HSs and towards cold-spots.

A second aspect of TR focusing is the choice of the polarization axis for the virtual source. The polarization problem in vector field focusing has been addressed in several ways in the literature, also recently (Iero et al 2016, Battaglia et al 2020). In the case of a cylindrical applicator, the antennas can be aligned to the symmetry axis (figure 2(a)): the virtual source should then also be aligned with this principal axis. In more complex array topologies such as a semi-spherical helmet, the E-field generated by each antenna can vary significantly in polarization, depending on the location of the focal spot (figure 2(b)). The resulting interference is no longer characterized by a single main polarization, and the concept of virtual source becomes inadequate. This concept should then be dropped in favor of the more generalized concept of complex power contribution at the focal spot. We know that the SAR at a location ζ resulting from the superposition of all channel fields at a particular frequency f is given by:

$$\begin{eqnarray}&&{{\rm{SAR}}}_{f}(\zeta )=\displaystyle \frac{1}{2}\displaystyle \frac{{\sigma }_{f}(\zeta )}{\rho (\zeta )}{{\boldsymbol{p}}}_{f}^{* }{}^{t}\langle {{\boldsymbol{E}}}_{f}^{* }(\zeta ),{{\boldsymbol{E}}}_{f}^{t}(\zeta )\rangle {{\boldsymbol{p}}}_{f},\end{eqnarray} \tag{ 3 }$$

where 〈 · , · 〉 denotes the scalar product, while σ_f and ρ are the local material conductivity (frequency-dependent) and density, respectively. The outer product, called SAR matrix and denoted Q _f , contains information about the complex SAR contributions from each channel irrespective of their polarization directions:

$$\begin{eqnarray}\begin{array}{rcl} & & {{\boldsymbol{Q}}}_{f}=\displaystyle \frac{1}{2}\displaystyle \frac{{\sigma }_{f}}{\rho }\langle {{\boldsymbol{E}}}_{f}^{* },{{\boldsymbol{E}}}_{f}^{t}\rangle \\ & & \langle {{\boldsymbol{E}}}_{f}^{* },{{\boldsymbol{E}}}_{f}^{t}\rangle ={\left({{\boldsymbol{e}}}_{f}^{X}\right)}^{* }{\left({{\boldsymbol{e}}}_{f}^{X}\right)}^{t}+{\left({{\boldsymbol{e}}}_{f}^{Y}\right)}^{* }{\left({{\boldsymbol{e}}}_{f}^{Y}\right)}^{t}+{\left({{\boldsymbol{e}}}_{f}^{Z}\right)}^{* }{\left({{\boldsymbol{e}}}_{f}^{Z}\right)}^{t},\end{array}\end{eqnarray} \tag{ 4 }$$

where a Cartesian coordinate system has been assumed. By construction, Q _f is Hermitian and positive semi-definite. This property allows us to perform an inverse operation and decompose Q _f into an outer product, such that:

$$\begin{eqnarray}&&{{\boldsymbol{Q}}}_{f}={\hat{{\boldsymbol{p}}}}_{f}^{* }{\hat{{\boldsymbol{p}}}}_{f}^{t},\end{eqnarray} \tag{ 5 }$$

where ${\hat{{\boldsymbol{p}}}}_{f}$ is a column vector containing the complex contribution of every channel at ζ. The solution ${\hat{{\boldsymbol{p}}}}_{f}$ of the decomposition can be obtained as the only non-singular eigenvector of Q _f . The sought set of steering parameters for focusing at ζ in a generalized TR perspective is the complex conjugate of this solution:

$$\begin{eqnarray}&&{p}_{c,f}={\hat{p}}_{c,f}^{\,* }.\end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

The proposed iterative TR (i-TR) focusing technique takes into account all the points discussed above to improve an initial TR solution for HTP. The schematic of the procedure is shown in figure 3. The initial solution p is found at the center of mass of the tumor using equation (6). The resulting SAR distribution is analyzed and the most prominent cold-spot is identified as the center of mass of the contiguous sub-volume containing the lowest 1-percentile of SAR within the target volume. The TR solution at the cold-spot is calculated and the iterations begin for one frequency at a time, shifting the phase of the current parameter set towards that of the cold-spot solution p ^C .

Zoom In Zoom Out Reset image size

The iterative section aims at minimizing a novel cost function specifically tailored for i-TR, which we will call the hot-to-cold spot quotient (HCQ), defined as:

$$\begin{eqnarray}&&{\rm{HCQ}}=\displaystyle \frac{{\overline{{\rm{SAR}}}}_{{\rm{R1}}}}{{\overline{{\rm{SAR}}}}_{{\rm{T1}}}},\end{eqnarray} \tag{ 7 }$$

where ${\overline{{\rm{SAR}}}}_{{\rm{R1}}}$ is the average SAR in the sub-volume of remaining HSs containing the highest 1-percentile of SAR values, while ${\overline{{\rm{SAR}}}}_{{\rm{T1}}}$ is the equivalent for the lowest 1-percentile in the target volume. The individual frequency iteration is carried out stepping by a geometrical factor α, with 0 < α < 1:

$$\begin{eqnarray}&&{{\boldsymbol{p}}}_{f}^{{\prime} }={{\boldsymbol{p}}}_{f}^{{\prime\prime} }\cdot {\rm{\Delta }}{{\boldsymbol{p}}}_{f}^{\alpha },\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{{\boldsymbol{p}}}_{f}=\left\{\begin{array}{l}({e}^{j\angle {{\boldsymbol{p}}}_{f}^{C}})/{{\boldsymbol{p}}}_{f}^{{\prime\prime} },\,\mathrm{for}\,\mathrm{cold} \mbox{-} \mathrm{spot}\,\mathrm{iteration}\\ (| 1/{{\boldsymbol{p}}}_{f}^{H}| )/{{\boldsymbol{p}}}_{f}^{{\prime\prime} },\,\mathrm{for}\,\mathrm{hot} \mbox{-} \mathrm{spot}\,\mathrm{iteration}\end{array}\right.,\end{eqnarray} \tag{ 9 }$$

where j is the imaginary unit, and all operations are intended to be element-wise. The buffer variable ${{\boldsymbol{p}}}_{f}^{{\prime\prime} }$, initially set as p _f , is replaced by the search direction ${{\boldsymbol{p}}}_{f}^{{\prime} }$ only if the stepping was successful, i.e. if the resulting ${\rm{HCQ}}({\boldsymbol{p}}^{\prime} )$ is better than the current HCQ( p ''). Note that, while only the subset of steering parameters relative to one frequency is iterated at a time, the assessment of HSs, cold-spots and HCQ at each iteration is always carried out over the total SAR distribution:

$$\begin{eqnarray}&&{\rm{SAR}}={{\rm{SAR}}}_{1}+{{\rm{SAR}}}_{2}\,+...+\,{{\rm{SAR}}}_{{n}_{f}}.\end{eqnarray} \tag{ 10 }$$

If the stepping was not successfull, the step factor α is geometrically reduced and another attempt is made. If the norm of the step ∣Δ p ^α ∣ falls below a threshold value Δp₀ without any improvement on the HCQ, the next frequency is considered and iterated. Once all frequencies have been iterated, only the subset of parameters p _f relative to the frequency that yielded the best improvement in HCQ is updated with ${{\boldsymbol{p}}}_{f}^{{\prime\prime} }$.

A similar procedure is then carried out for the HS. The updated SAR distribution is analyzed and the most prominent HS is identified as the highest 1-percentile of SAR within the HSs. However this time the amplitude of the parameters is shifted away from the HS solution p ^H . If at least one spot iteration was successful in improving the HCQ, a new cold-spot is considered and iterated, and so on. Otherwise, the algorithm terminates.

The parameter Δp₀ is directly related to the desired precision on the steering parameters. We select a precision of 1% (Δp₀ = 0.01) since RF systems usually do not provide steering accuracy higher than ±5 % (Wust et al 1998, Paulides et al 2007). The parameter α₀, on the other hand, decides the length of the first iterative step and its subsequent geometrical progression. As shown later in section 4.3, the value of α₀ affects mainly the algorithm's convergence speed and number of iterations, and has almost no effect on the overall accuracy. Its optimal value α₀ = 0.1 is determined by means of a sensitivity analysis, and this value is kept throughout all executions reported here.

The HTP algorithms are tested on two array topologies: cylindrical and semi-spherical. The cylindrical applicator (figure 4(a)) is meant for treatment of tumors in the neck region. It consists of 10 radiators arranged on one ring of 20 cm diameter. The radiators are self-grounded bow-tie antennas (SGBT), which work across the range 400–800 MHz and have been adapted specifically for ultra-wideband MW-HT (Takook et al 2017). A cylindrical water bolus, with average thickness 4 cm, is included between the antennas and the patient to improve matching. The SGBT antennas are further immersed each into a protruding cylindrical extension of the water bolus, to enhance their directivity and reduce cross-coupling between array elements.

Zoom In Zoom Out Reset image size

The semi-spherical applicator (figure 4(b)) represents a simplified version of a helmet applicator for brain tumors. It consists of 10 SGBT antennas radially arranged from the top of the scalp and following the skull's curvature. The applicator covers only a small part of the head and, as such, can only be meaningfully used for targeting volumes close to the the skull. This applicator is not optimized for an effective treatment of deep-seated brain tumors, but it is intended for the investigation of the i-TR method's robustness to mixed polarization axes. The average water bolus thickness for this applicator is 1.5 cm.

The Duke human voxel model from the IT'IS Foundation (Gosselin et al 2014) is used as virtual patient for both regional treatment plannings. The original resolution of 1 × 1 × 1 mm is down-sampled to 2 × 2 × 2 mm to spare computational resources. The constant (density) and dispersive (permittivity, conductivity) tissue properties are obtained from the IT'IS Foundation database (Hasgall et al 2012). A 230 × 260 × 390 mm subset of the model is selected to represent the neck region, figure 4(a). The resulting voxel matrix includes 47 segmented tissues. For the brain applicator, a smaller subset is selected, 230 × 260 × 138 mm, figure 4(b). The resulting voxel matrix includes 31 tissues.

For the preliminary SAR investigation, three spherical target volumes are defined in MATLAB^® after the E-field simulations. As the original dielectric properties within the target volumes are kept, these targets are heterogeneous. In addition to these, two realistic test cases are considered: one in the larynx for the neck applicator and one in the meninges for the brain applicator. These targets are defined before the E-field simulations, and their volumes are filled with a homogeneous material whose properties are determined by taking the weighed average of the materials originally inside. This results in a conservative approach where the tumor does not exhibit increased conductivity. The volume and composition of the 5 targets are described in table 1 and shown in figures 5 and 6. Resulting tumor properties for the realistic cases are reported in table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Properties of the homogeneous materials filling the target volumes in the realistic cases, obtained as the weighed average of the original materials within the volume. Some properties are adjusted for hyperthermic conditions.

Tumor	[∼]	σ (S m−1)	ρ (kg m−3)	k (W Km−1)	c (K K−1 Kg−1)	ωb (kW Km−3)	Qm (kW m−3)
Larynx	47.8	0.674	1088	0.423	3.285	22.033	5.456
Meninges	57.3	1.372	1063	0.535	3.754	17.048	5.478

In all cases, the target volume is a contiguous shape, as multi-target treatment planning is beyond the scope of the present work. The artificial targets are included to evaluate the i-TR scheme's ability to achieve target coverage and HS suppression for different tumor sizes. The larynx case is included to assess the HTP's ability to target concave shapes with both superficial and deep areas. Note that the trachea lumen creates a cavity in the target volume. The meningioma is included to investigate the case where multiple polarization directions are interfering in the target volume.

The E-field distributions of each antenna within the applicator are obtained numerically using the commercial software CST Microwave Studio^® (Dassault Systèmes SE, Vélizy-Villacoublay, France 2019). The software implements an FDTD scheme on a hexahedral mesh with possibility for locally denser sampling grids. The complex geometrical curves of the SGBT antennas are sampled at $0.5/\sqrt{3}\,{\rm{mm}}$ to correctly represent at any angle the 0.5 mm thick copper plate with which they are realized. The mesh resolution within the patient varies between this value and a maximum of 2 mm. Together with the water bolus, this results in ≈125/220 million cells for the collar/helmet setups, respectively. The almost doubled amount of mesh cells for the helmet setup is due to the antenna arrangement occupying more effectively all 3 dimensions in space. Absorbing boundary conditions (PML) are defined for all sides of the computational domain.

Upon import for post-processing in MATLAB^®, the E-field distributions are re-sampled to a constant grid resolution of 2 mm. Out of these distributions, the local SAR matrices Q_f (x, y, z) are constructed according to equation (4). This results in 1.1 million matrices for the neck model and 0.3 million for the brain one (considering the patient only). Each element q_i,j (x, y, z) is further spatially averaged according to a 1g SAR mass-averaging scheme. The averaging scheme is similar to the CST Legacy one as described in (Qi et al 2008): we treat surface voxels by expanding a spherical kernel until the mass of patient tissues within reaches 1g.

To validate the HTP plans, steady-state thermal simulations are performed with CST Microwave Studio^®. We consider only the realistic scenarios and only the frequency combination that yields the lowest HTQ after i-TR optimization. This combination would be the one selected for treatment after quickly surveying all possible combinations with i-TR. The thermal properties of the HSs are also taken from the IT'IS database. However, to mimic the hyperthermic stress condition of the body, the muscle's blood flow perfusion is increased by a factor 4 (Lang et al 1999, De Greef et al 2010). Moreover, the thermal conductivity of the cerebrospinal fluid is increased by a factor 10 to emulate the convective transport of heat (Schooneveldt et al 2019). The tumor thermal properties are calculated as a weighed average of the tissues originally composing its volume, in the same way as for the electromagnetic properties. The tumor blood perfusion rate is further diminished by a factor 0.7 to account for the chaotic vascular structure that characterizes neoplasms (Song 1984, Lang et al 1999). The tumor properties are summarized in table 2.

The set of possible operating frequencies is constructed by stepping 100 MHz within the antenna working range 400–800 MHz, yielding the following set:

$$\begin{eqnarray}&&F={[400,500,600,700,800]}^{t}\,{\rm{MHz}}.\end{eqnarray} \tag{ 13 }$$

Out of this set, all the single and double frequency combinations are evaluated, resulting in 15 operating frequency settings. This is to highlight potential frequency-dependent phenomena within the patient and to evaluate the i-TR scheme's robustness to different combinations of operating frequencies, while at the same time keeping the amount of test cases at a minimum.

The EV implementation is similar to the one used in Guérin et al (2018). It determines via a direct solution the set of steering parameters that maximizes the SAR amplification factor (SAF), defined as the following quadratic ratio:

$$\begin{eqnarray}&&{\rm{SAF}}=\displaystyle \frac{{\overline{{\rm{SAR}}}}_{{\rm{T}}}}{{\overline{{\rm{SAR}}}}_{{\rm{R}}}},\end{eqnarray} \tag{ 14 }$$

where ${\overline{{\rm{SAR}}}}_{{\rm{T}}}$ is the average SAR in the target, and ${\overline{{\rm{SAR}}}}_{{\rm{R}}}$ is the average SAR in the rest. An important consequence of how equation (14) is constructed is that only one frequency at a time from any working set can maximize such ratio (Zanoli and Trefná 2020). The amplitude of the other frequencies must be set to zero, and the only active frequency is the one whose interference pattern best fits the target shape and location. This approach differs substantially from the continuous FIR solution described in Converse et al (2006), which, to the best of the authors' knowledge, is the only published study on wide-band EV beam-forming. Their approach is not applicable here since the spectrum is discrete. Naturally, maximizing the SAF is not the same as minimizing nonlinear indicators such as HTQ, which usually correlate better with the resulting thermal dose.

The PS implementation used for this study is the one readily available in MATLAB^® (particleswarm). The cost function is the HTQ as described at the beginning of this section. In order for PS to properly converge regardless of the problem's complexity (Piotrowski et al 2020), some of its settings are made proportional to the number of variables, n_v = 2 · n_c · n_f , i.e. phase and amplitude for each channel at each frequency. In particular, the number of particles is set to 10 · n_v , and the optimization is halted when n_v consecutive iterations have produced a relative change in the objective function smaller than 0.1. Other settings have been kept to their default values. When the particleswarm algorithm has completed, the optimization is handed over to a local optimizer (fmincon) and the solution is further refined until the relative parameter change falls below 0.01, thus reaching the desired precision of 1%.

SAR-based HTP quality indicators for the artifical test cases are reported in figure 8. The lowest HTQ for any given problem and frequency combination is achieved by PS, being configured to minimize this indicator. However, this does not always translate into better coverage. The i-TR returns sensibly higher TC₂₅ values in the large and medium cases, at the expense of an increase in HTQ with respect to the initial TR solution. While some solutions become impractical due to this increase, there is always at least one optimal frequency combination where i-TR achieves HTQ ≈ 1 and TC₂₅ ≈ 100%. In the small target case, i-TR also lowers the HTQ provided by classic TR while achieving remarkably higher TC₅₀ values. For this target, which is less heterogeneous than T₁ and T₂ (see table 1), the solutions provided by PS and i-TR are comparable in quality. EV exhibits consistently high HTQ and poor coverage in all cases.

Zoom In Zoom Out Reset image size

SAR-based HTP indicators for the realistic targets are reported in figure 9. Once more, PS returns the lowest HTQ for each case. All algorithms yield comparable HTQ values below 1 in the larynx, while only PS manages to reach values below 1 in the meninges. EV and classic TR perform similarly in the larynx case, achieving acceptable HTQ, but suboptimal coverage. EV fails to provide viable solutions for the meninges, with HTQ always above 1.5, and TC₂₅ above threshold only at 500 MHz and 700 MHz. In contrast, i-TR yields HTQ values much closer to the global PS optimum for this case. More importantly, it does so with simultaneous extensive target coverage. In the larynx model, TC₂₅ is often higher for i-TR than for PS when more than one frequency is considered. The SAR results for these two realistic cases are further summarized into two combined HCQ–TC₂₅ plots in figure 10. In both plots, the dual-frequency solutions provided by i-TR are seen to be located near the TC₂₅ optimum (100%) and close to PS's global HTQ optimum.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Speed performance metrics for the artificial test cases are reported in figure 11. Since EV and TR are direct solution methods, they take less than a second to complete. About 6000 cost function evaluations are needed for PS to return a single-frequency solution, and up to 30 000 in the case of two frequencies. This translates into long execution times, with a typical run taking about 10 minutes for one frequency and 50 minutes for two frequencies. Conversely, i-TR manages to return a single-frequency solution within at most 50 evaluations, and a double-frequency solution in 100 evaluations. This reflects into at most 10 s running time in the single-frequency case, and on average 50 seconds for a double-frequency problem. Performance metrics for the realistic cases are not reported since they indicate similar trends.

Zoom In Zoom Out Reset image size

The results of the sensitivity analysis with respect to the initial step factor α₀ are reported in figure 12. Both HTP quality indicators exhibit negligible variations across the sweep, meaning that the solution provided by i-TR is stable in this respect. A closer look suggests that the optimum is located at the lowest end of the α₀ range. This can be intuitively explained by the fact that smaller stepping allows the algorithm to get closer to the HCQ optimum. High α₀ values exponentially increase the convergence time and computational burden of the algorithm, indicating that the initial classic TR solution is already close to the optimum and only small steps are needed to improve its HCQ value. On the other hand, both t_r and n_e tend to increase again towards α₀ = 0, as a consequence of the step being too small to quickly reach the optimum HCQ. In conclusion, α₀ = 0.1 seems to be an appropriate value for this parameter.

Zoom In Zoom Out Reset image size

Figures 13 and 14 show the thermal distributions for the realistic cases. It can be seen that both direct solvers, EV and TR, fail to reach therapeutic temperatures in the deeper parts of the tumors. In the larynx case, the trachea acts as a heat barrier between opposing sides of the target volume, whose morphology is toroidal. PS and i-TR manage to overcome this barrier by extending the SAR coverage to the inner side, at the expense of higher average temperatures in the healthy tissues. This higher power deposition, however, does not result in additional HSs. i-TR performs better in this sense than PS by gaining half a degree in median temperature, see table 3. For the meningioma, both EV and TR end up maximizing the peak SAR in the target which results in a focal spot close to the skull. This spot is however mitigated by the cooling effect of the bolus, so that neither high peak temperatures nor sufficient coverage are achieved within the target. Both PS and i-TR shift this main SAR peak deeper in the brain, reaching high median temperatures. i-TR closely follows PS by only half a degree lower T₅₀. These results are further illustrated in the cumulative histograms, figure 15.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size