IEEE TRANSACTIONS ON MAGNETICS, VOL. 51, NO. 3, MARCH 2015
7001204
Temperature Field Optimization and Magnetic Nanoparticles Optimal Approximation of MFH for Cancer Therapy Guanzhong Hu, Yuling Li, Shiyou Yang, Yanan Bai, and Jin Huang College of Electrical Engineering, Zhejiang University, Hangzhou 310027, China The magnetic nanoparticles (MNPs) in ac alternating magnetic fields will produce a sufficient amount of heats owing to the Néel and Brown relaxations. Magnetic fluid hyperthermia (MFH), based on this mechanism, is a new promising approach for tumor treatments. The temperature field distribution in the cancer and its neighbor regions have a significant effect on the therapeutic effect of MFH. As a result, it is generally required to maintain the temperature in the cancer in the range of 42 °C–90 °C while guaranteeing a sharp temperature gradient in the neighbor regions of the cancer and healthy tissues. This paper provides an automated shape and MNP volume fraction solid design optimization methodology with finite element analysis, coupled field analysis, and an enhanced multiobjective quantum particle swarm optimization, to realize the aforementioned two ultimate goals. Moreover, an optimal radial basis function approximation technique is proposed to approximate the distribution of the optimized nanoparticles volume fraction solids to give a smooth and implementable nanoparticles distribution. Index Terms— Magnetic fluid hyperthermia, multiobjective optimization, quantum behaved particle swarm optimization (QPSO), radial basis function (RBF).
I. I NTRODUCTION AGNETIC nanoparticles (MNPs) have been used in biomedical applications and in vitro diagnostic protocols during the last 50 years. Hyperthermia is a promising therapy which heats cancer tissues to a certain temperature to moderate cellular inactivation or directly induce the localized cancer cell to death in a dose-dependent manner [1], [2]. In the therapy, MNPs are the heat generator and placed in the lesions of the human body by using a targeting technique. Temperature inside the treatment region is generally limited between 42 °C and 90 °C, because the temperature beyond this range has either less treatment effects (90 °C) [3], [4]. A completely new iron-free magnetic field exciter system is proposed and analyzed based on the uniform distribution of NPs for achieving a more uniform magnetic field in [5]–[8]. In [2] and [9], the NPs multiinjections procedure is synthesized to increase the temperature up to a certain therapeutic values by assuming a normal dispersion of NPs. However, the number of injections is not easily determined. In this paper, an automated magnetic-thermal field analysis is conducted with finite element analysis (FEA) and coupled field analysis (CFA), an automated design by using an enhanced multiobjective QPSO optimization algorithm, and an optimal tensor product radial basis function (RBF) approximation of the volume fraction, considering the shape of magnetic generating source and the distribution of NPs in cancer regions, is proposed for directly homogenizing the temperature fields.
M
II. FEA M ODEL AND CFA A. Magnetic Model of the Magnetic Fluid Hyperthermia Device The prototype device is originally designed to produce uniform magnetic fields [5]. The system consisted of two Manuscript received May 23, 2014; revised August 25, 2014; accepted October 7, 2014. Date of current version April 22, 2015. Corresponding author: S. Yang (e-mail:
[email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TMAG.2014.2363108
Fig. 1.
3-D model of the MFH device.
Fig. 2.
2-D simplified model of the MFH device.
induction coils, the main and correction coils, and the directions of currents of the two coils are opposite. The magnetic field excited by the main coils is used to sustain the suitable thermal field distribution while the magnetic field produced by the correction coils is to ensure a higher temperature gradient along the border between the healthy and cancer lesions. The 3-D configuration of the device is shown in Fig. 1. Because of the axisymmetry, the system is simplified to a 2-D model, as shown in Fig. 2. B. Thermal Field Analysis In an alternating magnetic field, the MNPs will produce a large amount of heats owing to the effect of Néel and Brown relaxations [10]. Using the first law of thermodynamics and with the human body as an incompressible fluid system,
0018-9464 © 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
7001204
IEEE TRANSACTIONS ON MAGNETICS, VOL. 51, NO. 3, MARCH 2015
Fig. 3. MNPs distribution. i is the ith volume fraction solids which varies from 0.05 to 0.1.
the power dissipation of MNPs is P = πμ0 χ0 H02 f × 2π f τ /(1 + (2π f τ )2 )
(1)
where H0 is the amplitude of the magnetic field intensity, χ0 is the chord susceptibility corresponding to the Langevin equation, f (450 kHz) is the frequency of the alternating magnetic field, μ0 is the permeability of the vacuum, and τ is the relaxation time. If the nanoparticles distribution be axisymmetric about the Z -axis, then the thermal distribution should also be axisymmetric. Therefore, the thermal field can also be simplified to be 2-D. The function of the heat transfer is cγ ∂ T /∂t = κ∇ 2 T + P(x, y, z, H, φ)
(2)
where P(x, y, z, H , ) is the power dissipation as given in (1), κ is the thermal conductivity of the tissue, c is the specific heat, γ is the density of the tissue, is the volume fraction solids, and T (K ) is the thermodynamic temperature. III. M ULTIOBJECTIVE A LGORITHM AND RBF A PPROXIMATION A. Optimal Model of MFH Device As described previously, the temperature in the treatment region should be as uniform as possible, and the differences between the temperatures of the treatment and normal tissues should be as large as possible. As a result, a two objective function problem is optimized for a magnetic fluid hyperthermia (MFH) device, that is, f1 to measure the uniformity of the temperature field while f 2 to define the temperature gradient between the boundary B and the average temperature of the tumor region. Moreover f 1 = (Tmax − Tmin )/Tmean (3) f 2 = 1/Tadj − Tmean where Tmax , Tmin , and Tmean are, respectively, the maximum, minimum, and average temperatures in the tumor region and Tadj is the average temperature of the nodes on boundary E adjacent to the boundary B, as shown in Fig. 2. For medical purpose, the temperature is limited to ⎧N ⎪ ⎪ φi − N · φav = 0 ⎨ i=1 (4) ⎪ ≥ 42 °C T ⎪ ⎩ mean Tmax ≤ 90 °C. The decision variables are including the independent geometric parameters, as shown in Fig. 2, and the distributions of the MNPs in each element of a 4 × 4 subdomain, as shown in Fig. 3; totaling to 26 ones in a compact form of x = (x 1 , x 2 , . . . , x 26 ).
Fig. 4.
Flow chart of the proposed MOQPSO-DE.
To be fairly compared, the total amount of MNPs in the treatment region is kept identical, and is limited in [0.05, 0.1]. B. Multiobjective QPSO Algorithm 1) Algorithm: QPSO is a population-based optimal algorithm which is inspired by the classical Particle Swarm Optimization (PSO) and the quantum delta potential well model. The state of a particle in QPSO is controlled by a wave function which describes the particle’s appearing in position x in a certain probability. QPSO is suitable for complex optimal problems because the algorithm has a strong global searching ability and high convergence speed [12], [13]. The position of a particle in QPSO is updated from x t +1 = pt ± β · |Mbest − x t | · ln(1/u t )
(5)
where u t is the random number uniformly distributed in [0, 1], β is the contraction–expansion coefficient which can be tuned to control the convergence speed, and Pt is the local attractor and defined as pt = ϕ · Pbi + (1 − ϕ) · Gbi
(6)
where Pbi and Gbi are the optimal location according to the personal and the global optimal values of the last i th iteration, respectively. The parameter Mbest can be the personal optimal point Pbi , the mainstream thought point Pms , or the median point Pmp according to different problems [12]. To combine with nondominated crowding distance sorting and Differential Evolutionary (DE) strategy [13], [14], a MultiObjective (MO) algorithm, called MOQPSO-DE, as shown in Fig. 4, has been implemented to improve the global convergent speed and increase the diversity of the population. The algorithm introduces two external files: one for storing the solutions with Pareto dominated Rank 1 (Archive) and the other for the local minima (Pbest ). Pbi is selected randomly from Pbest and Gbi in the Archive. To guarantee the stability of the convergences, parameter Pmp is usually selected as Mbest . 2) Validation: To analyze the strength and weakness of the proposed MOQPSO-DE, it is numerically experienced on two test problems: Schaffer’s function (SCH) and 15D-ZDT3 [13], and compared with the fast elitist multiobjective genetic algorithm: NSGA II [13]. To give a fair comparison, every algorithm is run for 30 times with a
HU et al.: TEMPERATURE FIELD OPTIMIZATION AND MNPs OPTIMAL APPROXIMATION
7001204
TABLE I C( A) AND S M ETRIC OF P ROPOSED MOQPSO-DE AND NSGA II
Fig. 6. Searched Pareto Front of MOQPSO-DE algorithm. The objective functions f 1 and f 2 are all normalized to [0, 1] by using antitangent function. TABLE II P ERFORMANCE PARAMETERS OF THE S EARCHED PARETO S OLUTIONS Fig. 5. True Pareto and searched Pareto solutions searched by MOQPSO-DE and NSGA II algorithm for the different standard problems.
maximum of 10 000 function evaluations by using the same parameter settings. To evaluate the quality of the final solution of an optimizer, the average convergence metric C(A) [13] and spacing metric S [15], as proposed by Deb and Schott, to compare the diversity and uniformity of Pareto solutions are used, and compared in Table I for the two algorithms. Fig. 5 shows the final solutions of the two algorithms in a typical run. Moreover, the averaged CPU time for 20 independent runs of the proposed algorithm is only about 15% of that of NSGA II for solving ZDT3. Obviously, the proposed algorithm is definitely superior to NSGA II for solving SCH, and not inferior to NSGA II for solving ZDT3 but using a very small number of CPU times. C. RBF Approximation The final distribution of MNPs in the neighbor elements obtained using an optimizer is generally noncontinuous, and infeasible in an application view. In this regard, some technique should be designed to address this problem. Therefore, the RBF is proposed to act as an approximation technique to smoothly interpolate the optimized result to produce a feasible MNPs distribution. Radial function is defined as H (x, x j ) = H (r ) = H (x − x j )
(8)
and the interpolation of a function is f i (x) =
N
c j H j (x i − x j )
X = C F.
(10)
Moreover, the tensor product MQ function is adopted
0.5
0.5 · (y − yi )2 + h . (11) H (r ) = (x − x i )2 + h
(7)
where x j ( j =1, 2, …, N) is the coordinate of the i th sampling point. Because of the accuracy, stability, and efficiency in the scattered data interpolation, the most commonly used RBF is multiquadrics (MQ) function and used. Moreover H (r ) = (r 2 + h)0.5
Let i varies from 1 to N, C = [c1 c2 . . . cN ]T , F = [ f 1 f 2 . . . f N ]T , and X is N × N matrix, one has
(9)
j =1
where the shape parameter h affects the shape of the interpolation and the value must be close to or less than the average distance among the interpolated points.
IV. N UMERICAL R ESULTS A prototype device is optimized using the proposed model and methodology. The searched Pareto front using the proposed MOQPSO is shown in Fig. 6. Table II lists some typical ones. As shown in Table II, n tot is the total number of nodes in the tumor region whose temperature is higher than a predefined value, such as 41.5 °C or 42 °C, and Tx (x is minimum, maximum, or average) is the minimum (maximum or average) temperature of the tumor region. For performance comparisons, the last two columns give the results consistent with [7] and [8]. Compared with [7] and [8], it is obvious that an improved performance of MFH device, especially the solution of No. 5, is observed for the proposed algorithm. More specially,
7001204
IEEE TRANSACTIONS ON MAGNETICS, VOL. 51, NO. 3, MARCH 2015
RBF-approximated smoothly distributed MNPs, and Fig. 7(c) shows the corresponding temperature fields. Obviously, the approximated MNPs distribution is a reliable and practical one. Moreover, the focus of this paper is on the optimization of the NPs distribution. However, the optimal injection points for achieving such distributions are being investigated by us, and will be presented in another paper. V. C ONCLUSION This paper proposes an automated optimization methodology, by using the numerical magnetic-thermal CFA, MOQPSO-DE, and RBF interpolation method, for the design optimization of an MFH device. The salient feature of MOQPSO-DE is that it has strong global search ability even on disconnect nonconvex problems, and obtain more reliable and practical results of an MFH device. ACKNOWLEDGMENT This work was supported in part by the 973 Program of China under Grant 2013CB035604, in part by the National Natural Science Foundation of China (NSFC) under Grant 51377139, and in part by the Specialized Research Fund for the Doctoral Program of Higher Education of China under Grant 20120101110110. R EFERENCES
Fig. 7. Distribution of temperature and MNPs of No. 5 solution. (a) Optimized MNPs distribution. (b) RBF interpolated MNPs distribution. (c) Temperature distribution.
the n tot /N(%) of No. 5 solution has obtained higher temperature percentages. It reveals that an effective hyperthermia pattern can be refined by the tuning of the distribution of MNPs . Fig. 7 shows the distribution of temperature and for No. 5 solution. Fig. 7(a) shows the MNPs distribution of a directly optimized results of No. 5 solution, and Fig. 7(b) shows the
[1] P. Di Barba, F. Dughiero, and E. Sieni, “Synthesizing distributions of magnetic nanoparticles for clinical hyperthermia,” IEEE Trans. Magn., vol. 48, no. 2, pp. 263–266, Feb. 2008. [2] A. Jordan, R. Scholz, P. Wust, H. Fähling, and R. Felix, “Magnetic fluid hyperthermia (MFH): Cancer treatment with AC magnetic field induced excitation of biocompatible superparamagnetic nanoparticles,” J. Magn. Magn. Mater., vol. 201, nos. 1–3, pp. 413–419, 1999. [3] P. Wust et al., “Hyperthermia in combined treatment of cancer,” Lancet Oncol., vol. 3, no. 8, pp. 487–497, 2002. [4] A. Candeo and F. Dughiero, “Numerical FEM models for the planning of magnetic induction hyperthermia treatments with nanoparticles,” IEEE Trans. Magn., vol. 45, no. 3, pp. 1658–1661, Mar. 2009. [5] P. Di Barba, F. Dughiero, and E. Sieni, “Magnetic field synthesis in the design of inductors for magnetic fluid hyperthermia,” IEEE Trans. Magn., vol. 46, no. 8, pp. 2931–2934, Aug. 2010. [6] P. Di Barba, F. Dughiero, and F. Trevisan, “Optimization of the Loney’s solenoid through quasi-analytical strategies: A benchmark problem reconsidered,” IEEE Trans. Magn., vol. 33, no. 2, pp. 1864–1867, Mar. 1997. [7] P. Di Barba, “Evolutionary multiobjective optimization methods for the shape design of industrial electromagnetic devices,” IEEE Trans. Magn., vol. 45, no. 3, pp. 1436–1441, Mar. 2009. [8] G. Hu, Y. Li, S. Yang, J. Chu, Y. Bai, and J. Huang, “Temperature field distribution optimization of MFH for cancer therapy,” Int. J. Appl. Electromagn. Mech., vol. 45, nos. 1–4, pp. 697–702, 2014. [9] H. Guanzhong, Y. Shiyou, and L. Yuling, “Temperature field optimization of nano-magnetic fluid hyperthermia device,” Trans. Chin. Electrotech. Soc., vol. 28, no. 2, pp. 194–199, 2013. [10] R. E. Rosensweig, “Heating magnetic fluid with alternating magnetic field,” J. Magn. Magn. Mater., vol. 252, pp. 370–374, Nov. 2002. [11] W. Fang, J. Sun, Y. Ding, X. Wu, and W. Xu, “A review of quantumbehaved particle swarm optimization,” IETE Tech. Rev., vol. 27, no. 4, pp. 336–348, Jul. 2010. [12] S. N. Omkar, R. Khandelwal, T. V. S. Ananth, G. N. Naik, and S. Gopalakrishnan, “Quantum behaved particle swarm optimization (QPSO) for multi-objective design optimization of composite structures,” Expert Syst. Appl., vol. 36, no. 8, pp. 11312–11322, 2009. [13] K. Deb, “A fast and elitist multiobjective genetic algorithm: NSGA-II,” IEEE Trans. Evol. Comput., vol. 6, no. 2, pp. 182–197, Apr. 2002. [14] R. Storn and K. Price, “Differential evolution—A simple and efficient heuristic for global optimization over continuous spaces,” J. Global Optim., vol. 11, no. 4, pp. 341–359, 1997. [15] J. R. Schott, “Fault tolerant design using single and multicriteria genetic algorithm optimization,” M.S thesis, Massachusetts Inst. Technol., Cambridge, MA, USA, 1995.
