4.3.2. In-Vivo Predictions

The computational model is a multiphysics FEM model, thus, it took into account optical and thermal effects. Light distribution was based on the diffusion approximation of the transport theory [29] and temperature distribution by Pennes bio-heat transfer equation [30], which takes into account the effect of cell death on blood perfusion and the dependence of cell death and properties of the tissue. Cell death was determined by an Arrhenius based integral injury model [53].

*Light Distribution.* The optical diffusion approximation of the transport theory [29] was used to describe light distribution due to the dominance of scattering over absorption in biological tissues. It is defined by:

$$\frac{1}{c\_{\text{fl}}} \frac{\partial}{\partial t} \boldsymbol{\varrho} = D \nabla^2 \boldsymbol{\varrho} - \mu\_a \boldsymbol{\varrho} + \boldsymbol{S} \tag{6}$$

*cn* (m s<sup>−</sup>1) is the speed of light in a medium, *ϕ* (W m<sup>−</sup>2) is the fluence rate, *μ<sup>a</sup>* (m<sup>−</sup>1) is the absorption coefficient, *S* (W m−3) is the light source term and *D* = *μa*/*μ*<sup>2</sup> eff (m) is the diffusion coefficient. *μ*eff = 3*μa*(*μ<sup>a</sup>* + *μ*- *<sup>s</sup>*) (m−1) is the effective attenuation coefficient and *μ*- *<sup>s</sup>* (m−1) is the reduced scattering coefficient. Assuming that the light source was a continuous wave Gaussian NIR laser beam that was incident onto the breast model, the *ϕ* can be defined by

$$\phi(\vec{r}) = \frac{P\_0 \exp(-\mu\_{\rm eff} \vec{r} \cdot \hat{n})}{4\pi Dr} \tag{7}$$

where *P*<sup>0</sup> is the laser power and *n*ˆ is the direction of the beam. A summary of the values of the optical properties of the tissue used in the simulation is presented in Table 2.


**Table 2.** Optical properties of the biological domains that were used in the simulations. The values were obtained from Refs. [34–36].

*Temperature Distribution.* The Pennes bio-heat transfer equation [30] was used to estimate the temperature distribution. An additional term was added to account for the external heat source. The resulting equation is given by:

$$
\rho c\_{\rm P} \frac{\partial T}{\partial t} = \lambda(T) \nabla^2 T + \rho\_{\rm b} c\_{\rm b} \omega\_{\rm b}(\Omega) (T\_{\rm b} - T) + Q\_{\rm met} + Q \tag{8}
$$

where *ρ* (kg m<sup>−</sup>3) is the density, *c*<sup>p</sup> (J kg−<sup>1</sup> K<sup>−</sup>1) is the specific heat capacity at constant pressure. *λ*(*T*) (W m−<sup>1</sup> K−1) is the temperature dependent thermal conductivity, which is assumed to be a linear function defined by [54]:

$$
\lambda(T) = \lambda\_{(\text{\textdegree C} \text{\textdegree C})} [1 + 0.0028(T - 293.15\text{K})] \tag{9}
$$

where *T* (K) and *T*<sup>b</sup> (K) are the normal body and arbitrary temperatures, respectively. *ρ*<sup>b</sup> is the density of blood, *c*b, the specific heat capacity of blood and *ω*b(Ω) is the coefficient of blood perfusion assumed to be dependent on the cell damage, Ω, and defined by [33,38]:

$$
\omega\_b(\Omega) = \begin{cases}
\omega\_b^0 & \text{if } \Omega = 0 \\
(1 + 25\Omega - 260\Omega^2)\omega\_{b'}^0 & \text{if } 0 < \Omega \le 0.1 \\
(1 - \Omega)\omega\_{b'}^0 & \text{if } 0.1 < \Omega \le 1 \\
0, & \text{if } \Omega > 1
\end{cases} \tag{10}
$$

*ω*0 *<sup>b</sup>* (*s*<sup>−</sup>1) is the baseline coefficient of blood perfusion. *<sup>Q</sup>*met (W m<sup>−</sup>3) is the metabolic heat. *<sup>Q</sup>* accounts for external heat sources, which varies for the different domains of geometry. The heat generated after the absorption of NIR light is defined as *μaϕ*(*r*) (W m−3) and *Nσ*a*ϕ*(*r*) (W m−3) for the tissue and tumor domains respectively. *N* is the number volume of Fe3O4 NPs and *σ*<sup>a</sup> (m2) is the absorption cross-section of nanoparticles. Table 3 presents a summary of the values of the thermo-physical properties that were used in the simulation.

**Table 3.** Thermo-physical properties of the biological domains that were used in the simulation. The values were obtained from Refs. [31,33]


*Thermal Damage.* The Arrhenius injury model was used to estimate tissue destruction. The model, which relates temporal temperature to cell death, is defined by [53]:

$$\Omega(\tau) = A \int\_0^\tau \exp\left(\frac{-E\_a}{RT(t)}\right) dt\tag{11}$$

where *E*<sup>a</sup> (J mol<sup>−</sup>1) is the activation energy, *A* (s<sup>−</sup>1) is a scaling factor and *R* = 8.3 (J mol<sup>−</sup>1K−1) is the gas constant. The values for *<sup>E</sup>*<sup>a</sup> and *<sup>A</sup>* were obtained from Ref. [33] as 302 kJ mol−<sup>1</sup> and 1.18 × 1044 <sup>s</sup>−<sup>1</sup> respectively. Ω = 1 corresponds to the 100% irreversible cell damage.

*Model Implementation.* This FEM model was developed with the COMSOL Multiphysics 5.2 software package (Comsol, Inc. Burlington MA, USA). All properties and dimensions were added explicitly to the FEM model as parameters and variables under the "Global Definition" and "Model" nodes, respectively. Equations (9) and (10) were added as analytic functions under the "Global Definition" node. The 2D axisymmetrical model was used to reduce simulation time.

The light distribution was achieved by implementing Equation (5) as an analytic function. The temperature distribution was achieved using the bio-heat heat transfer application mode. Each tissue was represented with a separate "biological tissue" node. The boundary and initial conditions were specified as follows: a Dirichlet condition, T = 37 ◦C, at Γ1; a Neumann condition, **<sup>n</sup>**·(*λ*∇*T*) <sup>=</sup> *<sup>h</sup>*·(*T*ext − *<sup>T</sup>*) at <sup>Γ</sup><sup>2</sup> where the heat transfer coefficient, *<sup>h</sup>* was equal to 13.5 Wm−2K−<sup>1</sup> and *T*ext = 25 ◦C and continuity, **n**·(*λ*1∇ T1 − *λ*2∇ T2) = 0 at all interior boundaries. A temperature of 37 ◦C (for the normal body) was used as the initial temperatures in all domains of the model. The heat source was added to the bio-heat transfer application mode as a user-defined heat source.

The cell death model was implemented with the "Coefficient Form PDE" application mode. To achieve a time integration, the coefficients *da* and *f* were set to 1 and Equation (11), respectively. All other coefficients were set to zero. The initial conditions were: *S* = *∂S*/*∂t* = 0.

In order to enhance the accuracy of results, we resolved the model with successively smaller element sizes and compared results, until an asymptotic behavior of the solution emerged. The comparison was done by analyzing the temperature at the interface between tumor and the tissue. The choice of 2D-axisymmetric model allowed for the use of a physics-controlled mesh with the triangular element with sizes: maximum = 0.24 cm and minimum = 0.0024 cm for the tumor region and element sizes: maximum = 0.42 cm and minimum = 0.018 cm for the other regions. This resulted in 2656 domain elements and 277 boundary elements. The numerical solutions were obtained using the time-dependent solver "GMRES" with its default settings. The simulations were run on a mid-range workstation with Intel(R) Xeon(R) E5-1620 CPU and 8 GB of RAM.
