Hyperthermia phased arrays pre-treatment evaluation

Abstract

Purpose: In the hyperthermia treatment of deep-seated tumours by a phased array of radiofrequency (RF) antennas, heatability will be investigated in terms of power-to-tumour and other figures-of-merit of hyperthermia treatments to be optimised. The assumption is that each source is individually constrained to not exceed a maximal nominal power. The nominal power may differ from a source to another as a physical limit or an operative modality.

Method: Under such constraint, new procedures for the maximisation of (i) power-to-tumour, (ii) heating efficiency and, in general, (iii) power ratios as tumour-heating selectivity are proposed. (iv) The problem whether a tumour is equally heatable after turning off some antennas is addressed as array thinning.

Case study: An array of eight dipoles arranged on two lines around a head/neck is introduced to perform a numerical analysis. The achievable power-to-tumour according to the new optimizations and other performance indices adopted from the literature is tested against values of power that can be found to be sufficient for heating tumours to clinical temperatures. New solutions to data rendering in hyperthermia heating are proposed.

Keywords:

• Hyperthermia phased-array
• deep-seated tumors
• pre-treatment evaluation
• radiofrequency/microwave

Introduction

Clinical trials have shown that the addition of hyperthermia to radiotherapy or chemotherapy results in improved clinical outcome [1,2]. There is general consensus that effectiveness of hyperthermia technology is a prerequisite for clinical effectiveness [3,4]. An ideal hyperthermia treatment would yield a desired thermal dose distribution in the tumour [5,6] without heating normal tissues excessively. Solutions to this problem in terms of theoretical power patterns have been investigated [7].

To improve hyperthermia treatment quality, optimisation of power or temperature distribution in the tumour by hyperthermia treatment planning (HTP) is considered before or during treatment. Exhaustive reviews of the available techniques for HTP and specific simulation tools are in [8,9]. The optimisation procedures reported in the literature determine the source parameters for either temperature or specific absorption rate (SAR) pattern. Main motivations for the former start from the consideration that the human body is a complex thermodynamic system, characterised by blood perfusion and thermal conduction, whose behaviour may change during treatment as a consequence of thermoregulatory processes. As a result, even if power is focussed on the tumour, a preferential temperature rise is not guaranteed to occur there [10]. Furthermore, clinical response correlates with temperature elevation achieved, not power deposited [3,11]. These facts have motivated the development of treatment planning models and systems that directly optimise the temperature distribution inside the patient [12,13] or the thermal dose [14]. Online treatment-focussing algorithms, capable of dynamically adapting the model to uncertainties and changes in model parameters during treatment, would greatly enhance the clinical outcomes of hyperthermia. Feedback methods triggered by measurements of either the electrical field [15] or the temperature [16] at points within the patient body have been considered. Adaptive methods based on temperature feedback allowing more accurate settings received particular attention in conjunction with the proposal of non-invasive techniques for temperature measurement since the early times of hyperthermia and currently by magnetic resonance [17,18]. Using a feedback control loop to drive the phased array towards the optimal adjustment during the therapy circumvents the need for an exhaustive physical model by directly controlling the temperature during the therapy. However, system identification from measured data is necessary [19] and may take some time. Fast operation that is needed to compensate for sudden temperature changes in the patient motivated the elaboration of dedicated strategies [20,21]. Advancements in this field are reported in reference [22]. Alternatively, the SAR distribution, being proportional to the temperature increase at treatment start, was selected for optimisation [23,24]. As the robustness of electromagnetic solvers has increased significantly in recent years, the presently available codes provide excellent opportunities to perform SAR characterisation as a starting point for SAR optimisation [25]. In this study, heatability will be investigated in terms of power-to-tumour and other power-related figures-of-merit to be optimised.

Radiating systems that are currently used to heat deep-seated tumours include phased arrays of elementary radiators for which SAR-optimisation-based HTP is pursued by feed amplitude and phase setting. Finding the excitation setting for optimal tumour heating remains a difficult task because of the degrees of freedom available and because a large temperature elevation may occur at healthy tissue sites (hot spots) resulting in unwanted side effects, e.g. pain or healthy tissue damage [26]. In regional hyperthermia, the majority of treatment sessions is limited due to localised pain caused by hot spots [27]. Objective functions have been proposed for SAR optimisation according to various quality indicators, a comprehensive list of which is in Canters et al. [25]. They include maximum absolute SAR and SAR homogeneity in the target both to be optimised as well as hot spots to be minimised. The power deposited in the tumour as a fraction of the power released by the amplifiers (heating efficiency, η_H) or delivered to the patient’s healthy tissues (selectivity of tumour heating, η_Ω) are further indicators [23,28,29].

For given array, patient and tumour, assessing whether or not the power forwarded by the amplifier set suffices to heat the tumour is an important task taking into account (i) the maximal achievable power-to-tumour; (ii) various treatment-quality performance indicators; and (iii) the main causes for power loss instead of reaching the tumour, that are body and bolus absorption as well as mismatch, mutual coupling between antennas, and radiation leakage. A solution to this task, which is somehow different from HTP, logically precedes treatment planning. Only in the case of positive answer will the search for optimal feed settings be meaningful even in the case of adaptive optimisation.

To our knowledge, however, the whole task has not yet been undertaken along the general definition above. Moreover, single issues have not yet been addressed under the constraint that each amplifier does not overpass a maximal nominal power, which may differ from channel to channel, as is the case for modern phased arrays [30]. In this study, heatability will be theoretically investigated under this constraint. In particular, (i) a proof is introduced that the maximal power-to-tumour is obtained when each generator provides all its nominal power; (ii) a procedure for its numerical computation is given; (iii) the optimisation of η_H, η_Ω and other ratios between powers structured as quadratic forms through the solution of suitable eigenproblems [31,32] will be extended to include the constraints on the nominal powers; (iv) whether heatability can be preserved when one or more amplifiers are turned off according to some rules will be quantitatively investigated as array thinning. The next section on Theory is devoted to these contributions.

An application of these methods is proposed in the numerical analysis structured as a Case Study in this article. SAR optimisation schemes integrate the volumetric power density from antennas over domains for which the average SAR is determined, e.g. a tumour to be heated or an organ to be protected. In turn, hot spots and SAR homogeneity in the tumour are better evaluated from SAR distributions over the irradiated portion of the patient’s body and the tumour, respectively. In this study, a voxel-based analysis is used [33,34] as a trade-off between the former integral approach and the latter point-wise one. The analysis is performed for eight dipoles on two partially overlapping rings, positioned around a head/neck volume within a distilled-water bolus. Powers and performance indicators are computed and compared following three optimisation schemes. Cumulative frequency diagrams are proposed. The results are stored on a voxel basis for the head/neck allowing an exploration by slicing. Then, results for an extended target, i.e. a 60 ml sphere in central neck, will be carried out also in the presence of array thinning. Additional quality parameters related to tumour heating uniformity and hot spots are introduced and evaluated.

Symbolism

 real scalar

m, n, k integer spanning a set

 real or complex vector in a 3-dimensional (3D) physical space

[Q] M × M matrix whose complex entries are denoted either Q_mn or [Q]_mn, for m, n = 1, …, M 

 M – tuple of scalars a_m arranged as a column vector, its transpose conjugate

P a real power (W), further specified by subscript and/or superscript

η an efficiency, i.e. a quotient of powers, further specified by subscript and/or superscript

Nomenclature

j imaginary unit

 vector of an observation point

d_m distance between the centre of radiator m and a point of τ specified in the text

x, y, z Cartesian frame

τ, τ_i computation volume, sub-volume

Ω patient’s body portion belonging to τ 

T_k kth voxel ∈ Ω 

S tissue sphere (64 ml), whose centre is specified in the text

δ, δ³ grid spacing, volume of a voxel

δ_mn Kronecker symbol

M, N, K number of elements in an array, sub-array or volume

f, ω operation frequency, angular frequency

a_m , b_m complex amplitude of the forward and reflected wave at port m (W^1/2)

A_m, ϕ_m amplitude and phase of the forward wave at port m, i.e.  

[Q_Ω] influence matrix for a general target, kth target T_k, sphere S, and body Ω

ψ_mn phase of Q_mn

Q_{BP, BP} Q_mm for m = BP, i.e. the best performing radiator

 electric field at due to a_m = 1 (W^1/2), the other generators being kept off

 total electric field at

σ_e effective electrical conductivity, including ionic currents and dielectric losses

ε_r,ε₀ relative permittivity real part, permittivity of the vacuum

 power delivered to T_k in coherent (incoherent) modality

P_G, P_Gm net power from generators, power carried by mth-generator forward wave

P_N, P_Nm nominal power, nominal power of the mth-generator

η_H, η_Ω heating efficiency, tumour heating selectivity

Theory

The model

The computation volume, τ, includes the body portion Ω to be heated (the head/neck region in this study), the antennas, the surrounding bolus and free air. Ω is segmented into voxels. T_k denotes the kth voxel when k spans the set of voxels. T_k is heated by an array of M antennas that are accommodated around the body, as a consequence of power deposition by the radiated electromagnetic field. Location and orientation of each antenna are kept constant during a treatment session. Each antenna is connected to an amplifier through a line, typically a coaxial cable. Amplifiers are endowed with independent level and phase regulations. Modern amplifiers exhibit an impedance to the line that does not change with amplitude and phase setting. We assume that this property is owned by the system we deal with and that the amplifiers are matched to their lines. Therefore, from a circuital point of view an amplifier is modelled as a matched generator. Generator, line and antenna constitute a radiating channel.

The computation volume τ is modelled as a microwave junction whose M ports are the antenna terminals [35,36]. a_m will denote the complex amplitude of a direct wave at port m, where m spans the M channels. The available power at port m is given by , where and the asterisk is for conjugate. Amplitude A_m and phase ϕ_m can be adjusted for SAR optimisation. Finding the most effective phase-amplitude setting is a purpose of HTP.

P_Nm will denote the maximal available power (nominal power) at port m, therefore (1)

P_Nm is a physical limit of the amplifiers or an operative modality [30]. In this study, P_Nm may differ with the channel, so generalising a constraint on the maximal individual power that was introduced by Das et al. [37]. In [32] a steepest descent optimiser is used while the power per antenna can be a predefined maximum percentage of the total power. In this study, a power-to-tumour resulting from an optimisation process is said to be absolute if it is constrained as shown in Equation (1)(1) .

To shorten the forthcoming equations we shall use matrix notations so that the feed vector a is an M-dimensional column vector with complex elements a_m (forward wave amplitude in W^1/2). Let be the electric field due to the mth generator at , a point of T_k, when the other generators are off. Assuming linearity in the response, the total field at , when all generators are switched on, is obtained by superposition [38] as for general a. The power absorbed by T_k is [39] (2) where the dimensionless influence matrix [Q] has entries [31] . The dot is for the usual scalar product. is contributed by the mth generator when its direct wave carries P_Gm = 1 W with phase ϕ_m = 0, and the other generators are off. [Q] is Hermitian and positive definite. A similar start-up has been adopted in several studies, e.g. in references [33,40].

For later use, we define the Best Performing (BP) radiator as the one which, if individually excited, contributes the largest amount of power to the target. The BP radiator is characterised by the largest diagonal entry Q_mm, for 1 ≤ m ≤ M.

Phase-only optimisation

An optimisation problem aims at maximising power to target, , over only phases ϕ_m by letting all feed amplitudes A_m be invariant. The choice for any m will provide the largest P_k complying with (1). A proof of this statement, based on the positiveness of [Q], is in the Appendix. This choice will be referred to as Saturated Amplifier (SA) optimisation.

Optimised phases can be found by zeroing the gradient of P_k. Explicit Q_mn in modulus and argument as . The following system is obtained (3) where . This formulation extends a result in [34]. The system will be solved by a procedure making only use of trigonometric direct and inverse functions. To do this, the system is rewritten as (4) with and D_m = . We recognise Equation (4)(4) as fixed-point iterations. This means that, at the th iteration, is updated as a function of the unknowns that, in turn, received values at the kth iteration. This scheme gives rise to the sequences , , for any m, starting from a guess . An interesting feature of the algorithm is that the residual of (4) evaluated for the current phases coincides with the difference (absolute accuracy) between them and updates.

The maximal absolute power to the target is (5)

is the largest power to the target when the channel powers are individually constrained by (1) and the phases are the only operator’s choice.

The heating efficiency is the ratio of the power-to-target and the net available power from the generators given by . Its equation is (6)

Following the above SA scheme the heating efficiency is . The heating efficiency of the BP radiator, if individually excited, is .

Linear algebra optimisation

The maximisation of η_H and other quotients between powers, with relevance as figures-of-merit in hyperthermia, has been largely studied [25]. Here we are interested in the maximisation of η_H under the additional constraint (1). To our knowledge the following derivation is new. It includes the possibility of generalisation to the optimisation of other quotients. We take advantage from the spectral decomposition [41, p.253] of the Hermitian matrix [Q] which gives , where [D] is a diagonal matrix of eigenvalues and [V] is a unitary matrix whose columns are the corresponding eigenvectors. Such decomposition is a one-line instruction in MATLAB. The columns V_m of [V] are orthonormal so that , where δ_mn is the Kronecker symbol. Let λ and Vλ denote the largest eigenvalue and the corresponding eigenvector, respectively. The Rayleigh Ritz theorem [41, p. 267] states that the choice a = V_λ optimises P_k under the constraint . If we select a = sV_λ, with s a real constant, and substitute a into (2), we obtain . This is the largest power-to-tumour under the constraint P_G = s². We select s in such a way to have the generators not exceeding their individual nominal power. If V_{λ, m} denotes the mth component of Vλ, s is found as the following minimum (7)

The minimum is taken for over the ratios obtained by spanning m from 1 to M. In conclusion, (8) is the largest absolute power that the array constrained as shown by (1) can deliver to T_k under heating efficiency optimisation. This result extends a result in [37] to possibly different individual constraints. The power from the generators is also scaled by the same factor, i.e. . Therefore, the heating efficiency is unchanged by scaling and is equal to λ. The above scheme will be referred to in the following text as linear algebra optimisation (LAO) of η_H.

Note that maximises the power-to-target (5) for the phase-setting solution to (3), i.e. . In turn, maximises the heating efficiency for a_m = sV_{λ, m}, for m = 1, …, M, i.e. a different feed setting. In other words, no combination of amplitudes and phases in agreement with (1) can deliver more power to the tumour, in case of , or attain a larger heating efficiency, in case of . Consequently, , while , i.e. the power delivered to a tumour under phase-only optimisation is larger than that of linear algebra optimisation at the cost, however, of worsening the heating efficiency.

From Equation (7)(7) , the smaller is, the larger the absolute power to the target. Indeed, the array effect is more important if the channels contribute comparable powers to the target. This is an important result. In the numerical analysis we shall see that the voxels located centrally in the head/neck benefit from this remark, while the voxels close to the body surface do not, being mainly heated by the dipole just in front of them. Consequently, a low scale factor, s, is obtained for superficial voxels. This behaviour is well represented by an array factor [42].

Many authors have formulated the optimisation problem as extremization of a performance index structured as ratio between two Hermitian quadratic forms (Generalized Rayleigh Quotient, GRQ). A general form of this ratio has been proposed by Paulides et al. [43] as shown in the next equation with the symbols used in this study: (9)

If the computation volume τ is partitioned into subvolumes τ_k, then is the power delivered to τ_k that is characterised by the influence matrix [Q_k]. The numerator is a weighted sum of powers with weights w_k. The divisor has a similar structure with different weights w_i and matrices [Q_i]. An important application of (9) is the selectivity of tumour heating given by the ratio (10) between power-to-tumour (2) and power-to-healthy tissue . The optimisation of η_Ω over amplitudes and phases has been investigated by e.g. [33,44]. In Kowalsky et al. [10], only the phases are allowed to vary and the maximum is found by a conjugate gradient algorithm.

The maximisation of η_Ω over the set of feed vectors is equivalent to finding the largest eigenvalue μ of ([41] p. 285). An extension of LAO to η_Ω and other efficiencies amenable to a GRQ takes care that is not Hermitian. Therefore, a new normalisation of the eigenvectors is introduced as (11)

A call to the routine “eig” of MATLAB provides a direct solution to this problem. Denote V_μ a normalised eigenvector corresponding to the largest eigenvalue μ. If we select , the power-to-tumour is obtained as from (10), (11). The scale factor s is chosen as shown by (7), with as divisor.

For brevity, the objective function is maximised in the numerical case study with an equal constraint on the individual source power, i.e. P_Nm = P_N for any generator m.

Antenna ordering and array thinning

Array thinning, i.e. partial array activation, is the practice of keeping off some generators due to the assumption that a subset of the available radiators suffices to heat a target [45–47]. In [48], turning off unneeded elements during treatment is included in adjusting the element amplitudes to achieve the desired temperature distribution.

Methodologically, off generators are matched loads for the feeding lines. Off dipoles behave like parasitic antennas that partially re-irradiate the field they receive. In this study, array thinning is considered in the attempt to answer the following questions that are naturally posed: (i) whether a subset of the available antennas individually constrained by (1) can deliver enough power to the target, (ii) which is a suitable choice for this subset, (iii) whether there is a price to pay in terms of quality of heating, and (iv) in case of positive answer to the last question, what is the price.

Different rules for antenna turning off can be devised. The method followed in this study starts with a rule for antenna sorting, from the best antenna to the worst one according to a performance parameter. Equivalently, a hierarchy of antennas is introduced which may change according to the parameter adopted for ordering. Then, the antennas are turned off step-by-step starting from the last one and some evaluations are performed. Three rules for sorting will be compared in the case study: (i) the dipoles are sorted in descending order according to the amplitude of the on-diagonal entries of the influence matrix; (ii) an intuitive rule is based on the distances between the antenna centres and a target; (iii) a more elaborate rule finds the subset of N radiators maximising the heating efficiency for any fixed N < M.

A case study

A numerical analysis was performed as a benchmark for the optimisation schemes introduced in the section on theory. The section is organised in subsections. The first one presents the configuration of the case study. The other subsections are devoted to independent numerical experiments, i.e. a voxel-based analysis, an extended target analysis, and finally to an array thinning one.

Configuration

Heatability was evaluated in the context of a simulated head/neck treatment by a phased array. The frequency of operation was set to 434 MHz. M = 8 antennas, made by a half-wave dipole backed by a rectangular metallic reflector (Figure 1), were positioned around the head (Figure 2). A deionised-water bolus surrounds the head allowing skin cooling and wave-impedance matching to tissues. The bolus has a cylindrical shape that does not include nose or eyes. The bolus cross-section is a circle cut by a chord on the nape side (Figure 3).

Figure 1. Two-element Yagi-Uda antenna immersed in water, operating at 434 MHz

Figure 2. Four antennas located on each side of the head. Adult male phantom by CST, water bolus, a Cartesian frame and volume τ also shown.

Figure 3. Bolus cross section.

Computations were performed with CST Microwave Studio^® (Darmstadt, Germany) suite, which uses Finite Integration Technique (FIT) and includes the anthropomorphic phantom “Gustav”. Bolus, surrounding air and human tissues are segmented by the solver into 2.08 × 2.08 × 2 mm³ elements. A total number of about 3.7 million of elements were used. The complex permittivity is constant within each element and takes values from a database also provided by the solver. ε₀ is vacuum permittivity and ω = 2π f, with f operation frequency. A list of head/neck tissues likely involved in the simulations is in Table 1. The values reported in Table 1 match those in the literature [49]. In turn, the values for the deionised bolus-filling water match those in [50] at 20 °C.

Table 1. ε_r and σ_e for head/neck tissues and distilled water.

	εr	σe (S/m)
Bone	13.07	0.095
Brain	49.24	0.603
Eye bulb	69.00	1.533
Eye lens	47.96	0.675
Fat	5.57	0.042
Lung	23.58	0.380
Muscle	58.87	0.805
Oesophagus	68.19	1.013
Skin	46.06	0.702
Spinal cord	35.04	0.456
Teeth	13.07	0.094
Thyroid	61.33	0.886
Trachea	43.93	0.644
Water	80.00	0.040

The time required to solve for each excitation is about 7 min on a computer running a 12 thread CPU (INTEL X5650 at 2.67 GHz) and 1 GPU (NVIDIA Tesla 20). To speed up the algebraic computations that follow, the results were archived for Gustav’s intersection Ω with a smaller box τ shown in Figure 2. A Cartesian frame also shown in Figure 2 was introduced in τ. The electric field vector and the power of the Joule effect were saved at the nodes of a regular 3D grid having side δ = 2 mm. 102 × 121 × 76 nodes were located along x, y, z, respectively. In the following analysis the voxels are cubes having side δ with centres at the nodes and sides parallel to the axes. MATLAB^TM (Natick, MA) was used for algebraic operations on the field data provided by the solver. The antenna labels are shown in Figure 4.

Figure 4. Antenna labelling.

Voxel-based analysis

A voxel-based analysis (VBA) is here preferred to a point-wise one. With the performance indices being defined with reference to integrals over volumes, a voxel is suitable for their evaluation. In a voxel-based analysis, the target T_k is assumed with the size of one voxel and spanned within the head/neck segment Ω, intersection of the box τ and Gustav. The graphic rendering problem for 3D results was dealt with by providing head/neck cuts that are parallel to the box faces. Numerical results will be presented for 1W of source nominal power.

Some thresholds were proposed for an optimal treatment. It was found that a SAR in the range 20–50 W/kg is needed in the target region to reach therapeutic temperatures [26,51]. To take these thresholds into account when a voxel-based analysis of heatability is performed, the power delivered by the array to a δ³ = 8mm³ voxel will be compared with P_l = 1.6 10 ^− 6 W and P_h=4 10 ^− 6 W. These values for P_l and P_h were obtained considering an average tissue mass density of 1 g/cm³ while the SAR of 20–50 W/kg was referred to a presumed nominal power of 100 W and suitably scaled to P_N = 1W.

The computation of the influence matrix [Q_k] was repeated for every voxel T_k belonging to Ω. Preliminarily, heatability was analysed by computing incoherent array heating modelled by . A cumulative frequency (CF) diagram of is shown in Figure 5 (dotted line). The other diagrams in Figure 5 refer to deliverable absolute powers following LA and SA optimizations, respectively. They will be discussed later. A CF diagram as those in Figure 5 shows on the vertical axis the number K of voxels in τ (alternatively, their net volume Kδ³) for which the power released by the array is ≤P. The condition P_Q > P_l (P_Q > P_h) is unfulfilled by 5.7% (33.1%) of the voxels, respectively, i.e. by 158 ml (925 ml) of tissue. If the computation box is enlarged these volumes increase because τ will include additional poorly irradiated voxels, and viceversa.

Figure 5. Volume (ml) for which P_k < P vs. P where P_k is maximal absolute power to voxel T_k. P_N = 1W nominal power. Thresholds P_l and P_h are also shown. Incoherent heating (dotted line), LAO of η_Ω and η_H (black and grey solid lines, respectively), and SAO (dashed line).

CF diagrams of are reported in Figure 5 (solid line) for the feed settings obtained from LA optimisation of η_H, η_Ω. They show expected improved absolute power to the targets in comparison with incoherent heating. Indeed, the condition is unfulfilled by 1.9% (2.2%) of voxels, i.e. by 52 ml (62 ml) of tissue, in the case of LAO of η_H (η_Ω). In turn, the stronger condition is unfulfilled by 12.4% (14%) of voxels, i.e. by 346 ml (389 ml) of tissue, for η_H (η_Ω) optimisation. If the source nominal power is doubled the poor-heated percents decrease to 0.25 (0.32) and 3.0 (3.5), respectively.

The above steps were repeated for the phase-only optimisation. The fixed-point procedure was started by the phases of the LAO of η_H as guess, and was stopped after a number of iterations that were sufficient to attain an absolute accuracy better than ± 10⁻⁵. On average, 8–9 iterations were sufficient. The number of voxels in the head/neck being 348 522, the total time required by SAO for all the voxels was 949s, and the mean time for a single voxel optimisation was 2.7 ms.

Absolute power-to-target, , is obtained from Equation (5)(5) . A CF diagram of is reported in Figure 5 (dashed line) showing expected improved absolute power to the targets in comparison with linear algebra optimisation and, a fortiori, incoherent heating. The condition () is unfulfilled by 0.07% (1.6%) of voxels.

The above optimizations can be compared on the basis of maps of performance indices such as P_k, η_H, η_Ω, and an array factor, AF, that was introduced in [42] to evaluate the array effect. AF is defined as the following quotient of heating efficiencies (12) where is the heating efficiency of the one-element array, i.e. the BP one. The power delivered to the target by the phased array can be equally delivered by the BP channel alone if, however, its power is raised to AF times the power P_G of the array.

In order to map the results of the voxel-based analysis, for each T_k we store the feed vector that optimises P_k (or η_H or η_Ω) to that voxel. Then is substituted into the equations for the various parameters (power, efficiencies and array factor), generating scalar functions of the voxels. This procedure generates virtual solids having the head/neck shape, whose voxels store the values of the parameters. After establishing a correspondence between scalars and colours, 2-D colour maps are finally obtained as an intersection of these solids by orthogonal planes passing for a head/neck central point (principal planes). For brevity, only one example relating to η_Ω optimisation is presented in Figure 6, where colour maps of P_k are shown. They show the largest power deliverable to each voxel when optimality is searched specifically for that voxel. A logarithmic scale is used to smooth the spread of P_k over the head/neck. Due to a larger distance from the dipoles, the voxels that can be heated less are mainly located in the upper head. Due to lower values of the conductivity σ_e, they can also be found in trachea, spine, maxillary bones/teeth, and air ducts.

Figure 6. in μW, on the sagittal plane x = 55δ (left) and on the transversal plane z = 25δ (right), as result of LAO of η_Ω. P_N = 1W. Axis unit is δ = 2 mm.

Extended target analysis

During a second numerical experiment the target was a sphere, S, having a volume of 65.74 ml. The sphere consists of all voxels having their own centre at a distance no greater than 2.5 cm from the point 55δ, 53δ, 25δ. The sphere is strongly inhomogeneous in its dielectric properties. On the basis of the correspondences in Table 1 between tissues and σ_e, the most common tissues in S can be inferred as muscle, tongue, trachea and glands. The volume filled by air inside the sphere is 4.9 ml. The influence matrix for the whole sphere is denoted [Q_S]. It is obtained as , where [Q_k] is the influence matrix of the kth voxel belonging to the sphere and the sum is over all such voxels. If a target is the sum of parts, the heating efficiency by LAO of the whole target is upper bounded (i.e. ≤) by the sum of the heating efficiencies of its parts. This is a consequence of the property that the largest eigenvalue of the sum of Hermitian matrices is not larger than the sum of the respective largest eigenvalues ([41], p. 274). So some decrease in the maximal power deliverable to an extended target is expected in comparison with the sum of the maximal powers deliverable to its parts. The BP radiator is 3. Results are reported for nominal power P_Nm = P_N = 1W for all channels.

The feed vector obtained by LAO of η_H has the following modules and phases

In turn, the phases obtained by SAO with guess the above LAO phase vector are:

The fixed-point algorithm was used. It was stopped after 8 iterations, that were sufficient to attain an absolute accuracy better than ± 10⁻⁵. The largest difference between the phases following LAO and SAO is less than 14° (<8.7○ averagely). The closeness of the solutions can explain the good performance of LAO as guess for SAO. A hybrid feed is that when the amplifiers are saturated while the phases are those of LAO instead of those of SAO. Interestingly, this choice does not penalise the power-to-tumour to a large extent. Indeed, the relative difference in the power delivered to the sphere-target is <3%.

Finally, the LAO of η_Ω gave the following results

In a comparison with the feeds optimised for η_H, they show major changes in the phases while all the amplitudes of the non-saturated amplifiers more or less decrease.

Further results are reported in Table 2. To complete the analysis of heatability, two parameters to evaluate intratumoral heating uniformity and hot spot bearing are given in addition to those already introduced. As far as intratumoral uniformity is concerned, some results of a statistical analysis of the measured SAR at points within the tumours were shown by a histogram of the SAR distribution [25]. A non-normal distribution was found. Intratumoral temperature uniformity was evaluated in [52] through CF diagrams of the difference T₁₀ − T₉₀, where T₉₀ (or T₁₀) is the temperature above baseline temperature achieved by at least 90% (or 10%) of the voxels. In this study a similar definition is used in terms of power. P₉₀ (or P₁₀) is the power at least achieved by 90% (or 10%) of the voxels. The ratio (P₁₀ − P₉₀)/P₉₀ is an estimate of heating uniformity. The smaller the ratio is the greater is the uniformity.

Table 2. Sphere optimisation parameters (best result per row in bold).

	Symbol (unit)	LAO of ηH	LAO of ηΩ	SAO
Absolute power to the sphere	P^LA, P^SA(W)	0.063 PN	0.046 PN	0.122PN
Efficiency	ηH (mW/W)	17.4	16.0	15.3
	ηΩ(mW/W)	37.4	40.7	30.4
Array factor	AF	2.04	2.81	1.8
Power from generators	PG(W)	3.64 PN	2.84 PN	8 PN
Percentile ratio	(P10 − P90)/P90	3.4	3.9	1.8
Hot spot volume 3x	HSV 3x (ml)	39.6	31.3	49.5

Hot spots are treatment limiting factors. A discussion about hot spots and their delineation is not a purpose of this study. Following Wiersma et al. [33], a voxel belongs to a hot spot if the power delivered to it overtakes n times the power to the voxel target. The total volume of such voxels constitutes the Hot Spot Volume (HSV). The voxels close to the body surface, i.e. at depth <1 cm, are not included in HSV in the assumption that they are well managed by surface cooling. Here the hot spot volume for n = 3 is a further parameter for comparisons of different heating modalities.

In Table 2, the best result for each parameter is put in bold. In a comparison of the optimisation modalities, SAO provides 2–3 times the absolute power of LAO, at the cost of slightly less heating efficiency and larger hot spot volumes. The heating uniformity evaluated by the percentile ratio is better for SAO. Averaging on the whole sphere, the heatability condition SAR >50 W/kg gives at least 3 W for S. This condition is fulfilled by a set of ∼25, 50, or 65 W amplifiers in the case of SAO, η_H − LAO, or η_Ω − LAO.

The powers delivered following SA and LA optimizations are shown in some detail in Table 3 for P_N = 1W. The net power radiated by the antennas, P_V, differs from that from sources for the power returned to the generators by mismatch or mutual coupling. It is obtained from P_G through a feed efficiency [42]. P_V is split as P_V = P_Ω + P_R, where P_R is the sum of power to bolus and that radiated to environment. P_Ω e P_R are comparable. The powers dissipated in the dipoles, connectors and cables are not included in the balance. A fair sensitivity to the optimisation to which the data in Table 2 refer can be appreciated. All the powers decrease passing from LAO of η_H to that of η_Ω. A decrease in P_Ω is accompanied by a corresponding decrease in the power to the sphere. The power of 3W, that has been assumed to heat adequately the target S, entail 75.7 or 102.1 W to the body according to the LAO of η_Ω or SAO, respectively. In comparison with the latter, the former optimisation allows delivering about 25% less power to healthy tissue.

Table 3. Powers (in W) resulting from central sphere optimisation and P_N = 1W.

	PG (sources)	PV(antennas)	PR (bolus & leak)	PΩ (body)	PS (sphere)
LAO of ηH	3.64	3.21	1.46	1.75	0.063
LAO of ηΩ	2.84	2.34	1.18	1.16	0.046
SAO	8.00	7.30	3.15	4.15	0.122

Array thinning

Heatability will now be discussed in the presence of array thinning. In the first instance, the dipoles are sorted in descending order according to the amplitude of the on-diagonal entries of [Q_S]. The resulting order is 3 6 1 7 4 2 8 5 for the central sphere, where the correspondence between numbers and dipoles is shown in Figure 2. The thinning is obtained by cancelling the last digit at each step. Thus, the SA and LA optimizations are performed for N lighted dipoles. For brevity, only the LAO of heating efficiency is discussed. Some results are summarised in Table 4. The BP dipole is 3, to which the last column refers when it is the only one fed. The absolute power to the sphere as well as and AF^LA decrease monotonically with N. Being [Q_S] a Hermitian matrix this behaviour can be predicted ([41], p. 269). Also and, trivially, decrease monotonically with N. The other parameters also basically decrease showing, however, small local maxima with the exception of the hot spot volumes which first oscillate and then take large values for N ≤ 2. The last result is evidence that a phased array advantage compared to a single radiator is found in hot spot reduction. Moreover, having available 100 W amplifiers, the power to the sphere overtakes the minimal one, i.e. 3W, for N ≥ 4 or N ≥ 3 in the case of LA or SA optimisation, respectively. High values of the ratios (P₁₀ − P₉₀)/P₉₀ have been introduced as indicators of heating inhomogeneity. These ratios increase when the number N of lighted dipoles decreases, evidencing a worsening of homogeneity. Moreover, , so that SA optimisation seems preferable from this point of view. We conclude this example observing that the central sphere can still be heated when only the first 3–5 dipoles in the sequence of the best performing ones are lighted even if almost all quality parameters worsen with N decreasing.

Table 4. Results of thinning for central sphere and LA of η_H & SA optimizations.

Lighted dipoles	N	8	7	6	5	4	3	2	1
Absolute power to the sphere × PN (W)		0.063	0.060	0.049	0.046	0.039	0.026	0.020	8.5× 10^−3
		0.122	0.112	0.096	0.076	0.053	0.037	0.020
Heating efficiency (mW/W)		17.3	17.1	16.8	15.7	13.6	12.4	10.1	8.5
		15.3	16.0	16.0	15.2	13.3	12.3	10.0
Array factor	AF^LA	2.04	2.01	1.97	1.84	1.60	1.46	1.18	1
	AF^SA	1.80	1.89	1.87	1.79	1.56	1.44	1.18
Percentile ratio (P10 − P90)/P90	by LAO	3.39	3.70	5.60	4.73	2.27	10.9	7.97	8.8
	by SAO	1.79	1.89	2.87	2.5	2.12	7.01	7.7
Hot spot volume 3x (ml)	HSV^LA	39.6	39.3	43.1	40.0	44.0	37.0	67.7	121
	HSV^SA	49.5	45.9	49.2	44.1	48.2	35.5	67.7

A rule to define a minimal number of channels to keep lighted may have interest in applications. The array factor is an estimate of the advantage of a phased array over BP antenna performance in terms of the ratio of heating efficiencies. From Table 4, AF^SA does not change significantly for N ≥ 5. From this remark it follows that a rule based on less change in array factor with N increasing is suitable.

Other rules for thinning can be considered and examples are now given. A more elaborate rule finds the subset of N radiators maximising a performance index for any fixed N. If heating efficiency is such an index, the largest eigenvalue is computed for every principal N × N submatrix of [Q_S]. A principal submatrix corresponds to a subset of N lighted radiators. Any subset of lighted radiators has its own principal submatrix, for which the largest eigenvalue and the heating efficiency can be computed with the procedure used for the whole [Q_S]. Then, a heating efficiency is associated to each subset, allowing for a comparison among subsets and the selection of the optimal subset of radiators. According to this rule the resulting order for the central sphere is 3 4 1 2 6 7 8 5.

An intuitive rule for dipole sorting is based on the distances d_m (m = 1, …, M) between the centres of the dipoles and S. The distances are as follow: d_m = 74.7 81.7 81.7 102.5 112.7 93.2 92.8 85.2 (mm). The resulting order is 1 2 3 8 7 6 4 5. The sequences of lighted dipoles obtained by the three thinning rules differ by, at most, 2 elements. The heating efficiencies are presented in Table 5. The number of lighted dipoles being equal, the difference in heating efficiency is small between the two first rules, increasing when the distance-based rule is adopted for thinning. In this example, a sorting based on the influence matrix diagonal entries is a trade-off between efficiency and simplicity.

Table 5. Heating efficiency (mW/W) for different rules for array thinning.

Rule	N = 7	N = 6	N = 5	N = 4	N = 3	N = 2
Best performing	17.15	16.79	15.70	13.64	12.40	10.09
Best principal submatrix	17.15	16.79	16.04	15.26	13.71	11.14
Decreasing distance	17.15	15.10	12.77	12.37	12.37	7.11

The effects due to thinning were investigated with reference to other spherical targets having same size but differently centred in the head/neck. Similar results were obtained, not reported here for brevity.

Discussion

In this study, the heatability of a tumour by a phased array means ability to deliver enough power to the tumour, avoiding power excess to healthy tissues. An investigation of heatability requires the computation of the maximal power that the system can deliver to the tumour taking into account the physical limits of the system or operative modalities, and other constraints. Healthy tissue overheating is opposed by reducing total heat to the body (systemic heating) and contrasting hot spots individually. In the literature, a trade-off between power from generators to tumour, on one hand, and to healthy tissue, on the other, is commonly reached through the maximisation of some quotients structured as RRQ or GRQ. Heating efficiency η_H (Equation (6)(6) ) and selectivity of tumour heating η_Ω (Equation (10)(10) ) are largely used for the purpose.

Patient’s body, tumour, phased-array system (including bolus and operation frequency), and relative position of antennas and body are already set when this investigation on heatability starts. It is further assumed that some preliminary steps of HTP have been performed. In particular, a dielectric map of the tissues for a body segment including the tumour to be heated is available. The influence matrices [Q] have been computed for the volumes to be heated and for those to be protected, according to (9) [43]. It should be noted that other start-ups have been explored [51] that do not exploit influence matrices. Using [Q] allows a separation of electromagnetic computations from those strictly related to optimisation. Moreover, some properties of achievable optimizations can be directly evaluated on the entries of [Q] [42]. A hierarchy of antennas is possible by ordering the antennas according to Q_mm decreasing. In particular, a best performing antenna can be singled out. Each target to be heated determines its hierarchy, which may change if the relative position of patient and array changes.

The only degrees-of-freedom considered in this study are feed parameters, i.e. amplitude and phase for each channel. Three optimisation strategies have been explored under the common constraint of a maximal output power per channel, P_Nm, which may differ from an amplifier to another. (i) The first strategy optimises the power-to-target subject to this only constraint (SAO); (ii) the second strategy maximises the power-to-tumour ratio to the total array power (LAO of η_H); (iii) the third one maximises the power-to-tumour ratio to the power delivered to another volume, e.g. healthy tissue (LAO of η_Ω). In the literature, problem (ii) has been addressed with an equal constraint on the nominal powers for all the channels [37]. Problem (iii) was solved by search algorithms for fixed values of the power per channel [10,32].

A first nontrivial result of SAO is that the largest power-to-tumour is obtained when all the amplifiers forward their nominal power, thus reducing the problem of parameter setting to that on the phases only. The latter is formalised as finding the roots of a system of nonlinear equations. This is performed iteratively starting from a guess by a fixed-point algorithm introduced in this study. The electromagnetic power delivered to the tumour according to SAO is the largest one for given patient, tumour, array and maximal power per channel. The maximal absolute power following an SA optimisation is larger (typically twice or more) than that resulting from an LA one. The price to pay is a worsening of the heating efficiency. Proofs of this statement have been given.

A target can be heated by a phased array if the maximal absolute power obtained by SAO and/or LAO overtakes selected comparison powers. At a second stage, different heating modalities that comply with this condition are compared on the basis of other tests related to hot spot generation and heating homogeneity. In the case study, related methods are taken or adapted from the literature [33,52]. Hot-spot relevance in a treatment is estimated by the body segment volume where power dissipation overtakes n times the average power to the target. Heating homogeneity is evaluated through the difference P₁₀ − P₉₀ of the percentile powers received by the most and least irradiated portions of the target normalised to P₉₀. The complex of these tests is a possible definition of heatability.

If a tumour has been found heatable by a phased array of M radiators, the question whether it is still heatable by a subset of N < M antennas is naturally posed. This question introduces the matter of array thinning that has been addressed in this study by numerical evaluation and comparison of possible rules. Array thinning has been investigated by (i) organising a hierarchy from the best performing radiator to the worst performing one; (ii) turning off the worst performing at each step of an iterative procedure; (iii) repeating the tests of heatability on the subset of lighted radiators. In the numerical analysis we have observed a progressive decay with thinning of both the power to tumour and the other quality parameters. The next question is whether a parameter may suggest the minimal number of antennas to keep lighted. The results of the numerical analysis suggest the number N for which the array factor AF^SA reaches a nearly stationary value as suitable to this purpose. Additional numerical investigation is necessary to confirm this conclusion.

A review of possible applications follows. When patient, tumour and phased array are given, the proposed methods can be exploited before HTP and treatment, the computational effort in addition to that required by HTP being minimal. They provide information on the largest power deliverable to the tumour and achievable performance indices, such as efficiencies, SAR homogeneity and localised SAR hot spots. A comparison with archived data may have a predictive value. In the case of adaptive hyperthermia, the phase/amplitude settings from SAO or LAO may be used when a treatment starts. Additional degrees of freedom, like adjustments in the relative position of patient and array, can be explored based on the influence matrix pertaining to each new arrangement.

A different application of the above methods is when the heating potential of a given phased array is investigated for a specific organ/tumour-location without referring to a particular patient or lesion. The question may be whether the phased array is able to heat a tumour localised somewhere in that organ. An actual lesion being the union of voxels belonging to the organ and neighbouring tissue, a voxel-based analysis of heatability appears to be the only viable one. The VBA rationale is that an extended target is the union of a subset of voxels so that its heatability is strictly related to the heatability of the composing voxels. In addition, the number of voxels being large enough in an extended target, some statistics can be drawn. However, the transfer of results from composing voxels to a tumour is trivial only if each voxel is not heatable. The opposite is not true a priori. Indeed, the feed setting that optimises the power to a voxel is only suboptimal for a close voxel. Consequently the average SAR to a subset of voxels is less than the average of the SAR values separately optimised for the composing voxels. In the case study an analysis based on percentiles P₁₀, P₉₀ is proposed to achieve some conclusion. The matter deserves further investigation in future work.

Conclusion

Methods that can be used in hyperthermia pre-planning for heatability evaluation have been addressed. The problem of finding the maximal power that a phased array can deliver to a tumour under a constraint on the maximal available power per channel has been investigated for comparison with assumed therapeutic doses.

Disclosure statement

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

Appendix

From Equation (2)(2) in the main text, the power to the kth voxel is given by . We let be the mn-entry of a Hermitian matrix [R]. The available power of each generator is constrained by (1) to be not larger than P_Nn. We search for a stationary point of P_k in the set for any n. Zeroing the first derivatives with respect to A_n we get for the A_n of a stationary point. The above homogeneous system of M equations in the unknowns A_n admits nontrivial solutions, i.e. different from the all-zero one, only if . However, since is positive definite, its determinant is strictly > 0, and we conclude that extremal values of P_k can only be found on the boundary, i.e. for , for any n.

is positive definite if [R] is such. It remains to show that [R] is positive definite, equivalently that its smallest eigenvalue is strictly positive. Let , for n = 1, …, M, be the entries of a vector ϑ, and define , so that , where ○ is for the Hadamard product. [Φ] is a correlation Hermitian positive-semidefinite matrix ([41], p. 243). Now, the positive definiteness of [R] is a consequence of that of [Q]. Equivalently, a theorem at p. 274 of [41] can be used.