Temperature-controlled power modulation compensates for heterogeneous nanoparticle distributions: a computational optimization analysis for magnetic hyperthermia

Purpose: To study, with computational models, the utility of power modulation to reduce tissue temperature heterogeneity for variable nanoparticle distributions in magnetic nanoparticle hyperthermia. Methods: Tumour and surrounding tissue were modeled by elliptical two- and three-dimensional computational phantoms having six different nanoparticle distributions. Nanoparticles were modeled as point heat sources having amplitude-dependent loss power. The total number of nanoparticles was fixed, and their spatial distribution and heat output were varied. Heat transfer was computed by solving the Pennes’ bioheat equation using finite element methods (FEM) with temperature-dependent blood perfusion. Local temperature was regulated using a proportional-integral-derivative (PID) controller. Tissue temperature, thermal dose and tissue damage were calculated. The required minimum thermal dose delivered to the tumor was kept constant, and heating power was adjusted for comparison of both the heating methods. Results: Modulated power heating produced lower and more homogeneous temperature distributions than did constant power heating for all studied nanoparticle distributions. For a concentrated nanoparticle distribution, located off-center within the tumor, the maximum temperatures inside the tumor were 16% lower for modulated power heating when compared to constant power heating. This resulted in less damage to surrounding normal tissue. Modulated power heating reached target thermal doses up to nine-fold more rapidly when compared to constant power heating. Conclusions: Controlling the temperature at the tumor-healthy tissue boundary by modulating the heating power of magnetic nanoparticles demonstrably compensates for a variable nanoparticle distribution to deliver effective treatment.


Introduction
Hyperthermia involves raising and sustaining the temperature of malignant tumors and adjacent tissues to about 41-46 °C for ~0.5 to 2 h to achieve a therapeutic effect [1,2]. The principal challenges encountered with all thermal therapies are to effectively deliver and control energy deposition. Magnetic nanoparticle hyperthermia (mNPH) has emerged as a treatment modality that offers benefits because the heat sources, i.e., nanoparticles, can be embedded within the target tissue [3,4]. The region containing the magnetic nanoparticles is then exposed to an alternating magnetic field (AMF) which generates heat via magnetic hysteresis loss [3,[5][6][7][8][9]. Heat generated by the nanoparticles is transferred by convective and conductive processes to the tumor, while offering the potential to minimize energy deposition outside the treatment margins. Thus, mNPH can overcome some limitations for hyperthermia that exist with other heating technologies because the energy is generated and controlled within the tumor.
While offering advantages, mNPH also presents challenges to clinical implementation. Nanoparticle delivery can be variable, and the resulting nanoparticle distribution within tumors is often heterogeneous [10][11][12][13][14]. The heterogeneity of nanoparticle distributions thus presents new challenges for precise energy delivery which can lead to off-target heating of healthy tissues or under-heating of regions within the tumor [14]. It has been recognized that clinical acceptance of mNPH has been slowed, in part, by the imprecise spatial control of the nanoparticles and energy within the tumor [10][11][12][13].
progress has been made to model various complex biological and physical factors relevant to mNPH. Most models (1) assume homogeneous or regular nanoparticle distributions within the model tissue to reduce computational complexity; and, (2) neglect complexities of variable capillary blood perfusion, an important physiologic response to temperature that dissipates heat energy [15,17,[24][25][26][27][28][29]. Typically, the focus of computational studies of mNPH is to determine optimal magnetic field conditions or magnetic nanoparticle parameters, i.e., specific loss power (SLP), to achieve a target temperature. Capillary blood perfusion, when included in computational studies for thermal therapy applications is typically assumed a non-physiologically constant value.
Effective hyperthermia is achieved when heat energy is controlled to minimize damaging energy deposition outside the target tissue while also achieving and maintaining an elevated temperature within the tumor. Controlling the temperature at the tumor-healthy tissue boundary offers potential to preserve normal tissue while simultaneously enabling therapeutic heating within tumors. However, controlling heat deposition with temperaturefeedback control at the tumor-tissue boundary may limit therapeutic heating throughout the tumor when the nanoparticle distribution is inhomogeneous [14]. Furthermore, locally high nanoparticle concentrations deposit sufficient energy to generate 'runaway' heating with local ablation and excess energy deposition that threatens to overheat adjacent normal tissues. Thermal dose regulation requires a responsive control system to rapidly modulate energy deposition with appropriate temperature feedback.
Proportional-integral-derivative (PID) control systems are widely used in industry and in medicine to regulate input power for heating [30][31][32][33][34]. Salomir et al. used a PID temperature controller to regulate the temperature at a focal point during MR-guided focused ultrasound therapy [30]. They obtained temperature control at nearly the precision of the temperature measurement, in vitro and in vivo. Mougenot et al. further developed three-dimensional spatial and temporal control of temperature during MR thermometry-guided focused ultrasound therapy [31]. Haemmerich et al. implemented an automatic PI closed-loop controller system and optimized the controller parameters for temperature controlled RF ablation [32].
Here, we describe a strategy to modulate AMF amplitude to regulate nanoparticle SLP using single-point temperature feedback with a PID controller to compensate for variable nanoparticle distributions. Our computational 2D and 3D phantoms incorporated a temperature-dependent capillary blood perfusion model that reduces with significant thermal dose. We considered six nanoparticle distributions and we compared time-dependent changes of temperature and thermal dose, the latter defined by cumulative equivalent minutes at 43 °C (CEM43). In our simulated models, we study the use of PID-controlled power modulation to reduce temperature heterogeneity in the tumor; and, to achieve a more rapid approach to the therapeutic thermal dose compared to constant power heating.

Nanoparticle distributions
In the present study, the magnetic nanoparticles were modeled as point heat sources [20]. Three distributions -E1, E2, and E3 in Figure 1(a) -were generated by digitizing (processed with MATLAB) images previously published from tissue sections obtained from human prostate tumor xenografts grown in mice which were injected with BNF-Starch (micro-mod Partikeltechnologie, GmbH, Rostock, Germany) nanoparticles [35]. An additional three were mathematical functions representing 'idealized' distributions -M1, M2, and M3 shown in Figure 1(b), and described below. The heating power of the phantom nanoparticles was modeled using previously reported specific loss power (SLP) values for magnetic iron oxide nanoparticles (see Figure S1) [36]. The corresponding alternating magnetic field amplitude can be estimated using the continuous polynomial approximation of SLP vs H given by Soetaert et al. [29].
To compare heating effects among the six tumor models, the total number of nanoparticles N (dark pixels) in a tumor was fixed at N = 1460 ±5 for each model [29]. For the present computational analysis, the area of a phantom 'nanoparticle' was fixed at ~35 pixels/mm 2 in the 2D models, as analyzed from imaging data obtained from a previous study and corresponding approximately to the known injected concentration of nanoparticles into human tumor xenografts (= 5 mg Fe/cm 3 of tumor) as reported in that study [14]. This means that each phantom tumor received a total heating power comparable to the others, allowing a quantitative comparison of temperatures and thermal dose. If the individual heating power of a phantom nanoparticle is Q i (=384.8 W-nanoparticle −1 ·m −3 ), and the total number of nanoparticles N (= 1460) distributed in the tumor is known (see Supplementary Materials for details of discretization of heating power and Figure S1), the total heating power Q NP power generated by the nanoparticles deposited in the tumor is then:

Idealized nanoparticle distributions
The three mathematical nanoparticle distributions were uniform (M1) [16], centrallyconcentrated (M2) [14], and Gaussian (M3) [15,28,29] (Figure 1(b)). The M1 distribution contained nanoparticles uniformly distributed throughout the tumor area. For the M2 distribution, nanoparticles were arranged uniformly in an area occupying only the central 40% region of the tumor. For M3, the nanoparticle distribution was governed by a Gaussian probability function originating at the tumor center and extending to the tumor-tissue boundary. The Gaussian probability function in Cartesian coordinates is where σ x = 2 (~44% of tumor major semi axis) and σ y = 0.5 (~17% of tumor minor semi axis) are arbitrarily selected standard deviations of x and y, respectively. The generated Gaussian distribution of points was truncated to fit within the tumor.

2D tumor model
We chose to develop, test and compare our strategy with constant power heating using a 2D tumor model. This approach helps us to efficiently test our strategy before translating our results to 3D and is recommended by American Society of Mechanical Engineers (ASME) [37][38][39]. Our chosen 2D model comprised two concentric elliptical regions, shown in Figure  1(c), representing the tumor and surrounding healthy tissue. The inner elliptical region (major axis l t = 9 mm, minor axis w t = 6 mm) represented the tumor seeded with nanoparticles. The outer elliptical region (major axis l = 19 mm, minor axis w = 16 mm) corresponded to modeled healthy tissue. We assumed that the outer boundary of the healthy tissue was in contact with the core body temperature (infinite heat reservoir) maintained at 37 °C. Heat transfer in the healthy tissue and tumor can be described by the Pennes' bioheat equation as [40] ρ n c n ∂T n (x, y, t) where the subscript n accounts for the tissue layer (n = 1 for the tumor and n = 2 for healthy tissue) and the subscript b for blood, respectively. ρ n is the respective tissue density (kg/m 3 ), c n the respective tissue specific heat (J/kg K), T n (x, y, t) the local tissue temperature (K), k n its thermal conductivity (W/m-K), Qm,n the metabolic heat generation rate in either the healthy tissue or tumor (W/m 3 ), and Q P is the heat rate per unit volume of tumor (W/m 3 ) that arises from nanoparticles if present in that volume (= 0 if no phantom nanoparticle is present). The subscript m in Equation (3) refers to metabolic heat generation in tissue layer n. ρ b , c b , ω b , n (x, y, T n ) and T b denote density, specific heat, perfusion rate, and temperature of blood, respectively. Table 1 summarizes the thermophysical properties used in this study [41][42][43][44][45][46]. We assumed no power deposition due to eddy currents by AMF interactions with tissue.
At the interface between the tumor and healthy tissue, continuity in temperature and heat flux is given by (see Figure 1(c)): In Equation (4), j denotes the surface normal to the elliptical boundary between tumor and healthy tissue. Initial temperature of tumor and healthy tissue was set to the core body temperature, i.e., T tissue boundary = 37°C.

Perfusion models
Effects of temperature on perfusion were tested using three models: (i) constant perfusion (i.e., 'classic' Pennes' equation) as described by Equation (3), where ω b (T) = ω bi (1/s); (ii) Arrhenius perfusion; and, (iii) modified Arrhenius perfusion, as described previously [47,48]. In the Arrhenius model, blood perfusion is described as a function of temperature and time, correlated to the degree of microvascular stasis (DS). DS is a dimensionless value expressed as In Equations (6) and (7), A is the frequency or pre-exponential factor (1/s), E a the activation energy (J/mol), R the universal gas constant (J/K·mol), T(τ) is the absolute tissue temperature as a function of time and ω bi the initial blood perfusion rate (1/sec). The degree of vascular stasis DS varies between 0 (no vascular damage) and 1 (complete vascular shutdown). Values for A and E a are shown in Table 1.
The foundation of the modified Arrhenius perfusion model is the first-order kinetic Arrhenius model of vascular stasis (Equations (6) and (7)), expanded by an additional term to increase the relative perfusion when vascular stasis is low. This latter model is derived empirically from data determined experimentally by He et al. [47]. Schutt When the degree of stasis is low (DS ≤ 0.02), there is an increase in perfusion. As the degree of stasis increases (0.02 < DS ≤ 0.08), perfusion first decreases moderately, and then more rapidly (0.08<DS ≤ 0.97), before eventually becoming zero (0.97<DS ≤ 1.0).

Constant power heating
For constant power simulations, the total heating power of nanoparticles, Q NP , was fixed. Temperature distribution T(x, y, t) and thermal dose were evaluated as a function of time for each simulation. The total thermal dose achieved in the tumor is expressed as [50]: In Equation (9), B is a constant with value 0.5 for T> 43 °C, and 0.25 for T≤ 43 °C [1].
CEM43 T90 is defined as the total isoeffect thermal dose expressed in cumulative equivalent minutes at 43 °C achieved or exceeded in 90% of the tumor area [1,51], with treatment considered effective when achieved for 60 min or longer, i.e., CEM43 T90 ≥ 60 [1,51]. The power required to achieve CEM43 T90 ≥ 60 within 20 min was computed. The heating power needed to achieve the target dose was defined as the isoeffect heating power, Q iso . It was computed by iteratively performing simulations with constant power settings until CEM43 T90 ≥ 60 was achieved.

Modulated power heating
For modulated-power heating, nanoparticle heating was varied as a function of the temperature computed at the tumor-tissue boundary [19,24]. Attaluri et al. have previously described methods to heat tumors using temperature-control at the tumor-tissue boundary to control energy deposition for a spherical tumor model with uniform distribution of nanoparticles [14]. Presently, we consider a 2D elliptical model, which accounts for variability in two directions, and different nanoparticle distributions. These modifications demand a more precise location of the temperature probe. To determine this, eight probe locations (P1-P8, Figure 2(a)) on the tumor-tissue boundary were considered. For all six models, the heating power was varied between two levels: a higher heating power when the probe temperature, T probe <43.5 °C, and a lower heating power when T probe >43.5 °C [14].
The heating algorithm is described mathematically as, The higher heating power was chosen to be greater than Q iso to achieve a rapid rise to target temperature in the tumor, which is clinically desired. A probe location that achieved a thermal dose of CEM43 ≥ 60 in 90 ±5% of tumor area and simultaneously limited the thermal dose in surrounding normal tissue to ≤5% was chosen for further temperature feedback control optimization.
For optimization with a controller, we aimed to achieve a target temperature of 43.5 °C at the probe location rapidly, and maintain it at that temperature for the remainder of treatment. The block diagram for the feedback loop of the PID (proportional-integral-derivative) temperature controller is shown in Figure 2. The temperature computed at the probe location T probe was input into the controller with the reference temperature T ref . The difference between T probe and T ref input into the controller and the controller output were then substituted into the model. The closed-loop transfer function (Tr) for the closed-loop system ( Figure 2) is given by where P is the plant transfer function and C is the controller function. For the controller, we used the general PID temperature control [34].
A proportional controller enables increasing the heating power proportionally to the difference between the computed temperature at the probe location and the target temperature [52]. In general, an integral controller ensures that the control output agrees with the target set-point temperature in steady state [34,52,53]. While 90% of controllers used in industry employ a PI controller, adding derivative control ensures closed-loop stability and enables more rapid rise to target temperature at the target location while minimizing temperature overshoot [34,52,53]. The PID controller is described by [52] (16) In Equation (15), θ(x,y,t) is the difference between the computed temperature at the probe location, T probe (x,y,t), and the target temperature, T ref . The gains K p (1/K), K i (1/(s·K)), K d (s/K) are the proportional, integral and derivative gains of the PID controller, respectively. In Equation (16), u ctrl adjusts the heat generation of nanoparticles, Q np (T), by modulating the heating power applied based on the difference in measured temperature at the probe location and target temperature. Q max is the maximum power that can be generated by the nanoparticles in the model, which depends on the SLP of the nanoparticles and number of nanoparticles. The maximum SLP of 537 W/g Fe (peak field strength of 47kA/m at frequency 150 kHz) for BNF nanoparticles reported by Bordelon, et al. was used to calculate Q max (see Figure S1) [36].
A first order low pass filter was added to the derivative control to limit high-frequency noise amplification. The transfer function for the controller C(s) is given by, where τ d is the filter time constant [33,52] given by where f is the closed-loop frequency and ζ is the damping ratio. The cut-off frequency (1/τ d ) for the low pass filter was chosen to be 4 rad/s [33,52].
The PID controller gains K p , K i , K d were determined by using model control design methods based on Youla parametrization (Q-parametrization) [33,34,52,54]. By Qparametrization theory, for a single-input single-output stable plant transfer function P(s), the family of all stabilizing controllers C(s) can be expressed as, where Q(s) is any stable transfer function. We approximated our model system by a second order plant transfer function P(s) with three parameters, given in Laplace domain by, where g is the static gain for step input, and τ 1 and τ 2 are time constants. The plant transfer function parameters were assessed by an open loop step response. The power was stepped to 30% of maximum power and the temperature response was recorded at the probe location. The static gain, g, is given by the ratio of temperature gain achieved with the step in control input g = Δ T u ctrl (21) The time constant, τ 1 , was calculated as the delay between the maximum rate of change in temperature and the temperature response at the probe location. The time constant, τ 2 , is the difference between the time taken for the temperature response to reach 63% of total temperature gain and the time constant τ 1 [33,52]. The relationship among these parameters and criteria for their selection has been well described by Ebert et al. [33] and visually characterized in Reference [52].
A second-order transfer function, Tr, was chosen, given in the Laplace domain, Using Equations (14), (19) and (22), the controller function C(s) was determined as, The individual gains can then be calculated as, After determining the PID controller gains, the power was modulated using temperature feedback from the probe. The computed temperature at the probe location, T probe , was compared with a reference temperature signal T ref and the difference was input into the PID controller. In the study, the reference temperature signal T ref was chosen to be ramp signal, beginning at 37°C and t = 0s, and achieving the target temperature at t = 30 s. This time point was chosen because it is clinically desirable to reach the target temperature rapidly. The output from the PID controller u ctrl modulated the nanoparticle heat output.

3D models
The healthy tissue was represented by a cube of length 5 cm (see Figures 3(a,b)). Tumors were represented by 3D ellipsoids of a total volume of 150 mm 3 (see Figure 3(c)). Three 3D nanoparticle distribution models were considered. The T1 (uniform) distribution model contained homogeneously distributed nanoparticles throughout the tumor. In the T2 (Gaussian-centered) model, nanoparticles were normally distributed about the tumor center. For the T3 (3 pt-Gaussian), nanoparticles were normally distributed about three centers (indicating points of injection) [14]. The same thermal model that was implemented for the 2D models was used.

Solution method
COMSOL Multiphysics 5.2a (COMSOL Multiphysics, Natick, MA), a commercial finite element solver, was used to solve numerically the governing differential equations (Equation (3)) subject to the boundary conditions described by Equations (4)  was used for solving the partial differential equations and a gen-eralized-α solver [56] was used for the transient solution.

Temperature distributions in 2D models
Temperatures achieved in the 2D models of tumor and healthy tissues after 20 min of constant-power heating and constant perfusion are shown in Figure 4, and plotted graphically in Figure S2 The power required to achieve equivalent thermal dose, i.e., isoeffect power with M1 as the reference was greater in all other models, see Table 2.
When temperature-dependent perfusion was considered, the range of tumor temperatures increased for all nanoparticle distributions compared to constant perfusion, and there was little difference between the two temperature-dependent perfusion models (see Figures S2 (b) and (c)). Conversely, by including temperature-dependent perfusion less power was required to achieve the target thermal (isoeffect) dose, and thermal dose to normal tissues was decreased for the off-center, concentrated nanoparticle distribution (E3), but changed little for other nanoparticle distributions (see Table 2).

Modulated power heating
Single-point thermometry at the simulated tumor-tissue boundary was used to provide data to modulate power during heating. The probe location for further optimization simulation with PID control was selected by simulating heating in each model using a two-step power function (see Equation (10) at each location. The probe location that predicted a thermal dose conforming to our chosen criteria was selected (see Table 3). Probe positions for distribution models M1 and M2 were identical due to symmetry. However, for the off-center concentrated nanoparticle distribution, E3, no location satisfied both criteria simultaneously. In this case, the location P6 (see Figure 2 and The time-dependent power application differed between dispersed and concentrated nanoparticle distributions (see Figure S4). For dispersed nanoparticle configurations, M1 and E1, initial power was at maximum and then decreased to a lower value for the remainder of heating. For distributions having concentrated regions of nanoparticles (i.e., E2, E3 and M2, M3), the applied power initially increased to maximum for a period of time, and then oscillated between a selected maximum and minimum for the duration of heating.
The results of simulated tissue heating are displayed graphically in Figure S5, which also contains tissue damage (a) calculations. High nanoparticle concentrations, located off-center within the tumor and proximal to a normal tissue boundary represent the least ideal treatment scenario ( Figure S4). It is worth noting that in all cases the ideal probe location was farthest from the nanoparticle clusters.

Comparing constant power with PID controlled heating
In Figures 5(a,b) we show temperature distributions and provide tabulated metrics comparing the two heating methods, while Figures S6 show temperature data obtained from each of the models after heating. For all simulated nanoparticle distributions, modulated power heating reduced both overall temperature and temperature heterogeneity. The latter was assessed by comparing T 10 and T 90 , (temperature achieved with at least 10% and 90% of tumor area, respectively) and the dimensionless heterogeneity coefficient HC defined as HC = (T 10 -T 90 )/(T 90 -T core ) [57]. For all models except E2, values of T 10 , T 90 , and HC decreased, signifying a more uniform temperature was achieved with power modulated heating.
The deposited thermal dose was comparable for both heating methods as illustrated in Figure 6(a). On the other hand, modulated power significantly decreased the time needed to heat >50% tumor area regardless of nanoparticle distribution, see Figures 6(b,c).

PID controlled modulated power heating in 3D models
Extending the model to three dimensions (see Figures 3) provides additional insights into power-modulated heating, and better reflects scenarios likely to be encountered in either preclinical or clinical settings. In Figures 7(a-c), we display the time variation of applied power, and in Figure 7(d) the total tumor volume having achieved the target thermal dose. For the uniform nanoparticle distribution, T1, model, the applied heating power peaked and then reduced to a lower power (see Figure 7(a)), similar to its 2D counterpart. For the Gaussiancentered, T2, nanoparticle distribution model, PID control initiated application of maximum power followed by a brief period of damped oscillation stabilizing at a lower power ( Figure  7(b)). Interestingly, the applied heating power oscillated continuously between maximum and minimum for the duration of heating for the 3-point Gaussian distribution, or T3, model (Figure 7(c)). All models exceeded the minimum 90% tumor volume CEM43 ≥ 60 min by the end of heating (Figure 7(d)).
Temperature distributions achieved in the tumor and healthy tissue along the XY, YZ, and ZX planes at the center of the tumor for the three 3D nanoparticle distribution models are shown in Figure 8. Regions with or adjacent to concentrations of nanoparticles achieved higher temperatures than did regions distant from the nanoparticles. The highest maximum intra-tumor temperature (91 °C) was predicted in the Gaussian-centered, or T2, model whereas the T1 (uniform) model achieved the lowest maximum temperature of 52 °C. Qualitatively, these results are similar to those observed with 2D distributions.
In Figure 9, we display predictions of damage to healthy tissue, given by the degree of stasis DS, for the 3D nanoparticle distribution models. It can be seen that the variation of damage in the healthy tissue depends on the proximity of the nanoparticles to the tumor-healthy tissue boundary.

Discussion
The effectiveness of magnetic nanoparticle hyperthermia is largely determined by both nanoparticle distribution in tissues and heat transfer. Heat transfer in living tissues is a complex combination of heat conduction and capillary blood perfusion. Nanoparticle distributions depend upon tumor physical properties; and, significant variation is encountered even with convection-enhanced percutaneous delivery [14,58]. Thus, it is impossible to precisely control nanoparticle distributions for mNPH, but energy deposition can be controlled by adjusting the AMF. Computational modeling can be an efficient tool to explore various optimization strategies. In this study, we sought to explore various power modulation scenarios to identify promising strategies that can compensate for the variable and heterogeneous nanoparticle tissue distribution that accompanies mNPH.
We began with a comparison of six different nanoparticle distributions using 2D computational phantoms to enhance computational efficiency. Local concentrations of nanoparticles, positioned towards the tumor center and away from the tumor-tissue boundary, offer benefit to deposit significant energy to the tumor in a short time (see Figure  4). An idealized uniform nanoparticle distribution tends to generate the lowest local maximum tumor temperature because the nanoparticles are distributed throughout the tumor and heat generated by isolated nanoparticles rapidly transfers throughout the tissue yielding smaller gradients. Consequently, and perhaps not surprising -an idealized uniform nanoparticle distribution will require longer heating times to achieve a specific treatment goal. The benefit, however, is that risks to normal tissue are reduced with such a distribution because high-temperature gradients are less likely. Not surprising, constant power application with a constant perfusion model generated temperature profiles within tumors that closely approximated nanoparticle distributions, with 'hot-spots' manifesting where nanoparticle concentrations were highest (see Figure 4). With constant power, we compared the effects of constant perfusion with two temperature-dependent perfusion variations. Both of the variable-temperature perfusion models predicted higher tissue temperatures than did the constant perfusion model (see Figure S2). As perfusion, and associated heat removal decreased with increasing thermal damage more heat was retained in the tumor, leading to higher local temperatures and thermal dose.
A uniform temperature distribution in the tumor is often considered desirable; however, achieving a minimum effective temperature in 90% of tumor with sufficient control to preserve surrounding normal tissues is more realistic and is demonstrably effective [1,59]. Achieving a minimum clinically relevant thermal dose (CEM43 ≥ 60min) in a larger volume of tumor is considered a better measure of treatment outcome [1,10,51]. Equally important is minimizing the damage to surrounding normal tissues. The latter is often ignored, particularly in computational studies of mNPH. A reference isoeffect heating power, Q iso , was chosen as the minimum power required to realize a thermal dose in a maximum tumor area. The constant perfusion model underestimates tumor temperature(s), and consequently overestimates the power required to achieve a target thermal dose ( Table 2) compared to temperature-dependent perfusion for different nanoparticle distributions We combined non-linear temperature-dependent perfusion with power modulation of (nonlinear) nanoparticle heating properties (see Figure S1) to reduce the impact of heterogeneous nanoparticle distributions in tissues. Attaluri et al. [14] previously showed with a 2D circular tumor model that modulating power using temperature feedback from a sensor at the tumortissue boundary can be used to achieve a minimal effective thermal dose (CEM43 T90). A non-spherical tumor shape and asymmetric nanoparticle distribution require more precise placement of temperature sensor. The criterion for choosing a probe location was to ensure that the thermal dose of CEM43 ≥ 60min was achieved in 90% of tumor area and in less than 5% of healthy tissue. It was possible to identify an ideal probe location for all nanoparticle distribution models, except E3 (Table 3). The concentration of nanoparticles close to the tumor-tissue boundary leads a predicted ineffective treatment and undesirable damage to healthy tissue.
We implemented a PID temperature controller (Figure 2(b)) to regulate power using temperature feedback from a probe located at the tumor-tissue boundary of our phantom models. The higher proportional (and consequently higher derivative) gains generated power oscillation or 'cycling' for all models having concentrated nanoparticle distributions (see Figure S3). For 'uniform' distributions M1 and E1, the applied power reached a maximum quickly then reduced to a lower power for the remainder of the treatment. For all other distributions, the applied power ramped quickly to maximum, and oscillated between maximum and minimum thereafter. For M3, E2, and E3, the peak of this oscillation was the maximum power possible Q max , while for M2 it was ~75% of Q max . The distance of the probe location from the nanoparticle heat sources for these distributions produced a delay between the control action and the temperature output. This required the controller to actively control the temperature output producing oscillatory input power. Arora et al. [53] reported similar oscillating control inputs in with their ultrasound treatments. Further tuning of the PID controller can improve its performance for more stable power application. Overall, our models predicted that PID controlled power modulation can improve tumor heating, regardless of nanoparticle distribution (see Figures S4 and S5), while simultaneously minimizing healthy tissue damage ( Figure S4(b)).
A direct comparison of constant power and PID modulated power heating highlights the advantages of the latter. Overall, lower and more homogeneous temperatures were achieved for all models with power-modulated heating compared to constant power heating ( Figure  5(a,b)). As expected, the uniform model (M1) showed no significant difference (~1%) between constant and power-modulated heating. The maximum benefit of modulated heating was observed for off-centered nanoparticle distribution (E3), where the maximum tumor temperature decreased by 16%. This can be potentially clinically significant, as previous studies have shown [10,60] that nanoparticle leakage and distribution occurs along the needle track, resulting in nanoparticle deposition close to tumor-healthy tissue boundary. Thus, even in cases where undesirable nanoparticle distributions occur, power modulation provides a better alternative over constant power heating to deliver effective thermal dose and in minimizing healthy tissue damage. With power-modulated heating, the time required to achieve therapeutic thermal dose in a large area of tumor (> 50%) is both considerably (~9x) shorter and less dependent on nanoparticle distribution than with constant power heating. Pulsed-power application may be the preferred method to perform mNPH because the time to reach therapeutic temperatures is short increasing the potential for effective therapy [61,62].
Following the optimization studies in 2D, we implemented the modulated power heating approach in 3D models. Results obtained demonstrate significant similarities to the 2D results. The PID controller generated a damped oscillating power during the first ~2min of heating for the Gaussian-centered nanoparticle concentration (i.e., T2, see Figure 3). This contrasts with trends observed in 2D results (M3, Figure S3). It is likely these differences occur because of different heat transfer between 2D and 3D models. An additional consideration, not assessed in the present work is the sensitivity of 3D models to variations of temperature probe placement. To reduce computational expense, 3D probe placement was chosen from our analysis of 2D scenarios which indicated the ideal placement of temperature probes was distal to nanoparticle localization. Our analysis of 2D models revealed greater sensitivity for inhomogeneous (i.e., concentrated, off-center) nanoparticle distributions. We expect that for 3D geometries such sensitivity may be heightened, warranting a further analysis for likely preclinical and clinical scenarios that have associated significant probe localization uncertainties. Nevertheless, we identify that a critical component to reduce such uncertainty is precise knowledge of nanoparticle localization.
Knowledge of the precise nanoparticle distribution, e.g., from imaging, would facilitate patient-specific optimized treatment planning by enabling computational modeling that incorporates active power modulation, and could be used to reduce uncertainties of temperature probe placement. Stated another way, the outcome of therapy can depend less on achieving a specific nanoparticle distribution within the tissue. A limitation of the current study was our use of a single temperature probe. Additional probes can potentially enhance power modulation to further optimize the simulated temperature distribution throughout the tumor; however, the invasive nature of probe-based thermometry may limit the number of probes used in clinical applications. Another limitation of our study was the exclusion of eddy current generation in tissues leading to off-target heating. A patient torso of ~30cm diameter exposed to AMF will display significant temperature rises from induced eddy currents, absent nanoparticles and depending upon the AMF conditions [63,64]. Additional study is warranted to explore further optimization by including these.

Conclusions
To address one of the principal challenges of magnetic nanoparticle hyperthermia, we proposed modulated heating of magnetic nanoparticles using temperature feedback from the tumor-healthy tissue boundary, in situ, as a viable means to compensate for variable nanoparticle distributions in tissues. With both 2D and 3D models, we demonstrate it is possible to achieve effective tumor treatment with modulated power heating that realizes a more rapid rise to therapeutic temperature in a larger volume of tumor, a more homogeneous tumor temperature distribution, and with overall improved control of temperature in normal tissues. Further experimental investigation is warranted to develop this approach.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.         Healthy tissue damage is higher for regions closer to nanoparticle distributions. Degree of thermal damage (DS) achieved inside the healthy tissue, (a) XY plane, (b) YZ plane, (c) ZX plane, for the three 3D nanoparticle distribution models T1, T2 and T3, after 20 min of heating by power modulation with PID control based on temperature feedback from probe at tumor-tissue boundary with the modified Arrhenius perfusion model. Thermophysical properties used in simulations for the tumor, surrounding healthy tissue, and blood.
Percent area of tumor and healthy tissue having received a thermal dose CEM43 ≥60 min after 20 min of heating with two-step power modulation (Equation (10) based on temperature feedback from eight probe locations (P1-P8, see Figure 2).  Table 4.
Open loop response and PID control parameters.