Ultimate hyperthermia: Computation of the best achievable radio-frequency hyperthermia treatments in non-uniform body models
Bastien Guerin1,2, Jorge F. Villena3, Athanasios G. Polimeridis4, Elfar Adalsteinsson5,6,7, Luca Daniel5, Jacob K. White5, Bruce R. Rosen1,2,6, and Lawrence L. Wald1,2,6

1A. A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Charlestown, MA, United States, 2Harvard Medical School, Boston, MA, United States, 3Cadence Design Systems, Feldkirchen, Germany, 4Skolkovo Institute of Science and Technology, Moscow, Russian Federation, 5Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA, United States, 6Harvard-MIT Division of Health Sciences Technology, Cambridge, MA, United States, 7Institute for Medical Engineering and Science, Massachusetts Institute of Technology, Cambridge, MA, United States


We propose a framework for the computation of the ultimate hyperthermia, which is the best possible hyperthermia treatment for a given frequency and non-uniform body model achievable by any multi-channel hyperthermia coil. We compute the ultimate hyperthermia treatment of two shallow (close to skull) and deep (close to ventricle) brain tumors in the realistic “Duke” body model and for treatment frequencies ranging from 64 MHz to 600 MHz. We characterize the convergence to the ultimate SAR pattern as well as temperature increase associated with the ultimate SAR distribution in the presence of non-uniform perfusion effects.

Target audience

RF engineers, oncologists, hyperthermia physicists.


To compute the best achievable RF hyperthermia treatment performance using sources outside the head (“ultimate hyperthermia”) in realistic body models with non-uniform electrical and thermal properties and non-uniform tissue perfusion.


Electromagnetic basis set: Our first step in computing the optimal external RF hyperthermia treatment is to obtain a basis set for electromagnetic (EM) fields in a realistic non-uniform body model. We have shown previously how to compute such a numerical basis [1] and used this approach to compute the ultimate signal-to-noise ratio in MRI [2]. First, we place a large number of dipoles (>105) around the body model at a distance no less than D=3 cm (Fig. 1). Second, we compute the incident E and H fields (i.e., in free space without the body model) created by random excitations of the dipole cloud. This is done by (i) exciting all the dipoles simultaneously using random excitations and (ii) computing the resulting E and H fields using the free-space Maxwell Green function. Finally, we compute the scattered E and H fields in the body model using the fast volume integral equation EM solver “MARIE” (https://github.com/thanospol/MARIE) [3]. Although bases of the solutions of Maxwell equations are strictly speaking infinite, we truncate the infinite basis to a finite number of basis vectors which converges to the ultimate with accuracy better than 1%. We model the Virtual Family “Duke” body model (head only, the body was cut off below the neck) [4] with a resolution of 3 mm isotropic (63×80×85 pixels), to which we added two tumors (one shallow, i.e. close to the skull, and one deep in the brain close to a ventricle). We compute EM basis sets in Duke at 64, 128, 297 and 600 MHz. Treatment optimization: We optimize the weights assigned to each basis vector (equivalent to voltage levels in a multi-channel hyperthermia coil) by maximization of the average specific absorption rate (SAR) inside the tumor while enforcing that the average SAR outside the tumor is below 3 W/kg. This is a convex constrained optimization problem that we solve using Matlab’s function “fmincon” using analytical expressions of the first and second derivatives of the objective and constraint functions for robust convergence. Temperature simulation: We solve Pennes bio-heat equation using the non-uniform thermal conductivity, heat capacity and perfusion maps of “Duke” using a Crank-Nicholson finite-difference time domain scheme with a time step computed using the Von Neuman stability criterion (code available at www.martinos.org/~guerin) [5]. Thermal parameter maps were obtained from the Gabriel database [6]. Perfusion maps were obtained from the Virtual Family database, except for tumors which were each decomposed in their necrotic core (perfusion 10% that of grey matter) and well-perfused penumbra (perfusion 200% that of grey matter) as shown in Fig. 3.


Fig. 2 shows that convergence of the hyperthermia treatment efficacy (“SAR amplification ratio”) is achieved using ~500 basis vectors for all frequencies and both tumor depths. Roughly 230 and 100 basis vectors are needed to achieve 80% of the ultimate for the shallow and deep tumors, respectively. This suggests that, in practice, multi-channel hyperthermia coils with as many channels are needed to achieve results similar to the ultimate. The treatment frequency seems to affect the specificity of the ultimate treatment for the deep tumor only. This is likely because complex E fields interferences are required to “focus” the energy deposited in the deep tumor, which is more easily achieved at high frequency (=small wavelength). We point out that it is likely that the optimal treatment efficacy drops above a certain frequency however because of the shallower skin depth with increasing frequency. Fig. 3 shows SAR maps of the ultimate treatment at 600 MHz. Although the treatment energy is correctly focused in the tumor, there is considerable “spillover” to neighboring tissues. Note that we constrained the SAR outside the tumor to be below 3 W/kg, however the peak local SAR can – and obviously is – much greater in healthy tissues. For the deep tumor, the horn of the ventricle was significantly heated because of its proximity to the tumor and the high electrical conductivity of CSF. This suggests that an explicit SAR constraint corresponding to the ventricle volume should be added to the optimization algorithm. Temperature simulations in Fig. 4 reflect the SAR results of Fig. 3 and show a temperature amplification factor of x18 and x1.7 in the shallow and deep tumors, respectively, compared to the average healthy tissues.


R01EB006847, P41EB015896, K99EB019482


[1] Villena JF (2014). ISMRM 22:623. [2] Guerin B (2014). ISMRM 22:5125. [3] Polimeridis A. (2013). JCP 269:280-296. [4] Christ A. (2010). PMB 55(2): N23. [5] Karaa S (2005). MCS 68(4): 375-388. [6] Gabriel C (1996). DTIC Document AL/OE-TR-1996-0037.


Fig. 1. Computational setup for the ultimate basis. The non-uniform body model (“Duke”) is shown in red. Every blue dot indicates a dipole location (there are six dipoles at every blue location: Three X, Y and Z electric dipoles and three X, Y and Z magnetic dipoles).

Fig. 2. a: SAR amplification factor for the shallow tumor as a function of the number of basis vectors and treatment frequency. b: Same thing for the deep tumor. The red dash lines indicate points in the convergence curves corresponding to 80% of the ultimate amplification factor.

Fig. 3. a: Linear and log maps of the SAR of the ultimate hyperthermia treatment of the shallow tumor at 600 MHz in transverse and sagittal views. Also shown is a map of the anatomy (tumor mask overlaid on the electrical conductivity map of Duke) showing the location of the tumor as well as its necrotic core (“C”) and well-perfused penumbra (“P”). b: Same thing for the deep tumor.

Fig. 4. a: Pennes bio-heat temperature simulation of the ultimate hyperthermia treatment of the shallow tumor at 600 MHz. The left plot shows temperature variations during the treatment time (=10 minutes) in the core and penumbra of the tumor as well as in the entire head (average). The rightmost maps are transverse and sagittal views of the temperature distribution at the end of the treatment through the tumor. b: Same thing for the deep tumor.

Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)