Numerical Model for Magnetic Fluid Hyperthermia in a Realistic Breast Phantom: Calorimetric Calibration and Treatment Planning

This paper aims to apply a proposed, based on calorimetric measurements, a reliable numerical model for magnetic fluid hyperthermia (MFH) treatment planning of breast cancer. Furthermore, we perform a comparative analysis of magnetic nanoparticles (MNPs) and tumour tissue interactions by means of the magnetic-field-dependent Néel and Brownian relaxation times. The analysis was based on an anatomically correct breast model (developed in-house) and a modified linear response theory, which was applied to investigate the heat dissipation from the magnetic nanoparticles dispersed in the breast tumour. The calculations of the single-domain magnetic power losses were conducted for a case where the magnetic field value and the applied frequency were known, but also for the different concentrations of the MNPs in the tumour. Two scenarios were considered: The MNPs mobilised and immobilised in the tumour. In parallel, the eddy currents effect, together with the related temperature distributions, were considered in order to analyse safety issues. By changing the MNP concentration in the tumour, the corresponding temperature distributions were calculated. The eddy current effect, together with the related temperature distribution, were considered in order to analyse safety issues. Varying the MNP concentration in the tumour, the corresponding temperature distribution was calculated. Moreover, the cumulative equivalent minutes at 43 ℃ were analysed. In the anatomically correct breast phantoms, the tissue location can lead to “hot spots” due to the eddy current effect and subsequently to the high gradients of the temperature. That is why the analysis of safety issues related to the overheating side effect should be taken into consideration during the treatment planning of magnetic fluid hyperthermia. The phenomenon of heat dissipation from MNPs is very sophisticated and depends on their concentration, the distribution and the relaxation mechanism in the tumour, together with magnetic field strength and frequency. Furthermore, we inferred that the phenomenon of heat dissipation from MNPs equally depends on MNP-tissue interactions, and it can lead to 30% differences in the power assessment. Nevertheless, the aforementioned factors should be considered in parallel using anatomical, volume-dependent models to enhance the efficiency of in vivo treatment.


Introduction
Novel therapy methods without negative impacts are needed to reduce mortality rates in patients with breast cancer. One such method is magnetic fluid hyperthermia (MFH) in which magnetic nanoparticles (MNPs) dissipate heat during radiofrequency (RF) application [1][2][3][4][5], activating biochemical pathways leading to necrosis or apoptosis [6][7][8][9]. In the past decade, attention has been given to MFH treatment for breast cancer either on its own or in combination with chemotherapy and radiotherapy [10][11][12][13][14]. However, dexterous techniques of modelling and validation will help to improve and standardise this potential clinical treatment. Effective MFH is accomplished when maximum cytotoxicity with low or no side effects are attained with a limited dosage of Alternating Magnetic Field (AMF) susceptible MNPs, so suitable heat characterisation of MNPs and treatment planning should be carried out prior to therapy. Besides, ethical factors, cost-effective clinical trials and infrastructure requirements restrict the number of teams performing in vivo studies, hampering significant progress in understanding MFH for the therapy of breast cancer. Furthermore, constraints such as those for the pain threshold of 4.85 × 108 A/m/s for induced eddy currents in patients during MFH treatment [9] and the field limit being sufficient to saturate all the particles within the distribution at a particular frequency, limit MFH treatment. However, it is possible to use virtual breast phantom models to evaluate the efficacy of MFH and the risk assessment of the radiofrequency.
Three dimensional (3-D) models were used in the past to investigate the effects of MNP volume fraction, nanocomposite geometry, and treatment parameters on thermal profiles regarding female breast cancer [15][16][17]. However, the proposed models were simplified to the semi-ellipsoidal block, the semi-sphere or to the anatomically correct shape without naturalistic breast tissues. This study aims to develop an improved numerical breast phantom model derived from human magnetic resonance images [18] to obtain a lifelike model. Moreover, the literature lacks numerical models, which include 3-D temperature distribution due to eddy currents and single-domain magnetic power losses in anatomically correct phantoms. Furthermore, the numerical investigation conducted in this work was based on the results received from calorimetric measurements of the ferrofluid using the magneTherm TM system [19,20], i.e., the numerical model developed for heat generation terms. Specific loss powers (SLP) were validated with calorimetric measurements and are now used to analyse breast cancer hyperthermia treatment mediated by magnetic fluid.
It is clear that the reported SLPs values [20] do not reflect the situation one might expect to find in human tissues where, for example, blood perfusion, spatial tissues location and magnetic nanoparticles-tissue interactions affect heat dissipation. That is why, in this work, the reported SLPs values together with the numerical model, have been transposed to realistic breast cancer.
Our previous studies presented mathematical models, i.e., the single-domain magnetic power losses model and the eddy current model. Both of them were used later in the modified Pennes equation to evaluate a temperature distribution. The current work details the anatomically correct breast model used in the simulations and its parameters, followed by several analyses of volumetric power losses, which determine spatial temperature distribution as a function of time for various concentrations of magnetic fluid preparations. Furthermore, mobilised and immobilised MNPs in the tumour were considered.

Results and Discussion
The mathematical model described in the previous section was used to investigate the temperature distribution during MFH. Different magnetic fluid concentrations were applied in the tumour to control the temperature in the target tissue, then the tumour position in the breast model was determined by maximum electrical loss. These power density values are referred to as "hot spots". In addition to the power losses due to the eddy current effect, a power dissipation from magnetic nanoparticles should be considered in parallel when dealing with MFH treatment planning. Moreover, particular interest was paid to the comparative analysis of MNPs-tissue interactions, as in MFH, and the heating rates of MNPs which were determined using Néel and Brownian losses, which can be inappropriate for MFH treatment. For a reliable model of heating, two cases were considered, i.e., the MNPs free to move in the tumour (Néel and Brownian relaxations included in the tumour) and the MNPs immobilised in the tumour (in which only Néel relaxation takes place).

Breast Phantom Model
In our previous study, the breast model containing only skin, breast fat and tumour tissues was evaluated based on field application by a pancake coil [17]. In the present study, more tissues were added to achieve a very detailed anatomical structure. The breast phantom model was taken from [18], and it was assigned to the Class 3 category-heterogeneously dense (51-75% glandular). The following breast tissues were included in the model: Breast gland (blue), breast fat (yellow), fat (green), skin (red) and muscle (orange) as shown in Figure 1. Moreover, the physical parameters of the tumour were considered as the same as the muscle tissues ones but with perfusion. The physical parameters of tissues were calculated for the frequency f = 171 kHz based on the Information Technologies in Society (IT'IS) database [21]. The parameters are shown in Table 1, where HTR = ρ b c b ρω is the heat transfer rate. In our previous work [17], the tissues included in the breast phantom models were averaged due to their volume fraction together with their physical parameters. Such approximation can be useful while planning the treatment when the applicator position pattern is required. In the current work, the model was improved, and the effect of the eddy current and "hot spots" (the places where high values of power losses may occur due to the anatomical configuration of breast tissues) were examined. The model cross-section is shown in Figure 1.

94
In our previous study, the breast model containing only skin, breast fat and tumour tissues was 95 evaluated based on field application by a pancake coil [17]. In the present study, more tissues were 96 added to achieve a very detailed anatomical structure. The breast phantom model was taken from 97 [18], and it was assigned to the Class 3 category-heterogeneously dense (51-75% glandular). The

98
following breast tissues were included in the model: Breast gland (blue), breast fat (yellow), fat (green), 99 skin (red) and muscle (orange) as shown in Figure 1. Moreover, the physical parameters of the tumour 100 were considered as the same as the muscle tissues ones but with perfusion. The physical parameters 101 of tissues were calculated for the frequency f = 171 kHz based on the Information Technologies in 102 Society (IT'IS) database [21]. The parameters are shown in Table 1, where HTR= is the heat transfer rate. In our previous work [17], the tissues included in the breast phantom models were averaged

Eddy Currents Effect on Temperature Distribution
The authors would like to emphasise that the aim of this part of the study is not the design a breast cancer applicator, but to fulfil the condition for the magnetic field strength (H max = 10.3kA/m in the middle of the applicator) and the frequency, as was performed during the calorimetric experiments using the magneTherm TM system. The location of the spherical tumour (2 mL volume) in the breast model was determined by the maximum volumetric eddy power density (P eddy ) to consider the worst-case scenario.
The volumetric power density due to the eddy currents was evaluated for frequency f 1 = 171 kHz, and is shown in Figure 2. It is obvious that the highest value of eddy power density occurred on the skin tissue, reaching 70 kW/m 3 . That is why the data was rescaled to the maximum value (35.5 kW/m 3 ), i.e., to one of the "hot spots", which were placed in the model excluding the skin layer to improve visualisation. It is important to note that the presence of "hot spots" should be considered during treatment planning as they are hard to predict and they can limit the amount of energy that can be deposited in a patient's body. Moreover, one should remember that maximal electrical power losses will occur on the surface of the tumour as the eddy currents have to circulate around the tumour centre where they are the smallest. This sometimes can be used as an advantage because high-conductivity tumours surrounded by low-conductivity tissues will have a local eddy current flowing around the approximate centre of the tumour leading to increased heating of deep-seated tumours.

118
The authors would like to emphasise that the aim of this part of the study is not the design a 119 breast cancer applicator, but to fulfil the condition for the magnetic field strength ( = 10.3 kA m ⁄ 120 in the middle of the applicator) and the frequency, as was performed during the calorimetric 121 experiments using the magneTherm TM system. The location of the spherical tumour (2 mL volume) 122 in the breast model was determined by the maximum volumetric eddy power density ( ) to 123 consider the worst-case scenario.
The volumetric power density due to the eddy currents was evaluated for frequency = 171 kHz,

125
and is shown in Figure 2. It is obvious that the highest value of eddy power density occurred on the 126 skin tissue, reaching 70 kW/m 3 . That is why the data was rescaled to the maximum value (35.5 127 kW/m 3 ), i.e., to one of the "hot spots", which were placed in the model excluding the skin layer to 128 improve visualisation. It is important to note that the presence of "hot spots" should be considered 129 during treatment planning as they are hard to predict and they can limit the amount of energy that  In our case, to investigate the worst-case scenario, the tumour was positioned in one of the "hot spots", and the power density was calculated (see Figure 2-right). One can see that the tumour position has changed the power distribution pattern and the "hot spot" shifted due to the physical properties of the tumour. However, taking the temperature distribution into account (see Figure 3), the area where the temperature rises above 39 • C (which is relatively high), means that overheating of tissues can occur primarily on the skin.  Table 2. 154 Table 2. Data used to calculate single-domain magnetic power losses [20].

156
The values of maximum magnetic power losses for the different relative concentrations are given

157
in Table 3 and also presented in Figure 4.
158 Table 3. Single-domain (S-D) magnetic power losses with regard to the relative concentration ( ) of 159 ferrofluid in the tumour [22].

Single-Domain Magnetic Power Loss Analysis
As the first step, the single-domain magnetic power losses for both the mobilised and the immobilised MNPs were calculated. The ferrofluid concentration was the same as what is was during the calorimetry experiments (ϕ = 5.0 mg/mL), and it was used as a reference. The parameters which have been used in Formula (1) are provided in Table 2. Next, the ferrofluid concentration was changed (φ r ) to control the temperature in the tumour. The values of maximum magnetic power losses for the different relative concentrations φ r are given in Table 3 and also presented in Figure 4. Table 3. Single-domain (S-D) magnetic power losses with regard to the relative concentration (φ r ) of ferrofluid in the tumour [22].     Taking into account the single-domain magnetic power losses, one can see that the powers had decreased by 30% when MNPs were immobilised in the tumour. It is also easy to linearly approximate both cases as follows: P mobilzed = 244.4φ r − 0.7826, P immoblized = 170.3φ r − 0.2539 in order to predict the magnetic power losses for other magnetic fluid concentrations.
For both cases (the mobilised MNPs on the left and the immobilised MNPs on the right) in Figure 5, the temperature distribution patterns are shown in the breast phantom for different concentrations of the ferrofluid after 1800 secs of heating. The black area in Figure 5A-C indicates the temperature of 40 • C for different φ r = [1, 0.5, 0.35]. Moreover, the temperatures over time in the middle of the tumour are presented in Figure 6 for the mobilised case.     Table 3). Finally, Figure 7 presents the relative concentration

180
In the case of the immobilised MNPs in the tumour, it can be seen that both the maximum power 181 dissipation (see Table 3), as well as the temperature, has reduced (see Figure 6). For the reference  Table 3). Finally, Figure 7 presents the relative concentration 190 as the function of the temperature difference between the two cases together with its linear 191 approximation. In the case of the immobilised MNPs in the tumour, it can be seen that both the maximum power dissipation (see Table 3), as well as the temperature, has reduced (see Figure 6). For the reference concentration, φ r = 1.0, the power decreased by 30%, while the temperature difference, ∆T, was approximately 5.4°C when comparing the two cases of MNPs in the tumour. Nevertheless, the relationship between MNPs concentration and the temperature increase is not straightforward, for example, for φ r = 0.5 the temperature difference is ∆T = 2.1°C and for φ r = 0.35 it is ∆T = 1.2°C. However, in order to compare the aforementioned cases, the concentration, φ r , was changed for the corresponding values of the power dissipation (for example, to receive P moblil = 121.0 kW/m 3 , the concentration should be φ r = 0.5 for the mobilised MNPs, while for immobilised MNPs concentration should be φ r = 0.71, see Table 3). Finally, Figure 7 presents the relative concentration as the function of the temperature difference between the two cases together with its linear approximation.
According to the data in Table 3 and Figure 5, it is possible to control the temperature in the tumour and its vicinity by changing the magnetic fluid concentration. The volumetric power density, required to raise the temperature to 43°C in the tumour was at least 84.2 kW/m 3 , as depicted in Figure 8, where the isosurface 43°C can be seen inside the tumour volume. According to the data in Table 3 and Figure 5, it is possible to control the temperature in the

195
According to the data in Table 3 and Figure 5, it is possible to control the temperature in the  However, in order to reach 43°C in the whole tumour volume 121.0 kW/m 3 was required but, in this case, the temperature inside the tumour was 45°C. On the other hand, it should be underlined that the heating power may have increased due to a possible agglomeration of MNPs in the tumour. Moreover, the heterogeneous temperature distribution in the tumour vicinity, which can be observed in Figures 4, 6 and 9 may lead to high-temperature gradients even if the eddy currents power losses are minimal when compared with the magnetic power ones. This phenomenon should be considered in patient treatment planning.

246
A model dedicated to the numerical analysis of heat transfer in human tissues was proposed by

247
Pennes [24]. In this model, blood perfusion is assumed to be uniform throughout the tissue and all 248 heat leaving the artery is absorbed by the local tissue, with no venous rewarming. Moreover, the

249
Pennes model is limited to the perfusion source term, i.e., the arterial temperature is assumed to be 250 equal to the body core temperature.

251
To account for the strong temperature dependence due to bio-regulatory processes [25], the 252 breast model tissues perfusion rates were assumed to be constant ( ( ) = . ) while the tumour 253 perfusion was assumed nonlinear according to the following formula [26]: which indicates the mass flow rate as presented in Figure 9.

Single-Domain Magnetic Power Losses Model
When polydisperse magnetic nanoparticles in the superparamagnetic regime are exposed to an AMF with the given magnetic field strength, H max , and frequency, f, the magnetisation lags behind the external field. According to linear response theory (LRT), the heat dissipation in W/kg can be expressed as the specific loss power (SLP) [23]: where g(µ, s, V M ) is the volume-weighted distribution (log-normal distribution), ρ and V M are the density and the magnetic volume of MNP, respectively, and χ is the average out-of-phase component of the susceptibility given by: where k B is the Boltzmann constant, m = M s V M is the magnetic moment where M s is the magnetic saturation and τ is the Néel-Brown relaxation time. The Néel (τ N ) and Brownian (τ B ) metrics, as dependent on an external magnetic field (H max ), are given by [23]: where τ B = 3ηV H k B T and η is the viscosity, h = H max H k and µ 0 H k = 2K M s represents the anisotropy field and f 0 = 10 9 s −1 is the attempt frequency. Finally, the existence of two relaxation times leads to the effective relaxation time given by τ(H max ) −1 = τ N (H max ) −1 + τ B (H max ) −1 . In this case, it is worth mentioning, that the above model contains few simplifications such as follows: All the MNPs have spherical shapes, MNPs do not create clusters or chains, and therefore there are not any inter-particle interactions and finally-the hydrodynamic particle size distribution is uniform. One should realise that if the concentration of the MNPs in the ferrofluid is relatively small and, for example, does not exceed 5 mg/mL, the effects caused by dipole-dipole interactions between the nanoparticles are considered insignificant at room temperature. Thus, in our case, the nanoparticles were considered to be magnetically isolated from each other.

Eddy Current Effect Model
When human tissues are exposed to an alternating electromagnetic field, the eddy current effect is observed due to non-zero conductivity of the tissues, which finally leads to their heating.
The magnetic field strength (H) generated by an applicator can be solved using the magnetic vector potential (A): where J s is the current density vector (A/m 2 ) and H = µ −1 0 ∇ × A. In order to include a volumetric power density (P eddy ) produced by an eddy current (J eddy ) effect, a magneto-quasi static algorithm was applied. Considering the scalar potential (u) together with a magnetic vector potential (A), one can formulate the issue as: where σ is the conductivity (S/m) and ω = 2πf is the angular frequency. Finally, P eddy can be expressed as: Discussing the eddy current or Foucault current is essential due to their effect on different types of tissues when exposed to a time-varying magnetic field. The eddy current effect can lead to unwanted non-specific heating of healthy tissues, so it must be considered in the applicator design and during treatment planning.

Pennes Equation
A model dedicated to the numerical analysis of heat transfer in human tissues was proposed by Pennes [24]. In this model, blood perfusion is assumed to be uniform throughout the tissue and all heat leaving the artery is absorbed by the local tissue, with no venous rewarming. Moreover, the Pennes model is limited to the perfusion source term, i.e., the arterial temperature is assumed to be equal to the body core temperature.
To account for the strong temperature dependence due to bio-regulatory processes [25], the breast model tissues perfusion rates were assumed to be constant (ω(T) = const.) while the tumour perfusion was assumed nonlinear according to the following formula [26]: which indicates the mass flow rate as presented in Figure 9. Hence, the temperature distribution in the breast model was investigated using the modified bio-heat transfer equation: ρc ∂T ∂t = ∇ · (κ∇T) + ρ b c b ρω(T)(T b − T) + ρHGR + ρSLP LRT + P eddy W m 3 (9) where ρ is the tissue density, c is the tissue-specific heat, ρ b is the density of blood, c b is blood specific heat, κ is the tissue thermal conductivity, ω(T) is the perfusion rate, T b is the arterial blood temperature [C], T is the local temperature, HGR is the metabolic heat generated rate and SLP LRT is the external power losses due to magnetic nanoparticles. The convection boundary condition (Robins type) was used: where h is the heat transfer coefficient, and κ is the thermal conductivity of the skin layer.

Conclusions
The modified Pennes heat equation with Robin's boundary condition to compute the temperature distribution within the realistic breast phantom was investigated. The single-domain magnetic and eddy current mediated power loss densities were applied as the parallel external heat sources. It was shown that due to the anatomical structure of the breast tissues, the "hot spots" could occur together with high-temperature gradients. That is why this phenomenon should be considered during treatment planning.
More importantly, the distribution of the magnetic power losses along with the temperature was not homogeneous in the cancer area. It was also demonstrated that the viscosity played a crucial role, and the magnetic fluid concentration was critical for effective MFH treatment. The behaviour of MNPs within the tumour volume was also analysed as their lack of mobilisation leads to lower power dissipation. Consequently, it is vital to know if the MNPs are free to move in the interstitial fluid or whether they are immobilised in the tumour tissue. Thus, magnetic nanoparticles-tissue interactions hold the equal potential to affect treatment.
To conclude, from the MFH-treatment planning point of view, it is crucial to prevent healthy surrounding tissues from overheating, so particular attention should be given to the temperature distribution in the tumour vicinity. The accuracy of magnetic fluid thermal therapy determines its role in prospective treatment planning for the individual patient because inadequate thermal doses to the tumour can cause a failure of treatment. Hence, the proposed numerical model will help to calibrate magnetic fluid temperature measurements when exposed to radiofrequency, and realistic models like the one presented here will enable in silico evaluation of MFH efficiency prior to clinical treatment.

Conflicts of Interest:
The authors declare no conflict of interest.