Evaluation of magnetic nanoparticles for magnetic fluid hyperthermia.

Abstract Background: Magnetic nanoparticles (MNPs) generate heat when exposed to an alternating magnetic field. Consequently, MNPs are used for magnetic fluid hyperthermia (MFH) for cancer treatment, and have been shown to increase the efficacy of chemotherapy and/or radiation treatment in clinical trials. A downfall of current MFH treatment is the inability to deliver sufficient heat to the tumor due to: insufficient amounts of MNPs, unequal distribution of MNPs throughout the tumor, or heat loss to the surrounding environment. Objective: In this study, the objective was to identify MNPs with high heating efficiencies quantified by their specific absorption rate (SAR). Methods: A panel of 31 commercially available MNPs were evaluated for SAR in two different AMFs. Additionally, particle properties including iron content, hydrodynamic diameter, core diameter, magnetic diameter, magnetically dead layer thickness, and saturation mass magnetization were investigated. Results: High SAR MNPs were identified. For SAR calculations, the initial slope, corrected slope, and Box–Lucas methods were used and validated using a graphical residual analysis, and the Box–Lucas method was shown to be the most accurate. Other particle properties were identified and examined for correlations with SAR values. Positive correlations of particle properties with SAR were found, including a strong correlation for the magnetically dead layer thickness. Conclusions: This work identified high SAR MNPs for hyperthermia, and provides insight into properties which correlate with SAR which will be valuable for synthesis of next-generation MNPs. SAR calculation methods must be standardized, and this work provides an in-depth analysis of common calculation methods.


Introduction
Hyperthermia has played an increasingly crucial role in treating cancer since ancient times, first having been mentioned in 1700 BC when Egyptians used hot rods to burn breast cancer cells [1]. However, its application was rekindled only in the 1960s when American surgeon George W. Crile Jr. discovered that under long-lasting heat stress tumors would shrink without significantly damaging the surrounding tissue [1,2]. Hyperthermia is based on the principle that heating malignant cells to around 40-44 C can induce their necrosis, apoptosis, or enhanced susceptibility to chemo-and radiotherapy [3]. Current hyperthermia treatments can be delivered locally, regionally, or to the whole body using instruments such as microwaves, radio waves, ultrasound, hot water perfusion, resisting wire embeds, or infrared radiators [4,5]. However, these methods of thermotherapy have some disadvantages. For example, local hyperthermia can be used only for superficially located tumors, while whole body hyperthermia can cause difficulties such as thermal pressure in the heart [5]. Magnetic fluid hyperthermia (MFH) can overcome these limitations and offer solutions to these problems.
MFH was first proposed in 1957 by Gilchrist et al. [6], and works by coupling magnetic nanoparticles (MNPs) to an alternating magnetic field (AMF) at frequencies between 100 kHz and 1 MHz in order to generate heat at a specific site within the body [7]. The foundation of MFH states that MNPs delivered to a tumor can specifically target the heat to the malignant cells and minimize the effect on surrounding healthy tissues. In an applied field, the atomic magnetic moment of MNPs will align with the direction of the field, where the degree of alignment depends upon the properties of the MNPs and the field. As the frequency of an applied AMF increases, the rotation of the magnetic moment begins to lag behind the field (N eel relaxation) resulting in magnetic hysteresis and heating. For multi-domain MNPs, the magnetic hysteresis arises due to domain wall reorganization and physical rotation of the particles (Brownian relaxation) [8]. Further explanation of the physics of MNP heating is discussed by Dennis et al. [8], Pankhurst et al. [9], and Carrey et al. [10].
MFH offers the possibility to locally heat tumors which can be found deeper within the body such as glioblastoma [5]. Despite clinical successes [11], pitfalls of current MFH treatment include the inability to deliver sufficient heat to the tumor due to insufficient MNP tumor loading, heterogeneous distribution of MNPs throughout the tumor, and loss of heat to the surrounding tissues [12,13]. Additionally, although heat generation increases with field strength and frequency, these parameters are limited in the clinic as harmful eddy currents may be generated at higher frequencies.
Heat generation by MNPs can be characterized by their specific absorption rate (SAR) [9,10,[14][15][16], and identification of MNPs with high SAR values will be invaluable to the success of MFH clinical treatments. In the context of MFH, SAR value is the power dissipation per unit mass of magnetic nanoparticles (or iron mass) in an applied AMF, and it is also known in the literature [17,18] as specific loss power (SLP).
High-SAR MNPs would be beneficial for research in many other areas utilizing magnetic heating, for example, magnetically triggered drug or biologic release [19][20][21], activation of latent growth factors [22], catalysis reactions in microfluidic reactors [23], or control of cell functions [24]. Commercially available MNPs are frequently used for a variety of research applications [25][26][27][28][29][30] and thus it would be useful for researchers to know their properties, particularly their magnetic characteristics. Unfortunately, commercially available MNPs are often not fully characterized or characterization data are limited. The SAR values for MNPs remain unreported in many cases. Therefore, detailed information on properties of commercially available MNPs will enable researchers to make informed decisions on which MNP will be best suited for their research. Three previous studies have sought to evaluate the heating potential of commercially available MNPs, but the variety and number of the MNPs evaluated in these studies are limited [28,31,32]. The current investigation will build upon the previous work to include a larger variety of companies and products, as well as comprehensive study of MNP properties in addition to SAR values. Since SAR values have high variation between labs due to poor standardization of measurement conditions [12,33], it is crucial for a wide variety of particles to be measured in one lab setting to provide an accurate direct comparison.
For this study, the heating potential of 31 commercial MNPs was measured in two different magnetic fields and quantified based on calculated SAR values and intrinsic loss parameters (ILP) [12]. Additionally, MNP properties including iron content, hydrodynamic diameter, magnetic diameter, core diameter, magnetically dead layer thickness, and saturation mass magnetization were evaluated. Correlations between these parameters and SAR values were analyzed using Pearson's product moment and Spearman's rank correlation coefficients. Identification of any specific properties of the MNPs which result in high SAR values will have a pivotal impact on the development of the next generation of MNPs for hyperthermia. To standardize the SAR values, we calculated SAR as energy per MNP mass and as energy per iron mass, as both of these annotations are commonly found in the literature. For SAR calculation, the initial slope method (ISM), corrected slope method (CSM), and Box-Lucas method (BLM) were used and validated. Evaluation of the SAR calculation methods and units will give researchers an idea of the variation of results that can be obtained due to different data analysis methods, and highlight the importance of specifying calculation method and units.

SAR measurements
AMFs were applied to samples using an Easyheat TM induction heater (Ambrell, Rochester, NY). The experimental setup is shown in Figure S1. Two different types of water-cooled (set-point 18 C) copper solenoid coils were utilized to test the heating potential of MNPs at two different field strengths and frequencies. The first coil (vertical; length ¼ 0.072 m; diameter ¼ 0.032 m; n-turns ¼ 8) operates at 401.1 A and 345 kHz, with a measured field strength of 59 kA/m, while the second (horizontal; length ¼ 0.135 m; diameter ¼ 0. 045 m; n-turns ¼ 14) operates at 270.9 A and 223 kHz, with a measured field strength of 41 kA/m. Field strengths B a (T) were calculated according to the formula: where U a is the root mean square voltage measured using an AMF probe (Life Systems, Everett, WA) (V); S a is axial sensitivity factor of the probe; f is the frequency (kHz). Aqueous dispersions of MNPs (200 mL, 1 mg MNP/mL) were measured in the center of the coil at least in triplicate. MNP concentration was expressed in mg MNP/mL rather than in mg Fe/mL due to two reasons: (1) most suppliers provide MNP dispersion concentration as mg MNP/mL, and (2) MNP samples with low iron content will have high nanoparticle concentration which can potentially cause sample aggregation. Custom 3D printed polymer holders were used to precisely place the samples in the center of the solenoid. Distilled water was used as a control sample for both fields. Fiber optic temperature sensors (NeOptix; OD ¼ 1 mm) were used to measure the sample and coil temperature for 180 s total in 1 s intervals: 30 s with no field application, 120 s with field application, followed by 30 s of no field application. Pyrex vials (ID ¼ 6.0 mm; OD ¼ 8.0 mm; height ¼ 45.0 mm, V ¼ 1.24 ml) were used for the 59 kA/m coil, while polystyrene 8 Stripwell TM plates (PN 9102, Corning; ID ¼ 6.4 mm; OD ¼ 8.5 mm; height ¼ 12.0 mm, V ¼ 0.36 ml) were used for the 41 kA/m. The coil was cooled to 18 C between measurements.
The coil system was arranged within a Styrofoam TM foam insulted box ( Figure S1) to decrease temperature fluctuation. For the coil operating at 345 kHz/59 kA/m, Styrofoam foam was also placed around the sample holder to ensure further insulation of the sample, while for the 223 kHz/41 kA/m coil, the 3D printed plastic holder for the sample required sample placement in the middle of the coil and insulation. Due to the different arrangements, heat loss may differ slightly at the two field strengths; however, all samples were treated the same way at each field strength for internal consistency. Parafilm V R film was used to seal the vials with the MNP samples. Three temperature sensors were inserted through the Parafilm V R film and positioned in the middle of each sample, while a fourth sensor was used to record the coil temperature.
SAR calculations were performed using the ISM, the CSM, and the BLM. The formula for SAR calculation using the ISM is [34]: where DT=Dt is the slope of the first 60 s of AMF exposure (K/s); C is the heat capacity of fluid (J/(KÁg)); q is the density of the fluid (g/mL); m Fe=MNPs is the concentration of iron or MNPs in magnetic fluid (mg/mL). The ISM assumes that the sample temperature increases linearly during heating and that heat losses are negligible during the time interval used. Given that some researchers report SAR values in W per mass of particle [35][36][37] instead of per mass of Fe [38][39][40], the calculations in this paper were performed using both units.
The CSM considers linear losses of the system due to a non-adiabatic set-up. The formula below was used to calculate SAR [12]: where L is the linear loss parameter (mW/K). The L was calculated in Excel V R software (Microsoft Office) using the open-source corrected slope data analysis algorithm [41]. This algorithm estimates L by finding the value that results in the smallest standard deviation between SAR values on the same heating curve (30 s intervals, 5 fits).
In the phenomenological BLM, the heating curve is fitted using the following expression: SAR values were calculated using the fitting parameters a and b as: The linear loss parameter L (W/K) for an individual sample was calculated as: The intrinsic loss power parameter (ILP, nHÁm 2 /kg) was proposed for normalization of SAR values measured at different magnetic field amplitudes/frequencies [12]: where H is the field strength (kA/m); f is the frequency (kHz). It should be noted that the ILP is limited to measurements for MNP samples with a dispersity greater than 0.1 in applied fields below the saturation field for the MNPs, and field frequencies between 10 5 and 10 6 Hz [32]. Additionally, the thermodynamic losses should not be greater than the power input from the field [42].

Iron quantification
The iron content in each sample was quantified in triplicate using a colorimetric method adapted from the American Society of Testing and Materials (ASTM) standard method E394-15 for detection of trace amounts of iron [43]. Ten to 50 mL of each sample (depending on sample concentration) was digested for 12 h in 1 ml of 70% HNO 3 at 101 C to remove any coatings that may have been present. Ten to 30 mL of the acid solution with digested particles was evaporated for an hour at 115 C, followed by the addition of 46 mL of deionized water and 30 mL of hydroxylamine hydrochloride (8.06 M) to each sample at room temperature. The solution was vortexed and left to react for 1 h to reduce the Fe(III) to Fe(II). After the reduction step, 49 mL of 1.22 M sodium acetate were added to each sample to buffer the solutions. Fe (II) cations were chelated by adding 75 mL of 13 mM 1,10-phenanthroline to form a colored complex (tris(1,10-phenanthroline) iron (II)) which have a maximum absorbance at 508 nm (e ¼ 14,332 ± 4594 L/(molÁcm) at pH 4.8) (see Figure S2). The absorbance of the solutions was measured in 96 well plates (PN 9017, Corning, Fairport, NY) using a SpectraMax V R M5 microplate reader (Molecular Devices, Sunnyvale, CA) in PathCheck V R mode. Iron concentration of samples was calculated using a standard calibration curve. Standard solutions of known iron concentration were prepared using an ICP iron standard in the concentration range from 0 to 5 mg/mL.

Dynamic light scattering (DLS)
Ultraviolet-visible (UV-vis) spectroscopy was used to determine the optimal concentration of MNP dispersions for DLS measurements. A SpectraMax V R M5 (Molecular Devices, Sunnyvale, CA) microplate reader equipped with a 24 well microplate and a 0.5 mm cover glass (Molecular Devices, SpectraDrop TM micro-volume starter kit, PN 0200-6262) was used to record UV-vis spectra. A calibration curve for each sample at 658 nm (the wavelength of the DLS laser) was made and used to estimate the optimal concentration of MNPs for DLS. A practical criterion is that the absorbance for a 1 cm cuvette must be smaller than 0.04 to exclude the effect of multiple scattering [44]. Using this criterion, the upper concentration threshold was determined for each type of MNP from the calibration curves (see Table S1), and this concentration was used for DLS measurements.
ZetaPALS TM analyzer (Brookhaven Instruments Corporation, Austin, TX) equipped with a 658 nm laser (35 mW) was used to determine particle size distribution (PSD, intensity-weighted) and hydrodynamic diameters of aqueous MNP dispersions. Each sample (100 lL) was measured in a plastic cuvette at 90 in triplicate (10 runs, 30 s) with the dust cutoff filter set to 40.00.
Intensity-weighted size distribution obtained from combined runs, or its first peak in the case of multimodal distribution was fitted to a lognormal function: where A is the area under the curve; d s is the scale parameter; r is the shape parameter. The arithmetic mean (d arithm:mean ), median (d median ), and mode (d mode ) diameters; variance (varðdÞ); and dispersity ( -D) were calculated from parameters of lognormal function according to the formulas below [45]: Given the non-symmetrical nature of pore size distribution (skewness), the mode diameter (d mode ) was used for the size characterization of the MNPs.

Nanoparticle tracking analysis (NTA)
A NanoSight TM LM14 instrument (Malvern Instruments Ltd., Worcestershire, UK) equipped with a 405 nm laser (max. power 5 W) was used to determine the number-weighted particle size distributions, hydrodynamic diameters, and estimated nanoparticle concentrations (particles/mL) of MNPs dispersions. Samples at a concentration of 10 mg MNP/mL in distilled water were measured at least in triplicate. For each type of MNP, three video files (30 fps) with a duration of 60 s were captured and analyzed using NanoSight NTA 2.3 analytical software. The data obtained using the NTA software was analyzed using the same methodology as described in DLS (Section 2.4).

Vibrating sample magnetometry (VSM)
For measurements of magnetic properties, 0.1 mg of MNPs were immobilized into agarose films. Briefly, 2% (w/v) agarose solution was prepared by dissolving agarose in distilled water upon heating with simultaneous mixing. After being completely dissolved, 500 mL of agarose solution was cooled to approximately 50 C and vortexed with 0.1 mg of MNPs. The agarose solution with MNPs was then transferred into a Costar V R 48 well-plate (PN 3548, Costar, Washington, DC) and dried overnight at 50 C.
For VSM measurements, films were fixed on 8 mm Pyrex V R sample holder (PN 1600-SHP, MicroSense) with double-sided carbon mounting tape. Magnetization of MNPs was measured in magnetic fields with strength from À1 T to 1 T (GMW Associates, Magnet Systems Model 3473-70 Electromagnet). Data were recorded using EasyVSM TM software (ADE Technologies, Campbell, CA). Measurements were performed once per film, while each magnetic field point was measured 4-times per measurement. The magnetization of blank agarose gel samples was measured in triplicate and used for the slope correction by its subtraction from the magnetization sample curves. Obtained data was plotted as mass magnetization (emu/(g MNP)) versus magnetic field (T) using OriginPro TM software (OriginLab, Northampton, MA). For comparison, mass magnetization of MNPs at 1 T (an average of four points) was used.

Magnetic diameter analysis
The magnetic diameter was calculated by fitting the obtained magnetization versus field curves (VSM data) to the Langevin-Chantrell function as previously described [46,47] using Igor Pro TM 6.37 software (WaveMetrics, Tigard, OR). The domain magnetization for magnetite was assumed to be equal to that of bulk magnetite (446 kA/m) [48]. The arithmetic mean for the volume-weighted magnetic diameter was found (D arithm:mean mv ), and was used later in calculations for the magnetically dead layer.

Transmission electron microscopy (TEM)
For TEM analysis, particle samples were placed on a copper grid with a carbon supporting film. A FEI Tecnai F20 S/TEM microscope with a Schottky field emission ZrO/W gun source was used to image the following samples: Micromod 79-01-102 and Sigma-Aldrich 747335 (lots MKBW1867V and MKCD1194), 900028, and 900042. ImageJ software (version 1.51j, NIH Image J system, Bethesda, MD) was used to measure the core diameter of the particles (n ¼ 100) from TEM images. For the rest of the MNPs, the TEM diameter provided by the manufacturer was used for calculation of the magnetically dead layer thickness (D mdl ). The magnetically dead layer is a surface layer of the magnetic particle in which magnetic dipoles are disordered and magnetic properties are degraded [46,49,50]. The thickness of the magnetically dead layer was calculated as where D TEM is the core diameter obtained from TEM images.

Statistical data analysis
The data are presented as an arithmetic mean ± standard uncertainty of the mean (sðYÞ) [51]. The sðYÞ was estimated as where s is the sample standard deviation; N is the number of parallel experiments; t p ðmÞ is the critical value for p percent of two-tailed t-distribution with defined degrees of freedom (m); Dy is the error of indirect measurements, which was calculated using the numerical differentiation method [52]: where f is a function of directly measured variables x i ; Dx i is the error of the direct measurements. The systematic error was considered negligible for all measurements.
Differences between groups were examined for statistical significance using two-tailed paired t-test or one-way ANOVA for a significance level a ¼ 0:05: The paired t-test was used to compare differences between ILP values for particles obtained at 345 kHz/59 kA/m and 223 kHz/41 kA/m. The oneway ANOVA was used to compare the three calculation methods (ISM, CSM, BLM). As per recommendation by National Institute of Standards and Technology (NIST), both the Bonferroni and Tukey's tests were used as post hoc tests for one-way ANOVA, and the smallest confidence band was selected [53].
To evaluate correlations between SAR (calculated by the BLM at 345 kHz/59 kA/m) and other parameters (iron content, hydrodynamic diameter, magnetic diameter, core diameter, magnetically dead layer thickness, and saturation mass magnetization), Pearson's product-moment correlation coefficient and Spearman's rank correlation coefficient were calculated. The latter one is the non-parametric version of the Pearson product-moment correlation. Both coefficients measure the strength and direction of association between two ranked variables. However, Spearman's correlation determines the strength and direction of the monotonic relationship between the two variables rather than the strength and direction of the linear relationship between the variables, which is the purpose of Pearson's correlation.
All calculations were performed using Microsoft Excel 2016 software, while Minitab 18 V R and OriginPro 2018 (b9.5.0.193) software were used for statistical analysis, nonlinear and linear fitting, plotting graphs, and graphical residual analysis.

SAR measurements and calculation
The temperature was measured in two different AMFs over time for each sample. The average and maximum temperature change are listed in Table SII. The average temperature difference for the water control was 0.6 ± 0.1 K for the 223 kHz/41 kA/m coil and 2.1 ± 0.5 K for the 345 kHz/ 59 kA/m coil. The same temperature heating curves were used to calculate SAR by three methods (ISM, CSM, and BLM) described previously in Section 2.2. The SAR values of MNPs in W/(g MNP) and W/(g Fe) calculated at 345 kHz/59 kA/m and 223 kHz/41 kA/m are listed in Tables 1 and 2, respectively. In Table 2, negative SAR values for Sigma-Aldrich 747440, 747459, 747467, and 900146 MNPs which do not have any physical meaning were replaced with zero. The SAR value for the water control samples was also calculated. In Tables 1 and 2, samples which have a SAR value similar (±3 SD) to the SAR value of water or whose average temperature change during heating is within three standard deviations of the mean provide no magnetic heating, and are highlighted in light gray and dark gray, respectively. The following MNPs: Micromod 79-01-102, Ocean Nanotech SEI-10-05 and SPP-05-02, and Sigma-Aldrich 544884, 637149, 747254, 747440, 747459, and 747467 have a high standard uncertainty due to low colloidal stability of sample (based on visually observed aggregation and sedimentation of MNPs), especially upon heating.
For visual comparison, the SAR values of MNPs in W/(g MNP) calculated at 345 kHz/59 kA/m and at 223 kHz/41 kA/m by three different methods (ISM, CSM, and BLM) are depicted in Figure 1.
Based on the results of the one-way ANOVA test, there were several MNP samples that had significantly different SAR values depending on the SAR calculation method (see Table SIII). For example, at 345 kHz SAR values for 10 samples are significantly different at a ¼ 0.05. According to the posthoc tests, the majority of cases correspond to the difference between ISM and BLM.
Both methods, ISM and CSM, are approximations of the phenomenological BLM method which provides the most accurate calculation of SAR among all three methods. Yet BLM is not a complete model for describing heating of MNPs in an applied AMF as concluded based on the results of graphical residual analysis ( Figure 2).
To perform the regression analysis, the data were fitted to the three different models (ISM, CSM, and BLM). The quality of the fitting was judged based on the values of adjusted R 2 , while the validity of the models was assessed using graphical residual analysis (see Figure 2). A visual examination of the residuals was conducted to check for obvious deviations from randomness or the possible presence of heteroscedasticity. The following were used in graphical residual analysis to assess the validity of linear and non-linear fits: (1) scatter plot of the response and predicted values versus the independent variable; (2) scatter plot of the residuals versus the independent variable; (3) scatter plot of the residuals versus the predicted values; (4) histogram of the residuals; and (5) normal probability plot of the residuals [54].
As an example of residual analysis, the plots for HyperMAG B MNPs are depicted in Figure 2. As can be seen from the graphs corresponding to BLM, the residuals are normally distributed around zero (Figure 2(e)), and a normal probability plot of the residuals is approximately linear and has a characteristic slightly s-shaped appearance (Figure 2(d)). There are no outliers, and the fit looks better using non-linear BLM (Figure 2(a)) rather than its linear approximation ISM (Figure 2(f)); though, for both models a non-random structure (curved pattern) in the residuals can be observed (Figure 2(b, c, g, h)).

High-SAR MNPS
The top-6 MNPs with the highest SAR values are highlighted in each column in red in Tables 1 and 2. The top 25% MNPs based on average temperature change are listed in Table SII and highlighted in yellow. The ranking of MNPs based on their SAR values marginally depends on magnetic field parameters and SAR units. Sigma-Aldrich 900042 MNPs have the highest BLM-SAR value: 2158.4 ± 110.1 W/(g MNP) at 345 kHz/59 kA/m. These MNPs have an iron content of 977. 8 ± 55.4 lg Fe/(mg MNP), a magnetic diameter of 12.8 ± 5. 8 nm, a mass magnetization value of 90.7 ± 4.5 emu/(g MNP) at 1 T, a core diameter of 22.9 ± 0.2 nm, and a magnetically dead layer thickness equal to 5.1 ± 2.9 nm. They can be heated by 36.9 ± 2.0 K (see Table SII) in 120 s under AMF application (345 kHz/59 kA/m) at a concentration of 1 mg ; values highlighted in dark grey mark samples whose temperature change during heating is within three standard deviations of the mean (Table S3); values highlighted in light grey have SAR values below SAR value of water (mean þ 3r). For each product, the determined iron content per particle is listed, as well as the SAR calculated by three methods: initial slope method (W/g MNP and W/g Fe), corrected slope method (W/g MNP and W/g Fe), and the Box-Lucas method (W/g MNP and W/g Fe). For notation, see notes for Table 1. MNP/ml, indicating that they can be an optimal candidate for MFH treatment.

ILP values
In addition to the differences in the final results due to the calculation methods, SAR values increase with increasing magnetic field frequency/strength. Thus, the ILP (nHÁm 2 /kg) was calculated to normalize the BLM-SAR values W/(g MNP) measured at two different magnetic fields. The corresponding data can be found in Table SIV and in Figure S3. The ILP values at different field strengths and frequencies are in good agreement within the standard deviation of the mean. Nevertheless, 10 samples were found to have significantly different ILP values between the two field strengths and frequencies (highlighted in blue in Table SIV).

Other MNP properties and their correlation with SAR values
To better understand heating characteristics of MNPs, the physical and magnetic properties were evaluated and compared with information on MNP properties provided by manufacturer (Table 3). Table 3 contains diameters provided by companies, hydrodynamic diameters measured by DLS and NTA, and the magnetic diameter (Fe 3 O 4 ) found by fitting magnetization curves to the Langevin-Chantrell function. Particle coatings and functional groups, which might affect the MNP properties, are also listed in Table 3. Additionally, the mass magnetization (r) at 1 T was determined using VSM and presented in Table 3 as well. TEM images and magnetization versus field (M versus H) curves for selected MNPs are depicted in Figure 3. All MNP samples imaged have similar hydrodynamic diameters (DLS) ranging from 45.7 ± 0.6 nm to 57.1 ± 12.3 nm and magnetic diameters (VSM fitting) ranging from 9.9 ± 4.4 to 12.8 ± 5.8. However, the SAR values for these particles change from 95.5 ± 13.8 W/(g MNP) to 2158.4 W/(g MNP) based on the BLM at 345 kHz. The SAR value and core diameter calculated from TEM images are listed in the bottom right corner of the TEM images.
M versus H curves shows that most particles are superparamagnetic in DC fields. Low field plots of MNPs from Figure 3 can be found in Figure S4 and may suggest that a small amount of nanoparticle population can be magnetically blocked. This also includes Sigma-Aldrich MNPs 773352, while 725331, 747440, 747459, and 747467 nanoparticles appear to be paramagnetic. Based on the discrepancy between physical core diameters from TEM and the magnetic diameter calculated from VSM curve fitting with the Langevin-Chantrell function, the magnetically dead layer thickness for each sample was calculated and listed in Table SV.
The correlation coefficients and corresponding p values between BLM-SAR values at 345 kHz/59 kA/m and iron content, hydrodynamic diameter, magnetic diameter, mass magnetization, core diameter, and magnetically dead layer thickness can be found in Table 4. Scatter plots showing the directions and strengths of the correlations can be found in Figure S5.
The iron content, mass magnetization, and magnetically dead layer thickness were found to have significant and positive correlations with SAR values based on Pearson's and Spearman's correlation coefficients. The magnetic diameter and the core diameter were found to have a significant and positive correlation with SAR values based on Spearman's correlation coefficient alone. Based on Spearman's correlation coefficients, in most cases the correlations are moderate (0.50-0.69 range) and moderate to low (0.30-0.49 range), with the strong correlation between SAR and magnetically dead layer thickness (0.83).

SAR values and calculation methods (ISM, CSM, and BLM)
Determined SAR values for MNPs often vary between labs due to poor standardization of methodology and measurement conditions. Therefore, the goal of this study was to compare a wide variety of MNPs under the same conditions. In addition, companies often do not provide detailed information on physicochemical, colloidal, and magnetic properties of MNPs, with SAR, in particular, being rarely reported. Thus, an extensive evaluation and characterization of 31 commercial MNPs including SAR (measured at two different AMFs), ILP, and physical and magnetic properties was performed to help to guide researchers in identifying the optimal MNPs for their studies. As a result, the top-6 MNPs highlighted in red in Tables 1 and 2    Two-tailed test of significance has been used; parameters significantly (at the 0.05 level) and positively correlated are highlighted in blue; correlation coefficients were rounded to the nearest hundred.
which includes a  [57] and 600 W/g (g of material and calculation method not specified) at 11 kA/m and 410 kHz [58], corresponding to ILP values of 3.8 nHm 2 /kg and 11.7 nHm 2 /kg, respectively [32]. In addition to MFH, photothermal therapy using nanoparticles is a promising approach for hyperthermia studies. One study used a panel of non-commercial nanoparticles to compare MFH and photothermal therapy and showed a maximum SAR value of $1000 W/(g Fe) for MFH using iron oxide nanocubes at 900 kHz/19 kA/m, and 1.4 mg Fe/mL and $12,000 W/(g Au) for photothermal therapy using gold nanorods at 808 nm, 1 W/cm 2 , and 0.15 mg Au/mL [59]. This study [59] also compared heating properties of nanoparticles within an aqueous environment, inside cells, and within tumors, and showed that MNPs lose heating power when internalized by cells. Despite similarities to ILP and SAR values from the literature, the experimental set-ups for heating measurements can have potential sources of error [18,60,61], with sample volume being one of them. In a few studies [18,60] it was shown that decreasing sample volumes results in lower SAR values. However, it remains unclear if the sample volume effect is statistically significant given that the SAR values were reported without standard uncertainty. Furthermore, another study showed that smaller sample volumes (0.25 mL) have lower rates of heat loss compared to larger volumes (1.0 mL) [12]. Apart from that, conduction losses and energy storage within the sample holders and different sample geometries were discussed as a potential source of error as well [12,60]. Additionally, it has been shown that the location of the temperature probe can have an effect on the SAR value [12,18,60]. Standardization and calibration of experimental set-up for SAR measurements using non-adiabatic systems is further discussed in references [12,62].
In our study, consistent volumes, sample vials, and temperature probe placement were used to minimize these potential sources of error. Additionally, the sample vials were placed in exactly the same location (center) within the solenoid coil using custom 3D printed sample holders. It is important to note that water control samples in our study showed insignificant heating (Table SII). The standard uncertainty for SAR measurements was less than 20% and also accounted for the error propagation during SAR measurements. Higher uncertainty in some samples can be explained by poor colloidal stability of MNP dispersions, especially upon heating, which resulted in poor reproducibility of heating curves for a few samples. Some MNPs quickly aggregated during heating (e.g. Sigma-Aldrich 544884 and 747459), whereas most MNPs remained stable during the measurements.
The SAR calculation methods were found to have a significant effect on its value based on the results of the oneway ANOVA test (Table SIII). Thus, it is crucial for the calculation method of SAR to be reported by other scientists. Between the three methods (ISM, CSM, and BLM) used for SAR calculation, the BLM provides the most accurate SAR estimation for the non-adiabatic system based on the results of graphical residual analysis. The ISM assumes only linear increases in sample temperature and ignores heat losses, while the CSM and BLM both include linear losses due to a non-adiabatic experimental set-up. The first two methods, ISM and CSM, are linear approximations of the phenomenological BLM and can be used in certain situations. ISM would be accurate for adiabatic measurements with no chemical or physical changes within the samples [61].
It should also be noted that different time intervals for ISM-SAR calculation have been used in literature and can affect its accuracy [18,61]. In this study, we chose 60 s time interval since the initial temperature rise was linear within this region for most MNPs. A study by Bordelon et al. [56] has established a protocol for empirically determining appropriate fitting intervals for the ISM to achieve results comparable to the BLM; however, it was only performed for three MNP types and would need to be studied in more detail to ensure it will be accurate for variety of MNPs. This method was further developed by Soetart et al. [61] to identify time ranges of heating which satisfy quasi-adiabatic conditions and it was discussed that an arbitrary choice of a single time interval will likely result in unreliable results. SAR calculation methods are also evaluated by Wildeboer et al. [12] and they found BLM and CSM to give more accurate SAR values than the ISM. However, this study [12] suggested that CSM gives the most accurate results because CSM-SAR values were higher than BLM-SAR value for the selected sample (28 mg/ mL Ferucarbotran) even though there is no direct connection between absolute values and model validity.
Herein, we used the graphical residual analysis (see Figure  2) described in Section 3.1 to validate ISM and BLM. However, as one can find in literature, many authors comparing various models often claim a model to be better/worse than others based only on marginal differences in correlation coefficients. However, using R 2 or adjusted R 2 is not sufficient and Anscombe's quartet is a good illustration for this point. Besides, the value of R 2 depends on the degrees of freedom. It means that by simply adding variables to a data set, one can increase the value of R 2 , but this does not improve a fit.
To exclude this effect, adjusted R 2 can be used instead. Further, R 2 , measures only the percentage of variance explained by the data and does not tell us anything about "strength of relationship" or "goodness of a model." Even though it does measure goodness of a fit, R 2 does not determine the goodness of a model overall. Therefore, it is necessary to diagnose the regression results and validate the tested models using the residual analysis, either graphical or quantitative.
The adjusted R 2 values for both models (ISM and BLM) shown in Figure 2 are very close to 1 and suggest a good fit for both models. Even though the BLM is the more accurate model compared to ISM as judged by visual examination of residuals, the observed trend in the scatter plot of the residuals versus independent/predicted variables indicates that the BLM is not complete. This leaves room for the development of more comprehensive models for SAR calculations based, for example, on non-linear heat losses [63].
Overall, the results of this study underline the importance of standardizing SAR calculation method. For the case of non-adiabatic systems, which are predominantly used in the literature [61], it is recommended that the BLM to be the method of choice for SAR calculation. Besides the calculation methods, the units for reporting SAR play an important role as well and the SAR values expressed as W/(g MNP) differ from values expressed as W/(g Fe). Until the method for SAR calculations and units is standardized, it is recommended that scientists provide both values in W/(g MNP) and W/(g Fe). The SAR value in W/(g MNP) is more indicative of the heat generation by the sample and it is known in most cases; thus, it should be the chosen annotation as it more accurately represents the heating that will be achieved in real-world applications. SAR values calculated as W/(g MNP) were generally lower than the SAR values calculated as W/(g Fe). Due to low iron content (approximately less than 50 mg Fe/(mg MNP)) and error propagation during calculation, SAR values reported in W/(g Fe) sometimes can have high standard uncertainty which makes such SAR values ambiguous and less precise. Thus, one can obtain a misleading result: samples which generate almost no heat will have high SAR values in W/(g Fe) if they have a low iron content (e.g. Sigma-Aldrich 747440, 747459, and 747467).

ILP values
The ILP, which corrects SAR for extrinsic parameters, field strength and frequency, was calculated for two measured field strengths and frequencies for each MNP. In general, the ILP ( Figure S3) did a sufficient job at normalizing the SAR values, even though significant differences were found for 10 of the evaluated MNPs.
The ILP values must be taken with reservations due to the narrow domain of validity for the relationship SAR a f Á H 2 0 [8,33,56,64]. Carrey et al. [10] examined the domain of validity for the linear response theory that this relationship is based upon and showed that it is mainly useful for highly anisotropic MNPs. As mentioned in the methods section, the ILP parameter is only valid for frequencies up to several hundred MHz, samples with a dispersity of more than 0.1, and field strengths well below the saturation field of the MNPs [32]. In addition, to compare results between different systems, similar environmental thermodynamic losses must be involved [32]. Further, the ILP parameter is only valid for single domain nanoparticles and cannot be used to compare single domain versus multidomain particles.

Other MNP properties and their correlation with SAR
MNP properties including iron content, hydrodynamic diameters, magnetic diameters, mass magnetization, core diameters, and thickness of the magnetically dead layer were determined for all samples in order to find correlations with SAR. Identification of any specific properties of the MNPs with high SAR values may have a pivotal impact on the development of the next generation of MNPs for hyperthermia. Importantly, correlations between SAR and iron content (moderate), magnetic diameter (moderate), mass magnetization (moderate), core diameter (moderate to low), and thickness of the magnetically dead layer (strong) were found. According to the existing models for MNP heating, differences in heat production arise for different MNPs due to their size and anisotropy, and the interplay between size and anisotropy is crucial [8]. Anisotropy can arise from numerous aspects of the MNPs including magnetocrystalline, shape, or colloidal anisotropy [8]. The optimal value/range for MNP size and anisotropy is dependent on field amplitude and frequency; for a given frequency and size, heat generation will increase with anisotropy until the anisotropy field is greater than the magnetic field [8]. Thus, there is a "sweet spot" for the ratio between the anisotropy and magnetic diameter for each magnetic field/frequency combination. This agrees with the findings that correlated magnetic diameter and core diameter with SAR. Other parameters for which correlations with SAR were found including iron content, mass magnetization, and thickness of the magnetically dead layer will have important effects on the anisotropy of the MNPs. Anisotropy can arise from numerous aspects of the MNPs including magnetocrystalline, shape, or colloidal anisotropy [8].
Microscopy (TEM) was used to elucidate whether the core shape and size is responsible for the discrepancies in SAR values. As one can see in Figure 3(a), Micromod particles 79-01-102 with SAR value of 95.5 ± 13.8 W/(g MNP) have an irregular shape which would lead to high shape anisotropy. For a given frequency/field strength, if the anisotropy field is higher than the magnetic field, the MNPs will not produce heat [8]. While it was difficult to discern the exact core size for the Micromod particles 79-01-102, it is clearly smaller than the core size for the rest of MNPs imaged with TEM, which ranged from 22.9 ± 0.2 nm to 28.8 ± 0.3 nm.
All other MNPs depicted in Figure 3 have relatively high SAR values and appear euhedral. The facets on the MNPs could contribute to anisotropy effects or to enhanced chain aggregation effects (i.e. colloidal anisotropy). Work by Martinez-Boubeta et al. has shown that single domain iron oxide nanocubes provide improved magnetic heating efficiency compared to spherical MNPs of analogous sizes, which was attributed to the increased anisotropy of the cubic MNPs and their tendency to aggregate into chains due to the cubic shape [65]. It remains unclear whether the euhedral shape found in the commercial particles is beneficial to SAR or not as it was found in all imaged samples and needs to be studied further. Further, a recent publication by Dennis et al. has shown that differences in heat generation from MNPs with similar size and bulk magnetic properties can be due to variations in the internal magnetic structure that influence coupling of magnetic domains within MNP cores [66].
The magnetically dead layer thickness was calculated by subtracting the magnetic diameter found from fitting of VSM curves to the Langevin function from the core diameter determined from TEM images [46]. A strong positive correlation between SAR values and the magnetically dead layer thickness was observed. Thus, the thickness of the magnetically dead layer might contribute to the anisotropy of the MNPs, even though other studies have stated that a small magnetically dead layer will result in improvement of magnetic performance of the nanoparticles such as saturation magnetization [46,67]. Therefore, this finding needs to be further studied in order to justify the strong correlation found between increased magnetically dead layer thickness and increased SAR value.

Conclusions
A comprehensive study of physical and magnetic properties of 31 commercial MNPs including calculation of SAR and ILP was performed. High-SAR MNPs such as Sigma-Aldrich 900042, 900028, 747335, 900062, and nanoTherics HyperMAG B were identified and three SAR calculation methods (ISM, CSM, and BLM) were validated using a graphical residual analysis. Based on the results of model validation, it is suggested that other researchers use the non-linear BLM for SAR calculation and express BLM-SAR values in both units: W/(g MNP) and W/(g Fe). Positive correlations (from moderate to strong) between SAR values and other MNP properties were found and discussed in detail. Studied MNP properties, especially the magnetically dead layer thickness, might contribute to the nanoparticle anisotropy and could be responsible for efficient heating of MNPs in an applied AMF.

Future directions
Two different batches of Sigma Aldrich 747335 MNPs were evaluated, but the data was not enough to make any comments on the statistical relevance as only two batches were assessed. However, scientists should be aware that batch to batch variability may exist, and should be studied further.
Another limitation of the current study is the frequency/field strength products evaluated are outside of the Atkinson's limit. Atkinson [68] identified that the biological range of AMF applications for humans should have a product of field strength and frequency less than 485 kA/mÁkHz. However, this is currently an area of debate as the value can be adjusted for coil diameter, tissue diameter, and electrical conductivity of the tissue [69][70][71]. As described by Carrião et al. [70], the free current loss is proportional to the square of the distance from the coil axis and thus a better estimation for the critical field for biomedical applications can be found using the equation Hf < ð0:150=rÞ Á 485kA/mÁkHz. Additionally, field conditions viable for human use could be expanded in the future as coil design is an area of active development.
Moreover, MNP heating studies performed over a range of magnetic field strengths and frequencies can provide important information about the samples [56,61]. Therefore, in future studies, select MNPs will be evaluated at numerous field strengths/frequencies. Further research for the selected MNPs will be performed at lower field strengths to assess their practical use in biomedical applications in addition to MNP testing in cellular and tumor environments as it was done in previous studies [59,[72][73][74]. Preliminary data collected by our lab using commercial MNPs immobilized in composite polymeric microparticles has shown that they continue to produce high amounts of heat when MNP incorporation is less than 1% of the microparticle volume. In addition, preliminary studies using HeLa cells with these commercial MNPs at a concentration of 1 mg MNP/mL showed high heating rates and subsequent cell death. Thus, it is anticipated that the MNPs with the highest SAR values identified in this study will preserve their heating characteristics once internalized by cells. Other parameters that were not studied in this work such as the crystal structure, crystal defects, and purity might also affect MNP heating efficiency and will require further investigation.