Heating Protocol Design Affected by Nanoparticle Redistribution and Thermal Damage Model in Magnetic Nanoparticle Hyperthermia for Cancer Treatment

Magnetic nanoparticle hyperthermia has attracted a lot of attention due to its ability to confine heating to targeted tumors [1–10]. In magnetic nanoparticle hyperthermia, injected magnetic nanoparticles in tumors can generate heat when subjected to an alternating magnetic field, and the volumetric heat generation rate in tumors is proportional to the local nanoparticle concentration. Thereby, the required thermal dosage for tumor destruction greatly depends on nanoparticle distribution in tumors. In recent years, with the advancements in computational resources, improved strategies such as using image scans to generate realistic tumor models have been proposed to develop patient- or tumor-specific approaches that deliver effective and safe thermal therapy [11–13]. Micro-CT scans have been demonstrated to visualize and quantify nanoparticle distribution in opaque tumors [11,14–20]. Although micro-CT image-assisting simulation mimicking realistic tumor geometry has improved accuracy in preclinical treatment design, those previous theoretical simulations of magnetic nanoparticle hyperthermia did not consider the possibility of nanoparticle migration during heating.

Micro-CT scans of resected PC3 (prostate cancer) tumors have been used by our group to visualize distribution of magnetic nanoparticles dispersed in tumors. It is shown that most of the nanoparticles injected at the tumor center are located in the vicinity of the injection site, and the brightness in the micro-CT images is an index of local nanoparticle concentration. Via a previously calibrated relationship between micro-CT Hounsfield unit values and local volumetric heat generation rate $QMNH′′′$⁠, one can simulate the temperature elevations in the tumors exposed to an alternating magnetic field, leading to heating protocol design to ensure thermal damage to the entire tumor [18]. However, a recent study [20] shows possible nanoparticle redistribution in PC3 tumors due to local heating. In that study, micro-CT scans were performed on resected tumors in a control group without heating and in an experimental group with 25 min local heating. The nanoparticle distribution volume in the tumors resected after local heating is 40% larger than that in the tumors without local heating. It suggests that nanoparticles may spread from the originally concentrated central region to tumor periphery, leading to a bigger nanoparticle distribution volume in the tumors in the heating group. Although the previous experiments did not directly shed light on possible transport mechanisms responsible for the redistribution of nanoparticles, we speculate that thermal damage to the local tumor cells may lead to an increase in tumor porosity, thus, increasing the diffusion coefficient of nanoparticles in interstitial fluid [21,22].

One limitation of the previous theoretical simulations of tumor temperature fields is that the possibility of nanoparticle migration during heating was not included in the model. If magnetic nanoparticles migrate during local heating, it would alter the volumetric heat generation rate distribution in tumors. Depending on how the nanoparticle distribution affects the temperature field, the dynamic behavior may necessitate different heating protocol designs. It is unknown to us the extent to which the dynamic nanoparticle redistribution affects the local temperature field.

It is well known that local blood perfusion acts as a heat sink carrying heat away from targeted tumor during hyperthermia treatment. Temperature elevations in tumors would be influenced significantly by the blood perfusion rate in tumors. In addition, blood perfusion in tumors may also change during heating due to both thermoregulation in the body and thermal damage to the vasculature in tumors. One of the limitations of most previous theoretical simulations used blood perfusion values available in online resources. Unfortunately, blood perfusion rates may be affected by tumor types, growing stages, and whether the tumor is embedded inside an organ or implanted on superficial skin layer. In PC3 tumors alone, we found that the local blood perfusion rate in most of the previous theoretical simulation studies cited data from Lang et al. [23]. It would be an improvement to acquire the blood perfusion rate based on the experimental data of temperatures on PC3 tumors implanted on flanks of the same nude mice. Further, dynamic responses of the vasculature in tumors to heating are unknown. Typically, the extent of increase in blood perfusion rate in tumors during heating is minor, due to the lack of smooth muscle in tumor vasculature that prevents compensatory vasodilation. During heating, blood vessels in the tumor also suffer permanent thermal damage, leading to a decrease in the average blood perfusion rate. The assumption that the blood perfusion rate decreases to 50% of its baseline value after heating in Lang et al. [23] may not be an accurate description of the dynamic responses. Coupling the local blood perfusion rate directly with the local thermal damage would improve the accuracy of theoretical simulations.

Typically, theoretical simulation of a transient temperature field in tumors is coupled with a thermal damage model to quantitatively assess thermal damage to tumor cells, leading to a heating protocol design. Among many cell survival models that have been developed and tested, the original Arrhenius integral developed by Henriques and Moritz [24] has been used extensively in hyperthermia treatment to quantify the extent of thermal damage in tumor regions. The original Arrhenius model used a linear relationship to curve fit the experimental data of cell death rate affected by heating time. Unfortunately, the linear relationship did not capture the shoulder region representing no cell death initially after heating starts, and the curve fitting was poor to the experimental data. Recently, the original Arrhenius model was modified by Dr. Pearce to include a time delay accounting for the initial cellular activities before cell apoptosis can occur [25,26]. With the time delay, fitting of the rest of the experimental data using a linear relationship is improved significantly. This time delay is a function of local temperature and is also tumor cell specific. In principle, the modified thermal damage model would change the heating time when used in heating protocol design. To the best of our knowledge, the newly developed thermal damage model has only been implemented in theoretical simulation in a few studies to improve the accuracy of heating protocol design [26].

In this study, a theoretical simulation was performed to evaluate the extent to which the nanoparticle redistribution affects the temperature elevation and heating time required to cause permanent thermal damage to PC3 tumors. The previously obtained micro-CT scans of PC3 tumors were used to generate physical models of the tumors, which later were imported into comsol software for temperature elevation simulation [17,18,20]. The baseline blood perfusion rates in the mouse body and PC3 tumor were first extracted based on the experimental data of mouse surface temperatures captured by an infrared camera. The control group consists of four PC3 tumors with nanoparticles distribution resected after an intratumoral injection of the ferrofluid, while the experimental group consists of another four PC3 tumors with a ferrofluid injection resected after 25 min of local heating. Further, we assessed the heating time needed to cause essentially 100% permanent thermal damage (Arrhenius integral Ω = 4) to the entire tumor, based on the obtained nanoparticle distribution in each tumor group, as the first step to evaluate how significant the nanoparticle distribution affects the treatment design. We implemented both the traditional Arrhenius model and the Pearce thermal damage model [25,26] in designing heating protocols for magnetic nanoparticle hyperthermia. It is expected that different heating durations are required to damage the tumors due to different nanoparticle spreading patterns in both tumor groups, and the modified thermal damage model predicts a shorter heating time than the traditional Arrhenius damage model.

Figure 1 shows maximum intensity projection (MIP) images of tumors obtained from our previous study [20]. 0.1 cc of a commercially available ferrofluid was injected directly to the center of PC3 tumors implanted on mice. The tumor was later resected and scanned by micro-CT. The top panel is for a tumor in the control group without heating, and the bottom panel is for a tumor in the experimental group, resected after local heating. Both tumor groups have similar size (control: 1261±177 mm³ versus experimental: 1304±317 mm³). The white clouds in the MIP images show the locations of nanoparticles in the tumor. In the resected tumors after local heating, the nanoparticles spread more from the injection site at the tumor center than that in the tumors of the control group [20].

The micro-CT scan of the tumor was then imported to generate a tumor physical model of a homogeneous material. The tumor was then attached to a previously generated mouse body, to mimic the situation in the original animal experiments [18,20], as shown in Fig. 2. The volumetric heat generation rate distribution $QMNH″′$ in each tumor was obtained from the study of Gu et al. [20]. Both tumor groups have similar energy deposition rate, approximately 0.37 W in the tumor under an alternating magnetic field of 5 kA/m at a frequency of 200 kHz.

where T is the temperature, t is the time, k is the thermal conductivity, ρ is the density, c is the specific heat, ω is the local blood perfusion rate, T_blood is the arterial blood temperature prescribed as 37 °C, $Qm′′′$ is the volumetric heat generation rate by metabolism, and $QMNH″′$ is the volumetric heat generation rate by nanoparticles. The subscripts tissue and tumor represent mouse tissue and PC3 tumor, respectively. It is assumed that the thermal properties remained constant and isotropic within each domain, as listed in Table 1. Similar to our previous study [18], the boundary conditions of the simulation are illustrated in Fig. 2. Specifically, the bottom surface of the mouse is a prescribed 37 °C, due to the use of a heating pad set at 37 °C to keep the mouse warm in the previous experiments [20]. The rest of the mouse and tumor surfaces are subject to a convection environment (h = 5.13 W/m² K and T_air = 25 °C). The initial condition is prescribed as the steady-state temperature field before the heating.

The normal blood perfusion rates in the tumor (ω_tumor,0) and the mouse body (ω_tissue,0) will be extracted via measured surface temperatures before the heating using an infrared camera. The volumetric heat generation rate by metabolism before the heating is assumed coupled with the local blood perfusion rate. The details of the blood perfusion rate determination are described in Sec. 2.2. In Eq. (2), Ω is the Arrhenius integral assessing thermal damage, described later. In this study, we assume that local heating damages local vasculature in tumors, therefore, both the blood perfusion rate and the volumetric heat generation rate by metabolism decrease with the progression of thermal damage.

where A is the frequency factor, E_a is the activation energy, i.e., a barrier that reactants must overcome for irreversible damage reaction, R_u is the universal gas constant, T_tumor is the tumor temperature in Kelvin, t is the time, and τ is the heating duration. The thermal damage threshold of Ω = 4 is used as an indicator of treatment end point for identifying the heating duration needed when the tumor is exposed to the specific alternating magnetic field. The values of the Arrhenius parameters shown in Table 2 are obtained from Dr. Pearce's papers [25,26].

where the coefficients in Eq. (4) are specific to the PC3 cells. Note that the tumor temperature in Eq. (4) is in Celsius.

Different tumor locations would be damaged at different paces. A high temperature in cells would lead to accelerated cell death if the exposure time is the same. Inside a tissue with nonuniform temperature elevations, the tumor cell locations with higher temperatures would reach the threshold of irreversible thermal damage first. Thermal damage then spreads to cells with lower temperatures, since it takes longer heating exposure to cause the same degree of damage. The heating treatment design is to identify a heating time to ensure all the tumor region reaches the threshold of irreversible damage. Therefore, the tumor location with the lowest temperature would require the longest time to be damage. Implementing Eqs. (3) and (4) in a transient simulation is difficult since the time delay depends on local tumor temperature that keeps changing with time. In this study, we only need to assess the thermal damage at the tumor location with the lowest temperature. We first simulate the steady temperature field to identify the tumor location with the lowest value. Then, we perform the simulation of the transient heat transfer process to determine the temperature rising curve versus time at the tumor location. The temperature versus time relationship described in Eq. (4) is also plotted on the same chart. It yields an intersection between those two curves, which gives the time delay τ_d for this location. This allows the assessment of thermal damage using Eq. (3) to determine heating duration for damaging the entire tumor. Since the temperature and thermal damage are coupled, several iterations are needed to determine the time delay with less than 1% variation in any two sequential values of τ_d in iterations.

In a recent study by our group (unpublished), an infrared camera was used to record the surface temperatures of both a mouse body and a PC3 tumor implanted on the left flank of the mouse. The six mice used in that study were the same kind of nude mice. Experimental procedures of tumor implantation were similar to that described previously [20]. Briefly, PC3 xenograft tumors were implanted in 6 Balb/c Nu/Nu male mice (The Jackson Lab, Bar Harbor, ME), with one tumor in each mouse. PC3 cell solution containing 5 × 10⁶ cells was then injected into the left flank of the mouse using a 27 gage needle (Tuberculin Syringe w/Needle by BD, Fischer Scientific, Springfield, NJ). The tumor size was estimated using a Vernier caliper, and the growth of the tumor was monitored three times a week until it reached 10 mm in the transverse diameter, which usually took 4–6 weeks. The mouse was brought to the lab and anesthetized via sodium pentobarbital solution injection (40 mg/kg, intra-peritoneal). The mouse was then placed on a heating pad with circulating warm water of 37 °C to maintain a normal body temperature of the mouse. An infrared camera (TIS45, Fluke Corporation, Everett, WA) was used to take the infrared image of the surface temperature field of the mouse body. In each image, a circular tumor region was identified and the average temperature of the tumor surface was determined. Similarly, a circular region of the same size was selected on the mouse body, and its average temperature was calculated. Based on the six infrared images of the six mice, the temperature values were further averaged.

Equations (1) and (2) can be modified slightly to simulate the steady-state temperature field of the mouse body and implanted tumor before the heating. For simplicity, we only simulate the steady-state temperature field of one of the generated tumor models attached to a mouse body. In Eq. 2, $QMNH′′′$ is set as zero due to no nanofluid injection in the tumor, and the Arrhenius integral Ω was also zero before any heating. $Qm,0′′′$ in Eqs. (1) and (2) is assumed proportional to the local blood perfusion rate. In the simulation, the adjustable parameters are the baseline blood perfusion rates in the tumor (ω_tumor,0) and the mouse body (ω_tissue,0). Their values are adjusted until the simulated average temperatures of the tumor surface and mouse body surface match that measured in the experiments using the infrared camera within ±0.01 °C. The extracted values of ω_tumor,0 and ω_tissue,0 are used in the transient temperature field simulations described in Sec. 2.1.

The numerical model has been solved using a commercially available finite element method based solver, comsolmultiphysics 5.2 (Comsol, Inc., Stockholm, Sweden). Extremely fine tetrahedral mesh elements were implemented to discretize the physical domain and quadratic Lagrangian elements were used to discretize the solution space. The numerical convergence has been found to take place with a prespecified relative tolerance of 0.001. Thermal damage assessment was calculated via converting Eq. (3) to a partial differential equation for solving for Ω. A direct linear solver, MUMPS (Multifrontal Massively Parallel sparse direct Solver), with a default pre-ordering algorithm was used for solving the temperature fields. The backward difference formula was used for time marching, and care has been taken to store the data at each step by assigning strict steps for the solvers. The Newton–Raphson algorithm in the direct solver was used for iteration. The temporal resolution was selected as 0.01 s. Lowering the temporal resolution by half resulted in a less than 0.01% change in the average temperature distribution. The total number of mesh elements in all the models in the study varied from 400,000 to 800,000. A fourfold increase in the mesh elements resulted in a difference of less than 0.1 °C in the minimal and maximal temperatures in the tumors.

The experimental measurements of the tumor and mouse body surface temperatures illustrated that the tumor surface temperatures were much lower than that of the mouse body temperatures. The circular region selected for the tumor and mouse body was 12 mm in diameter. The average values of the six tumor surface temperatures and the mouse body surface temperatures were determined as 33.38 ± 0.51 °C and 35.70 ± 0.64 °C, respectively. Simulations of the steady-state temperature fields in both the tumor and the mouse body were performed. Both the baseline local blood perfusion rates and the metabolic heat generation rates were extracted to agree with the average values of the surface temperatures with a deviation less than ±0.01 °C. The extracted values are given in Table 1.

Figure 3 gives the central slice of the three-dimensional steady-state temperature field inside the tumors and the surrounding mouse tissue. Temperature elevations are largely confined to the tumor, however, collateral thermal damage of the surrounding mouse body may occur. The top panel is for a tumor in the control group, where the maximal temperature reaches 79 °C and the minimal temperature of approximately 52 °C occurs at the tumor mouse body interface. As shown in the bottom panel for a tumor in the experimental group, the overall temperature field in this tumor in the experimental group varies from 48 °C to 78 °C, lower than that of the tumor in the control group. This may be due to the fact that in the experimental group, the nanoparticles are more spreading from the tumor center to the tumor periphery (average temperature in the control group: 70.6 ± 4.7 °C versus experimental group: 69.2 ± 2.6 °C). In the control group, the minimal temperature is also significantly higher than that in the experimental group (control group: 53.3 ± 3.5 °C versus experimental group: 50.0 ± 2.5 °C).

It would have been ideal to incorporate dynamic nanoparticle migration into the heat transfer model. In this study, as the first step, we simulated the transient temperature elevations based on the obtained two sets of nanoparticle distribution from the two tumor groups. Assuming that the nanoparticle distribution does not change during the simulation in each tumor group, the transient temperature distribution was simulated and thermal damage was assessed to design the heating protocol. Two sets of heating protocols should be determined from the two tumor groups. In principle, if one implements the dynamic responses of nanoparticle migration in the heat transfer simulation, the actual heating protocol should fall between the two sets of heating protocol designs in this study.

Figure 4 illustrates the transient temperature curves at two locations, one is the maximal temperature location and the other is the minimal temperature location in a tumor. For the tumor in the control group and the tumor in the experimental group shown in Fig. 4, it takes more than 40 min to establish its steady-state.

The traditional Arrhenius integral or the time-delay thermal damage model is then coupled with the transient temperature elevations to determine the distribution of the Arrhenius integral Ω. The top panel in Fig. 5 shows the damages on one cross-sectional plane of a tumor in the control group after 5, 10, 15, 20, and 25 min of heating, respectively, using the traditional Arrhenius model. The dark area represents the damaged tumor region when Ω ≥ 4. It can be clearly seen that as the heating time increases and heat conducts throughout the tissue, more and more damage in the tumor is observed. Our results have shown that it takes approximately 25 min to completely damage this tumor in the control group. The damage propagation in another tumor in the experimental group is given in the bottom panels in Fig. 5, again based on the traditional Arrhenius model. It takes much longer time (40 min) to damage the entire tumor, due to the lower temperature at the interface between the tumor and the mouse body, where the lowest temperature in the tumor occurs.

For each tumor, simulated steady-state temperature field showed that the minimal temperature elevation within the entire tumor usually occurred at the interface between the tumor and the mouse body. As described in the Methods section, the required heating time is determined by the temperature elevation history of this location. As shown in Eq. (4), the time delay should be determined first to implement the modified thermal damage model. Figure 6 illustrates how the time delay for this location in one tumor is determined. As shown in Fig. 6, at the early stage of the heating, the tumor temperature is low, resulting in a time delay τ_d that is bigger than the actual time t. As t increases, tumor temperature rises, leading to a decrease in the time delay. The temperature rising curve at the tumor location and the decaying curve of Eq. (4) would eventually cross paths, and the intersection gives the time instance when the actual time t and the time delay τ_d are exactly the same. The intersection of the temperature rising curve and that of Eq. (4) in Fig. 6 then determine the time delay for using the Pearce thermal damage model. τ_d is calculated as 540 s in Fig. 6. For the tumors in the control group, the time delays are calculated similarly and it varies from 396 to 632 s (498 ± 108 s). The values of the time delay in the experimental group are 16% longer by average, varying between 499 and 640 s (576 ± 59 s).

The effect of the Pearce model on the temperature elevations is shown in Table 3. Since the temperature and thermal damages are coupled, the influence of the Pearce model on the temperature elevations should be more evident in the earlier heating stage with larger local blood perfusion rate. When one examines the temperature elevations at the lowest temperature locations, the modified Pearce thermal damage model leads to almost the same temperature elevations in the later stage of the heating with deviations less than 0.2 °C, shown in Table 3.

Figure 7 illustrates the accumulation of thermal damage versus heating time at the tumor location with the lowest temperature, using either the traditional Arrhenius model or the Pearce thermal damage model. In the traditional Arrhenius model, the value of Ω starts from zero, then increases slowly at the beginning of the heating, and finally, accelerates toward the later duration of the heating due to large temperature elevations. The Pearce model assumes an initial delay time (τ_d = 540 s) within which no accumulation of thermal damage occurs. However, due to the thermal damage parameters that promote quicker thermal damage in the Pearce model, the damage curve accelerates after the initial time delay duration, and eventually, it requires a shorter heating time to achieve the same thermal damage threshold of Ω = 4 (Pearce model: 1370 s versus traditional Arrhenius model: 1645 s) for the same tumor.

Simulations are performed for all the tumors in both groups. As shown in the solid columns in Fig. 8, the required average heating time for the control group is 27 min using the traditional Arrhenius damage model. However, the more spreading of nanoparticles to tumor periphery in the experimental group results in a much longer average heating time of 39 min, which is 44% longer than that in the control group, to induce permanent thermal damage to the entire tumor.

The patterned columns in Fig. 8 illustrate the required heating time assessed using the Pearce thermal damage model. The Pearce thermal damage model leads to a designed heating time of 19 and 29 min for the control and experimental groups, respectively, at least 24% shorter than the predicted values by the traditional Arrhenius model. Since the heating time is determined based on an assumption of two nanoparticle deposition patterns, we speculate that the heating time needed when considering nanoparticle migration during heating is likely to vary between 19 and 29 min, based on the more accurate Pearce thermal damage model. The variation of the designed heating time in each group is probably due to difference of tumor shapes and relative locations of nanoparticles inside the tumors.

This study is focused on determining the extent to which the designed heating protocol in magnetic nanoparticle hyperthermia is affected by various factors such as the shape and size of the tumors, the thermal damage assessment model, as well as dynamic responses of parameters during heating. The required heating time is primarily determined by the tumor location with the lowest temperature value in the tumor. Therefore, factors affecting the value of the lowest temperature at that location can exert significant influence on the heating protocol design. From the point of view of traditional heat conduction in a solid, the shape and size of the tumor, the boundary conditions, the volumetric heat generation rate distribution due to injected magnetic nanoparticles, and dynamic responses of parameters during heating can all play important roles in the design of the heating protocol.

One improvement of the study is to use an experimentally extracted baseline blood perfusion rates in both the mouse body and the implanted tumor. As shown in Table 1, the extracted blood perfusion rate in the tumors is smaller than that in the mouse body, contributing to the lower temperature values recorded by an infrared camera. The extracted baseline blood perfusion rates are in the same order of magnitude as that reported by Lang et al. [23]. It is possible that the difference is due to different growth stages of the tumors and mice. In addition, the local blood perfusion rate in tumors was assumed to decrease to 50% of its baseline value for a temperature range higher than 47 °C in Lang et al. [23]. This may not be realistic since tumor vasculature may be completely damaged with tissue damage. In this study, we propose that the local blood perfusion rate is coupled with the thermal damage. With those improvements, we are confident that the simulation results are more accurate than the previous simulation data based on blood perfusion rates provided in Lang et al. [23]. The extracted blood perfusion rates from this study can also be used by other researchers in future theoretical simulation when they use similar mice and types of tumors.

Nanoparticles are more concentrated in the vicinity of the injection site in the control group [17]. As a result, there exists a much higher volumetric heat generation rate in the tumor central region which would result in a high temperature gradient in the region by conduction. Therefore, it is not a surprise to observe a higher maximal temperature in the control group than that in the experimental group. On the other hand, the value of the minimal temperature may be affected by various factors. As shown in the model, a majority of the tumor surface is exposed to the cold ambient environment, while a small portion loses heat by conduction to the mouse body. Via examining the heat loss rates from those two boundaries, one finds that the wider spreading of nanoparticles in the experimental group increases temperatures at the tumor surface exposed to air, thus, leading to a bigger convection heat loss rate to the ambient environment than that in the control group. Since both groups have similar energy deposition rates, the extra loss of heat to the ambient environment in the experimental group suggests less heat conduction to the mouse body at the tumor–mouse interface. This leads to a lower temperature at the tumor/mouse body interface in the experimental group than that in the control group. Our simulation results agree well with a previous temperature elevation simulation [18], when a more spreading of nanoparticles caused by different nanofluid infusion strategies also resulted in a lower temperature at the tumor–mouse interface.

The results have suggested that the shape and size of the tumor are very important factors in determining the temperature elevations when the tumor is subjected to an alternating magnetic field. The large standard deviation in the designed heating time in each group implies that the shape and size of an actual tumor may alter the location of the minimal temperature in the tumor and the steady-state temperature value at that location. Truly individualized heating protocol designs have to consider the shape and size of actually tumors, generated by image scans of the tumor. This information with the obtained volumetric heat generation rate distribution due to nanoparticle deposition would greatly improve the accuracy of heating protocol design specific for the tumor.

Micro-CT imaging scans have allowed the quantification of volumetric heat generation rate in tumors based on values of grayscale or Hounsfield unit of micro-CT scans. If one assumes that the accuracy of previously determined relationship between the Hounsfield unit value and the volumetric heat generation rate is acceptable, in theory, the obtained realistic three-dimensional $QMNH″′$ should improve the accuracy of the theoretical simulation for design. In this study, steady-state temperature simulation demonstrates that the maximal and minimal temperatures in the tumors of the control group can be several degrees Celsius higher than that of the experimental group, when the nanoparticles are spreading more over the tumor. Later, a transient heat transfer simulation is carried out based on the assumption that the nanoparticle distribution in each group is unchanged. This is certainly not an ideal approach. In fact, a better way to formulate the problem is to couple a mass transport model of nanoparticle diffusion with the heat transfer model to simulate real-time migration of nanoparticles, and interaction between the dynamic responses of nanoparticles and resulted heat transfer process. Unfortunately, it is not clear to us which transport mechanisms are dominant in determining nanoparticle migration during heating. We speculate that porosity change of interstitial fluid space may be one of the factors facilitating the nanoparticles to diffuse away from their original locations. However, more experiments using better experimental tools are warranted to illustrate the roles played by various transport mechanisms. Although our approach has limitations, it demonstrates the importance of including nanoparticle migration in the simulation. The heating protocol design for maximally damaging the entire tumor shows 44% (traditional thermal damage model) or 54% (Pearce thermal damage model) longer heating duration needed in the experimental group when the nanoparticles are more spreading, a result further motivates one to improve the model. Again, like any theoretical simulation, future experimental studies with imaging tools to provide real-time monitoring of nanoparticle distribution, local blood perfusion rates, and temperature elevations are needed to validate the model. After the heating, histological analyses of tumors should also be used to confirm nanoparticle migration and permanent tumor cell damage. One of the experimental studies by our group implemented a heating time of 25 min to the same kind of PC3 tumors, and found all tumors in the group disappeared three days after the heating and the tumors did not grow back within 8 weeks postheating. If one implements the dynamic responses of nanoparticle migration in the heat transfer simulation, the accurate required heating time is probably between 19 min and 29 min when using the more accurate Pearce's model. The theoretically predicted range of heating time agrees well with our previous experimental observations.

Cellular behavior in tissue during heating is complicated. When the temperature elevation from its baseline is relatively small, recoverable damages, such as intracellular edema, thermal inactivation of specific enzymes, as well as cellular membrane rupture may occur and it is possible that cells may tolerate and survive the heating [28]. However, high temperature elevations often result in necrosis (trauma induced cell death) and apoptosis (programmed cell death) when cells are unlikely to recover, since repair mechanisms or their mediators (deoxyribonucleic acid and ribonucleic acid enzymes) are affected due to the imposed thermal stresses [29]. To complicate the process further, after a heating treatment, the release of heat shock proteins can also be regarded as an indication of protective responses of cells from further stresses [30]. Mathematical models simulating cellular responses to heating have been proposed to aid the design of heating protocols, often incorporated with theoretical simulation of temperature elevations in tumors [24,31–33]. The traditional Arrhenius integral, based on the first-order chemical kinetics, has been used to fit experimental measurements of cell survival rate to determine the frequency factor and the activation energy in the integral. The Pearce model is to address the relatively poor fit of the first-order chemical kinetics to the experimental data at the early stage of heating. As shown in the Pearce model [25], the fitting can be significantly improved by including a time delay that is dependent on local temperature and cell type. Using the parameters provided for PC3 tumors, the initial time delay is calculated as approximately 8.5 min by average in the thermal damage assessment in this study. However, the initial time delay does not prolong the required heating duration. We compare the results of the heating duration using either the traditional Arrhenius model or the Pearce model. The actual heating time using the Pearce model is at least 24% shorter than that predicted by the traditional Arrhenius model, despite the initial delay. This is understandable since the accumulation of thermal damage grows exponentially with tumor temperatures. Our predictions agree with that in experimental results presented in Dr. Pearce's paper. Figure 1 in Pearce [25] shows the experimental data of PC3 tumor survival rates at 44 °C after specific heating durations. One can see a longer heating time needed if one uses the Pearce model when the survival rate is larger than 30% (or Ω < 1.2), however, the heating time is shorter than that predicted by the traditional thermal damage model, when the threshold survival rate is set smaller than 30% (or Ω ≥ 1.2). This agrees well with this study when Ω = 4 is selected as the threshold for irreversible thermal damage. One limitation of the Pearce model is the upper limit of the temperature elevation when using Eq. (4). More experimental studies are needed to examine the extent to which the modified model by Pearce influences the heating protocol design when tumor temperatures are elevated above 54.5 °C.

There are several limitations with computational model developed in this study. Although we found in general that the temperature elevations in tumors are affected by nanoparticle spreading, due to the small sample size in each group, the results may not show statistical difference between the two groups. In the future, tumor groups with a large sample size should be used to test whether statistical significance would be achieved when the sample size is getting bigger.

We use a previously calibrated relationship between the micro-CT Hounsfield unit and the volumetric heat generation rate when the magnetic nanoparticles are subject to an alternating magnetic field of 5 kA/m. The calibration was performed based on the assumption of no nanoparticle agglomeration and independence of the volumetric heat generation rate on the medium (tissue, tumor, or phantom gel). However, as shown in recent studies [34–37], the linear relationship may not be valid, especially in the high nanoparticle concentration range. In addition, the effect of possible nonuniform magnetic field within the volume of the tumor is unknown [38]. Development of a uniform alternating magnetic field to cover a large treatment volume is important for reliable thermal treatment design.

As shown in the temperature contours in the tumor and mouse body, significant heat loss to the mouse body occurs via heat conduction, leading to significant temperature elevations in the healthy tissue at the interface. It is well known that the local blood perfusion rate in healthy tissue normally increases with elevated temperatures at the early stage of heating when thermal damage is not severe [39,40]. As shown in a previous study [39], the local blood perfusion rate later decreases due to progressing vascular damage in the tissue region. In theory, the dynamic response of local blood perfusion may affect the temperature elevation and treatment design. When combining the initial increase and later decrease in local blood perfusion, we speculate that error due to our assumption of an unchanged local blood perfusion rate in the mouse body during the heating might be minor. However, it is a limitation of this study.

Although this study aims at improving simulation accuracy, real-time monitoring of tumor blood perfusion rate during would have been helpful to understand the dynamic responses. Implementation of blood perfusion rate monitoring using imaging techniques such as magnetic resonance imaging is not very easy due to poor spatial and temporal resolutions, as well as possible interference with the alternating magnetic field during heating and image scans. Simultaneous measurements of temperatures at various tumor locations could also be used to validate theoretical simulation results. Unfortunately, it is often not feasible due to low resolution of temperature measurement technique, small sizes of tumors, and physical intrusion of temperature sensors that may alter the original nanoparticle distribution in tumors. Future experimental studies are needed to continue to provide precise thermal/physiological properties and their dynamic responses, and simultaneous measurements of temperature elevations and other parameters to validate theoretical models.

In summary, in this study, a theoretical simulation is performed to evaluate the extent to which nanoparticle redistribution affects the temperature elevations and thermal dosage required to cause permanent thermal damage to PC3 tumors using magnetic nanoparticle hyperthermia. The micro-CT generated tumor model is attached to a mouse body model. The blood perfusion rates in the mice and embedded PC3 tumors were extracted based on the experimental data of recorded surface temperatures using an infrared camera. A previously determined relationship between the nanoparticle concentration and the volumetric heat generation rate is implemented in the theoretical simulation. Our simulation results show that the average steady-state temperature elevation in the tumors in the control group is approximately 2–3 °C higher than that in the experimental group when the nanoparticles are more spreading from the tumor center to the tumor periphery. Further, we assess heating time needed to cause permanent thermal damage to the entire tumor, based on the nanoparticle distribution in each tumor. The heating time in the experimental group is at least 44% longer than that in the control group. The modified damage model by Pearce excluding the thermal damage in an initial time delay leads to a shorter heating time than that predicted by the traditional thermal damage model. The results from this study suggest that the heating time needed when considering the influence of nanoparticle migration during heating is probably between 19 and 29 min based on the more accurate Pearce thermal damage model. In conclusion, the study demonstrates the importance to include modeling of nanoparticle spreading during heating and accurate thermal damage model into theoretical simulation of temperature elevations in tumors to determine thermal dosage needed in magnetic nanoparticle hyperthermia design.

The research was performed in partial fulfillment of the requirements for the Ph.D. Degree in Mechanical Engineering by Manpreet Singh from the University of Maryland Baltimore County, USA.

• National Science Foundation (Grant No. CBET-1705538; Funder ID: 10.13039/100000001).