In-depth magnetometry and EPR analysis of the spin structure of human-liver ferritin: from DC to 9 GHz

Ferritin, the major iron storage protein in organisms, stores iron in the form of iron oxyhydroxide most likely involving phosphorous as a constituent, the mineral form of which is not well understood. Therefore, the question of how the ca. 2000 iron atoms in the ferritin core are magnetically coupled is still largely open. The ferritin core, with a diameter of 5–8 nm, is encapsulated in a protein shell that also catalyzes the uptake of iron and protects the core from outside interactions. Neurodegenerative disease is associated with iron imbalance, generating specific interest in the magnetic properties of ferritin. Here we present 9 GHz continuous wave EPR and a comprehensive set of magnetometry techniques including isothermal remanent magnetization (IRM) and AC susceptibility to elucidate the magnetic properties of the core of human liver ferritin. For the analysis of the magnetometry data, a new microscopic model of the ferritin-core spin structure is derived, showing that magnetic moment is generated by surface-spin canting, rather than defects. The analysis explicitly includes the distribution of magnetic parameters, such as the distribution of the magnetic moment. This microscopic model explains some of the inconsistencies resulting from previous analysis approaches. The main findings are a mean magnetic moment of 337μB with a standard deviation of 0.947μB. In contrast to previous reports, only a relatively small contribution of paramagnetic and ferrimagnetic phases is found, in the order of maximally 3%. For EPR, the over 30 mT wide signal of the ferritin core is analyzed using the model of the giant spin system [Fittipaldi et al., Phys. Chem. Chem. Phys., 2016, 18, 3591–3597]. Two components are needed minimally, and the broadening of these components suggests a broad distribution of the magnetic resonance parameters, the zero-field splitting, D, and the spin quantum number, S. We compare parameters from EPR and magnetometry and find that EPR is particularly sensitive to the surface spins of the core, revealing the potential to use EPR as a diagnostic for surface-spin disorder.


TEM description
Regionprops is a function implemented in MATLAB.In particular, we used the regionprops function implemented for MATLAB2019a.This function detects connected components in an image (objects).For each object, the total area, centroids, and many other properties can be extracted.For detecting the ferritin cores, we first binarized the images to enhance the contrast between the ferritin cores and the image background, then the center and radii of all the circular objects that are detected in the image are calculated.To separate the ferritin cores from other possible circular objects a threshold parameter (in the binarization preprocess) is needed to prevent the detection of false positives and minimize the false negatives (the undetected cores).In the image below, encircled in red are all the individual ferritin cores that were detected with the algorithm.Encircled in light blue is an example of a false negative.7.7% of ferritin cores are undetected.The diameter estimation has a standard deviation of ± 0.35 nm, which is the uncertainty.A limitation of this algorithm is the assumption that the ferritin cores are spherical, which is not necessarily the case.Additionally, ferritin cores that are formed by two sub-nanoparticles inside one ferritin will be detected as two different ferritin cores. .

Ferritin purity assessment
NativeMark Unstained Protein Standard (Cat.No. LC0725) was obtained from Life Technologies Corporation.The purity and the homogeneity of the protein were assessed by a 7.5% non-denaturing polyacrylamide gels electrophoresis.Gel electrophoresis was run at 4°C and 90 V until samples reached the running gel, then the voltage was increased to 160 V for 3 hours.After electrophoresis, Coomassie blue staining was used to stain proteins.For iron staining, the gel was immersed in Prussian blue staining solution, which was freshly prepared by mixing 2% K 4 Fe(CN) 6 and 2% HCl (1g Potassium ferrocyanide + 47,3 ml H 2 O + 2,7 ml HCl) for 1h at RT, then 1h in water, changing water every 15-20 min.Finally, the gel was incubated in a solution containing 0.025% DAB (3,3'-Diaminobenzidine-Sigma) and 0.05% H2O2 in 1X TBS [12,5mg DAB resuspended in 500 µl DMSO, added to 50 ml TBS 1x, 75-180 µl H 2 O 2 added just before incubation] for 15-30 minutes at RT to enhance the signal.

EPR baseline correction
The baseline correction of the spectra was performed by a home-made script implemented in Matlab 2019a, which automatically looks for a linear transformation (Fig. S4a), such that the slope of the baseline correction is the optimal value by which the high magnetic field values of the second integral of the spectra reach a plateau-like region (Fig. S4c, which is equivalent to having a zero first integral value at the highest field values (Fig. S4b).

Materials:
Nitric acid (65%, Suprapur®, Merck) was used in the digestion process, and further diluted 1% nitric acid was used as a carrying solution throughout the ICP measurements.National Institute of Standards and Technology (NIST)-traceable 1000 mg/L elemental standards were used (TraceCERT®, Fluka) for the preparation of calibration and internal standards.
Approximately 18 MΩ cm −1 water (MiliQ) was used in all sample preparation and analysis steps.

Instrumentation:
Calibration standards were prepared in a Secuflow fumehood (SCALA) to prevent contamination by atmospheric particulates.The standards and samples were analysed for trace elements using the NexION® 2000 (PerkinElmer) ICP-MS equipped with a concentric glass nebulizer and peltier-cooled glass spray chamber.An SC2 DX autosampler (PerkinElmer) was connected to the ICP-MS for sample introduction.SyngistixTM Software for ICP-MS (v.2.5, PerkinElmer) was used for all data recording and processing.

Elemental Analysis:
Five trace elemental calibration standards for ICP-MS analysis were prepared using NISTtraceable 1000 mg/L Fe standard: 0, 1, 5, 20, 100 ug/L.10 ug/L Rh and In were used as internal standard.Fe was analyzed in kinetic energy discrimination (KED) mode with 10% He gas to minimize polyatomic interferences.To check the calibration, samples were analyzed consisting of a blank measurement and a repeat measurement of one of the calibration standards.For the calibration curve only to be accepted as correlation coefficient (Cor.Coeff) was found higher than 0.999.EPR simulations, from 20 K to 210 K  S1.The simulation parameters used for all experimental spectra are shown in Table S1.The total spin for component 1 (E1) and component 2 (E2) is 10. Figure S8 shows the system-parameters temperature dependence for E1 and E2.A gentle slope is observed for the D values of both components, with a total decrease of 44% for both of them (Fig. S8a) as temperature increases.While the weight of E1 increases with the increase of the temperature, the weight of E2 decreases, although it remains preponderant over all the temperature range (see Fig. S8b).While decreasing the temperature, a gentle slope is observed for the H strain (Gaussian broadening) of E1 and E2.However, around 50 K a sudden jump is observed in both cases.

EPR simulations parameters
The relative proportion of E1 decreases with temperature, from ∼93% at 20 K to ∼75% at 210 K (Fig. S8a).Finally, both components display similar temperature dependencies on D and H strain (Fig. S8).The real temperature dependencies of E1 and E2 depend on the scaling factor given by the ratio between the unknown S real , and S = 10 used for the model.For instance, S real = 100 yields n = 10, and, according to the scaling laws of eq. 10 (see main article), the real temperature in Fig. S8 is 10 times larger than the modeled values.The parameters are given in Table S1.

Interpretation of D and H strain
The temperature dependence of the H strain , the Gaussian broadening of both components, see Fig. S S8, decreases with increasing temperature, revealing an averaging process that was previously 1 interpreted as anisotropy melting.We also observe a slight, though barely significant trend of D towards smaller values with increasing temperatures.Such a trend could also be a reflection of a reduction of the anisotropy of the ferritin core by some thermal averaging process.That the trend in D is not very pronounced suggests that the main effect of the anisotropy melting reflects itself in the reduction of the broadening parameter of the line, the H strain parameter, in the simulation model we employ.S1.For the entire temperature range, the line shape presents negligible changes, confirming the scaling method as a suitable candidate to reproduce and study the EPR line shape of high-spin systems.The parameters S and D are inversely related.That means, an increase (decrease) of S of one of the components followed by a decrease (increase) of D of the same component by the same amount, results in a simulation with the same lineshape as the original simulation.It's important to mention that a correction in the weight of the component needs to be applied to preserve the same simulation amplitude.The correcting factor for components 1 (n1) and component 2 (n2) is defined by

EPR S-D inverse compensation
where SX orig refers to the S used in the original simulations and SX new to the change in S that is being applied and will result in the opposite change in D, by the same factor.X = 1, 2 refers to EPR components 1 and 2, respectively.

Ferritin core and mononuclear Fe(III) EPR intensities
Mononuclear Fe(III) gives a characteristic EPR signal, with a sharp resonance at g ′ =4.3, the so-called g ′ =4.3 signal.This signal is also observed in the ferritin EPR spectra and the temperature dependence of its intensity is that of a paramagnetic species.Here we compare the amount of these species with respect to the ferritin core signal, the broad g = 2 signal.The EPR active spin concentration ratio of Fe(III) (g ′ =4.3) and the ferritin core spin system can be determined from 2 where g is the g factor, S electron spin number of each species and A is the enclosed area (second integral) that is calculated from the blue (mononuclear Fe(III)) and red (for the core) curves shown in Fig. S12.'Scan' is the sweep range, G is the relative gain of the signal amplifier, B m is the modulation amplitude, R is a factor related to the hyperfine interactions, but for both, Fe m and core, the hyperfine interactions are neglected.Scan, G, B m , and R are the same for both Fe m and core, so they cancel out.Thus, equation 3 becomes and follows this expression is solely dependent on the total spin of the ferritin core signal.As an example, this ratio is 0.8 and 6.8 if the spin number S of the core is 100 and 300, respectively.
In order to compare the calculated EPR active spin concentration ratio [Fem] [core] to the number of paramagnetic Fe atoms, normalized by the number of superparamagnetic atoms, obtained by magnetometry (N p /N Fh ) we need to multiply [core] with the number of iron atoms contributing to the ferrihydrite-like ferritin-core-EPR signal (N fh ).The amount of iron determined by ICP-MS (N ICP ) includes both mononuclear (N Fem ) and ferritin-like iron ions (N fh ).However, the number of mononuclear Fe-ions that contribute to their EPR signal is small (2 nd integral: 4575) compared to the ones contributing to the ferritin EPR signal (6778579).Therefore, we approximate the value obtained by ICP-MS to the number of iron atoms that contribute to the ferritin EPR signal, i.e., N fh = 1967±78.The error of [Fem] [core]•N fh was calculated using the following error propagation expression where F is [Fem] [core]•N fh , M 1 is the second integral of the ferritin-core EPR signal, M 2 is the second integral of the mononuclear Fe EPR signal, M 3 is the number of iron atoms in the ferritin core N fh .The errors of M 1 and M 2 were estimated by considering mainly the associated error of the simulations, using the Montecarlo error propagation with a small sampling number.Thus, δM 1 is 310 and δM 2 is 41.Finally, the error of M 3 is 78.
Table S2 shows the [Fem] [core]•N fh estimated for both S core of 100 and 300.The largest number (3.6 • 10 −3 ) is obtained when considering S core = 300.This value is one order of magnitude smaller than the fraction N p /N Fh , found by magnetometry (0.0645).The relative errors in Table S2 round up to 10%.The difference in [Fem] [core]•N fh for the two different S of a factor of 10 suggests that the uncertainty of [Fem] [core]•N fh is the dominating source of uncertainty and could be in the range of an order of magnitude.
Table S2: EPR active spin concentration ratio per ferritin particle [Fem] [core]•N fh Table S3: Data used in the error propagation calculation of [Fem] [core]

EPR Model Simulations
In this section, some examples of the large set of spectra that constitute the EPR Model Simulations (80 K), are provided.The Model Simulations are used to constrain certain EPR parameters based on the trends seen.Figure S13 shows simulations of one EPR component with S = 10 by varying the D parameter from 100 to 5000 MHz and considering two different Gaussian broadening parameters (3000 and 8000 MHz).For the smallest broadening, D values larger than 500 MHz can not be used to account for the features observed in our experimental spectrum, this is shown in Fig. S13a where a sharp feature at low field is visible for the green curve (1000 MHz) and in Fig. S13c where the spectra become distorted.
For the largest broadening, D cannot be greater than 1000 MHz (see Fig. S13d) since the spectrum starts to get progressively sharp at the lowest fields and flat at higher fields.In this way, to simulate the experimental spectrum measured at 80 K, figures S13a,b are used as starting points to inspect all parameters for D, H strain and weights considering these two base components

EPR alternative fitting approach
Following a different approach, a systematical fit of the spectra recorded at different temperatures was performed.First, We determined the S, D, and weight values of each component by fitting the EPR spectra recorded at 80 K (close to the middle-temperature range).The fitting was performed with a Matlab (R2019a) script, using the EasySpin package (5.2.4), using the scaling approach described in the main text.
The optimized scaled S and D values found for both EPR components, g and their respective weights are given in Table S4.Then, to fit the spectra from 20 to 190 K we used the scaled S, D and weight of each component determined at 80 K and only the Gaussian broadening given by an H strain parameter was used as the fitting parameter for both components (see Fig. S16).The spectra between 5 and 15 K were omitted from the analysis because it was not possible to accurately discriminate the broad signal from the background noise.S4.The scaled H strain of both systems (fitting parameters) are depicted in figure S16.
The two individual components are shown in figure S15 at 20 K, 80 K, and 190 K. Component 2 (blue curve) is narrower and highly contributes to the total intensity of the spectra.Comparing the lineshape evolution at the different temperatures, it can be noticed that component 2 is responsible for the overall lineshape change.Whereas, component 1, being quite broad and lower in intensity, appears almost temperature independent.This is a sign of a poor fitting result since component 1 can be hiding many features.D parameter does not change the simulation shape.The H strain changes the lineshape, but in an unsystematic way.In spite of the good agreement of simulations with the experiment, simulation parameters do not give physical insight.Therefore, a different approach was used (see main text).

Other possible approaches to fitting EPR lineshapes and critical evaluation of the model used
Considering the quantum nature of small nano-particles, a quantum-mechanical model for the ferritin core is the method of choice.The Giant Spin model makes this computationally feasible and was shown to be valid also for large spin systems as the ferritin core. 35][6][7] These models are computationally less demanding, and, for example, distributions of magnetic parameters could be introduced more easily than in the Giant Spin model (see below).Usselman et al. 1 showed that the application of these models to the EPR spectra of magnetite nanoparticles lead to inconsistencies, either resulting in unphysical magnetic parameters or did not fit the lineshape in the low-temperature regime, making quantum mechanical models the more attractive choice.
To implement the distribution of magnetic parameters, suggested by the analysis of the magnetometry data, in the Giant Spin model, however, is not straightforward.In the present study, we use two EPR components, and a significant Gaussian broadening of each of these components to represent the distribution.This approach leads to good simulations of the experimental spectra, and, as shown in Table 2 (main text), overall the parameters derived from of this interpretation agree well with the magnetometry data.Yet, it cannot be excluded that other representations of magnetic properties could also work, given the broadness of the experimental spectra.

Equilibrium magnetization models
The evaluation of eq. 4 requires numerical integration.Standard numerical integration routines require excessive computation time to evaluate the five-fold integral, while Monte Carlo integration does not work well, due to the different evaluation points at the numerator and denominator.Therefore, integrals have been replaced with sums over discrete states for a finite number of particles, that is where δϕ = δλ = δψ = δθ = π/2N are exact dividers of π/2 and δε an exact divider of the maximum canting angle ε max .The maximum canting angle is defined as the mean root square of ϵ weighted by e −E/k B T and avoids the evaluation of strongly canted states with negligible occurrence.As ϵ increases with the applied field, ε max is evaluated iteratively starting for a series of field values starting from ε max = πB/2B E .The number N + 1 of samples in each dimension is chosen to be the largest possible for given memory and time constraints.The numerical procedure has been implemented in Wolfram Mathematica and compiled in C using the built-in compiler.A PC with 64 GB RAM supports calculations up to N = 27 for a total of 10 9 combinations of particle orientations and magnetic states, for an angular resolution of ∼ 3.3 • in ϕ, λ, ψ, and θ.

Spontaneous spin canting
Spontaneous canting is generated by an additional energy term proportional to c • (u a × u b ), with c being the unit vector perpendicular to the canting plane.In this case, the total magnetic energy per unit of volume becomes where the term E a accounts for the combined effect of surface and volume anisotropy.In the isotropic case (E a = 0), the total energy is insensitive to the direction n of canting, as long as it is perpendicular to c. Therefore, E a = 0 describes an azimuthal dependence of E with respect to c, which requires an additional degree of freedom with respect to the case of a defect moment.Given the minor role played by realistic values of the magnetocrystalline anisotropy in the defect moment simulations of Fig. 4, numerical calculations have been performed with E a = 0.

Figure S1 :
Figure S1: Transmission electron microscopy (TEM) characterization of human-liver ferritin.(a) TEM image of human-liver ferritin cores.(b) Stained-SEM image of human-liver ferritin particles allowing the protein shell visualization.The images were acquired with the JEOL JEM-1400 series 120 kV Transmission Electron Microscope.

Figure S3 :
Figure S3: Purity assessment of human-liver ferritin.(a,b) Dynamic light scattering showing the particle size (hydrodynamic diameter) distribution: one main Gaussian peak with µ=12.2 nm and σ=4.5 nm.(c) Non-denaturing polyacrylamide gel electrophoresis of human-liver ferritin (HSF) stained with Coomassie.The strong band of 500 kDa agrees with ferritin monomers (24-mer protein) while the weaker bands of 1.048 and 1.236 kDa correspond to ferritin multimers (dimers/trimers) or aggregates of ferritin protein.Markers refer to NativeMark Unstained Protein Standard (Cat.No. LC0725).

Figure S4 :
Figure S4: Automatic baseline correction algorithm.(a) The experimental spectrum (in blue), the optimal line (red dotted line), and the transformed spectrum (yellow curve).(b) The first integral of the baseline-corrected spectrum.Its high field values reach zero.(c) The second integral of the baseline-corrected spectrum.The high field values reach a plateau-like region.

Figure
Figure S5 shows an example of one of the best fits obtained by considering a simulation with only one EPR component.It shows the discrepancy between the data and the model, indicating the need of two or more EPR components.

Figure S5 :
Figure S5: Simulation (light-blue line) of EPR spectrum (black line) with single component only (80 K).The simulation is performed with one spin component of S =15, g = 2.01, D = −335 MHz, and H strain = 6500 MHz.

Figure S6 :
Figure S6: Simulations (red dotted lines) of EPR spectra (black lines) for all data recorded between 20 and 210 K.The simulations are performed with two spin systems, as described in the main text.The system parameters are shown in TableS1.

Figure S7 :
Figure S7: Total EPR simulations (red lines) performed at 30 K, 80 K, 130 K and 190 K.The simulations are performed with two spin systems, as described in the main text.The first system component is depicted in green.The second component is shown in blue.

Figure S9 :
Figure S9: Simulations (colored lines) of EPR spectra (black lines) recorded at 20 K, 30 K, 80 K, and 130 K.The simulations are performed with two spin systems using different scaling factors (n).The remaining parameters are shown in TableS1.For the entire temperature range, the line shape presents negligible changes, confirming the scaling method as a suitable candidate to reproduce and study the EPR line shape of high-spin systems.

Figure S10 :
Figure S10: Simulations of EPR spectra performed at 80 K.The original simulation is shown in red.The simulations obtained by considering both components with equal D, The yellow curve is obtained by simulating both components with equal D.

Figure S11 :
Figure S11: Simulations (colored lines) of EPR spectra (black lines) at 80 K performed by varying D2.An example of what is considered to be an acceptable simulation is depicted in red.Examples of rejected simulations are shown in light-blue and yellow.

Figure S14 :
Figure S14: Simulations (red dash lines) of EPR spectra (black lines) using the 'alternative fitting model' (see text) for EPR spectra between 20 and 190 K.The simulations are performed with the automatic fitting routine with two spin systems.The system parameters are shown in TableS4.The scaled H strain of both systems (fitting parameters) are depicted in figure S16.

Figure S15 :Figure S16 :
Figure S15: Total simulations (red lines) of EPR spectra (black lines) performed at 20 K, 80 K, and 190 K.The simulations are performed with two spin systems, as described above.The first system component is depicted in green.The second component is shown in blue.

Figure S17 :
Figure S17: (a) Initial magnetization curves measured at temperatures comprised between 3 K (blue circles) and 200 K (yellow circles), and global fit to the modified Langevin model with one superparamagnetic component without moment distribution.(b) Fit residuals as a function of the predicted magnetization.

Figure S18 :Figure S19 :Figure S20 :
Figure S18: Real µ ′ and Imaginary µ ′′ parts of the magnetic moment of the human-liver ferritin (raw data).The legend refers to the frequency of the oscillating B 1 field.

Figure S21 :Figure S22 :Figure S23 :
Figure S21: Comparison between a direct estimate of the energy barrier distribution, obtained from the blocking temperature distribution of the IRM Component 1, and the energy barrier distribution expected from E b = mB a /2 with the incorrect assumption that m and B a are independent variables.

E M 0
= B E u a • u b + B E tan(2ε s ) c • (u a × u b ) − B(u a + u b ) +E a (u a , u b , e),

Table S1 :
System parameters used for simulations of the EPR spectra.Simulation parameters: S = 10, D: axial zero-field splitting, H strain : Gaussian broadening parameter.The parameters at 30 K, 80 K, 130 K, and 190 K were determined first.