Coarse-grained modeling of the calcium, sodium, magnesium and potassium cations interacting with proteins

Metal ions play important biological roles, e.g., activation or deactivation of enzymatic reactions and signal transduction. Moreover, they can stabilize protein structure, or even be actively involved in the protein folding process. Therefore, accurate treatment of the ions is crucial to model and investigate biological phenomena properly. In this work the coarse-grained UNRES (UNited RESidue) force field was extended to include the interactions between proteins and four alkali or alkaline earth metal cations of biological significance, i.e., calcium, magnesium, sodium and potassium. Additionally, chloride anions were introduced as counter-ions. Parameters were derived from all-atom simulations and incorporate water in an implicit manner. The new force field was tested on the set of the proteins and was able to reproduce the ion-binding preferences.


Introduction
Approximately 40 % of the protein structures deposited in the Protein Data Bank (PDB) contain one or more ions. [12] In some cases these ions were introduced artificially to enable structure determination; nevertheless, most of them play important biological roles, e.g., structure stabilization, activation/deactivation of enzymatic reactions and signal transduction. [39] Therefore, accurate treatment of the ions is crucial to model biological phenomena properly.
The most abundant ions in proteins structures from the PDB are Zn 2+ , Mg 2+ and Ca 2+ . They have significant effect on the protein folding, stabilize the structure and are an integral part of numerous enzymes [9]. If the cell environment is taken into account Na + and K + are the metal ions with the highest concentration. They take part in signal transduction.
Alkali and alkaline earth metal ions have a relatively straightforward pattern of interactions, in contrast to, for example, transition-metal cations, that form coordinate chemical bonds with protein side-chain atoms. Therefore, Mg 2+ , Ca 2+ , Na + and K + are going to be the main subject of this work.
These ions have different interaction preferences for coordinating groups that can be summarized as follows: for Ca 2+ Asp, Glu, Asn, Gln and oxygen from the main chain, for Mg 2+ Asp, Glu and oxygen from the main chain, and for K + and Na + oxygen from the main chain. In case of the two last ions the interactions are almost entirely electrostatic but the coordination geometry is essential. [12] The average concentration of these ions in free state in blood plasma are: 145 mM [16], 4 mM [36], 1.5 mM [13] and 1.8 mM [3] for Na + , K + , Mg 2+ , Ca 2+ respectively. Disturbed homeostasis of these macroelements can lead to many disorders and diseases, such as insomnia, anxiety, urolithiasis or cardiac arrhythmia [46]. Moreover, there is evidence that altered Mg 2+ and Ca 2+ homeostasis may be related to the progression of neurodegenerative diseases associated with intrinsically disordered proteins (IDP): Alzheimer's disease and Parkinson's disease. [2] While investigating the protein-ion interactions, especially when one aims to obtain information about the dynamics, the use of the methods based on the all-atom force fields Anna M. Antoniak and Patryk Wesołowski contributed equally to this work. seems to be appropriate. The extensive review about metal ion modeling has been published, [19] it covers various strategies in quantum mechanics, classical and polarizable force field methods. Generally, three ion models are frequently used: nonbonded, bonded and with dummy atoms. Those representations are used in many common software like: AMBER [18,[20][21][22], CHARMM [27], and GROMOS [41]. Unfortunately, the use of all-atom (AA) methods require high computational power, especially because most ionbinding proteins are often large multichain subunits. In this regard, coarse-grained (CG) models can be a good alternative. However, introduction of metal ions is a challenging task and this issue is often dismissed in CG representations. In one of the most popular CG force field MARTINI ions are represented as a Q type particle (charged interaction site) [31]. In case of single-atom ions the first hydration shell is considered as a part of CG representation. The MARTINI ion model is only qualitatively accurate and the force field itself is unable to model protein folding events. [31,32] Another example is SIRAH force field [8,26] developed to treat proteins with explicit representation of water (WT4) composed by four linked beads and represents 11 tetrahedrally coordinated water molecules [7]. Hydrated ions: Na + , K + , Ca 2+ and Cl − are also considered explicitly, each extended ion particle comprising the actual ion and 6 water molecules corresponding to the first hydration shell. [7,26]. It has been used in the simulations of natively unfolded proteins and protein aggregation and has the ability to capture ionic strength effects.
The coarse-grained UNited RESidue (UNRES) [23,24,51] force field was already successfully extended to treat calcium-protein interactions [15] by deriving the interactions of the ion with Asp, Glu, Asn and Gln side-chains and the peptide group. For the side-chains that do not interact with calcium specifically, simple excluded-volume potentials were introduced. The parameters of the potentials were obtained from potential-energy surfaces of model system: Ca 2+ -acetate, Ca 2+ -propionate, Ca 2+ -acetamide, Ca 2+ -propionamide, and Ca 2+ -N-methylacetamide systems calculated by using ab inito molecular quantum mechanics at the Restricted Hartree-Fock (RHF) level with the 6-31G(d,p) basis set. It should be noted that those calculations were performed in vacuo and obtained potentials of mean force (PMFs) were divided by water dielectric constant. Thus, these potentials do not incorporate interactions with solvent and cover only the interactions with carbonyl groups.
In this article we extended the UNRES force field to investigate the interactions between proteins and four alkali and alkaline earth metals of biological significance, i.e., calcium, magnesium, sodium and potassium. Additionally the parameters for chloride ion were also included. The neutrality of the system is not necessary while performing simulations with the UNRES force field but in case when one wants to represent the environment more accurately, especially under high ionic strength, they can be required.
The method of deriving the potentials was modified with regard to the one previously used for the calcium ion. [15] The following two modifications were introduced. First, the new potentials were derived to incorporate interactions with solvent in an implicit manner, and therefore, the AA force field was used. Second, we considered all amino acid side-chains and peptide group. Hence, the corresponding modifications of the calcium ion-protein interaction center potentials were also recalculated. Introduced potentials have been tested based on the known structures of proteins with bound metal ions.

UNRES force field
In the coarse-grained UNRES model [23,24,51] a polypeptide chain is described by two geometry points per residue namely: C atom and the geometrical center of each sidechain (SC). Moreover, there are two types of interaction sites: peptide group (p), which is positioned halfway between consecutive C atoms and geometrical center of a side-chain. The peptide groups are modeled with spheres whereas sidechains are ellipsoids of revolution along the longer axis (Figure 1). It should be noted that the C atom is not a center of interaction and acts only as a geometric center. The effective energy function is represented by restricted free energy (RFE) or potential of mean force (PMF) of a given conformation ensemble restricted to coarse-grained conformation defined by C and SC atoms and is expressed by Eq. 1 [25]: where the U ′ s are energy terms, i is the backbone virtualbond angle between three consecutive C atoms, i is the backbone virtual-bond-dihedral angle between four consecutive C atoms, i and i are the angles defining the location of the center of the united side-chain of residue i ( Figure 1) with respect to C i−1 , C i and C i+1 plane, d i is the length of the ith virtual bond, which is either a C ⋯C virtual bond or C ⋯ SC virtual bond [15,24].
Each energy term is multiplied by an appropriate weight, w x , and the terms corresponding to factors of order higher than 1 are additionally multiplied by the respective temperature factors that were introduced in paper [23] and which reflect the dependence of the first generalized cumulant term in those factors on temperature, as discussed in refs [23] and [43]. The factors f n are defined by Eq. 2.
where T • = 300 K. The term U SC i SC j represents the mean free energy of the hydrophobic (hydrophilic) interactions between the side chains, which implicitly contains the contributions from the interactions of the side chain with the solvent. The term U SC i p j denotes the excluded-volume potential of the sidechain-peptide-group interactions. The peptide-group interaction potential is split into two parts: the Lennard-Jones interaction energy between peptide-group centers ( U VDW p i p j ) and the average electrostatic energy between peptide-group dipoles ( U el p i p j ); the second of these terms accounts for the tendency to form backbone hydrogen bonds between peptide groups p i and p j . The terms U tor , U b , U rot , and U bond are the virtual-bond-dihedral angle torsional terms, virtual-bond angle bending terms, side-chain rotamer, and virtual-bonddeformation terms, respectively; these terms account for the local properties of the polypeptide chain. The term U (3)

corr
(2) f n (T) = ln exp(1) + exp(−1) represents correlation or multibody contributions from the coupling between backbone-local and backbone-electrostatic interactions, and the term U (3) turn are correlation contributions involving 3 consecutive peptide groups; they are, therefore, termed turn contributions. The term U ion−ion represents the ion-ion interactions with parameters as in our previous work [15] (charge of each ion was adjusted appropriately). The term U ion−prot represents the ion-polypeptide-chain interactions (details of obtaining and verifying this term are described in the next parts of this paper).

Determination of the potentials of mean force for protein-ion interactions
To obtain the PMFs from which the effective energy terms accounting for protein-ion interactions ( U ion ) were derived, simulations of the pairs consisting of the ion and the model of the UNRES interaction center (Table 1) were performed with the all-atom ff14SB force field from AMBER16 package [4]. Each system was placed in a cuboid box containing TIP3P water molecules [14], and counter-ions in order to maintain electroneutrality were added if needed. For Ca 2+ and Mg 2+ the multisite ion model, where the total charges are distributed into n-dummy centers that are placed in the direction of the coordinating atoms (n=7 for Ca 2+ and n=6 for Mg 2+ ), were used. [40] The simulations of systems with Na + , K + and Cl − were performed with the nonbonded-ions models available in AMBER package. [22] The size of the box was set to ensure minimal distance of 12 Å between any atom in our system and the edge of the box. For each system Umbrella Sampling (US) simulations were performed with different harmonic-restraint potentials imposed on the distance between the centers of the mass of each particle. The distances varied from 2 Å to 12 Å, and were sampled every 1 Å.
The simulation procedure consisted of four steps: 1. minimization of the starting structure with restrain at 12 Å distance (for 4000 steps of steepest descent and 3500 steps of conjugate gradient); 2. heating of the system with restrain at 12 Å distance until it reached 300 K (NVT, 20 ps, dt=0.002); 3. equilibration of the system for each of the given harmonic restraints value (NPT, 1000 ps, dt=0.002); 4. molecular dynamic of the system at the given harmonic restraints value (NPT, 10000 ps, dt=0.002).
To determine the PMFs of the systems studied, we processed the results of all restrained MD simulations by using the weighted histogram analysis method (WHAM) [17] to remove the bias from restraints.

Analytical expressions for the energy of interaction of ions with polypeptide chain
The two types of the UNRES interaction sites cover the 19 amino acids (as Gly has no side-chain) and the peptide group. They can be divided into three categories, taking into account their physicochemical properties that are responsible for the nature of interaction with the ions: charged interaction sites (Asp, Glu, Arg, Lys), polar interaction sites (peptide group, Ser, Thr, Asn, Gln, Cys, His, Tyr) and hydrophobic sites (Ala, Val, Ile, Leu, Met, Phe, Trp, Pro). Hereby, the following three combinations of analytical expressions for energy interaction of ions with polypeptide chain were used. For the interactions of the ion with the charged sites, both negatively and positively, the following expression was used: where E GB is the van der Waals energy term for two ellipsoids expressed as the Gay-Berne-type potential, ΔF cav is the cavity term which accounts for the difference between the free energy of hydratation of the side-chain-ion pair and the sum of the free energies of the isolated components, E el is Coulombic term, and e gb is generalized Born potential.
The expression for the ion-polar sites interactions is defined by: where E pol is the polarization energy, and E cp is chargedpolar term.
For the hydrophobic group the expression is defined by: The general analytical expressions with which the effective energy of the interactions of the ions considered in this study are approximated are shown below. [28,29,52] The variables describing the geometry of an UNRES interaction site and an ion are illustrated in Figure 2. The Gay-Berne potential [11] is expressed by the equation: where r ij is the distance between the center of mass of united side-chain or peptide-group model and ion, ij is the van der Waals well-depth for a given orientation of the two interaction sites, ij is the distance corresponding to the zero value of E GB for arbitrary orientation of the particles, and 0 ij is the distance corresponding to the zero value of E GB for the orientation where ij = 90 • . The dependence of the ij and ij on the orientation of the interacting sites is given by Eqs. 7-9 where ij is the van der Waals interaction anisotropy, ′ ij is the anisotropy of the van der Waals depth, and 0 ij is the van der Waals depth of the potential energy minimum for the orientation where ij = 90 • . The fitting parameters are: 0 ij , 0 ij , ij , ′ ij . To account for the interactions with solvent, a cavity term was introduced [28,29]: with where (1) ij , (2) ij , (3) ij , and (4) ij are coefficients, 2 i and 2 j are the deviation of the Gaussian distribution of the solute Fig. 2 The model of the interaction of the UNRES interaction site (particle i) and ion (particle j) where: d head is the distance from the center of particle i and the interacting head of it, r ij is the distance between the center of the particle i and particle j, r ij is the unit vector along the r ij distance, û ij is the unit vector along the long axis of particle i, ij is the angle between vector u ij and vector r ij , r ′ ji is the distance between the interacting "head" of a side chain and ion density for particles i and j, respectively, and �� (1) ij is the anisotropy pertaining to the cavity term. The fitting parameters are: (1) ij The polarization energy originated in the deformation of the charge distribution of the side-chain or peptide group by the ion is expressed by the equation: with f GB being the generalized Born function [49]: where the coefficient 322 is there to express the energies in kcal/mol, in and out are the effective dielectric constant of the "inside" of the interacting particles and the solvent (equal to 80 for water), respectively, pol ji is the polarizability of the nonpolar part of the particle i, a i and a j are the Born radii of particles i and j, and r ′ ji is the distance between the interacting "head" of a side chain or peptide group and ion. The fitting parameters are: a i , in , The Coulombic term is expressed by the equation: where q 1 and q 2 are the charges of particles i and j, respectively. The fitting parameters are: in and d head .
The charged-polar term is expressed by equation [52]: where w ⟂ and w ∥ are fitting parameters, and r ′ ij is a distance between the ion and the center of the amphiphilic head group of side-chain or peptide-group. The fitting parameters are: w ⟂ , w ∥ and d head .
The approximate analytical expressions for ion-interaction-site energy functions were fitted to the PMF calculated numerically for each ion-interaction-site system by means of a nonlinear least-squares Marquardt method. The similar procedure was used before in introducing the potentials for protein-DNA interaction into UNRES force field [52] Subsequently, the obtained expressions were compatibly added to the already existing UNRES energy function.

Simulation procedures
The developed coarse-grained model of polypeptide-ion interaction was tested to predict the ion-binding sites on five proteins listed in Table 2. As a control experiment the all-atom simulations for these five proteins stating from the native structures were performed with the ff14SB force field from AMBER16 package [4]. Each system was placed in a cuboid box containing TIP3P water molecules [14], and counter-ions in order to maintain electroneutrality were added if needed. The size of the box was defined by the distance between the atoms of solute and the box wall and set to 20 Å. The same ion representation as the one used for deriving the potentials was used, being the multisite ion models for Ca 2+ and Mg 2+ [40] and the one-center nonbonded-ion models available in AMBER package [22] for Na + , K + and Cl − . Three 25 ns trajectories in NPT ensemble at 300 K and 1 atmosphere were obtained for each system.
A set of multiplexed replica exchange simulations (MREMD) [6] were performed in temperatures 260 K, 272 K, 279 K, 284 K, 288 K, 291 K, 294 K, 298 K, 308 K, 322 K, 341 K, 370 K. The temperature set was adjusted with the algorithm developed by Trebst et al. [50] which maximizes the walk of replicas in temperature. For each temperature 4 replicas were calculated. Every simulation took 10 000 000 steps with a time step equal to 4.89 fs. The periodic boundary conditions [44] were used with the cubic boxes edges length equal to 150 Å. The simulations were started from the experimental structures from the Protein Data Base (PDB) -including the position of ions. This starting point is sufficient because the MREMD allows to enhance probing of the conformation space and ions are able to dissociate and bound again to the protein. For analysis, only the second half of trajectories were considered. During the simulation temperatures were exchanged between trajectories every 10000 steps; therefore, to obtain the correct ensemble properties reweighing was performed with the bin-less Weighted Histogram Analysis Method (WHAM), which was implemented in UNRES and described in earlier work [23]. Subsequently, the five most populated clusters were obtained with the use of Ward's minimum variance method [47] in the temperature 280 K.
The simulations in higher ion concentrations were performed for the systems where ions dissociate permanently after starting from the PDB structure. It was the case for two proteins: 2Z2K and 5E56, and used concentrations were 55 mM and 70 mM respectively. For this test counterions were used but the rest of the procedure was the same as the one described above.
In the next test association binding constants of Ca 2+ , Mg 2+ , Na + and K + to the -lactalbumin were computed and compared to experimental values. For this purpose the molecular dynamic simulations were performed for the 1F6S [5] protein with the ion ( Ca 2+ , Mg 2+ , Na + or K + ) bound in the native position. The 200 independent trajectories were run at 293 K, each last for 5 000 000 steps with time step equal to 4.89 fs. The periodic boundary conditions were the same as in two other tests.
For the same protein (1F6S) we performed steered molecular dynamics (SMD) [45] with pulling velocity of 0.1 Å ∕ns and with spring constant k = 1kcal∕mol∕Å 2 . The force was computed with maximum distance 20 Å. The simulations for 4 metal ions were performed in order to calculate the protein-ion association strength. The langevin thermostat was used and was set to 293K as previously. The pulling anchors were Asp 88 residue and cation. The 20 independent 5 000 000 steps trajectories with time step equal to 4.89fs were performed. The stretch and force were computed every 1000 steps.

Analysis
The results for the AA simulations were presented as RMSD plots of the protein calculated for the positions of C atoms only and the plots of RMSD values of the ion position calculated when the protein structure was superposed.
For CG simulations starting from the experimental structure the following analysis was performed. First, the representative structure for each cluster (the conformation with the lowest RMSD to cluster average structure) was determined and the TM-score was calculated in order to check its quality.
Subsequently, the distances averaged over all structures for ion from its starting placement for every cluster were calculated. Normalized histograms of those distances were also prepared.
Moreover, the contact maps for ions and particular amino acid side-chains averaged over all structures were prepared for each cluster and compared with the corresponding experimental structure. The scoring function ( S i,j ) for the contact between any C atoms and ion ( d i,j ) is defined as follows [34]: Root-mean square fluctuation (RMSF) was used to analyze the fluctuation of C atoms in each cluster. As a reference structure the cluster representative was taken.
For simulations with higher ion concentration the cluster representatives were also verified by TM-score. Then the radial distribution function (RDF) for the distances between ions and protein were calculated.
In order to calculate the binding constants the concentration of unbound states as a function of time averaged over all 200 trajectories were calculated for each system. An unbound state was defined as a structure where ions were further than 8 Å from any C atom of the protein. Subsequently, function expressed by Eq. 20 was fitted with the use of the nonlinear least-squares (NLLS) Marquardt-Levenberg algorithm [30] implemented in the GNUPLOT program (http:// www. gnupl ot. info).

where E[t] is the concentration of unbound ion and W is expressed by:
Equation 20 is derived based on the method used in [48].

Obtained potentials of mean force and fitted functions
The PMFs, fitted analytical energy functions and parameters are presented in SI ( Figure SI 1-15 and Table SI 1-3). The strongest interactions appeared for the systems containing ions: Ca 2+ and Mg 2+ , and side-chain models for: Asp, Glu (Figure SI 3A-D and Figure SI 6A-D). This is reasonable as there is the strongest Coulombic interaction for oppositely charged molecules. Those two ions strongly interact also with Asn, Gln and the peptide-group models ( Figure  SI

Simulations on the test proteins set
The results for the AA simulations are in SI Figure 16 where the plots of RMSD values for proteins and ions are presented. The strong interaction with protein for all metal ions can be observed. In case of the 4RJ9 protein only one ion is bound to the protein for the entire simulation, while another one dissociates at some point but rebinds to the protein in one of the obtained trajectories. The dissociation of the Cl − is observed in case of the 2Z2K protein, which is expected behavior for a monovalent counter-ion in a MD simulation. Probabilities and population sizes of each cluster, TMscore for cluster representative and average ion displacement  The experimental structures (cartoon representation, gray) and the predicted coarse-grained models (sticks representation, red) for three out of five tested proteins. Ions are represented as spheres. The predicted models are the representatives of the best cluster. The distances between predicted ion placement and the one in the native structure are in the brackets. Results from the simulations started from the experimental structure are presented in Table 3. The probabilities and population sizes of each cluster represent satisfactory distribution of obtained structures. The TM-score shows that at least one cluster in each system has in about the same fold (the value is greater than ∼ 0.5).
Normalized histograms of distances between the predicted ion placement and the one in the native structure are in Figure 3. As assumed, even for the metal ions that strongly interact with the protein, the ion was able to dissociate from the molecule during the simulation (panels A and B). Whereas for the rest of the ions the broader distances are represented (panels C, D and E) indicating weaker interactions.
The contact maps for ion-protein interactions ( Figure 4) and structures ( Figure 5) are presented for proteins 6EKG, 4AQO and 4RJ9 because only for those systems contacts occur (Figure 3).
The best prediction of ion placement was determined for 4AQO protein with Ca 2+ , especially for the third and fifth clusters with average ion displacements 6.52 Å and 2.99 Å, respectively ( Table 3). The quality of predicted binding sites can be verified based on contact maps (Figure 4) where the prediction of the amino acid side-chains that interact with Ca 2+ is in very good agreement with the experimental data. Moreover, the structure of the representative of the cluster with the best prediction (cluster III, based on contact maps) compared with the experimental structure is presented in Figure 5.
Mg 2+ -binding protein (6EKG) strongly interacts with ion ( Figure 3) but based on average ion displacements (Table 3) and contact maps ( Figure 4) the predicted binding sites are 4AQO 4RJ9 B Fig. 7 Logarithm of experimentally determined association constant [35] against the UNRES determined average force for stretching an ion from protein in range from 3.8 Å to 20 Å in SMD simulation. The green solid line represents the linear fitted function with R 2 = 0.8 different from the experimental one. The best prediction is the structure form cluster IV and is presented in Figure 5.
High-quality prediction was also obtained for K + (4RJ9) but only in the least populated cluster. For this protein the maps for two K + are shown in Figure 4 where weak but correct contacts can be observed. The best prediction (the structure form cluster V) is presented in Figure 5.  The worst predictions were obtained for the interactions with Na + and Cl − (proteins 5E56 and 2Z2K, respectively), where ion dissociate from the protein and do not bind to the protein again ( Figure 3D and E).
RMSF plots are presented in SI ( Figure SI 17). Shaded areas indicate the amino acid residues that are in contact with the ion in the native structure. The lowest fluctuations are observed for 5E56 protein ( Figure SI  The results from simulations in high ion concentration are presented in Figure 6A. It can be seen that neither Cl − nor Na + interact with protein. For comparison RDF for the simulations that were started from the experimental structures (low ion concentrations) were calculated for Ca 2+ and K + ( Figure 6B). The strong interactions with both ions can be observed.
The plots of concentration of unbound states of bovine -lactalbumin as a function of time and curves (Eq. 20) fitted to those data are in Figure SI 18. The calculated and experimental equilibrium binding constant are in Table 4. Despite the fact that qualitatively binding affinities are not reproduced the tendency of binding affinity is in agreement for the simulations and experiment with the strongest binding for Ca 2+ then Mg 2+ , Na + and K + , respectively.
The results of the SMD simulations of the 1F6S protein are presented in Figure 7 where logarithm of experimentally determined association constant [35] against the UNRES determined average force for pulling the ion from the protein were plotted. The correlation coefficient for this value is R 2 = 0.8 which shows there is a good agreement between the UNRES force field predicted association strength and the experimental data.

Conclusions
The new approach to simulate the protein-ion interactions in a CG representation was reported in this work. We implemented Ca 2+ , Mg 2+ , Na + , K + and Cl − within the formalism of the physics-based UNRES model. Parameters were derived from AA simulations and incorporate water in an implicit manner.
In order to verify the choice of the AA force field as a method to obtain the new parameters the AA simulations of test proteins were performed. It confirmed the ability of the AA force field to bind the metal ions to the protein and showed that Cl − dissociate from the protein. It is important to remember that the results can be at best of the same quality as the ones obtained by the method used to derive the new potentials.
With the use of CG parameters we were able to predict the binding site of Ca 2+ , Mg 2+ and K + and the strongest interaction was observed for Ca 2+ and Mg 2+ . As for the simulations of Na + and Cl − no interactions with protein were observed. Based on the PDB analysis [12] Na + (and to some extend K + ) with water molecules around often appear in the bulk solvent in the crystal structure, or with a small number of protein oxygen atoms around them. The interactions are likely to be very labile and would not be expected to persist in the solution [12], and our simulations reflect that.
Comparing to results from the AA simulations, the same binding tendency was reproduced for Ca 2+ , Mg 2+ and Cl − . It should, however, be noted that the binding strength of Na + and to some extend K + is underestimated in the new UNRES force field.
Moreover, we were able to reproduce the experimental tendency of binding of the metal ions to bovine -lactalbumin that again verified the strongest affinity for Ca 2+ and Mg 2+ .
Despite the fact that the new UNRES force field is able to produce the good correlation of binding strength with experimental results for the tested systems those parameters should be handled with care when applied to other proteins and should always be validated with experimental results. Therefore, the new potentials allow to represent approximate interaction sites for Ca 2+ , Mg 2+ and K + . They also can be applied to reproduce the labile interaction pattern for Na + and K + . Author Contributions AGL gave the main idea, performed the AA simulation, fitted the analytical expressions for K + , performed and analyzed the test of the method, and prepared the manuscript. AMA fitted the analytical expressions for Cl − and Na + . PW fitted the analytical expressions for Ca 2+ an helped with the "Introduction" section of Table 4 The experimental apparent equilibrium constants for the binding of metal ions to bovine -lactalbumin at 20 °C ( K app ) and calculated association and dissociation rate constants ( k 1 , k −1 ) and binding equlibrium constant ( K a ) helped with the AA simulations. AKS helped with implementation into UNRES software, performed the analysis for kinetic constants, and helped with the rest of data analysis. All authors read and revised the final version of the manuscript.

Funding
The funding for this work was obtained from the National Science Centre (Poland): UMO-2016/21/N/ST4/03154 (to AGL).

Data availability
The datasets generated during and analyzed during the current study are available from the corresponding author on reasonable request.

Declarations
Competing interests The authors declare no competing interests.