**2. Results**

**Optical and structural characterization of Fe**3**O**<sup>4</sup> **NPs.** Structural characterization of the Fe3O4 NPs—purchased commercially—were done to verify the specification provided by the manufacturer and also predict the optical absorption coefficient. X-ray diffraction spectra of the Fe3O4 NPs revealed the presence of peaks at 2*θ* = 31.5◦, 35.8◦, 38.35◦, 42.75◦, 47.2◦, 54.04◦, 57.24◦, and 62.75◦ (Figure 1a). The observed peaks correspond to diffraction planes: (220), (311), (222), (400), (110), (422), (511), and (440), which have been attributable to cubic spinel phase of Fe3O4 (space group, *Fd-3m*, JCPDS-#19-0629). Since no other prominent phase was detected, the result implied that the NPs are essentially crystalline Fe3O4. Transmission electron method (TEM) confirmed the morphology of the NPs to be spherical (with agglomerations) and size distribution to be between 15 and 20 nm in diameter as indicated by the manufacturer (Figure 1b). The agglomeration revealed in the TEM image have been attributed to dipolar coupling between the NPs [22,23]. For any NP, its NIR photothermal effects are controlled by their NIR optical absorbance. UV-vis-NIR spectra of the NPs showed an extended optical absorption that slowly increased in the NIR region relative to the visible light region (Figure 1c). The absorbance intensity at 810 nm increased linearly with concentration, from 0.35 ([Fe3O4] = 6 mM) to 1.51 ([Fe3O4] = 24 mM) (Figure 1d). The absorbance band in the NIR region of UV-vis-NIR optical spectra is consistent with the results in the literature and has been attributed to multiple charge (electron) transfer [24]. Furthermore, the linear increase of absorbance for the range of concentration tested in this work has been previously reported elsewhere [19,25]. Shen and co-workers [25] showed that saturation starts occurring at high concentration (100 mM, absorbance values > 3 at 808 nm). In an effort to translate the experimentally measured photothermal heat generation capabilities of the Fe3O4 NPs tested in this study, we followed the flow-chart shown in Figure 1e to obtain the extinction cross sections of the MNPs, which was then used in Equation (2) to predict absorbance, *A*pred. The validity of *A*pred was tested by evaluating its agreement with the experimentally measured absorbance, *A*exp, for the different concentrations of Fe3O4 (6, 12, 24 mM). We observed that the predictions agreed reasonably well with experiments to within 2% for all concentrations when the sample size, *n*, in Equation (5) was equal to 5 (see Figure 1f).

**Figure 1.** Structural and optical characterization results. (**a**) X-ray diffraction spectra at a power of 45 kV × 40 mA. (**b**) Transmission electron microscopy at magnification of 0.5 mm, Scale bar: 50 nm. Absorbance as a function of (**c**) wavelength (*λ*) and (**d**) concentration ([Fe3O4]) at *λ* = 810 nm. (**e**) Flow-chart for the comparison of theoretical predictions and experiment measurements. (**f**) Comparison of the *A*pred and *A*exp for different [Fe3O4] (6, 12, 24 mM).

**Photothermal effects of Fe**3**O**<sup>4</sup> **NPs.** The influence of laser power (*P*<sup>0</sup> = 0.5 and 1.0 W) and NP concentration ([Fe3O4] = 0–24 mM) on photothermal effects was accessed in aqueous solution (deionized water) to quantify their heat generation capabilities under an irradiation duration of 5 min. Pure deionized water—containing no Fe3O4 nanoparticles—was used as a control. The rate of change of the temporal curves increased with concentration at 5 min independent of the laser power that was used (Figure 2a,b). For *P*<sup>0</sup> = 0.5 W, the temperature change, Δ*T*, increased approximately by 44.4% (from ≈9 to 13 ◦C) when concentration was increased from 0 to 24 mM (Figure 2c). When the power was increased to 1.0 W, Δ*T* increased by approximately 83.3% (from ≈12 to 22 ◦C) for the same concentration. Photothermal conversion efficiency, *η*exp, decreased with concentration and laser power (Figure 2d). For instance, *η*exp for the 6 mM solution decreased from approximately 66% to 51% when *P*<sup>0</sup> was increased from 0.5 to 1.0 W. Furthermore, when the concentration was increased from 6 to 24 mM, *η*exp decreased from 46% to 39% using the same power regimes. Generally, the trend of Δ*T* recorded in this study was in agreement with measured absorbance properties and also consistent with previously reported studies [18,19,26]. For small NPs (<30 nm) and low concentrations, absorption dominates scattering leading to high *η*exp. On the other hand, scattering dominates the extinction efficiency as nanoparticle size or concentration is increased. As [Fe3O4] increases, clusters are formed due to the high surface area to volume ratio of nanoparticles [27]. These clusters act as large particles to enhance scattering leading to the reduction in *η*exp [28]. Several approaches are available for the prevention of clusters.

**Figure 2.** Photothermal characterization results. Temporal response curves for different concentrations after 5 min of irradiation with laser powers: (**a**) *P*<sup>0</sup> = 0.5 W and (**b**) *P*<sup>0</sup> = 1.0 W. Comparison of the corresponding (**c**) temperature change (Δ*T*) and (**d**) experimental photothermal conversion efficiency (*η*exp) as a function of laser power. Error bars: s.d.

**Computational modeling of NP-mediated photothermal heating of breast tumor.** The use of computational model as quantitative frameworks enables assessment and customization of the treatment parameters (NP concentration, treatment duration and irradiation protocols: duration and laser power) to potentially enhance efficacy. Thus, FEM simulations were applied to approximate photothermal heating of a Fe3O4-containing tumor embedded within a female breast using the optical diffusion approximation of the transport theory [29] and the Pennes bioheat transfer equation [30].

Figure 3 shows a schematic of 2D representation of the axisymmetric geometry of the computational model. It was configured as a heterogeneously dense [31] multi-layer block of tissue with proportions assigned according to the Breast Imaging Reporting and Data System (BIRADS) developed by American Cancer Research [32]. It consisted of various layers of normal tissue with unequal thickness. The dimensions of the model were chosen to represent a "heterogeneously dense" breast model [31], which consists of 20% muscle layer, 60% glandular layer and 20% fat layer. Also, a tumor is located at 55 mm from the base. The laser source was assumed to be a diode laser 810 nm placed close to the top surface of the breast model. The inset is a fragment of geometry showing control points P1−P4, where temperatures were recorded. The assigned optical, thermal and physical properties of different tissue layers were approximate values obtained from the literature [31,33–36]. Nanoparticles were assumed to be intravenously injected and uniformly distributed throughout the tumor.

**Figure 3.** FEM geometry. Schematic of the photothermal therapy consisting of a normal multi-tissue breast domain with an embedded spherical tumor (blue sphere) and NIR (810 nm) laser source. Inset: Fragment of geometry showing controls point P1−P4, where temperature were recorded.

To characterize the temperature and thermal damage profiles, we simulated temperaturecontrolled heating at a maximum tumor temperature, *T*max = 85 ◦C, for *t* = 15 min. The radius of the tumor, *R*, and *P*0, were chosen to be 2.5 mm and 1 W respectively. The predicted temperature distribution (Figure 4a) was revealed to be non-uniform with the maximum temperature occurring within the tumor and decreasing radially outwards into the surrounding tissue. The latter suggests that the heat transfer was predominantly conductive. For the case of the predicted thermal damage shown in Figure 4b,c, it can be seen that the entire tumor area, plus margins of up to 1 mm around it, was completely destroyed (Ω = 100%). A comparison of temporal response curves for temperatures (Figure 4d) at different control points (Figure 3) within the tumor (P1)and at the tumor-gland boundaries (P2–P4) revealed that the temperature rise as well as the final value was higher at (P1) relative to the boundaries: P2 (top), P3 (bottom) and P4 (side). This phenomenon can be attributed to factors such as relatively low blood perfusion and high metabolic heat of the tumor leading to high retention of heat within the tumor [37]. However, at all the locations, the temperature plateaued after about 2–3 min. The consequence of the high temperature within the tumor is revealed in corresponding

predicted temporal curves for the thermal damage (Figure 4e), which shows that 100% thermal damage occurs faster in the innermost part of tumor (P1)—≈3 min—compared to the peripherals, which take up to about ≈10 min (P4). Consistent with the literature [33,38], the model predictions showed the dependence of thermal damage spatial profile on the temperature distribution, which decreased with distance away from center of the tumor (see Figure 4f).

**Figure 4.** Simulation results. Cross-sectional view of the (**a**) temperature distribution, (**b**) thermal damage, (**c**) thermal damage showing the lesion parameter. Temporal response curves for (**d**) temperature and (**e**) thermal damage at the control points (P1–P4, cf. Figure 3). (**f**) Temperature and thermal damage as a function distance from P1. Simulation settings: *P*<sup>0</sup> = 1 W, *t* = 15 min and *T*max = 85 ◦C.

Ablative temperatures between 60 and 100 ◦C cause irreversible damages to key cytosolic and mitochrondrial enzymes [39,40]. For any tumor ablation therapy to be considered successful and thus reduce the chance of recurrence, it is critical to ensure that the entire volume of the tumor reaches therapeutic temperatures that ensures complete thermal damage (Ω = 100%). Such a goal can be achieved through the use of an appropriate maximum temperature, which takes into consideration the tumor dimensions. For NP assisted photothermal therapies such as the one being proposed in this study, maximum ablative tumor temperatures, *T*max, can be controlled by varying parameters such as NP number density, *N* (or volume fraction, *φ*v), the laser power, and treatment duration, *t*. To demonstrate this, a parametric study was used to determine *N* required to achieve a given *T*max (70, 85, 100 ◦C) and the corresponding volume of the lesion *V*<sup>L</sup> for different tumor sizes, *R* (1, 2.5, 5 mm). *V*L, was assumed to be spherical [41,42]; its radius, *R*L, was calculated as half the axial length of the predicted cross-sectional area where Ω = 100% (see Figure 4c). A summary of the results is presented in Table 1. The simulations were run with *P*<sup>0</sup> = 1 W and *t* = 15 min. Generally, it can be observed that *T*max required to achieve complete thermal damage increased with size of the tumor. For instance, *T*max = 70 ◦C produced a lesion with *V*<sup>L</sup> = 2.95 mm3, which was insufficient

to completely ablate the entire volume of tumor with *R*<sup>T</sup> = 1 mm (*V*<sup>L</sup> = 4.19 mm3). On the other hand, *T*max = 85 produced a lesion with *V*<sup>L</sup> = 2.95 mm3, which was big ensure to ensure complete thermal damage. Since *P*<sup>0</sup> was held constant for all simulation, it meant that *N* had to be increased to achieve the given *T*max. The results reveal that *N* required to achieve *T*max = 70 ◦C decreased with tumor size. For instance, *<sup>N</sup>* required to achieve *<sup>T</sup>*max = <sup>70</sup> ◦C decreased from 112.37 × <sup>10</sup><sup>14</sup> mL−<sup>1</sup> to 5.54 × 1014 mL−<sup>1</sup> when *<sup>R</sup>*<sup>T</sup> was increased from 1 to 5 mm. Lastly, the nanoparticle concentrations that were required to achieve the different values of *T*max corresponded to volume fractions in the range between 0.004% and 10.6%. A review of the nanoparticle delivery to tumors in the literature between 2006 and 2016 by Wilhelm et al. [43] revealed that only approximately 1% of administered nanoparticle dose reached the tumor. Therefore, it is important that the *φ*<sup>v</sup> is kept at the low value for practical applications. This can be achieved by through several means such as increasing the laser power or exploiting the capability of the Fe3O4 NPs to generate synergistic heat during simultaneous exposure to NIR laser and alternating magnetic field as previously reported elsewhere [19].


**Table 1.** Comparison of volume, *V*L, of predicted lesions and the number density, *N*, of nanoparticles (or volume fraction, *φ*v) used to achieve maximum tumor temperatures, *T*max (70, 85, 100 ◦C) in different tumor sizes, *R* (1, 2.5, 5 mm). *R*<sup>L</sup> is the radius of the lesion.

These predictions are consistent with previously reported experimental and computational results in the literature. Kannadorai et al. [44], developed a treatment planning model for the optimization to parameters such as laser power density, nanoparticle concentration and exposure time in an effort aimed at potential enhancement of treatment outcome. Their predictions showed that any change made to any of the parameters can be compensated by altering the remaining parameters. Using an integrated strategy that combined x-ray computed tomography or ex-vivo with a 4-dimensional FEM model, Maltzahn and co-workers [20] simulated photothermal heating with polyethylene glycol PEGylated gold nanorods (PEG-NR) and used the results to guided pilot therapeutic studies on human xenograft tumors in mice. Their simulations revealed the extension of thermal flux vectors from the region where PEG-NRs were located as well as the expected thermal profile.
