Simulation of Colloidal Stability and Aggregation Tendency of Magnetic Nanoflowers in Biofluids

Abstract

A population balance model for the aggregation of iron oxide nanoflowers (IONfs) is presented. The model is based on the fixed pivot technique and is validated successfully for four kinds of aggregation kernels. The extended Derjaguin, Landau, Verwey, and Overbeek (xDLVO) theory is also employed for assessing the collision efficiency of the particles, which is pertinent to the total energy of the interaction. Colloidal stability experiments were conducted on IONfs for two dispersant cases—aqueous phosphate buffered saline solution (PBS) and simulated body fluid (SBF). Dynamic light scattering (DLS) measurements after 24-h of incubation show a significant size increase in plain PBS, whereas the presence of proteins in SBF prevents aggregation by protein corona formation on the IONfs. Subsequent simulations tend to overpredict the aggregation rate, and this can be attributed to the flower-like shape of IONfs, thus allowing patchiness on the surface of the particles that promotes an uneven energy potential and aggregation hindering. In silico parametric study on the effects of the ionic strength shows a prominent dependency of the aggregation rate on the salinity of the dispersant underlying the effect of repulsion forces, which are almost absent in the PBS case, promoting aggregation. In addition, the parametric study on the van der Waals potential energy effect—within common Hamaker-constant values for iron oxides—shows that this is almost absent for high salinity dispersants, whereas low salinity gives a wide range of results, thus underlying the high sensitivity of the model on the potential energy parameters.

1. Introduction

The stability of colloids is a very important aspect that may affect pertinent applications of the synthetic nanoparticles, and therefore an understanding of the aggregation mechanisms and subsequent time-frames of a specific suspension is crucial. In particular, inorganic nanoparticles exhibit a tendency towards aggregation when dispersed in saline solutions. In addition, with the ion’s specific effects, there is a significant influence of charge, shape, and surface functional groups on the aggregation of different colloidal particles [1,2,3]. Although experimental studies are usually used for defining the time-dependent particle size distribution (PSD) of nanoparticles in a particular dispersant, in silico methods are emerging as an easy-to-use and cost-effective alternative. Based on the Smoluchowski’s coagulation equation that describes the population balance of the particles in conjunction with the process and kinetics of aggregation, several methods have been developed for numerically simulating these mechanisms by discretizing the size distribution range into consecutive sections and solving the aforementioned equations successively for each section. An extensively used method is the one developed by Hounslow et al. [4], which discretizes the range into consecutive volume sections, each of which is double in size in terms of volume compared to the previous. The equations are expressed in such a way so as to ensure the correct prediction of the total number of particles and, in addition, a correction parameter was introduced so as to ensure conservation of mass. The limiting factor of having sections of predefined volume was overcome by Kumar and Ramkrishna [5], who introduced a method of discretizing the volume range in sections of an arbitrary size, thus providing the flexibility to focus on the behavior and outcome at certain volumes of the spectrum, which was also capable of correctly evaluating two different moments of the distribution. An extensive review of the solution methods of the population balance equations (PBE) is given by Vanni [6].

There has since been a wide utilization of the aforementioned methods for modeling aggregation in liquid solutions. Biggs et al. [7] applied the method developed by Hounslow et al. [4] for modelling activated sludge flocculation, providing a good approximation of the change in mean floc size with time. The flocculation of colloidal suspensions by salts and polymers was studied by Runkana et al. [8], also by applying the method by Hounslow et al. [4], this time in a more sophisticated model encompassing the Derjaguin, Landau, Verwey, and Overbeek (DLVO) theory. The model predictions are in close agreement with the experimental results for the flocculation of colloidal hematite suspensions in the presence of KCl and polyacrylic acid at different concentrations. Atmuri et al. [9] performed a combined experimental and modelling study to explore aggregation in concentrated attractive colloidal suspensions using a PBE model in conjunction with the DLVO theory for describing the interparticle potential. Although the model is reported to predict much faster aggregation rates than those observed experimentally, the authors argue that it is the first PBE model that can successfully predict stable aggregate size and aggregation regimes of charged colloidal particles over a range of salt concentrations and particle concentrations. Liang et al. [10] used the Kumar and Ramkrishna [5] method, also taking into account the DLVO theory, for studying the aggregation behavior of submicron-sized particles of praseodymium-doped zirconium silicate, a ceramic pigment, in an aqueous suspension. The model was modified to predict the stable state of the aggregation by introducing the volume mean size of the aggregate to the stability ratio. In this study, experiments were also performed for investigating the aggregation of the particles in an aqueous suspension in the presence of sodium dodecyl benzene sulfonate or potassium chloride, and the authors report that the predicted data were in reasonable agreement with the experimental results, with an exception regarding predictions for certain levels of suspension concentrations. A similar study employing the same model was performed by Huang et al. [11], who studied the dispersion/aggregation behavior of colloidal particles in the suspension under different conditions by a modified population balance model (PBM). The colloidal system consisted of a dispersion of praseodymium-doped zirconium silicate pigment in an aqueous suspension by modified hydroxyl copolymer and an experimental study was also performed, with the authors ultimately reporting that predicted data are similar to the experimental results.

Iron oxide nanoparticles (IONps) are a prominent class of nanostructures, with impactful medical applications owing to their intrinsic magnetic properties. IONps lay on the superparamagnetic regime, allowing them to be used as contrast agents in magnetic resonance imaging (MRI), as well as nanoheaters for the magnetic hyperthermia treatment of cancer [12]. The shape and size of IONps dictates their colloidal stability, magnetic properties, and the corresponding biomedical application [13,14]. Several synthetic methodologies have been reported for the preparation of IONps, the most common of which is the “co-precipitation method” due to its high yield and low cost. The main disadvantages of this method are increased polydispersity of the final ferrofluid, sedimentation after prolonged storage, and aggregation effects. To overcome these limitations, a shift of interest has been noted towards the “polyol method”, where IONps are prepared in hydrophylic polyol solvents that act both as reducing and capping agents [15]. Through this method, IONps with narrow size distributions, controllable shape, and exceptional colloidal stability in aqueous solutions can be formulated. The formation of citrate coated high crystallinity iron oxide nanoflowers (IONfs) can be achieved via polyol synthesis by controlling the ratio on the polyols used, as reported by Hugounenq et al. [16]. The multicore IONf shape yields superior magnetic properties than the plain spherical nanoparticles, rendering them suitable for biomedical applications [17]. The aim of the current study is three-fold: (i) to carry out case-specific experiments on the colloidal stability of IONfs and specify the effects of using citrate coating, (ii) to develop an in silico model for simulating colloidal stability taking into account the xDLVO theory, and (iii) to study by simulation the effects on the colloidal stability of parameters such as ionic strength and van der Waals forces.

2. Materials and Methods

2.1. Numerical Model

The population balance equation (PBE) [18] is a basic tool in the field of colloidal science for the description of interparticle behavior and for the subsequent prediction of the evolution particle aggregation and consequent size distribution with time. Based on this, a discretized form of the PBE for aggregation of particles is given as [10,11]:

$$\frac{d{N}_i}{dt}=-N_i{\displaystyle \sum }_{j=1}^{\infty }{\beta }_{i,j}{N}_j+\frac{1}{2}{\displaystyle \sum }_{j=1}^{i-1}{\beta }_{i-j,j}{N}_{i-j}{N}_{j }$$

(1)

where i represents a particle-size section; $N_i$ is the number concentration of particles in section $i$; and ${\beta }_{i,j}$ is the aggregation kernel, which is a combination of collision frequency k and collision efficiency λ, and is expressed as

$${\beta }_{i,j}=k_{i,j}·{\lambda }_{i,j}$$

(2)

The first term on the right-hand side of (1) describes the loss of particles of size i because of their aggregation with other particles, whereas the second term represents the addition of particles of size i by their formation due to the aggregation of smaller particles.

The part of the aggregation kernel that corresponds to collision frequency is pertinent to the Brownian motion of the particles, i.e., perikinetic aggregation, and is given by the equation developed by Smolukowski [18], assuming the collision of two particles of diameters$d_i\mathrm{and} d_j$:

$$k_{i,j }=\frac{2k_{B}T}{3\mu }\frac{{\left(r_i+r_j\right)}^2}{r_i{r}_j}$$

(3)

where $k_B$ is Boltzmann’s constant, $T$ is the absolute temperature, and $\mu$ is the fluid viscosity.

The part of the aggregation kernel that corresponds to collision efficiency is pertinent to the total energy of interaction, which is a function of the forces acting between surfaces, and includes van der Waals force, electrostatics, short range repulsion by hydration, long range hydrophobic attraction, magnetic forces, steric forces or bridging due to adsorbed polymer, etc. [19]. This term is generally assumed to be the inverse of the Fuchs’ stability ratio, W, and is given by the following equation [20]:

$$\frac{1}{{\lambda }_{i,j}}=W_{i,j}=\left(r_i+r_j\right){{\displaystyle \int }}_{r_i+r_j}^{\infty }\frac{exp\left(V_{tot}/k_{B}T\right)}{R^2}dR$$

(4)

where $V_{tot}$ is the total interaction energy and R is the center-to-center distance between the particles. Based on the classic and extended DLVO theory, the total interaction energy concerning the present case is given by the following:

$$V_{tot}=V_{vdW}+V_{edl}+V_{mag}$$

(5)

where $V_{vdW}$ is the energy due to van der Waals forces, $V_{edl}$ is due to the electrical double layer, and $V_{mag}$ is due to the magnetic forces. In the Hamaker theory, the van der Waals energy of attraction for two spheres is given by the following:

$$V_{vdW}\left(R\right)=-\frac{A}{6}\left\{\frac{2r_i{r}_j}{R^2-{\left(r_i+r_j\right)}^2}+\frac{2r_i{r}_j}{R^2-{\left(r_i-r_j\right)}^2}+\ln \left[\frac{R^2-{\left(r_i+r_j\right)}^2}{R^2-{\left(r_i-r_j\right)}^2}\right]\right\}$$

(6)

where A represents the Hamaker constant [10]. For estimating the electrical double layer repulsion $V_{edl}$ between two spherical particles having surface potentials ${\phi }_i$ and ${\phi }_j$, the following analytical expression used in the past [8,10] is employed:

$$V_{edl}=64\pi {\epsilon }_{\Gamma }{\epsilon }_0{\left(\frac{k_{B}T}{z_{C}e}\right)}^2\left(\frac{r_i{r}_j}{r_i+r_j}\right)\tanh \left(\frac{z_{C}e{\phi }_i}{4k_{B}T}\right)\tanh \left(\frac{z_{C}e{\phi }_j}{4k_{B}T}\right)exp\left(-\kappa h_0\right)$$

(7)

where $e$ is the elementary charge; $z_C$ is the valence of the counter ions; ${\epsilon }_0$ and ${\epsilon }_{\Gamma }$ are the dielectric constants of vacuum and the solvent, respectively; and $\phi$ represents the surface potential of the particles. The Debye−Hückel parameter $\kappa$ is a function of the salt concentration $c_i$, the valence $z_i$ of electrolyte ions, and temperature $T$:

$$\kappa ={\left(\frac{e^2{N}_A{{\displaystyle \sum }}_i{c}_i{z}_i^2}{{\epsilon }_0{\epsilon }_{\Gamma }{k}_{B}T}\right)}^{1/2}$$

(8)

where $N_A$ is the Avogadro number. The magnetic energy between two particles that behave as magnetic dipoles and are oriented in a head-to-tail configuration is given by the following [21]:

$$V_{mag}=-\frac{8\pi {\mu }_0\left(M_i{r}_i^3\right)\left(M_j{r}_j^3\right)}{9R^3}$$

(9)

where M is the magnetization of the material and ${\mu }_0$ is the permeability of the vacuum.

According to the method by Kumar and Ramkrishna [5], the size range can be divided into consecutive sections of arbitrary volume sizes (also referred to as size bins) from x₁ to x_K of arbitrary size (Figure 1), and (1) thereby becomes:

$$\frac{d{N}_i}{dt}=-N_i{\displaystyle \sum }_{k=1}^K{\beta }_{i,k}{N}_{k }+{\displaystyle \sum }_{x_{i-1}\le \left(x_j+x_k\right)\le x_{i+1}}^{j\ge k}\left(1-\frac{1}{2}{\delta }_{j,k}\right)\omega {\beta }_{j,k}{N}_j{N}_k$$

(10)

where ${\delta }_{j,k}$ is the Kronecker delta. The method indicates that when the collision two of particles of volume x_j and x_k form a particle of volume υ = $x_j+x_k$ between preset grid volumes x_i-1 and x_i (Figure 1), then this is arithmetically split between these neighboring grid volumes. For the preservation of numbers and mass within the whole spectrum, this split is defined by $\omega$ in (10) as follows:

$$\omega =\left\{\begin{array}{c}\frac{x_{i+1}-\upsilon }{x_{i+1}-x_i}, x_i\le \left(x_j+x_k\right)\le x_{i+1}\\ \frac{\upsilon -x_{i-1}}{x_i-x_{i-1}}, x_{i-1}\le \left(x_j+x_k\right)\le x_i\end{array}\right.$$

(11)

For the numerical solution of (10), a FORTRAN code is developed with exploitation of the ODEPACK solver [22].

2.2. Experimental Procedure

All of the reagents were purchased from Sigma-Aldrich Co. (St. Louis, MO, USA) and the solvents were purchased from Fischer Scientific (Waltham, MA, USA, HPLC grade). Fetal Bovine serum (FBS) was obtained from PAN Biotech GmbH (Aidenbach, Germany). All of the other chemicals and materials were commercially available and of high grade, and are described below. Citrate coated IONfs were synthesized through the polyol method according to Hugounenq et al. [16]. The final ferrofluid was an aqueous colloidal solution of 15 nm IONfs, measured by transmission electron microscopy (TEM) (Figure 2), with a maghemite crystal structure, assessed through X-ray diffraction spectroscopy (XRD). The volume magnetization (M) of the IONfs was determined at 0.735 A/m via SQUID-based magnetometry. To assess their colloidal stability in terms of aggregation, size measurements were conducted in a Malvern Zetasizer Nano ZS dynamic light scattering (DLS) instrument.

IONfs were dispersed in an aqueous phosphate buffered saline solution (PBS) at a 1000μM concentration of iron content. Separate measurements for the same iron concentration were conducted in PBS supplemented with 10% v/v fetal bovine serum (FBS) to use as the simulated body fluid (SBF). FBS is a biofluid containing over 3700 different proteins, among them albumin represents ~60% and globulins ~40%, while the rest is ~1% [23]. The samples were first treated at 37 °C (mean physiological temperature) and measured in 15 consecutive cycles with 30 s duration and 30 s time gap. Two time points were selected for the measurement, the first (t₀) was within 3 min after mixing the IONfs with PBS or SBFand the second (t₂₄) after 24-h of incubation at 37 °C.

3. Results and Discussion

3.1. Model Validation

Validation of the model was carried out based on analytical expressions for three kinds of collision frequency kernels and where the collision-efficiency effect was not considered, i.e., λ = 1. For this purpose, the PBE equation was expressed in non-dimensionalised form, as follows:

$$\frac{d{\tilde{N}}_i}{d\tau }=-{\tilde{N}}_i{\displaystyle \sum }_{j=1}^{\infty }{\tilde{\beta }}_{i,j}{\tilde{N}}_j+\frac{1}{2}{\displaystyle \sum }_{j=1}^{i-1}{\tilde{\beta }}_{i-j,j}{\tilde{N}}_{i-j}{\tilde{N}}_{j }$$

(12)

where ${\tilde{N}}_i=N_i/N_0$, ${\tilde{\beta }}_{i,j}={\beta }_{i,j}/{\beta }_0^{\prime }$, $\tau ={\beta }_0^{\prime }{N}_{0}t$ and $N_0$ is initial number of particles, i.e., the number of particles at t=0. In addition, ${\upsilon }_0$ is the volume used for non-dimensionalisation of the particle volumes and, in this case, it is assumed as the median of the size bins. The discretization was carried out assuming $x_{i+1}=2x_{i }$and 21 size bins. The benchmark analytical solutions from the literature are provided for the number−density function n, the relation of which with the number concentration function is given in discretized form as follows:

$$N_i={{\displaystyle \int }}_{\left(x_i+x_{i-1}\right)/2}^{\left(x_{i+1}+x_i\right)/2}{n}_{i}dx$$

(13)

and its non-dimensionalised form is ${\tilde{n}}_i=n_i/\left(N_0/{\upsilon }_0\right)$.

3.1.1. Constant Kernel

This kernel is given as ${\beta }_{i,j}={\beta }_0$ and therefore, ${\beta }_0^{\prime }={\beta }_0$ and ${\tilde{\beta }}_{i,j}=1$. The analytical solution for the initial condition

$${\tilde{n}}_0=e^{-\tilde{\upsilon }}$$

(14)

where $\tilde{\upsilon }=\upsilon /{\upsilon }_0$, is given by Gelbard and Seinfeld [24]:

$$\tilde{n}=\frac{4}{{\left(\tau +2\right)}^2}exp\left(-\frac{2\tilde{\upsilon }}{\tau +2}\right)$$

(15)

Using (13), we obtain the benchmark solution for the number concentration, and the comparison with the numerical solution is shown in Figure 3a.

3.1.2. Linear Kernel

This kernel is given as ${\beta }_{i,j}={\beta }_1\left(x_i+x_j\right)$, and therefore ${\beta }_0^{\prime }={\beta }_1{\upsilon }_0$ and ${\tilde{\beta }}_{i,j}={\tilde{x}}_i+{\tilde{x}}_j$ where ${\tilde{x}}_i=x_i/{\upsilon }_0$. The initial condition is again taken as in (14), the analytical solution of which is given by Gelbard and Seinfeld [24]:

$$\tilde{n}=\frac{1-{\rm T}}{\tilde{\upsilon }\sqrt{{\rm T}}}exp\left(-\left(1+{\rm T}\right)\tilde{\upsilon }\right){{\rm I}}_1\left(2\tilde{\upsilon }\sqrt{{\rm T}}\right)$$

(16)

where ${\rm T}=1-exp\left(-\tau \right)$ and I₁ is the modified Bessel function of first class and order one. Using (13), we obtain the benchmark solution for the number concentration, and the comparison with the numerical solution is shown in Figure 3b.

3.1.3. Quadratic Kernel

This kernel is given as ${\beta }_{i,j}={\beta }_2{x}_i{x}_j$, therefore ${\beta }_0^{\prime }={\beta }_2{\upsilon }_0^2$ and ${\tilde{\beta }}_{i,j}={\tilde{x}}_i{\tilde{x}}_j$. For this case, the initial condition ${\tilde{n}}_0=exp\left(-\tilde{\upsilon }\right)/\tilde{\upsilon }$ is used, the analytical solution of which is given by Ernst et al. [25]:

$$\tilde{n}=\frac{1}{{\tilde{\upsilon }}^2\sqrt{\tau }}exp\left(-\left(1+{\rm T}\right)\tilde{\upsilon }\right){{\rm I}}_1\left(2\tilde{\upsilon }\sqrt{\tau }\right)$$

(17)

where

$$1+{\rm T}=\left\{\begin{array}{c}1+\tau , \tau \le 1\\ 2\sqrt{\tau }, \tau >1\end{array}\right.$$

and I₁ is the modified Bessel function of first class and order one. It is important to point out that for this kernel, the mass is not conserved for $\tau >1$ due a loss of mass from finite size clusters (sol particles) to the infinite cluster (gel or super particle) [25]. Using (13), we obtain the benchmark solution for the number concentration, and the comparison with the numerical solution is shown in Figure 3c.

3.1.4. Brownian Kernel

Taking into account Equation (3), this kernel can be expressed as ${\beta }_{i,j}={\beta }_3{\left(d_i+d_j\right)}^2/d_i{d}_j$, therefore ${\beta }_0^{\prime }={\beta }_3$ and ${\tilde{\beta }}_{i,j}={\left({\tilde{d}}_i+{\tilde{d}}_j\right)}^2/{\tilde{d}}_i{\tilde{d}}_j$, where ${\tilde{d}}_i=\sqrt[3]{{\tilde{x}}_i}$. The initial condition employed is given in (14) and the results obtained when solving (13) are given in Figure 3d and are compared with past results [26].

It can be seen that there is a very good agreement between the numerical and analytical solutions for all four cases for all of the time instances considered. This yields even for $\tau =2$ in the case of the quadratic kernel, meaning that the model is also capable of taking into account the loss of mass.

3.2. Aggregation Experiments

The IONfs in an aqueous dilution right after synthesis demonstrated a hydrodynamic radius (by number) of 16.56 ± 4.53 nm, whereas after one year, no significant alterations were observed (17.75 ± 4.70 nm), indicating a colloidally stable solution free of sedimentation or aggregation. The storage condition was refrigeration at 2–4 °C.

The measurements corresponding to simulated body fluids in terms of salinity (PBS) and the presence of proteins (SBF) are discussed for t₀ and t₂₄. The average size value by number of 15 consecutive measurements of IONfs in PBS at t₀ and t₂₄ are 16.22 ± 8.034 nm and 79.60 ± 19.25 nm, respectively. A significant increase in size was reported after 24-h incubation in PBS, dictated by the presence of salts in the colloidal solution (Figure 4a). For the sample of IONfs dispersed in SBF at t₀, the average size value by number is 25.80 ± 9.01 nm, whereas at t₂₄ it is at 35.62 ± 11.71 nm. It is worth mentioning that the presence of FBS adds a second peak in these measurements, which was determined at 5.22 ± 1.43 nm, by measuring the plain SBF without IONfs. Again, there is an apparent increase in size that corresponds to the formation of protein corona (PC), but not that substantial to indicate aggregation. This can be confirmed by the ζ-potential decrease before and after treatment with 10% FBS in PBS from −27 mV to −19 mV. The latter can be attributed to the attractive interactions induced between the carboxylic acid of IONfs and the protein’s functional groups present in FBS (e.g., NH₂). PC forms when nanoparticles are in close contact with proteins, such as albumins, which adsorb on their surface, further stabilizing them in the medium [27]. According to our data, IONfs seem to aggregate more in plain PBS than in SBF; this behavior can be attributed to the formation of PC, which averts the adsorption of the salts due to the higher affinity for the proteins rather than the salt counter ions [28]. In low nanoparticle concentrations, like the one used in our research, aggregation phenomena stemming from the PC formation are minimal (Figure 4b). FBS is a biological fluid that contains more than 3700 different proteins, with 12 orders of magnitude differences in their relative concentrations. Among all of the proteins in the serum, albumin represents the most abundant fraction and accounts for ~60% mass of the total proteins in the solution (1 mM), with globulins being the second most abundant (~40%) and all the rest representing less than 1% [23]. The observation of aggregation blocking due to PC is consistent with previous studies [29].

3.3. Simulations

In the present case and according to the discrete measurements by the dynamic light scattering (DLS) analyzer, the computational grid included 28 size bins, i.e., K = 28 in (10), corresponding exactly to the aforementioned discretization obtained from the experiments. The values of the parameters involved in the equations are given in

Table 1. Parameters used in the simulations.

Parameter	Value
Μ 1	0.6864 cP
T 2	310 K
zc 3	1
εΓ 4	78.5
Φ 5	−0.027 V
M 6	0.7346 A/m
A 7	3.1 × 10−20 J
PBS 8	KCl 2.68 mM KH2PO4 1.47 mM NaCl 136.89 mM K2HPO4 8.10 mM

¹ viscosity of solvent; ² temperature; ³ valence of counter ion; ⁴ dielectric constant of the solvent; ⁵ surface potential (approx. zeta potential); ⁶ intensity of magnetization; ⁷ Hamaker constant; ⁸ phosphate buffered saline solution.

.

The PBS results highlight the progress regarding the percentage number concentration ${N_i}^{\prime }$, where ${N_i}^{\prime }=N_i/N_0^t\%$ with $N_0^t$ showing the total number of particles at the current time instant. Simulations show a fair agreement with the experimental results (Figure 5a). The observed deviation from the experiments can be attributed to the uncertainty regarding the Hamaker constant, but also to the balance between the potential energies, as it seems that when the ionic strength is low, then the der Waals forces are dominant and can affect the aggregation mechanism. With respect to the physics of aggregation, a higher aggregation rate in simulations compared to experiments has also been reported in the past [9]. The effect of patchiness of the particles can be significant in the case of nanoflowers, as these particles deviate markedly from perfect spheres and therefore there can be important differences in the surface potential, which cause the particles to have a reversible aggregation, i.e., they also undergo breakup. In addition, the fact that the particles do not coalesce but retain their shape forming fractal aggregates, the shape of which is far from a perfect sphere, may therefore cause the corresponding kernel to not exactly represent this differentiated phenomenology of aggregates. Therefore, this effect is difficult to be accounted in the PBE model. Next, the magnitude of potentials is studied. It should be noted that in the experiments, no magnetic field was applied, so there is no magnetic potential. However, this is taken into account in the simulations, so as to study its effect when a magnetic field is assumed and to confirm its negligibility [30]. The magnitude of the potentials as a function of the dimensionless distance between the particle surfaces is shown in Figure 5b. It can be seen that the magnetic potential is the least effective of the potentials considered, and therefore has the lowest contribution in the total potential energy even in its most prominent effect, which is the head-to-tail configuration of particles, as assumed in (9), and which would anyway not be possible in the present case where IONfs aggregate in near spherical clusters. In view of this, the effects of the ionic strength and van der Waals forces were further studied.

3.3.1. Effects of Ionic Strength

Various cases of ionic strength are considered, namely the aforementioned PBS case and cases that are usually used in experiments [31] involving 1:1 electrolyte (e.g., NaCl or KCl) with a molarity of 1 mM and 10 mM, whereas a case of no electric double layer (EDL) potential energy was also considered. The van der Walls potential energy corresponds to A = 3.1 × 10⁻²⁰ J. It can be seen that, as expected, the introduction of dispersants of a higher molarity weakens the effects of EDL, as the Debye length increases and consequently imposes exponential decreasing of the EDL potential. Therefore, attractive forces (Van der Waals) remain dominant, as in the case of PBS, where the total potential energy remains negative for all distance values (Figure 6b). Consequently, the aggregation rate is higher in the cases of a higher molarity, and the PBS case is almost identical to having no EDL potential, as also observed in the past [32]. In addition, molarities as low as 1 mM induce lower aggregation rates than in the experiment (Figure 6a).

3.3.2. Effects of van der Waals Forces

Van der Waals potential energy effect on the aggregation rate is calculated for two cases of ionic strengths. The Hamaker constant value was used as a parameter and the range of values considered were obtained from the literature. These values can be as low as n 2 × 10⁻²⁰ J for IONps [33], and reach values of up to 4.5 × 10⁻²⁰ J especially for maghemite NPs [34]. In the case of PBS, it seems that the changes within the range of Hamaker constants studied are marginal, i.e., there is no apparent change in the aggregation rate and the results after 24 h are almost identical (Figure 7a). In this case, the total potential energy is again dominated by van der Waals forces and shows no repulsive effects for any surface distance (Figure 7b).

Conversely, the total potential energy for lower molarities such as 10 mM solution induces dominant repulsive effects within the surface distance range (Figure 8b). In this case, there are differences in the effect of the Hamaker constant values on the aggregation rate within the studied range (Figure 8a), and this seems to be pertinent to the balance between repulsive and attractive potentials, which consequently makes the integral of the total potential with respect to the NP surface distance very sensible to variations. The rates are lower than the previous case and for A = 2, the result for a 24-h period is almost identical to the experimental observations.

4. Conclusions

IONfs were synthesized using the polyol method and their colloidal stability in terms of aggregation was assessed by DLS measurements in both a PBS and an SBF suspension. In the first case, a significant increase in the size is observed after 24-h incubation, whereas in the second case, the presence of proteins induces the formation of a protein corona in the IONfs, thereby preventing aggregation, as also confirmed by the DLS measurements. Simulations using a validated in silico model for solving the particle balance equation in conjunction with the xDLVO theory show a slight overprediction of the aggregation rate. This can be attributed to the patchiness on the NP surface or the fractal nature of aggregates that cannot be considered in the present model. Further study on the effects of the ionic strength shows a scarcity of repulsion forces with increasing salinity, thereby inducing a higher aggregation rate, whereas the PBS case is almost identical as having no EDL effects. In addition, the range of Hamaker constant values for IONs seem to have no effect on the aggregation rate in the case of the PBS, whereas for lower salinity, the experimental case can even be predicted exactly, thus underlying the high sensitivity of aggregation models on the potential energy parameters.