Multiscale Modelling of Milk Proteins Adsorption on Metallic Iron Surfaces

1. Introduction

Biological activity and biocompatibility of inorganic materials are controlled by interactions at the bionano interface – a nanoscale layer where the material gets in contact with biomolecules. Through biomolecule adsorption, change of conformation, and surface chemistry, the bionano interface plays a key role in medical applications and food processing. Understanding the underlying mechanism for these interactions may help to control bio film growing, fouling, or contamination by metabolites arsing from bacterial activity at the contact surfaces [1,2,3,4]. Therefore, a molecular level insight into the main contributions governing the adhesion of biomolecules and the structure of the bionano interface would be beneficial.

Among the materials used in industrial and medical applications one of the most common one is the stainless steel, predominately consisting of metallic iron alloyed with chromium and carbon. These materials are widely used due to a low cost, relative biological neutrality, and corrosion resistance at standard low heat conditions. However, this might be not true at the higher temperatures: e.g. milk fouling and milk contamination is a critical problem in the heating process for food industry [5]. Daily, about a billion liters of milk being being processed. This causes fouling of the surfaces of heat exchangers which in turn leads to decreased heat transfer efficiency and frequent pressure drop events. Efficient and fast cleaning of factory’s heat transfer equipment is a costly and challenging task. When done often it can result in reduced production volumes and economical losses [6]. When liquid milk is heated, several processes take place that result in the deposition of minerals and proteins on the surface of containers. Generally, a thin protein layer is formed on a material’s surface right after the liquid milk is brought in contact with a food processing unit [7]. Adsorbed proteins may contribute to the decrease of heating by forming an insulating layer between the heater and material, can affect milk pasteurization and sterilization, and may cause contamination [8]. Physicochemical characteristics of the heat exchanger surface as well as the conditions and biochemical composition of the milk determine the outcome of this interaction. [9,10].

Beside the food industry, biomolecular interactions at metallic surfaces are essential for implants and other medical devices. When a tissue makes contact with such surfaces, unfavorable reactions such as infections, fibrosis, thrombosis, and inflammation may result [11]. Recent research demonstrates that quickly after implantation a layer of proteins form interstitial fluids covers biomaterial implants. Thus, the nature of this protein layer and its properties play a key role in how cells responses. To improve the operation and safety of medical implants, it is imperative to develop materials that cause the necessary tissue reactions. In order to increase implant biocompatibility, the majority of research has concentrated on the engineering of surface features that can change the numbers and types of bound proteins as well as the conformation, orientation, or binding strength of the adsorbed proteins [11,12,13].

Finally, the interest to bionano interactions is also driven by the concerns of safety of nanoparticles (NP) for human and animal health. The NP toxicity correlates to chemical aggressiveness of the material and scales with its physicochemical properties such as surface area, charge, or reactivity [14]. Establishing the relationship between metal properties and the reaction with different biomolecules is therefore key to screening the materials for potential health risks. In practice, the safety assessment is often costly and time-consuming and includes animal studies. In silico modelling may help predicting the interactions of nanomaterials with living organisms and providing the required information in a humane and cost-effective way [15,16,17,18,19]. Statistical data-driven methods are used for this purpose, where sufficient data are available [20,21,22]. Recently, also physics-based models have addressed the bionano interface. In particular, mechanisms of formation of NP protein corona have been studied my multple labs [23,24,25,26]. It is expected that the composition and configuration of corona determine the biochemical reactivity and sensitivity of NPs as well as their cell uptake and systemic transfer [23,27]. Yet, to allow predictive modelling one needs more information about the interactions at the bionano interface and their relation to the material and protein properties.

In this work, we investigate interactions between milk proteins and three face-centered-cubic (fcc) slabs of zero-valent iron constituting a simplified model of the stainless steel. To accomplish this, we have selected six of the most abundant proteins from the major protein groups found in natural cow milk [28]. Our primary goal is to quantify the binding affinity of these proteins on iron surfaces by evaluating their adsorption energy at their different orientations in respect to the solid iron. To calculate these energies and to predict the composition of insulating protein layer on the metallic surfaces, we will invoke a three-stage multiscale computational method including all-atom [29,30] and coarse-grained (CG) United Atom (UA) [26,31] and Monte-Carlo [22] simulations. This method was previously applied to explain various examples of interfacial phenomena between biomolecules and inorganic materials, including nanotoxicity [32,33]. The remainder of the paper is organized as follows. In the "Materials and Methods" section (Section 2) we will provide a detailed explanation of the theoretical model built to study the protein-metal interaction and the rationale behind the model parametrization scheme. In Section 3, we will discuss simulation results, analyze individual adsorption affinities predicted for molecules representing the bio-part of the interface (amino acids, milk proteins and carbohydrates), and report their preferred orientations. We will also discuss kinetics of competitive adsorption for six most abundant milk proteins to understand the process of protein deposition on the metallic surfaces. Obtained results will be also compared with existing experimental data on protein adsorption. And finally, in Section 4, we will summarize the key insights gained from this study.

2. Materials and Methods

We aim to calculate the content of milk protein layer on iron surfaces using first principles multiscale simulations. To do so, we will use a CG kinetic Monte Carlo (KMC) method [22] to model competitive adsorption of the six most abundant milk proteins. The KMC simulation requires the knowledge of individual binding energies at different orientations (heatmaps) for each selected protein immobilised at each fcc configuration of metallic surface. Heatmaps for individual proteins can be obtained by UA simulations as described in a following sub-section. Although the UA method was already parameterised for the range of crystalline surfaces for noble metals (Ag, Au, Cu), oxides (TiO${}_2$, SiO${}_2$, Fe${}_2$O${}_3$),organic NPs (graphene, carbon nanotubes and carbon black), and semiconductors (CdSe) [34,35] the method is missing the set of short-range potentials essential for calculating milk protein-iron adsorption energies. These potentials in a form of tabulated potentials of mean force (PMF) will be calculated from explicit all-atom simulations through a previously introduced scheme [29].

2.1. Protein-solid surface interaction in UA

Generally, the interaction between a protein and a solid surface, such as an inorganic engineered material (EM), includes several contributions, both specific and non-specific, and depends on the chemical composition, size, shape, surface roughness, surface charge, surface functionalization, and hydrophobicity of the EM. All these aspects should be considered for building CG models of the bio-metal interfaces. Various examples of such models (incl. UA model) have been previously reported to study the competitive adsorption of proteins onto solid surfaces [26,36].

The UA model, previously developed in our lab [26,31], includes all major non-covalent contributions to the interactions between the NP and the protein in a simplified way. The solid surface can be represented by a rigid flat (slabs), spherical (NP), or cylindrical (nanotubes) shape. Proteins are considered as rigid-body structures consisting of 20 different amino acids (AA). Each AA is represented by one bead in the UA model and respective positions of these beads within the protein are fixed. The center of the AA bead is positioned on the AA’s C$\alpha$ atom. The CG model of the protein has the same three-dimensional structure and key details as the original all-atom model. However, fewer beads will decrease the dimensionality of the model and result in lower calculation times. This model is illustrated in Figure 1 along with the presentations on all-atom and CG UA. Further coarse-graining will be applied to reduce the dimension of modelling the kinetics of the competitive adsorption. The ultra-CG hard-sphere (HS) model for selected proteins will be parameterised based on the UA heatmaps.

The UA model divides the NP into a core and a surface segments according to the distances between them and the protein. Then, the interaction potential between each AA and the NP can be represented by a combination of a short-range surface van der Waals (vdW) potential ($U_s^{vdW}$), a long-range core vdW potential ($U_l^{vdW}$), and an electrostatic potential($U^{el}$). The interaction potential between NP and the entire protein ($U_{p-NP}$) is written in a pairwise-additive way via interaction potentials for individual AAs with the NP:

$$U_{p-NP}=\sum _{i=1}^{N_{AA}}{U}_i\left(d_i(\theta ,\varphi )\right)=\sum _{i=1}^{N_{AA}}{U}_i^{el}\left(d_i(\theta ,\varphi )\right)+\sum _{i=1}^{N_{AA}}{U}_i^{vdW}\left(d_i(\theta ,\varphi )\right)$$

(1)

This potential depends on the distance $d_i$ between centers of mass (COMs) of the NP and each AA in the protein. This distance is determined by the protein’s overall orientation relative to the surface of the NP, which is set by two rotational angles $(\theta ,\varphi )$ with respect to the protein’s initial orientation as defined in the PDB file. The electrostatic interaction between NP and AA is described by the screened Coulomb potential:

$$U_i^{el}\left(d_i(\theta ,\varphi )\right)=\sum _{j=1}^{N_e}{l}_B{k}_{B}T{q}_i{q}_j\frac{e^{-\kappa r_{ij}}}{r_{ij}}$$

(2)

where $r_{ij}$ is the distance between the residue of charge $q_i$ and the point charge $q_j$ on the NP surface in terms of the elementary charge $e_0$, $k_B$ is the Boltzmann constant, T is the temperature, $\kappa =\sqrt{4\pi l_{B}I}$ is Debye length, $l_B=\frac{e_0^2}{4\pi ϵ{ϵ}_0{k}_{B}T}$ is the Bjerrum length, and I is the ionic strength, $I=\frac{1}{2}{\sum }_i^{N_{ions}}{c}_i{z}_i^2$ with $c_i,z_i$ being the ion concentrations and valencies, respectively. The properties of the solvent are reflected in the dielectric constant $ϵ$. The VdW potential, representing combined dipole-dipole and dispersion interactions between i-th AA and the NP, includes short-range and long-range terms:

$$U_i^{vdW}\left(d_i(\theta ,\varphi )\right)=U_{i,s}^{vdW}\left(d_i(\theta ,\varphi )\right)+U_{i,l}^{vdW}\left(d_i(\theta ,\varphi )\right)$$

(3)

We extracted the short-range surface potentials from all-atom adaptive well-tempered Metadynamics (AWR-MetaD) simulations by following the procedure described in next section. The long-range term arising from the vdW forces acting through the water medium between the core of the NP and the i-th AA can be approximated by the Hamaker procedure:

$$\begin{array}{c}\hfill U_{i,l}^{vdW}(R_{NP},R_{AA},d(\theta ,\varphi ))=-\frac{A_{123}}{12k_{B}T}\left(\frac{4R_{NP}{R}_{AA}}{d_i{(\theta ,\varphi )}^2-{(R_{NP}+R_{AA})}^2}+\right.\\ \hfill \left.\frac{4R_{NP}{R}_{AA}}{d_i{(\theta ,\varphi )}^2-{(R_{NP}-R_{AA})}^2}+2ln\frac{d_i{(\theta ,\varphi )}^2-(R_{NP}-R_{AA}^2)}{d_i{(\theta ,\varphi )}^2-{(R_{NP}+R_{AA})}^2}\right)\end{array}$$

(4)

Here, the coefficient $A_{123}$ corresponds to the interaction of material 1 with material 3 through the medium 2, $R_{AA}$ and $R_{NP}$ are the radii of AA and NP, respectively.

We then calculate the interaction energy for all possible orientations of the protein, described by orientation angles $\theta$ and $\varphi$ with respect to the initial orientation (taken from the origin PDB file). The total potential energy is found as a function of distance of the protein center of mass from the NP surface.

Integration of the interaction potential over all possible orientations $({\theta }_k,{\varphi }_l)$ and corresponding distances $0\le z\le a({\theta }_k,{\varphi }_l)$ will give the mean interaction energy $E({\theta }_k,{\varphi }_l)$. For flat surfaces, the energy is calculated as:

$$E({\theta }_k,{\varphi }_l)=-k_{B}Tln\left[\frac{1}{a({\theta }_k,{\varphi }_l)}{\int }_0^{a({\theta }_k,{\varphi }_l)}exp\left[\frac{-U_{p-NP}(z,{\theta }_k,{\varphi }_l)}{k_{B}T}\right]dz\right]$$

(5)

For a protein interacting with a spherical NP,

$$E({\theta }_k,{\varphi }_l)=-k_{B}Tln\left[\frac{3}{{(R_{NP}+a({\theta }_k,{\varphi }_l))}^3-R_{NP}^3}{\int }_{R_{NP}}^{R_{NP}+a({\theta }_k,{\varphi }_l)}exp\left[\frac{-U_{p-NP}(r,{\theta }_k,{\varphi }_l)}{k_{B}T}\right]r^{2}dr\right]$$

(6)

The set of rotational configurations along with corresponding $E({\theta }_k,{\varphi }_l)$ are stored in the heatmap. The next step is to calculate the average adsorption energy by using the potential energy as a function of distance for each angle with Boltzmann averaging and weighting factors:

$$E_{ads}=\frac{{\sum }_k{\sum }_l{P}_{kl}E({\theta }_k,{\varphi }_l)}{{\sum }_k{\sum }_l{P}_{kl}},$$

(7)

$$P_{kl}=sin\left({\theta }_k\right)exp\left[-\frac{E({\theta }_k,{\varphi }_l)}{k_{B}T}\right],$$

(8)

where $P_{kl}$ is the Boltzmann weighting factor. Below, we will show how the maximum binding energy affecting the protein adsorption the most can be defined from the adsorption heat maps.

2.2. All-atom model for recovering short-range potentials

The interaction of atoms with a solid surface can be best done using enhanced sampling methods [37,38,39,40]. These techniques include estimating the free energy during simulation and feeding that information back into the dynamics as a statistical bias. The common goal is to spend as little time as possible sampling already sampled free energy zones [41]. Adaptive biasing force methods, which employ an approximation of the mean force to bias the dynamics, are one of the strategies that can speed up the simulation [42,43]. In this work, we use the Adaptive Well-Tempered Metadynamics method that has previously been described for measuring the adsorption of the biomolecules for TiO${}_2$ [29] and Ag [30]. The approach keeps the adsorbate away from previously visited regions of the collective variable space by introducing a time-dependent bias to the system’s potential energy. GROMACS software was used to carry out the simulations in this section of the work [44,45]. CHARMM-GUI/Nanomaterial Modeler was used to create three fcc surfaces of iron: (100), (110), and (111) [46].

The AA side chain analogues (SCAs are the compounds that result from cutting off the side chains of AAs at the protein backbone and substituting a hydrogen for the carbon) were placed in the center top of the slab. In this work, we calculated the energy for 22 SCAs, including two forms of Histidine, HID and HIE, depending on the location of the protonation on the nitrogen atom and two forms of Glutamic acid, GLU and GAN, negatively charged and neutral respectively (See the Supporting Materials for more details). The system was solvated using the original form of TIP3 water model [47]. It was then neutralized by NaCl regarding the charge of the whole system. In the MD calculations, the energy of the system has been minimized by Verlet particle-based cut-off scheme using charge groups and steepest gradient method for 1000 steps. Figure 2 shows the simulation box, the Fe-100 slab and the SCA after the minimization, the water molecules have been removed for better visualization. The system was equilibrated under constant pressure 1.0 bar and temperature $300K$ conditions (NPT ensemble) using Berendsen weak coupling [48] for 1.0 ns. Then, the system was pre-equilibrated for 10 ns in the NVT ensemble. The Nose-Hoover thermostat’s relaxation time constant for the NVT ensemble was 5 ps. The cut-off distance was set to 1.0 nm for the VdW short range interactions. In the metadynamics biased simulation, the surface separation distance (SSD) was used as the reaction coordinate. This SSD measures the minimum distance of the COMs of the SCA $R_{mol}$ and the surface atoms $r_i$ along the z axis,

$$SSD=min|R_{mol}-r_i|,$$

(9)

In the adaptive well-tempered metadynamics, the adsorption energy is measured as a function of SSD. These simulations were done in the NVT ensemble with PLUMED software plugin in GROMACS [44,45,49]. The simulation time was between 500-600 ns for each SCA to get the set short-range one-dimensional PMFs in a tabulated form. The bias factor and the time factor was set to 10 and 500 ps respectively. Every 0.5 ps the Gaussian was added and the initial height was $1.75$ kJ/mol and the temperature was kept at 300 K.

2.3. Preparation of starting coordinates for biomolecules and surfaces

We simulated three fcc surfaces of iron slabs (100, 110, 111) using GROMACS with CHARMM-GUI/Nanomaterial Modeler tool and force fields [46,50,51]. The slab thickness for all surfaces (100, 110, 110) was between 3.53 to 3.71 nm, elongated in z direction in vacuum and were in periodic boundary conditions along x and y coordinates. The final simulation box size were $3.572\times 3.572\times 9.127$ nm${}^3$ for Fe-100, $3.548\times 3.561\times 9.914$ nm${}^3$ for Fe-110 and $3.541\times 3.943\times 9.258$ nm${}^3$ for Fe-111. To calculate the PMFs at the atomistic level, 22 sets of SCAs were chosen which are sufficient to assess adsorption affinities for different kinds of proteins in the CG UA model.

In order to investigate the interaction of entire protein with the NP, the UA model reads the atomistic 3D structure of the protein. These structures might not be available from experiments such as X-ray crystallography, NMR spectroscopy and cryo-electron microscopy. Thus, we predicted 3D structures of the proteins using the sequence data. Here, we used I-TASSER (Iterative Threading ASSEmbly Refinement) software [52] to generate 820 milk protein structures [53,54]. It should be noted, that the majority of these proteins do not exhibit a well defined globular structure. Table 1 shows the list of 6 selected representative milk proteins, UniProt ID, molecular weight, charge and number of AAs in each molecule. The charge information was obtained from PROPKA calculation [55,56] at pH 7.0. All the proteins were subsequently equilibrated in water for 50 ns under NVT and NPT ensembles.

To model competitive protein adsorption from milk, we imitated the natural protein concentrations. The majority of the milk ingredients are water (86–88%), fat (3–6%), protein (3–4%), lactose (5%), and minerals (0.7%) [28]. Milk proteins have many physicochemical properties that allow them to be used in a wide range of applications, from nutritional to functional and biological functions [57]. Caseins and soluble (whey) proteins are classified into two broad categories of milk proteins. Caseins account for about 78 % of bovine milk proteins in total and whey proteins for 17% [58]. Caseins can be classified into four genetic groups based on the similarity of their primary AA sequences [28]. These four groups are $\alpha$s1-caseins, $\alpha$s2-caseins, $\beta$-caseins, and $\kappa$-caseins and correspondingly make up $38\%$, $10\%$, $36\%$, $13\%$ of the total casein [59]. In decreasing quantities, the whey proteins consist of $\beta$-lactoglobulins, $\alpha$ lactalbumins, immunoglobulin and serum albumin of the total $60\%$, $20\%$ and $10\%$ of the whey proteins [60,61]. In the competitive adsorption simulations, we set the protein concentrations of representative in solution according to their relative mass concentrations in milk. We estimated this concentration based on the the percentage of each protein and considering the fact that cow milk has 30-39 g/L protein in total. The molar mass of each protein was taken from AlphaFold database [62].

2.4. Competitive adsorption model

With extracted individual protein UA adsorption heatmaps, we can predict the composition of the adsorbed protein layer on the metallic surface. Proteins, composing this layer, can be fundamentally divided into two classes of hard and soft bound protein [63]. The term "hard" refers to proteins that are bonded directly to the surface of the material with strong affinity and "soft " labels a group of proteins that are weakly bound on top of the adsorbed hard corona via protein-protein interaction [22,63]. We will further study the hard bound layer. The composition of the hard bound layer evolves with time and, depends the morphology and composition of the adsorbent material. In this dynamical process, which can last several hours in real systems, proteins diffuse from the bulk of the milk and adsorb on solid surfaces before occasionally desorbing and being replaced by other proteins. Since the milk sample contains hundreds of proteins, all-atom and even UA simulations of simultaneous adsorption of all milk proteins would be very time consuming [64]. Here, we use an ultra-CG method to describe the actual adsorption-desorption dynamics representing each protein as a single ultra-fine CG bead. Since the adsorption rate of incoming proteins is dependent on previously adsorbed molecules, the whole surface area of the metal will not be accessible for a full coverage by proteins [65]. The HS model, which was initially developed for predicting adsorption onto NPs, considers proteins as rigid spheres with established positions, while the solid is modelled as a spherical NP. In this model, the protein is in a reversible physical contact with the NP’s surface and occupies space that becomes unavailable to other proteins. For each protein, we set a molar concentration in solution as well as rate constants for adsorption and desorption. Due to the slow rates of protein adsorption and desorption, we use reaction-diffusion equations and KMC method to simulate the process [22]. Figure 7 below illustrates the corona formed by milk proteins on a large iron NP using KMC calculations.

3. Results and Discussion

3.1. Short-range potentials for SCAs on iron surfaces

The short-range potentials (in $k_{B}T$) for 22 AA SCAs on Fe (100, 110, 111) obtained by atomistic MD with Metadynamics are shown in Figure 3. Water density profiles from MD simulations for the slab-water system indicate that two water layers with elevated density are formed near the surface at around 0.2 nm and 0.6 nm for all three fcc surfaces. We can identify three different regions: one close to the wall with strong repulsion, an intermediate region between the two layers of water, and the bulk region when the PMF curve comes to a plateau. Traditionally, AAs are classified into four groups as hydrophobic, charged, polar, and aromatic. According to our data, these groups demonstrate different binding affinities to the zero-valent iron surfaces. Hydrophobic AAs (such as ALA, VAL, LEU) bind weakly while aromatic and charged AAs (such as PRO, TYR, ARG and HIS) are more strongly attached. Analysis of the PMFs suggests that the deepest minimum for SCAs adsorption is located close to the first water adlayer which can facilitate adsorption of hydrophilic and charged aminoacids. Similar behaviour was reported for adsorption of AAs onto noble metal surfaces [30,66,67,68,69,70] with the only difference to a sulfur-containing AAs. For example, for the most studied AuNPs the binding affinity ranking was as following: S-containing > aromatic > amines > aliphatic > amides > hydroxylic > carboxylic. The energy of adsorption for individual SCAs were ranging from $-0.35k_B$T to $-39.69k_B$T. This predicted range of adsorption energies for the small molecules in our study is in line with the adsorption free energy of furfural onto nZVI reported in [71]. For ALA SCA (which model is chemically equivalent to methane) onto Fe-110 we found adsorption energy close to $0k_{B}T$, corresponding to the lack of binding. Similar adsorption pattern was reported in [72] where no stable binding complexes were predicted for CH${}_4$ by DFT calculations.

3.2. Protein adsorption energies, preferred orientations and heatmaps

To model milk protein adsoprtion at the metallic surfaces we have selected six most abundant proteins listed in Table 1. These proteins were a subject to UA model simulations to obtain their adsorption binding affinities and preferred orientations for the protein immobilised on the iron surfaces. For these calculations, we used $R=80$ nm and $\zeta$ potential $-5$ mV. We report a single value for each protein after averaging of the energy over all possible orientations. Obtained adsorption energy rankings based on the lowest energy values of the adsorption heatmaps are listed in Table 2. The $\alpha s1$ casein (AS1C) and $\alpha s2$ casein (AS2C) was predicted to be the most strongly bound proteins and $\beta$ lactoglobulin (BLAC) and $\alpha$ lactalbumin (ALAC) are the weakest bound proteins on all Fe (100, 110, 111) surfaces. The preferred orientation of the proteins would not change that much on different fcc faces of iron.

The preferred orientation of all 820 milk proteins based on the lowest energy from heatmaps are reported in the Table S2 of the Supplementary Material. We performed calculations with five Fe NPs of radii 5 nm, 20 nm, 50 nm, 80 nm, and 100 nm to study the size dependence of adsorption energies. The data are presented in the Supplementary Materials (Figure S2). Figure 5 shows the output of the UA model for six selected milk proteins on iron NPs. The heatmap contains the adsorption energies for all values of $\theta$ and $\varphi$. The blue spots (with lower energy) correspond to the more favourable orientations of immobilised proteins. 3D presentation of the protein on the NP surface is shown on top of each heatmap corresponding to the same protein. The closest AAs to the surface of the NP are marked in each case. Analysis of the configurations shows that LYS, GLU, PRO, ASN are most likely to be in contact with the metal surface. This observation is in agreement with the results from Ref. [73] that will be discussed in Section 3.3 section.

The protein adsorption energies predicted by the UA model (Table 2) are high and correspond to irreversible binding of these molecules. This is a result of a high activity of zero-valent Fe. For example for BSA protein the average adsorption energy for all fcc facets was -23.67 $k_B$T for Ag0 [30]and -21.68 $k_B$T for Au(0) [32,74] vs. -75.23 $k_B$T for zero-valent Fe.

3.3. Validation of UA model parameters

To validate our model, we further considered interaction between zero-valent FeNPs and human serum albumin (HSA) through both experimental and docking simulation methods. Analysis using fluorescence spectroscopy showed that the FeNPs formed a complex with HSA through hydrogen bonds and van der Waals interactions. Furthermore, circular dichroism spectroscopy showed that the secondary structure of HSA was not affected by the FeNP. An MD study indicated that the FeNPs interacted with polar residues on the surface of the HSA molecule. The docking study found that ASP, ARG, SER, LYS, GLU residues are most likely to be on the FeNP surface. The reported adsorption free energy G was $-204.80$ kJ/mol, which is comparable to our observations. Taking the same initial protein PDB structure and temperature and size for Fe-110 in UA model we measured the minimum free energy was equal to $-317.84$ kJ/mol at $\theta ={50}^{\circ }$ and $\varphi ={340}^{\circ }$. The nearest interacting AAs in our case are LYS-564, GLU-565, LYS-573, GLU-505, ASP-562, GLU-82, ARG-81, LEU-80, ASP-56, ASP-129.

In reality, the metallic iron surfaces are quickly oxidized and hydroxylized, which changes their adsorption affinity to water and reduces the binding strength for AAs and proteins. Another assumption that affects the range of the observed energies is that the orientation distribution for adsorbed proteins is equilibrium. This assumption is used implicitly in the ensemble averaging procedure represented by Eq. (7). This average energy is dominated by a single preferred orientation. If, however, we assume that proteins bind in random orientations, then the simple average over all angles would reflect the actual binding strength. Yet, in the previous study [30] we found that the Boltzmann average correlates better with the adsorption affinity ranking for proteins, so we used these figures below to analyse binding affinities of individual proteins and to model the competitive adsorption. Section 5 of the Supplementary Materials discusses another experimental approach to validate the binding properties of Hen egg white lysozyme (HEWL) on zero-valent FeNP.

3.4. Competitive adsorption and milk protein layer

We next determined the composition of the milk protein layer at the iron surfaces. As discussed above, in this method we assumed that the surface is represented by a spherical NP with a the protein layer adsorbed on the whole surface of the NP (protein corona) and the adsorption energy is determined by the UA method. The adsorption free energy controls the behaviour between the NP and the protein by using the HS KMC model. The higher the adsorption energy the more probable it is for the protein to be bound on the surface for long periods of time. Calculating the actual adsorption kinetics for different proteins is more challenging. After giving the system enough time the proteins will compete and form a stable protein corona around the NP. The concentrations of milk proteins to determine the corona has been chosen from literature [28,59]. The adsorption-desorption process will stabilise after certain time that depend on the NP radius.

Figure 7 shows a snapshot of 6 milk proteins adsorbing on the iron surface from solution imitating natural milk. Our estimates suggest that the layer of the milk proteins should come to the equilibrium within about 40 minutes due to the high individual adsorption affinities for selected proteins. The process stops when the surface reach full loading capacity with no more space left for further protein adsorption and the protein replacements cease. We can see from the Figure 7 and Table 3 that AS1C (orange colored beads) is adsorbed in largest numbers on the surface and BSA (brown colored beads) has the least abundance in the corona. The left panel of Figure 7 demonstrates the nontrivial kinetics of the corona formation: the concentrations of bound BC, BLAC and ALAC decrease with time after a very quick buildup, while the concentration of AS2C grows much slower only to exceed the amount of the former after about 0.3 hours, evidently replacing the originally bound proteins. The concentration of AS1C remains high after the initial buildup and the concentration of BSA stays low at all times.

Figure 7. Protein adsorption kinetics on a FeNP plot that shows the adsorption-desorption of each protein over time (left). A snapshot of a NP (grey) of radius size 20 nm with protein corona of milk proteins (coloured) adsorbing on the surface using KMC model (right)

Figure 7. Protein adsorption kinetics on a FeNP plot that shows the adsorption-desorption of each protein over time (left). A snapshot of a NP (grey) of radius size 20 nm with protein corona of milk proteins (coloured) adsorbing on the surface using KMC model (right)

We repeated the KMC simulations several times for each NP size and averaged the observable over at least three randomly selected runs. Results of these simulations are listed in Table 3. There, we present the average abundance on the unit area and percentage of relative mass abundance for each milk protein on each surface of fcc Fe. The results in the table show that AS1C is the most abundant both in terms of number of molecules and mass on all the iron surfaces while BSA has the least number and mass abundance. BLAC, and BC are almost equally present in the corona are the second and third most abundant proteins. ALAC and AS2C respectively are the forth and fifth most abundant proteins. This is explained by relatively low molar fraction of BSA in milk as compared to caseins and lactalbumins.

The corona results obtained on NPs of $R_{NP}=80$ nm essentially represent flat iron surfaces. In previous works [31], we have demonstrated that for majority of proteins, the adsorption energies and preferred orientations for most proteins do not change after NP radius of about 30 nm as the surface curvature becomes too large compared to the protein size. Therefore, we expect that the predicted abundances and amounts of deposited proteins should be valid for industrial applications involving steel or iron devices. In the future, we are planning to extend the analysis to competitive adsorption involving lactose and milk fat at different temperatures. Further extension of this work is possible with an inclusion of protein unfolding and denaturation at high temperatures and different pH values. This, however, requires much more complex simulations with flexible proteins and 3D structures adjusted to the specific conditions and may be currently unfeasible in the most detailed form.

4. Conclusions

In this work, we presented the results of modelling the interaction between iron surfaces and most abundant milk proteins in a multiscale scheme based on combination of CG UA model with KMC model for predicting the composition of protein in the deposited layer at iron surfaces. We have considered a simplified model of milk formed from the solution of six most abundant proteins occurring in natural cow milk. Our study ranked the proteins by the adsorption strength as follows: $\alpha$s1 caseins, $\alpha$s2 caseins, $\beta$ casein, bovine serum albumin, $\alpha$ lactalbumin and $\beta$ lactoglobulins, respectively, from strong to weak binding. We have found that the amount of the bound protein depends on its concentration in solution, so the KMC simulation of the adsorption kinetics ranked the proteins differently (in terms of mass fraction in the milk model solution): $\alpha$s1 caseins, $\beta$ lactoglobulins, $\beta$ casein, $\alpha$ lactalbumin, $\alpha$s2 caseins and, finally, bovine serum albumin. In the future, we are planning to extend the multiscale model of the milk adsorption by adding lactose and milk fat molecules.

Our multiscale model of protein corona formation on solid NPs and surfaces can be generalised to a large variety of systems. It essentially relies only on the existence of atomistic force field for the solid. It can be used for variety of applications such as controlling fouling, biofilm growth in food processing and packaging, and medical devices. We also plan to extend the boundaries of the multiscale model to a wider range of temperatures and pH values employed in the industrial settings.