Magnetic fluid hyperthermia simulations in evaluation of SAR calculation methods

Physica Medica 71 (2020) 39–52

Contents lists available at ScienceDirect

Physica Medica journal homepage: www.elsevier.com/locate/ejmp

Original paper

Magnetic fluid hyperthermia simulations in evaluation of SAR calculation methods

T

Costas Papadopoulosa, Eleni K. Efthimiadoub, Michael Pissasc, David Fuentesd, Nikolaos Boukosc, ⁎ Vassilis Psycharisc, George Kordasc, Vassilios C. Loukopoulose, George C. Kagadisa,d, a

Department of Medical Physics, School of Medicine, University of Patras, Rion, GR 26504, Greece Inorganic Chemistry Laboratory, Department of Chemistry, National and Kapodistrian University of Athens, Panepistimiopolis, Zografou, GR 15784, Greece c Institute of Nanoscience and Nanotechnology, NCSR Demokritos, Athens, GR 15310, Greece d Department of Imaging Physics, The University of Texas MD Anderson Cancer Center, Houston, TX 77030, USA e Department of Physics, University of Patras, Rion, GR 26504, Greece b

ARTICLE INFO

ABSTRACT

Keywords: Magnetic fluid hyperthermia Magnetic nanoparticles Specific absorption rate SAR ILP Simulations COMSOL Box-Lucas Initial slope method Corrected slope method

Purpose: The purpose of this study is to employ magnetic fluid hyperthermia simulations in the precise computation of Specific Absorption Rate functions -SAR(T)-, and in the evaluation of the predictive capacity of different SAR calculation methods. Methods: Magnetic fluid hyperthermia experiments were carried out using magnetite-based nanofluids. The respective SAR values were estimated through four different calculation methods including the initial slope method, the Box-Lucas method, the corrected slope method and the incremental analysis method (INCAM). A novel numerical model combining the heat transfer equations and the Navier-Stokes equations was developed to reproduce the experimental heating process. To address variations in heating efficiency with temperature, the expression of the power dissipation as a Gaussian function of temperature was introduced and the LevenbergMarquardt optimization algorithm was employed to compute the function parameters and determine the function’s effective branch within each measurement’s temperature range. The power dissipation function was then reduced to the respective SAR function. Results: The INCAM exhibited the lowest relative errors ranging between 0.62 and 15.03% with respect to the simulations. SAR(T) functions exhibited significant variations, up to 45%, within the MFH-relevant temperature range. Conclusions: The examined calculation methods are not suitable to accurately quantify the heating efficiency of a magnetic fluid. Numerical models can be exploited to effectively compute SAR(T) and contribute to the development of robust hyperthermia treatment planning applications.

1. Introduction Hyperthermia is a well-established therapeutic modality that utilizes heat to impair cancer cells with negligible damage to surrounding healthy tissues and minimal side effects. Depending on the temperature range of the treatment, hyperthermia can be distinguished as moderate hyperthermia (41–46 °C) or thermo-ablation (T > 46 °C). This classification is based on different biological mechanisms of cellular death present at each temperature range [1,2]. Several local hyperthermia methods [3,4], either external or interstitial, have been employed in clinical practice, including capacitance, ultrasound, microwave, radiofrequency, and laser-mediated hyperthermia. However, these methods may be invasive or limited in ⁎

terms of tumor targeting, energy localization, and temperature distribution homogeneity [1,5,6]. To transcend such limitations, novel hyperthermia techniques exploiting image-guidance and specialized energy deposition [7–9] have been introduced. Among these techniques magnetic fluid hyperthermia (MFH) has attracted considerable scientific interest. MFH is a minimally invasive therapy that uses a magnetic fluid as the heating source [10]. A magnetic fluid is a stable colloidal suspension [11] of magnetic nanoparticles (MNPs) dispersed in an aqueous medium. An analysis of intrinsic MNPs’ properties is beyond the aims of the current study, and is available in the literature [12–20]. The heating efficiency of magnetic colloids is typically quantified through the specific absorption rate (SAR) [2], also referred as the

Corresponding author at: Department of Medical Physics, School of Medicine, University of Patras, Rion, GR 26504, Greece. E-mail address: [email protected] (G.C. Kagadis).

https://doi.org/10.1016/j.ejmp.2020.02.011 Received 28 November 2019; Received in revised form 21 January 2020; Accepted 13 February 2020 1120-1797/ © 2020 Associazione Italiana di Fisica Medica. Published by Elsevier Ltd. All rights reserved.

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al.

specific loss power [18,21,22], which is defined as the rate of energy dissipation to the medium per unit mass of the solid fraction of the sample. As the solid fraction either the MNPs’ or the iron’s mass is considered. Thus, research groups should provide a relevant notation with their results. Another quantification index is the intrinsic loss power (ILP) [21,23], also reported as effective specific absorption rate [1], that normalises SAR over the external parameters that affect heating efficiency, i.e. the product of the AMF frequency times the square of the amplitude. However, this normalization emanates from the linear response theory [24], which has been shown to be valid only at low fields and non-interacting systems [25–28], thus it should be carefully used. Finally, the volumetric power dissipation (P) [24], is another performance metric, that treats the whole magnetic fluid’s volume as a homogenous heating source. The accurate determination of heating efficiency is imperative in optimizing MFH treatment planning. Different approaches have been reported in the literature, including calorimetric and theoretical methods. Theoretical models are interesting from the perspective of studying the magnetic and physical properties of nano-systems, and have been employed in previous studies to quantify heating efficiency, however their predictive power is still ambiguous [21]. A recent study [29] proposed the giant-spin non-linear model to address the issues of the linear response theory with respect to AMF amplitudes. However, the collectivity effects with increasing sample concentrations that have been shown to drastically affect magnetic properties [26,30–32] are still ignored. Another study [33] attempted to resolve the issue of local dipolar interactions through an adapted Coffey model. Yet the model is limited within the linear response region. More robust theoretical approaches include the Stochastic Landau-Lifshitz-Gilbert [34,35] and Kinetic Monte Carlo [25,36] models. However, these models remain bound to the validity of the MNPs’ properties used in the calculations. In addition, even magnetic fluids fabricated through standardized fabrication methods exhibit varying properties between batches [23,37]; thus, every batch should be characterized to recalibrate theoretical models. The recalibration should then be validated against experimental data. From this perspective a theoretical approach while still of interest, is not currently applicable in clinical practice. Studying a magnetic fluid’s heating efficiency via calorimetric measurements constitutes a more straightforward approach. Several methods have been developed by research groups to estimate SAR values directly from the heating curve. The pulse heating method under adiabatic conditions is reported as the most accurate calorimetric method and allows the determination of the temperature-dependent heating efficiency SAR(T) [38–40]. However, achieving adiabatic conditions is time-consuming and challenging; thus, calorimetric measurements are typically carried out in non-adiabatic conditions. The prominent calculation methods based on calorimetric data include the initial slope method (ISM) [41–43] and the Box-Lucas method (BLM) [23]. More methods have been introduced to optimize the ISM estimates, including the corrected slope method (CSM) [42] and the incremental analysis method (INCAM) [21,44]. The decay method is another approach that uses the steady state temperature and a cooling parameter to determine SAR in non-adiabatic conditions [42,45]. While these methods are relieved from the complexity of theoretical modelling, they are also bound to limitations that must be considered during their application. Particularly, the ISM and INCAM assume adiabatic conditions, while the BLM, CSM, and decay method assume linear losses. However, experimental conditions are typically non-adiabatic, and the linear loss region is limited to a specific temperature range that depends on these conditions and must be determined [42]. Considering that in real systems multiple loss mechanisms are evident – including conduction, convection and radiation - compressing their effect to a linear loss parameter is an oversimplification [38,42]. Furthermore, the experimental design across the literature is not standardized and several factors that may affect the SAR estimation

vary. Particularly, superficial insulation may or may not be applied around sample containers [22,42]. Different volumes of magnetic fluids are employed, while the containers vary in shapes and materials [22,42,43]. The water-cooled solenoids that are typically used create a microenvironment that depends on their geometrical characteristics and may result in a cooling or heating effect, depending on the current flow [21]. Inhomogeneities of the magnetic field across the magnetic fluid’s volume [22,43,46] constitute an additional source of error to be considered during experimental design. The optical fiber tip’s positioning also has been reported to affect the calculations as the magnetic fluid is not an isothermal domain [21,22,43]. Finally, all calculation methods ignore any variation of physical and chemical properties of the sample during the heating process [21] and are based on the fundamental acceptance that SAR is constant with temperature. However, several studies [39,40,47–50] demonstrate the dependence of SAR on temperature. To address these issues well-elaborated numerical models can be developed. In this work, we exploit detailed magnetic fluid hyperthermia simulations to evaluate the accuracy of calculation methods in SAR estimations, including ISM, BLM, INCAM, and CSM. We accomplish this venture by introducing a simple yet robust method to calculate in silico the exact power dissipation of the magnetic fluid as a function of temperature. 2. Materials and methods 2.1. Characterization of magnetite (Fe3O4) SPIONs In the present study, uncoated iron oxide nanoparticles were fabricated through the co-precipitation method [51]. The iron oxide core size of the MNPs was studied with transmission electron microscopy, using a 200 kV acceleration voltage TEM FEI CM20 Microscope. The observations indicated intense agglomerations and a wide size distribution of log-normal type with a median iron oxide core diameter DM = 10.7 nm and a dispersion parameter w = 0.35 (σ = 0.38 nm) (Fig. 1). The results were verified with x-ray diffraction [52], using a Siemens D-500, Siemens AG diffractometer. A median size of 10.6 nm was derived using the Scherer’s law [53] upon a Voight pseudo-function fitting on the major reflection of the pattern. A median hydrodynamic size DH = 372 nm was measured for a 0.1 mg∙mL−1 sample through dynamic light scattering using a multimodal Malvern Instruments Series, Nano-ZS device. Magnetization measurements under high field were carried out for a powder sample to determine the saturation magnetization MS = 58.1 emu∙g−1, using a 5.5 T Quantum Design PPMS DC SQUID Magnetometer. The absence of a hysteresis loop at room temperature, verified the superparamagnetic behavior of the sample [15] (Fig. 1). 2.2. Calorimetric measurements An aqua-based magnetic fluid was prepared and exposed to three different AMF properties to obtain experimental hyperthermia heating curves. In a borosilicate glass vial (V = 2 mL, ID = 13 mm, OD = 15 mm, height = 31 mm), 30 mg of uncoated magnetite MNPs was dispersed in 1 mL of deionized water. Before data acquisition the sample was exposed to ultrasonication for 20 min to enhance dispersion and colloidal stability. The AMF used in the experiment was generated by an eight-turn, water-cooled (Twater = 17 °C), copper solenoid (ID = 2.6 cm, L = 3.5 cm, pitch = 0.5 cm) connected to an Ambrell EASYHEAT system. The vial of magnetic fluid was placed on the sample holder within the coil, and the magnetic fluid’s volume was centered with respect to the coil’s z-axis, as well as the coil’s diameter. No superficial insulation was applied around the sample. A Yokogawa fiber-optic system was used to monitor the temperature 40

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al.

Fig. 1. Characterization of magnetite SPIONs. (a) Magnetization curve, (b) x-ray diffraction, (c) MNPs’ size distribution derived by TEM measurements, (d) TEM image.

at the central point of the magnetic fluid’s volume where the tip of a 1 mm optical temperature probe was immersed. The system was allowed an appropriate time for temperature stabilization before activating the AMF. The duration of the heating process was 1200 s with a sampling rate of 0.5 s. The interval of 1200 s was used to allow the calculation of P(T), SAR(T) within a wide MFH-relevant window and demonstrate how temperature affects SAR. The sample was heated under three different AMF properties, with field amplitudes of 10.1, 17.2, and 30.3 kA∙m−1 and frequencies of 291, 283, and 275 kHz, respectively. Upon deactivation of the AMF, a cooling curve was obtained for the measurement of the maximum temperature increase to determine the linear loss region. To evaluate the response with respect to the concentration, two more samples with concentrations of 1 and 10 mg∙mL−1 were prepared and heated under the highest AMF properties. Triple measurements were taken for every experimental setting. Control experiments were carried out using a deionized water sample of 1 mL to evaluate possible heating or cooling effects from the coil. The experiments were indicative of minimal heating effects that can be attributed to the development of a microenvironment within the coil. To implement this effect in the simulations, we acquired measurements of the air temperature in the center of the coil for each of the AMF properties (Fig. 2).

that at the start of the heating process thermal losses are negligible, SAR is calculated as [20]:

SARISM = cmf

mmf mMNPs

T t

, t

0

(1)

where cmf is the specific heat of the magnetic fluid, mmf the magnetic

2.3. SAR calculation methods 2.3.1. Initial slope method The ISM has known widespread use in SAR estimations. Assuming

Fig. 2. Coil microenvironment for different AMF properties. The air temperature at the center of the coil, was recorded and implemented in the simulations. 41

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al.

fluid’s mass, and mMNPs the mass of the dispersed MNPs. The magnetic fluid’s power dissipation, P, is accordingly derived as:

PISM =

mf cmf

T t

, t

2.3.3. Corrected slope method As the name of this method implies, the CSM [42] corrects the estimates of the ISM by considering the early linear losses that occur from the start of the heating process. Multiple linear fittings are applied over an interval within the linear loss region. The authors have developed and provide an open source script for MATLAB that can be utilized to apply the method [55]. The SAR and P values are derived from the following equations:

(2)

0

where ρmf is the magnetic fluid’s density. The P and SAR values were calculated through differentiation of a linear fitting on the first 20 s of each curve, upon removal of any intervals of delayed [42] heating.

SARCSM =

2.3.2. Box-Lucas method Assuming linear losses, a least-squares fit of the Box-Lucas phenomenological model to the experimental data describes the temperature elevation [54]:

T = a (1

e

PBLM =

mf

(4)

a cmf

(5)

The linear loss parameter L can also by calculated: (6)

L = mmf cmf

mf m

dT dt

i

+ L T)/ i mMNPs

(7)

MNPs

mmf

(8)

2.3.4. Incremental analysis method The incremental analysis method was introduced [44] to determine the optimal time interval for the application of the ISM through analysis of a plot of the temperature difference of successive increments Tn − Tn−1 against time. In a later work [21], it was suggested that the method should be restricted to the increments that fall within the socalled quasi-adiabatic phase of the process, which is limited to the 5% of the net temperature increase along the linear region of the heating curve. The INCAM is useful for conduction of intra-curve statistical analysis of all possible SAR values. The numerical temperature derivatives for different time increments that fulfill quasi-adiabatic conditions may be reflected in box and whisker plots to provide insights on the variance of SAR. We determined the quasi-adiabatic phase via graphical analysis for each of our measurements. Instead of using 5% of the net temperature increase, we further restricted the compliant interval to the region in which the numerical temperature derivative remains stable with respect to a least-squares type linear fitting (Fig. 4). Different intervals were evaluated through linear fittings to localize the derivative plateaus. We observed that the initial phases of heating met the quasiadiabatic criterion upon subtraction of any intervals of delayed heating. Eqs. (1) and (2) were used to calculate the SAR and P values for each increment. Instead of using linear fittings to derive the temperature elevation rates as in the ISM, direct divisions of ΔT/Δt were employed in the INCAM calculations.

a cmf mmf mMNPs

i

(mmf c mf

We carried out the calculations for each heating curve, using ten consecutive 5 s intervals along a 50 s timeframe, which exhibited the maximum efficiency.

where a, λ are the fitting parameters. From the time derivative of Eq. (3), SAR and P are calculated respectively as:

SARBLM =

N

PCSM = SARCSM

(3)

t ),

1 N

Proper application of this method preconditions that the fitting is maintained within the linear loss region. However, the selected time interval should provide enough saturation of the model along the heating curve, otherwise unacceptable standard errors arise in the calculation of the fitting parameters. Thus, having knowledge of the linear loss region is required. To access the linear loss region, a cooling experiment was conducted using the 30 mg∙mL−1 sample and the highest AMF properties. Upon deactivation of the AMF a cooling curve was obtained and upon comdT putation of the curve’s numerical derivative, the product Cmf dt was

( )

plotted versus T = Tmax T (t ) , where Cmf is the magnetic fluid’s heat capacity, Tmax the maximum temperature of the measurement and T(t) is the temperature decay throughout the cooling phase (Fig. 3). A linear loss region of ~20 °C was observed. The Box-Lucas fittings were then applied across an 80 s time interval for all measurements, except for the 1 mg∙mL−1 sample, in which an interval of 240 s was required to reduce the fitting parameter’s standard errors within reasonably low levels.

2.4. Statistical and residual analysis Statistical analysis was carried out to calculate the standard deviation s, of the mean value x , between Ν = 3 parallel experiments for each measurement:

s=

1 N

N

(x i i=1

x )2

(9)

Any systematic error was neglected as marginal. Residual analysis was also carried out according to [54,56] including: a) scatter plot of the response and predicted values versus the independent variable; b) scatter plot of the residuals versus the independent variable; c) scatter plot of the residuals versus the predicted values; d) normal probability plot of the residuals, and e) histogram of the residuals. The formation of trends was visually observed. 2.5. Modelling 2.5.1. Magnetic field simulations To compute the AMF amplitude (Fig. 5) and homogeneity within the magnetic fluid’s volume, a 2D-axisymmetric model was created in

Fig. 3. Determination of linear loss region via cooling experiment upon heating of a 30 mg∙mL−1 sample under a 30.3 kA∙m−1 and 275 kHz AMF. 42

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al.

(10.d)

×H=J

where J is the current density, E the electric field, D = ε0εrΕ the electric displacement field, ε0 the permittivity of free space, εr the relative permittivity of the material, ω the angular frequency, H the magnetic field, A the magnetic vector potential, Je the externally generated current density, and B the magnetic flux. The model was optimized according to the guidelines concluded by Kennedy et al. [58]. These guidelines include the selection of a sufficiently large air domain surrounding the coil; thus, an infinite elements domain was applied around the air domain to allow proper formation of the magnetic field. A perfect magnetic conductor n × H = 0 and a magnetic insulation n × A = 0 boundary conditions were set across the inner and outer boundaries, respectively, of the infinite elements’ domain. In addition, the use of current-driven, single-conductor coils to provide the correct magneto-motive force, as well as the use of boundary meshes on the spirals’ surfaces to account for the limited skin depth of the tube, are suggested in coil modelling. A direct, stationary, MUMPS [59] solver using linear equations with automatic rescaling was selected. The relative tolerance was set to 0.001. The complete mesh consisted of 75,430 domain elements plus 2130 boundary elements. The number of degrees of freedom rose to 165,938 and the solution time to 7–9 s. The inhomogeneity, UM, of the magnetic field within the magnetic fluid’s volume was evaluated using the expression [60]:

UM =

Hmax Hmin Hmean

(11)

2.5.2. Hyperthermia simulations A 2D-axisymmetric model was created in COMSOL Multiphysics 5.2a to emulate the heating process. Upon detailed measurements, the vial, the magnetic fluid’s volume and the sample holder geometries were imported into the model. The material properties were retrieved from the literature [24,61–64]. The contribution of the MNPs to the thermal properties of the solvent was considered through the following equations [65,66]:

kmf = k w (T )(1 + 10.5 )0.1051

(12)

cmf = (1

) c w (T ) + cMNPs

(13)

)

(14)

mf

COMSOL Multiphysics 5.2a [57], which uses finite element methods. The coil geometry was imported in the model upon detailed measurements. Frequency domain studies were carried out for all the AMF properties, using the Magnetic Fields interface. The model is governed by Ampere’s law assuming time-harmonic fields:

B=

E=

×A

j A

× B + Je

w (T )

+

MNPs,

where kmf , cmf , and mf are the magnetic fluid’s thermal conductivity, specific heat, and density, respectively. kw(T), cw(T) and ρw(Τ) are the respective temperature-dependent thermal properties of water; cMNPs = 670 J∙kg−1∙K−1, ρMNPs = 5180 kg∙m−3, respectively, are the specific heat and density of magnetite MNPs; and φ is the volume fraction of the MNPs. Eq. (12) is applicable for magnetite nanofluids with 0 < < 0.02 . The “Non-Isothermal Flow” COMSOL Multiphysics interface, a double coupling of the “Heat Transfer in Fluids” and the “Laminar Flow” modules, were used to simulate the heating and stirring of the magnetic fluid (Fig. 6). The governing equations for the thermal problem are the following:

Fig. 4. (a) Determination of quasi-adiabatic phase via graphical analysis of temperature derivatives. (b) Quasi-adiabatic phase superimposed on the net temperature increase, and (c) statistical analysis of possible SAR values.

J= E+j D+

= (1

mf cmf

q=

T + t kmf T,

mf cmf

u

T+

q = P + Qp + Qvd

(15.a) (15.b)

where u is the fluid velocity field, q the heat flux by conduction, P the volumetric power dissipation of the magnetic fluid, Qp the volumetric work rate of pressure changes, and Qvd the volumetric viscous dissipation. To solve the thermal problem within the solid domains of the model, a “Heat Transfer in Solids” node governed by the corresponding equation was added.

(10.a) (10.b) (10.c) 43

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al.

Fig. 5. (a) Eight-turn coil dimensions and magnetic field distribution; (b) magnetic field contours (H0 = 30.3 kA∙m−1).

Fig. 6. (a) Model geometry, boundary conditions and temperature distribution, and (b) velocity field for the 30 mg∙mL−1 sample under 30.3 kA∙m−1 and 275 kHz AMF.

The developed flow is described by the full Navier-Stokes equations: mf

u + t

mf

(u

)u =

pI + µ ( u + ( u )T )

2 µ( 3

gravity. An open boundary condition was employed to account for the interaction of air layers within and outside the vial:

u) I + F

pI + µ ( u + ( u )T )

(16.a)

T = Tamb,

t

+

(

mf

u) = 0,

(16.b)

n q = 0,

u) I n = 0

(17.a) (17.b)

n u<0 n u

2 µ( 3

0

(17.c)

Natural convection boundary conditions were set across the external boundaries of the geometry:

where μ is the dynamic viscosity, p the pressure of the fluid, exponential T indicates a tensor, and I is the identity matrix. F is the volumetric force with an algebraic value Fz = g mf along the z-axis. A hydrostatic pressure p = mf g (h z ) is assumed for the magnetic fluid, where h is the z-coordinate of the magnetic fluid’s surface and g the acceleration of

q0 = h (Tamb

T)

h = hair (L , pA , Tamb), 44

(18.a) (18.b)

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al.

Insulation” boundary condition was set at the bottom of the sample holder. We accounted for the main component of radiation losses by applying a diffusion surface boundary condition on the outer boundary of the vial. Assuming an ideal grey body, the radiation losses to each direction defined by the unity vector n may be quantified through the Stefan-Boltzmann equation:

n q=

Table 1 Element parameters for five different mesh resolutions. Min Element Size (m)

Max Element Size (m)

Number of Domain Elements

Number of Boundary Elements

Coarse Normal Fine Finer Extra fine Extremely fine

4e−5 3e−5 2e−5 1e−5 4e−6 1.5e−6

8.7e−4 6.7e−4 4.5e−4 3.5e−4 2.8e−4 1.3e−4

2163 3067 5177 12,947 32,459 53,344

253 311 379 691 1222 1399

T 4 ),

(19)

where ε = 0.87 [67] is the emissivity of borosilicate glass and σ the radiating surface. A homogenous volumetric power source was assigned to the domain of the magnetic fluid to emulate the magnetic fluid’s power dissipation P, upon activation of the AMF. The P value for each experiment was determined through an optimization study. A least-squares objective was defined at the central point of the magnetic fluid’s domain, where the optical fiber’s tip was positioned to compare the temperature calculated in the simulations with experimental measurements. The objective function was controlled by the power dissipation P of the magnetic fluid. A P estimate derived from the calculation methods based on the experimental data was set as the initial condition in the optimization study. The gradient based LevenbergMarquardt [68] algorithm, which is specifically designed to solve leastsquares-type problems, was used to minimize the objective throughout 1200 s time-dependent studies, with respect to the mean values of the experimental data, returning the optimal value of P. A fully coupled direct linear PARDISO [69] solver implementing a constant Newtonian non-linear method was used throughout the studies. The default options for Jacobian update, damping factor, tolerance factor, and maximum number of iterations, were selected. The optimality tolerance of the control variable was set to 0.001. The first series of the optimization studies, in which a constant P value was assumed, were indicative that a constant power dissipation is insufficient to describe the heating process (Fig. 7), except for measurements where insignificant net temperature changes were monitored. Thus, a function that adequately describes the power dissipation variation with temperature must be defined. Considering that the heating efficiency of a magnetic fluid may increase at the start of the heating process, reach a maximum value Pmax at a temperature TPmax, and then decrease with increasing thermal energy [47], we suggest the expression of power dissipation as a Gaussian function of temperature:

Fig. 7. Simulated heating curve assuming constant power dissipation (small dots), and variable power dissipation (dashed line), superimposed on the mean experimental data (large dots) of a 30 mg∙mL−1 sample, under 10.1 kA∙m−1 AMF. Error bars represent the standard deviation of the mean, derived from triple measurements.

Mesh Resolution

4 (T amb

where q0 is the convective energy loss, Tamb the ambient temperature, T the temperature at the boundary, and h the heat transfer coefficient. The software handles the dependency of the heat transfer coefficient hair as a function of the height L of the interacting vertical wall, the pressure pA of the air, and the ambient temperature Tamb. A “Thermal

Fig. 8. (a) Time-dependent studies for six different mesh element sizes. Temperature elevation at the center of the magnetic fluid's volume. (b) Temperature profiles along the radial direction from the central point of the magnetic fluid, towards the vial walls. 45

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al.

Table 2 P, SAR and ILP values and the respective standard deviations for five different experiments estimated through calculation methods and simulations. Experiment

ISM

BLM

CSM

H0 (kA∙m−1)

c (mg∙mL−1)

P (104∙W∙m−3)

SAR (W∙g−1)

ILP (nH∙m∙kg−1)

P (104∙W∙m−3)

SAR (W∙g−1)

ILP (nH∙m∙kg−1)

P (104∙W∙m−3)

30.3 30.3 30.3 17.2 10.1

1 10 30 30 30

3.25 ± 0.17 29.5 ± 1.17 77.89 ± 4.31 31.79 ± 1.64 12.32 ± 0.47

32.47 ± 1.72 29.26 ± 1.16 25.35 ± 1.40 10.35 ± 0.53 4.01 ± 0.15

0.129 0.116 0.100 0.124 0.143

2.86 ± 0.15 30.53 ± 1.12 80.43 ± 4.36 33.01 ± 1.67 12.72 ± 0.48

28.54 ± 1.48 30.29 ± 1.11 26.18 ± 1.42 10.75 ± 0.54 4.14 ± 0.16

0.113 ± 0.006 0.12 ± 0.004 0.104 ± 0.006 0.128 ± 0.006 0.148 ± 0.006

3.75 ± 0.2 32.93 ± 1.25 81.02 ± 4.46 32.75 ± 1.68 12.49 ± 0.57

Experiment H0 (kA∙m

−1

CSM

30.3 30.3 30.3 17.2 10.1

P (T ) = Pmax e

SAR (W∙g

)

37.61 ± 32.77 ± 26.45 ± 10.69 ± 4.080.19

1.98 1.25 1.46 0.55

ILP (nH∙m∙kg 0.149 0.130 0.105 0.128 0.145

± ± ± ± ±

−1

)

0.008 005 0.006 0.007 0.007

4

P (10 ∙W∙m

SIM −3

)

4.17 35.67 97.74 38.25 14.87

(T TPmax )2 2 2

−1

SAR (W∙g 41.67 35.38 31.81 12.45 4.84

mMNPs

(21)

30.3 30.3 30.3 17.2 10.1

1 10 30 30 30

SAR Relative Error % ISM

BLM

CSM

INCAM

−14.99 −20.63 −32.29 −29.08 −17.66

−25.28 −17.85 −30.08 −26.35 −14.98

−1.54 −11.12 −29.35 −26.73 −16.22

9.07 −4.03 −15.03 −14.67 −0.62

P (104∙W∙m−3)

SAR (W∙g−1)

ILP(nH∙m∙kg−1)

3.83 37.26 115.29 44.94 15.01

38.2 36.87 37.44 14.59 4.87

0.151 0.146 0.148 0.174 0.164

The P and SAR values acquired through each method and for each experiment are summarized in Table 2. The initial values P(T0) and SAR (T0) derived from the optimization simulations were used throughout the comparison section. The relative errors of each method with respect to the reference simulations are listed in Table 3. Our results indicate that the INCAM better approximates SAR with respect to our model. At the same time, the ISM, BLM, and CSM exhibit systematically higher relative errors. It is evident that in most of the examined cases, all four calculation methods underestimate SAR. In Fig. 9, the SAR values for different concentrations, and different AMF properties are illustrated. In Table 4 the effect of time interval selection on the ISM calculations is demonstrated. In Fig. 10 the CSM is illustrated. A decrease in SAR with increasing concentrations and an increasing SAR with increasing field is derived. The discrepancies in ILP values indicate that within the studied range of

Table 3 Relative errors for all four methods with respect to the simulations. C (mg∙mL−1)

)

3. Results

Prior to the optimization studies a mesh independence study was carried out for physics controlled meshes with the tested resolutions

H0 (kA∙m−1)

ILP (nH∙m∙kg

−1

ranging from coarse to extremely fine (Table 1), indicating no effect of the element size above the normal refinement (Fig. 8). Therefore, to reduce computational costs, a normal mesh was assigned throughout the geometry. The mesh consisted of 3067 domain elements plus 311 boundary elements. The control model included 5590 plus 486 internal degrees of freedom. The solution time of the optimization studies ranged 3–15 min depending on the initialization of the control variables.

(20)

mmf mf

)

0.165 0.140 0.126 0.149 0.173

Using this formalism, the distribution parameters are employed as the control variables of the optimization studies. The solution returns the distribution parameters that provide the best fitting of the simulated heating curve to the experimental one, and the effective branch of the Gaussian function, which quantifies the magnetic fluid’s power dissipation within the measurement’s temperature range. Temperaturedependent SAR, is derived by:

SAR (T ) = P (T )

0.007 0.005 0.006 0.006 0.005

INCAM −1

)

± ± ± ± ±

Fig. 9. (a) SAR vs concentration, and (b) SAR vs H0 for each calculation method. The error bars represent the standard deviation of the mean SAR values. 46

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al.

coil. According to our simulation, the solenoid used in the experiments generated a highly homogenous magnetic field. The optical fiber’s tip was immersed at the center of the magnetic fluid’s volume as suggested for cylindrical flat-based vials. Other studies recommend the positioning of the probe at the point of maximum temperature instead [22]. Yet in our approach, the experiment is backed by simulations that compute the temperature distribution; thus, the systematic positioning of the probe at the same point is what matters and not locating the maximum temperature region. The analysis of the obtained heating curves was carried out with respect to each of the examined calculation method’s restrictions. In the ISM calculations, different methods have been reported [21,38], including differentiation of linear or polynomial fittings, or direct calculations of ΔT/Δt, while the selected intervals range from tens to hundreds of seconds. Others select maximum or constant slopes. The ISM assumes an adiabatic process; therefore, by definition, it refers to a time limit near zero. For this condition to be fulfilled, a high sampling rate is ideal, and the fitting must be applied within the first seconds of the heating process, upon removal of possible delayed heating effects, to avoid underestimation in SAR calculations. Although a relatively small interval of 20 s was selected, considering the quasiadiabatic phase derived from the INCAM, the ISM significantly underestimated the heating efficiency. An increase in relative error with increasing heating efficiencies was observed. This behavior may be attributed to the more spontaneous development of energy loss mechanisms in more efficient systems. We suggest that the ISM should be carefully used, as it is highly sensitive to the interval selection as illustrated in Table 4. Despite the proper selection of time intervals for the Box-Lucas fittings, considering the linear loss region as well as the maintenance of reasonably low standard errors for the fitting parameters, the BLM exhibited higher relative errors compared with the CSM and INCAM, and similar with the ISM. However, a correlation between heating efficiency and relative error is not obvious. It is important to point out that the BLM is reported [54] to be a better estimate than the ISM and CSM, that are approximations of the phenomenological model. Still, the improper application of the BLM may lead to miscalculations in SAR values. Allowing the model to saturate beyond the linear loss region, evaluating small temperature gradients, or selecting limited time intervals [44] to obtain higher SAR estimates constitute possible misuses of the BLM. The fitting parameters must always be accompanied with a reasonably low standard error; otherwise, a fitting of the model cannot be considered successful.

Table 4 % Deviation from SAR values obtained for 20 s interval. H0 (kA∙m−1)

c (mg∙mL−1)

30 s interval

60 s interval

100 s interval

30.3 30.3 30.3 17.2 10.1

1 10 30 30 30

−10.81 −1.28 −3.10 −1.60 −1.09

−12.16 −5.55 −8.97 −6.42 −5.46

−12.16 −12.23 −12.81 −10.87 −9.93

MNPs’ concentrations and the used AMF properties, the LRT is invalid. In Fig. 11, the effect of interval selection on the BLM and ISM is illustrated: In Table 5 the optimal fitting parameters of the numerical models are summarized. In Fig. 12, the curves obtained assuming variable power dissipation are superimposed on the experimental heating curves. The effective branches of the respective Gaussian P(T) and SAR(T) functions are also illustrated. From the magnetic field distribution simulations, Eq. (11) yields a low inhomogeneity of 3.03%, within the magnetic fluid’s volume, which implies negligible effects of the magnetic field variation on the performance of the magnetic fluid. 4. Discussion In this study, we developed a well-elaborated model that considers most of the phenomena involved in the heating process, including thermal losses due to conduction, convection, and radiation. The double coupling of the heat transfer equations with the Navier-Stokes equations, the definition of dynamic boundary conditions, and the employment of a robust, fully coupled solver ensure a more detailed description of the process compared with simplified models or analysis limited to the heating curve. Another limitation of the calculation methods that we addressed through simulations is the dependence of the magnetic fluid’s thermal properties on temperature. In addition, the contribution of the MNPs to the solvent’s thermal properties were accounted for. In the experimental design, we followed best practices regarding the proper conduction of non-adiabatic measurements [42]. Particularly, we defined the linear loss region using cooling experiments, and the baseline temperature, considering the microenvironment within the

Fig. 10. Application of the CSM for five different experiments. (a) SAR values derived from 5 s intervals along 50 s timeframes of maximum efficiency, and (b) respective ILP values. 47

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al.

Fig. 11. Residual analysis showcasing (a–e) proper and (f–j) improper use of the BLM. Residual analysis showcasing (k–o) proper and (p–t) improper use of the ISM.

In a recent study [54], Lanier et al. compared the ISM, the BLM, and the CSM throughout a high volume of measurements and concluded that the BLM is the most accurate approach. The authors used residual analysis to evaluate each method’s goodness of fitting on experimental data. Accordingly, to assess the validity of the ISM and BLM, we carried out residual analysis for each measurement. Our results concurred that the BLM is more efficient from a statistical perspective. Interesting results were extracted by selecting different increments to conduct the analysis. As illustrated in Fig. 11(a–e), in the BLM the first 80 s seemed a proper interval selection with minimal evidence of residuals’ trends and reasonably low standard errors for the fitting parameters. Allowing further saturation of the model for 300 s (Fig. 11(f–j)), further reduced the standard errors of the fitting parameters but resulted in clear trend observation in the residual plots. In the ISM method, while for limited selections of 20 s no clear trends were evident (Fig. 11(k–o)), regular residuals exhibited trends for 60 s interval selections (Fig. 11(p–t)). Therefore, even an approximation of the BLM, the ISM may deliver

comparable precision if used adequately as reflected in our results. While residual analysis ensures statistical relevance, in this study we mainly focused in using simulations to compute reference SAR(T) functions and investigate if SAR calculation methods based on calorimetric data deliver any predictive capacity. Regarding the CSM method, selecting ten intervals of 5 s with a 0.5 s sampling rate along the 50 s timeframe of the maximum efficiency, seems a proper approach. Our results agree with the statement by Wildeboer et al. [42] that the CSM and BLM provide better estimates compared with the ISM. Yet, considering the limited number of the curves that were subjected to our analysis, the method of highest efficiency is still debatable. The CSM is also useful to evaluate the compliance of the experiments with linear response theory. In Fig. 10(b), the ILP plots indicate that the linear response theory limitations, regarding the magnetic field, and the inter-particle interactions, are violated. These findings were verified via simulations (Fig. 12(e)). For the INCAM, the statistical analysis led to the conclusion that selecting increments of 2 s yielded the highest median SAR (SAR50%) and third quartile SAR (SAR75%) values, with reasonably low variance as reflected in box and whisker plots (Fig. 4(c)). The SAR variance decreases as the increments’ duration increases [21]. However, because of non-negligible losses that the linear loss parameter and the quasiadiabatic criterion are not sensitive to, SAR50% and SAR75% decrease with increasing increments. The selection of a representative SAR value among the wide range provided from the statistical analysis is deterministic. For extremely small increments (0.5–1 s) and low heating efficiencies, negative numerical derivatives are apparent in the statistical analysis yielding the extraction of SAR values infeasible and physically invalid. The prominent contribution of the present study is not limited to the development of a detailed numerical model that delineates the non-

Table 5 Optimal Gaussian function parameters derived by simulations. H0 (kA∙m−1)

c (mg∙mL−1)

Pmax (104 W∙m−3)

TPmax (K)

σ (K)

30.3 30.3 30.3 17.2 10.1

1 10 30 30 30

5.047* 52.559 27.609 58.309 127.000

284.96* 310.61 303.32 310.46 312.85

12.046* 19.846 8.97 22.946 42.264

*The fitting parameters are affected by the initial conditions due to low and unstable temperature increase.

48

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al.

Fig. 12. (a, b) Simulated curves obtained from the optimization studies, superimposed on the experimental data. (c) P, (d) SAR, and (e) ILP plotted against temperature. (f) Normalized SAR for the measurement, which exhibited the highest SAR variation ~45% with temperature.

49

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al.

adiabatic nature of the heating process. We introduce the scenario of precise SAR(T) calculations in silico, through inverse optimization studies. The main advantage of our methodology over theoretical models is that no characterization of the MNPs is required. Our initial optimization studies indicated that a constant value of power dissipation is not suitable for quantifying the heating efficiency. This finding agrees with previous studies [39,40,47–50], in which the dependence of SAR on temperature was demonstrated, both theoretically and experimentally. Particularly, Natividad et al. [40] employed the pulse heating method in adiabatic conditions to demonstrate the significant discrepancy of SAR with temperature. The calorimetric results were verified via AC susceptibility measurements. In a later study, the authors [39] used the same calorimetric method to determine SAR (T) for self-regulated manganite nanoparticles. In recent studies, Garaio et al. used AC magnetometry to examine the SAR variation with temperature [48,50]. Their results demonstrated both ascending and descending trends in heating efficiency with temperature increase, and SAR variations reached 40% within MFH-relevant temperature windows. SAR variations that exceed 50% have also been reported by Regmi et al. [49]. Landi [47] discussed the impact of the thermal energy and sample’s polydispersity on SAR in a theoretical analysis based on the Neel relaxation. The unblocking process of MNPs [70,71] may explain observations of increasing SAR values at the start of the heating process. As smaller MNPs deposit energy to the medium, larger MNPs that are blocked at room temperature transition to the superparamagnetic state. Their contribution yields an increasing trend in the overall power dissipation. However, this trend diminishes as the unblocking process saturates and equilibrium between new contributions and net power dissipation depreciation is restored. At this point, a maximum value Pmax is reached, and eventually the heating efficiency of the system decreases. The process described above is governed by the size distribution parameters [24], the magnetic properties of the MNPs [49,70], and the AMF properties [39,50]. In any case, it can be well described by a function of temperature that facilitates an ascending and a descending part. Therefore, we suggest the expression of power dissipation as a Gaussian function of temperature. Normally, determining the function parameters as well as the effective branch within the measurements’ temperature range would be challenging. Using the optimization algorithms provided in the simulations’ software resolves this issue. The importance of considering the temperature dependence of SAR is reflected by our results that demonstrate significant variance in heating efficiency with temperature. For example, the 30 mg∙mL−1 sample heated under the highest AMF exhibits a 10% increase in SAR between the initial temperature and the therapeutic range. At 10.1 kA∙m−1 SAR (T) variations up to 45% were observed (Fig. 12(f)). Our results concur with previous findings, in which significant SAR variations within MFH-relevant temperature windows have been recorded [48–50]. However, our approach detours the complex and time-consuming adiabatic setups or AC susceptibility measurements and requires only a proper description of the experimental setup and conditions. Implementing SAR(T) in biological applications is important. Ultimately, in an ideal scenario in which control over the MNPs’ distribution within the ROI has been accomplished, ignoring the dependence of SAR on temperature leads to overestimation in prescribed doses, inducing increased toxicity levels and possible overheating that could be avoided. Expanding our calculations to a wider range of measurements could potentially provide indirect insights to the MNPs’ intrinsic properties through the magnetic fluid’s macroscopic response. Yet the method is bound to the modelling assumptions and requires detailed information on the experimental setup and conditions. The development of protocols that standardize the experimental process would address this issue and enable effective inter-laboratory analysis of hyperthermia data via protocol-specific simulations.

5. Conclusions In this work, a robust method is presented to derive the temperature-dependent heating efficiency SAR(T) of a magnetic fluid by fitting a sophisticated numerical model to the experimental heating curve. We generated reference data in silico to evaluate the SAR estimates of prominent calculation methods including ISM, BLM, CSM, and INCAM. The results of this study indicate that the INCAM exhibited the lowest relative errors with respect to the simulations’ results. We used residual analysis on the BLM and ISM to demonstrate the significance of the selected time intervals while fitting the models to experimental data. In terms of predictive power, all four SAR calculation methods were susceptible to errors introduced due to experimental sources, insensitivity to heat transfer mechanisms, and neglect of the SAR dependence on temperature. To avoid such errors, we recommend employment of simulations and development of universal measurement protocols that enable proper inter-laboratory analysis. Our results affirm that temperature-dependent SAR should be considered in treatment planning applications, as significant variations occur within the MFH-relevant temperature range. Our prospective studies include the complete power dissipation mapping of efficient magnetic fluids throughout a wide range of concentrations. In addition, the development of numerical models aimed for simulations of biological systems, combining mapped magnetic fluids with computational phantoms such as XCAT or MOBY [72–77]. Such models could be exploited to improve thermal dose delivery by optimizing the MNPs’ distribution within the ROI. The effect of the tumor’s site and the surrounding tissues’ thermal properties on the required MNPs’ concentrations to reach therapeutic levels, could also be investigated. Verifying such models through image-guided in vivo experiments constitutes an important step in progress. Such accomplishments would contribute in the development of precise hyperthermia treatment planning applications. Acknowledgments This research is co-financed by Greece and the European Union (European Social Fund- ESF) through the Operational Programme “Human Resources Development, Education and Lifelong Learning” in the context of the project “Strengthening Human Resources Research Potential via Doctorate Research” (MIS-5000432), implemented by the State Scholarships Foundation (ΙΚY). Special thanks to Dr. George Loudos, Assistant Professor at the University of West Attica, for providing scientific insights and fruitful discussions during this study. Editorial support was provided by Bryan Tutt in UT MD Anderson Cancer Center, Scientific Publications Services, Research Medical Library. Declaration of Competing Interest There is no conflict of interest to declare. References [1] Deatsch AE, Evans BA. Heating efficiency in magnetic nanoparticle hyperthermia. J Magn Magn Mater 2014;354:163–72. [2] Jordan A, Scholz R, Wust P, FaKhling H, Felix R. Magnetic fluid hyperthermia (MFH): cancer treatment with AC magnetic field induced excitation of biocompatible superparamagnetic nanoparticles. J Magn Magn Mater 1999;201:413–9. [3] Habash RWY. Therapeutic hyperthermia. Handb Clin Neurol 2018;157:853–68. [4] Chicheł A, Skowronek J, Kubaszewska M, Kanikowski M. Hyperthermia – description of a method and a review of clinical applications. Rep Pract Oncol Radiother 2007;12:267–75. [5] Dutz S, Hergt R. Magnetic particle hyperthermia–a promising tumour therapy? Nanotechnology 2014;25:452001. [6] Thiesen B, Jordan A. Clinical applications of magnetic nanoparticles for hyperthermia. Int J Hyperthermia 2008;24:467–74. [7] Dhavalikar R, Bohórquez AC, Rinaldi C. Image-guided thermal therapy using magnetic particle imaging and magnetic fluid hyperthermia. Nanomater Magn Opt Hyperthermia Appl 2019:265–86.

50

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al. [8] Du Y, Liu X, Liang Q, Liang XJ, Tian J. Optimization and design of magnetic ferrite nanoparticles with uniform tumor distribution for highly sensitive MRI/MPI performance and improved magnetic hyperthermia therapy. Nano Lett 2019;19:3618–26. [9] Santos MA, Wu SK, Li Z, Goertz DE, Hynynen K. Microbubble-assisted MRI-guided focused ultrasound for hyperthermia at reduced power levels. Int J Hyperthermia 2018;35:599–611. [10] Ling-Yun Z, Jia-Yi L, Wei-Wei O, Dan-Ye L, Li L, Li-Ya L, et al. Magnetic-mediated hyperthermia for cancer treatment: research progress and clinical trials. Chin Phys B 2013;22:108104. [11] Latorre M, Rinaldi C. Applications of magnetic nanoparticles in medicine: magnetic fluid hyperthermia. P R Health Sci J 2009;28:227–38. [12] Torres-Lugo M, Rinaldi C. Thermal potentiation of chemotherapy by magnetic nanoparticles. Nanomedicine 2013;8:1689–707. [13] Kappiyoor R, Liangruksa M, Ganguly R, Puri IK. The effects of magnetic nanoparticle properties on magnetic fluid hyperthermia. J Appl Phys 2010;108:094702. [14] Laurent S, Dutz S, Häfeli UO, Mahmoudi M. Magnetic fluid hyperthermia: focus on superparamagnetic iron oxide nanoparticles. Adv Colloid Interface Sci 2011;166:8–23. [15] Pankhurst QA, Connolly J, Jones SK, Dobson J. Applications of magnetic nanoparticles in biomedicine. J Phys D Appl Phys 2003;36:167–81. [16] Jeong U, Teng X, Wang Y, Yang H, Xia Y. Superparamagnetic colloids: controlled synthesis and niche applications. Adv Mater 2007;19:33–60. [17] Neuberger T, Schöpf B, Hofmann H, Hofmann M, von Rechenberg B. Superparamagnetic nanoparticles for biomedical applications: possibilities and limitations of a new drug delivery system. J Magn Magn Mater 2005;293:483–96. [18] Pankhurst QA, Thanh NTK, Jones SK, Dobson J. Progress in applications of magnetic nanoparticles in biomedicine. J Phys D Appl Phys 2009;42:224001. [19] Arruebo M, Fernández-Pacheco R, Ibarra MR, Santamaría J. Magnetic nanoparticles for drug delivery. Nano Today 2007;2:22–32. [20] Kozissnik B, Bohorquez AC, Dobson J, Rinaldi C. Magnetic fluid hyperthermia: advances, challenges, and opportunity. Int J Hyperthermia 2013;29:706–14. [21] Soetaert F, Kandala SK, Bakuzis A, Ivkov R. Experimental estimation and analysis of variance of the measured loss power of magnetic nanoparticles. Sci Rep 2017;7:6661. [22] Wang S-Y, Huang S, Borca-Tasciuc D-A. Potential sources of errors in measuring and evaluating the specific loss power of magnetic nanoparticles in an alternating magnetic field. IEEE Trans Magn 2013;49:255–62. [23] Kallumadil M, Tada M, Nakagawa T, Abe M, Southern P, Pankhurst QA. Suitability of commercial colloids for magnetic hyperthermia. J Magn Magn Mater 2009;321:1509–13. [24] Rosensweig RE. Heating magnetic fluid with alternating magnetic field. J Magn Magn Mater 2002;252:370–4. [25] Ruta S, Chantrell R, Hovorka O. Unified model of hyperthermia via hysteresis heating in systems of interacting magnetic nanoparticles. Sci Rep 2015;5:9090. [26] Carrey J, Mehdaoui B, Respaud M. Simple models for dynamic hysteresis loop calculations of magnetic single-domain nanoparticles: application to magnetic hyperthermia optimization. J Appl Phys 2011;109:083921. [27] Branquinho LC, Carriao MS, Costa AS, Zufelato N, Sousa MH, Miotto R, et al. Effect of magnetic dipolar interactions on nanoparticle heating efficiency: implications for cancer hyperthermia. Sci Rep 2013;3:2887. [28] Mehdaoui B, Meffre A, Carrey J, Lachaize S, Lacroix L-M, Gougeon M, et al. Optimal size of nanoparticles for magnetic hyperthermia: a combined theoretical and experimental study. Adv Funct Mater 2011;21:4573–81. [29] Carrião MS, Aquino VRR, Landi GT, Verde EL, Sousa MH, Bakuzis AF. Giant-spin nonlinear response theory of magnetic nanoparticle hyperthermia: a field dependence study. J Appl Phys 2017;121:173901. [30] Castillo VLC-Dd, Rinaldi C. Effect of sample concentration on the determination of the anisotropy constant of magnetic nanoparticles. IEEE Trans Magn 2010;46:852–9. [31] Aslibeiki B, Kameli P, Salamati H. The effect of dipole-dipole interactions on coercivity, anisotropy constant, and blocking temperature of MnFe2O4 nanoparticles. J Appl Phys 2016;119:063901. [32] Kolhatkar AG, Jamison AC, Litvinov D, Willson RC, Lee TR. Tuning the magnetic properties of nanoparticles. Int J Mol Sci 2013;14:15977–6009. [33] Osaci M, Cacciola M. An adapted Coffey model for studying susceptibility losses in interacting magnetic nanoparticles. Beilstein J Nanotechnol 2015;6:2173–82. [34] Cabrera D, Lak A, Yoshida T, Materia ME, Ortega D, Ludwig F, et al. Unraveling viscosity effects on the hysteresis losses of magnetic nanocubes. Nanoscale 2017;9:5094–101. [35] Alkadour B, Southern BW, Whitehead JP, van Lierop J. Triangular array of iron oxide nanoparticles: simulation study of intraparticle and interparticle magnetism. Phys Rev B 2019:100. [36] Tan RP, Carrey J, Respaud M. Magnetic hyperthermia properties of nanoparticles inside lysosomes using kinetic Monte Carlo simulations: influence of key parameters and dipolar interactions, and evidence for strong spatial variation of heating power. Phys Rev B 2014;90. [37] Taddei C, Sansone L, Ausanio G, Iannotti V, Pepe GP, Giordano M, et al. Fabrication of polystyrene-encapsulated magnetic iron oxide nanoparticles via batch and microfluidic-assisted production. Colloid Polym Sci 2019;297:861–70. [38] Andreu I, Natividad E. Accuracy of available methods for quantifying the heat power generation of nanoparticles for magnetic hyperthermia. Int J Hyperthermia 2013;29:739–51. [39] Natividad E, Castro M, Goglio G, Andreu I, Epherre R, Duguet E, et al. New insights into the heating mechanisms and self-regulating abilities of manganite perovskite nanoparticles suitable for magnetic fluid hyperthermia. Nanoscale 2012;4:3954–62.

[40] Natividad E, Castro M, Mediano A. Adiabatic magnetothermia makes possible the study of the temperature dependence of the heat dissipated by magnetic nanoparticles under alternating magnetic fields. Appl Phys Lett 2011;98:243119. [41] Natividad E, Castro M, Mediano A. Adiabatic vs. non-adiabatic determination of specific absorption rate of ferrofluids. J Magn Magn Mater 2009;321:1497–500. [42] Wildeboer RR, Southern P, Pankhurst QA. On the reliable measurement of specific absorption rates and intrinsic loss parameters in magnetic hyperthermia materials. J Phys D Appl Phys 2014;47:495003. [43] Huang S, Wang SY, Gupta A, Borca-Tasciuc DA, Salon SJ. On the measurement technique for specific absorption rate of nanoparticles in an alternating electromagnetic field. Meas Sci Technol 2012;23:035701. [44] Bordelon DE, Cornejo C, Grüttner C, Westphal F, DeWeese TL, Ivkov R. Magnetic nanoparticle heating efficiency reveals magneto-structural differences when characterized with wide ranging and high amplitude alternating magnetic fields. J Appl Phys 2011;109:124904. [45] Teran FJ, Casado C, Mikuszeit N, Salas G, Bollero A, Morales MP, et al. Accurate determination of the specific absorption rate in superparamagnetic nanoparticles under non-adiabatic conditions. Appl Phys Lett 2012;101:062413. [46] Attaluri A, Nusbaum C, Wabler M, Ivkov R. Calibration of a quasi-adiabatic magneto-thermal calorimeter used to characterize magnetic nanoparticle heating. J Nanotechnol Eng Med 2013;4. [47] Landi GT. Simple models for the heating curve in magnetic hyperthermia experiments. J Magn Magn Mater 2013;326:14–21. [48] Garaio E, Sandre O, Collantes JM, Garcia JA, Mornet S, Plazaola F. Specific absorption rate dependence on temperature in magnetic field hyperthermia measured by dynamic hysteresis losses (ac magnetometry). Nanotechnology 2015;26:015704. [49] Regmi R, Naik A, Thakur JS, Vaishnava PP, Lawes G. Temperature dependent dissipation in magnetic nanoparticles. J Appl Phys 2014;115:17B301. [50] Garaio E, Collantes J-M, Garcia JA, Plazaola F, Sandre O. Harmonic phases of the nanoparticle magnetization: an intrinsic temperature probe. Appl Phys Lett 2015;107:123103. [51] Massart R. Preparation of aqueous magnetic liquids in alkaline and acidic media. IEEE Trans Magn 1981;17:1247–8. [52] Cullity BD. Elements of X-ray diffraction. Addison-Wesley Publishing Company Inc.; 1956. [53] Scherrer P. Göttinger Nachrichten. Math. Phys. 1918;2:98–100. [54] Lanier OL, Korotych OI, Monsalve AG, Wable D, Savliwala S, Grooms NWF, et al. Evaluation of magnetic nanoparticles for magnetic fluid hyperthermia. Int J Hyperthermia 2019;36:687–701. [55] Corrected slope SAR and ILP analysis programmes, written in Microsoft Excel and Matlab formats. [56] Samchenko Y, Korotych O, Kernosenko L, Kryklia S, Litsis O, Skoryk M, et al. Stimuli-responsive hybrid porous polymers based on acetals of polyvinyl alcohol and acrylic hydrogels. Colloids Surf A 2018;544:91–104. [57] Kashevsky BE, Kashevsky SB, Terpinskaya TI, Ulashchik VS. Magnetic hyperthermia with hard-magnetic nanoparticles: in vivo feasibility of clinically relevant chemically enhanced tumor ablation. J Magn Magn Mater 2019;475:216–22. [58] Kennedy MW, Akhtar S, Bakken JA, Aune RE. Analytical and experimental validation of electromagnetic simulations using COMSOL®, re inductance, induction heating and magnetic fields. COMSOL Conference. Stuttgart2011. [59] Amestoy PR, Duff IS, L’Excellent J-Y, Koster J. A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM J Matrix Anal Appl 2001;23:15–41. [60] Di Barba P, Dughiero F, Sieni E, Candeo A. Coupled field synthesis in magnetic fluid hyperthermia. IEEE Trans Magn 2011;47:914–7. [61] Borosilicate Glass Properties. SCHOTT AG. [62] Borosilikatglas BOROFLOAT® – Thermische Produkteigenschaften | SCHOTT AG. SCHOTT AG. [63] Kol HŞ, Sefil Y. The thermal conductivity of fir and beech wood heat treated at 170, 180, 190, 200, and 212°C. J Appl Polym Sci 2011;121:2473–80. [64] Troppová E, Tippner J, Hrčka R, Halachan P. Quasi-stationary measurements of lignamon thermal properties. BioResources 2013;8:6288–96. [65] Wu L, Cheng J, Liu W, Chen X. Numerical analysis of electromagnetically induced heating and bioheat transfer for magnetic fluid hyperthermia. IEEE Trans Magn 2015;51:4600204. [66] Sundar LS, Singh MK, Sousa ACM. Investigation of thermal conductivity and viscosity of Fe3O4 nanofluid for heat transfer applications. Int Commun Heat Mass 2013;44:7–14. [67] Jones JM, Mason PE, Williams A. A compilation of data on the radiant emissivity of some materials at high temperatures. J Energy Inst 2019;92:523–34. [68] The JJM. Levenberg Marquardt algorithm: implementation and theory. Numer Anal 1978;630:105–16. [69] Schenk O, Gärtner K. Solving unsymmetric sparse systems of linear equations with PARDISO. Future Gener Comput Syst 2004;20:475–87. [70] Shim H, Dutta P, Seehra MS, Bonevich J. Size dependence of the blocking temperatures and electron magnetic resonance spectra in NiO nanoparticles. Solid State Commun 2008;145:192–6. [71] Mohapatra J, Zeng F, Elkins K, Xing M, Ghimire M, Yoon S, et al. Size-dependent magnetic and inductive heating properties of Fe3O4 nanoparticles: scaling laws across the superparamagnetic size. Phys Chem Chem Phys 2018;20:12879–87. [72] Segars WP, Norris H, Sturgeon GM, Zhang Y, Bond J, Minhas A, et al. The development of a population of 4D pediatric XCAT phantoms for imaging research and optimization. Med Phys 2015;42(8):4719–26. https://doi.org/10.1118/1.492684. [73] Norris H, Zhang Y, Bond J, Sturgeon GM, Minhas A, Tward DJ, et al. A set of 4D pediatric XCAT reference phantoms for multimodality research. Med Phys 2014;41(43):033701https://doi.org/10.1118/1.4864238.

51

Physica Medica 71 (2020) 39–52

C. Papadopoulos, et al. [74] Kostou T, Papadimitroulas P, Papaconstadopoulos P, Devic S, Seuntjens J, Kagadis GC. Size-specific dose estimations for pediatric chest, abdomen/pelvis and head CT scans with the use of GATE. Phys Med 2019 Sep;65:181–90. https://doi.org/10. 1016/j.ejmp.2019.08.020. [75] Papadimitroulas P, Erwin WD, Iliadou V, Kostou T, Loudos G, Kagadis GC. A personalized, Monte Carlo-based method for internal dosimetric evaluation of radiopharmaceuticals in children. Med Phys 2018. https://doi.org/10.1002/mp.13055.

[76] Kostou T, Papadimitroulas P, Loudos G, Kagadis GC. A preclinical simulated dataset of S-values and investigation of the impact of rescaled organ masses using the MOBY phantom. Phys Med Biol 2016;61(6):2333–55. [77] Segars WP, Tsui BM, Frey EC, Johnson GA, Berr SS. Development of a 4D digital mouse phantom for molecular imaging research. Mol Imaging Biol 2004;6:149–59.

52