Elise Duboué-Dijon

Pauline Delcroix

Hector Martinez-Seara

Jana Hladílková

Pavel Coufal

Tomáš Křížek

Pavel Jungwirth
email: pavel.jungwirth@uochb.cas.cz

Binding of Divalent Cations to Insulin: Capillary Electrophoresis and Molecular Simulations

In the present study, we characterize the binding of divalent cations to insulin in aqueous salt solutions by means of capillary electrophoresis and molecular dynamics simulations. The results show a strong pH dependence. At low pH, at which all the carboxylate groups are protonated and the protein has an overall positive charge, all the cations exhibit only weak and rather unspecific interactions with insulin. In contrast, at close to neutral pH, when all the carboxylate groups are deprotonated and negatively charged, the charge-neutralizing effect of magnesium, calcium, and zinc in particular, on the electrophoretic mobility of insulin is significant. This is also reflected in the results on molecular dynamics simulations showing accumulation of cations at the protein surface, which becomes smaller in magnitude upon effective inclusion of electronic polarization via charge rescaling.

Introduction

Insulin is a key hormone for regulating sugar levels in blood. [START_REF] Cohen | The Twentieth Century Struggle to Decipher Insulin Signalling[END_REF] It is enzymatically active in the bloodstream as a monomer, while being stored in pancreatic secretory granules in the form of inactive hexamers, which crystallize at low pH in the presence of Zn 2+ . [START_REF] Greider | Isolation and Properties of Secretory Granules from Rat Islets of Langerhans II. Ultrastructure of the Beta Granule[END_REF] The role of zinc, which is present in high concentrations (27 mmol/kg dry weight) in the secretory granules, [START_REF] Foster | Elemental Composition of Secretory Granules in Pancreatic Islets of Langerhans[END_REF] in enabling insulin hexamer formation has been known for a long time. [START_REF] Dodson | Structural Relationships in the Two-Zinc Insulin Hexamer[END_REF][START_REF] Dodson | The Role of Assembly in Insulin's Biosynthesis[END_REF][START_REF] Dunn | Zinc-Ligand Interactions Modulate Assembly and Stability of the Insulin Hexamer -a Review[END_REF] Strong binding sites exist for Zn 2+ at the site formed by three B10-chain histidines of three insulin monomers. [START_REF] Dodson | The Role of Assembly in Insulin's Biosynthesis[END_REF] In contrast, much less is known about the interactions of insulin with Ca 2+ , which is present in the granules at even higher concentrations than Zn 2+ . [START_REF] Foster | Elemental Composition of Secretory Granules in Pancreatic Islets of Langerhans[END_REF] Calcium has been implicated to be important in inducing insulin hexamerization by screening the repulsion between negatively charged groups in the hexamer. [START_REF] Coffman | Insulin-Metal Ion Interactions: The Binding of Divalent Cations to Insulin Hexamers and Tetramers and the Assembly of Insulin Hexamers[END_REF][START_REF] Howell | Role of Zinc and Calcium in the Formation and Storage of Insulin in the Pancreatic β-Cell[END_REF] It has been suggested to bind at some of the zinc binding sites. [START_REF] Coffman | Insulin-Metal Ion Interactions: The Binding of Divalent Cations to Insulin Hexamers and Tetramers and the Assembly of Insulin Hexamers[END_REF] Moreover, a binding site for Ca 2+ in the hexamer has been shown to exist in the central cavity, surrounded by six glutamate carboxylate groups. [START_REF] Dunn | Zinc-Ligand Interactions Modulate Assembly and Stability of the Insulin Hexamer -a Review[END_REF][START_REF] Hill | X-Ray Structure of an Unusual Calcium Site and the Roles of Zinc and Calcium in the Assembly, Stability, and Storage of the Insulin Hexamer[END_REF][START_REF] Sudmeier | Cadmium-113 Nuclear Magnetic Resonance Studies of Bovine Insulin: Two-Zinc Insulin Hexamer Specifically Binds Calcium[END_REF] Most of the data about Ca 2+ binding were obtained for the insulin hexamer, [START_REF] Coffman | Insulin-Metal Ion Interactions: The Binding of Divalent Cations to Insulin Hexamers and Tetramers and the Assembly of Insulin Hexamers[END_REF][START_REF] Sudmeier | Cadmium-113 Nuclear Magnetic Resonance Studies of Bovine Insulin: Two-Zinc Insulin Hexamer Specifically Binds Calcium[END_REF][START_REF] Storm | The Glu(B13) Carboxylates of the Insulin Hexamer Form a Cage for Cadmium and Calcium Ions[END_REF] with little concerning the active monomeric form. In this work we employ capillary electrophoresis together with molecular dynamics (MD) simulations to characterize the binding of calcium, magnesium, and zinc divalent cations to the insulin monomer, comparing it to the behavior of monovalent sodium or potassium. Electrophoresis represents a technique well suited for a qualitative comparison of binding of different ions to a protein solute. Such a binding results in measurable changes of the electrophoretic mobility, as demonstrated successfully for small peptides previously. [START_REF] Kubickova | Overcharging in Biological Systems : Reversal of Electrophoretic Mobility of Aqueous Polyaspartate by Multivalent Cations[END_REF][START_REF] Wernersson | Counterion Condensation in Short Cationic Peptides: Limiting Mobilities beyond the Onsager-Fuoss Theory[END_REF] Since this technique requires only very low (micromolar) solute concentrations, we can explore conditions where insulin is mainly monomeric. Combining electrophoresis with MD allows us to characterize the ion-insulin binding with atomistic detail. It is important to realize that insulin can be found in two conformations (T or R), differing in the secondary structure of the B1-B8 residues. [START_REF] Dunn | Zinc-Ligand Interactions Modulate Assembly and Stability of the Insulin Hexamer -a Review[END_REF][START_REF] Brader | Insulin Hexamers: New Conformations and Applications[END_REF] The T form is the most stable form of insulin in solutions in the absence of specific ligands such as phenol. [START_REF] Knegtel | The Solution Structure of Monomeric Insulin[END_REF] In the present simulations, we thus consider only this form.

Methods

Experimental details

Chemicals. Human insulin, recombinant, expressed in yeast, 3-(N-morpholino)propanesulfonic acid (MOPS), ≥ 99.5 % and 2-(N-morpholino)ethanesulfonic acid (MES), ≥ 99.5 % (Sigma, USA), sodium hydroxide, p.a. (Penta, Czech Republic), phosphoric acid, p.a., 85 % (Lach-Ner, Czech Republic), formic acid, p. a. (Merck, Germany), 4-(2-hydroxyethyl)-1piperazineethanesulfonic acid (HEPES), ≥ 99.5 % (Carl Roth, Germany), sodium nitrate, p. a., calcium nitrate tetrahydrate, p. a., magnesium nitrate hexahydrate, pure, zinc nitrate hexahydrate, p. a., zinc oxide, p. a. , calcium hydroxide, pure, magnesium hydroxide, pure, (Lachema, Czech Republic) were used as purchased. Deionized water prepared by a water purification system (Premier MFG'D, USA) was used for preparation of all aqueous solutions.

Electrophoretic conditions. Capillary electrophoresis experiments were carried out on an

Agilent 1600 CE instrument (Agilent Technologies, Germany) using unmodified fused-silica capillaries (Polymicro Technologies, USA), 50 µm I.D., 375 µm O.D.. A capillary of 63.0 cm total length and 54.5 cm effective length was used for measurements at pH 3.0 and 6.3. A capillary of 35.0 cm total length and 26.5 cm effective length was used for measurements at pH 7.5. Prior to the first measurement of the day, the capillary was activated by rinsing with a 1 M NaOH solution for 10 min, followed by water for 10 min using a pressure of 93 kPa.

Between individual runs, the capillary was flushed for 2 min with background electrolyte. The samples were introduced with a pressure of 5 kPa for 5 s. A diode array UV detector was used, measuring the absorbance at 200 nm. The temperature of the capillary cassette was maintained at 25 °C. In measurements at pH 3.0, the separation voltage was set to 30 kV and a pressure of 5 kPa was applied to the inlet buffer vial during the measurement to shorten its duration. In measurements at pH 6.3, the separation voltage was set to 30 kV, and no pressure was applied during the run. In measurements at pH 7.5, the separation voltage was set to 25 kV without applying any pressure.

Background electrolyte and sample solutions. Background electrolyte solutions were prepared by dissolving individual acids in water and adjusting the pH to the desired value with the appropriate salt hydroxide (in solid form or as 1 M aqueous solution). In measurements at pH 7.5, sodium nitrate, calcium nitrate, magnesium nitrate and zinc nitrate were added to the HEPES/NaOH buffer in form of 0.1 M aqueous solution. The ionic strength of the buffers was calculated using Peakmaster 5.2 software. [START_REF] Jaroš | Eigenmobilities in Background Electrolytes for Capillary Zone Electrophoresis: IV. Computer Program PeakMaster[END_REF] Insulin was dissolved in water to 1 mg mL -1 stock solution. From this solution, the sample was prepared by mixing appropriate amounts of insulin stock solution, aqueous solution of thiourea (1 mg mL -1 ), which serves as an electroosmotic flow marker, and background electrolyte.

Computational details

We performed force field molecular dynamics simulations of a monomeric human insulin in the T-conformation in different aqueous salt solutions. The starting protein configuration was taken from a crystallographic structure of the T6 insulin hexamer (PDB code 1MSO), where we selected the first monomer (chains A and B). To model a system at pH close to neutral, both histidine residues of the monomer were considered as deprotonated (HID).

All simulations were performed using the Gromacs 5.1.2 software. [START_REF] Van Der Spoel | GROMACS: Fast, Flexible, and Free[END_REF] The insulin monomer was initially solvated in a box of 5057 water molecules. The system was minimized for 5000 steps using the steepest descent algorithm, before a 10 ns equilibration in the NpT ensemble. Further water molecules and ions were then added to the system to reach a total of 8354 water molecules, 91 cations, and 180 (for CaCl 2 or MgCl 2 ) or 89 (for NaCl) chloride anions. The system was then re-minimized for 5000 steps and equilibrated for 5 ns in the isothermal-isobaric NpT ensemble.

The average box size was then used for each system for further equilibration in the canonical NVT ensemble. The backbone heavy atoms of the insulin monomer were kept fixed during this part of the equilibration procedure (using a harmonic force constant of k=10 5 kJ/mol/nm 2 ), so that the final production runs with different ions started with the same initial protein configuration. The length of the NVT equilibration was varied to ensure proper equilibration of ion distributions around the protein, slower for magnesium due to the slow exchange of water molecules in the hydration shell, [START_REF] Bleuzen | Water Exchange on Magnesium (II) in Aqueous Solution: A Variable Temperature and Pressure[END_REF] than for calcium and sodium (250 ns in NaCl, 500 ns in CaCl 2 , and 1.5 µs for MgCl 2 ). The final 1.6 µs-production runs were performed in the NVT ensemble (without any restrains on the protein backbone) using a 2 fs time step and periodic boundary conditions. Long-range electrostatic interactions were accounted for using the Particle Mesh Ewald method [START_REF] Darden | Particle Mesh Ewald: An Nlog(N) Method for Ewald Sums in Large Systems[END_REF] with a 12 Å cutoff for non-bonded interactions. Configurations were saved every 5 ps. In all simulations, the temperature was kept at 300 K using a velocity-rescaling thermostat with a stochastic term [START_REF] Bussi | Canonical Sampling through Velocity Rescaling[END_REF] and a time constant of 1 ps. Simulations in the NpT ensemble used a Parinello-Rahman barostat [START_REF] Parrinello | Polymorphic Transitions in Single Crystals : A New Molecular Dynamics Method[END_REF] with a 5 ps time constant. Bonds involving hydrogen atoms were constrained by the LINCS algorithm, [START_REF] Hess | A Linear Constraint Solver for Molecular Simulations[END_REF] and water molecules were kept rigid using the SETTLE algorithm. [START_REF] Miyamoto | An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models[END_REF] Water molecules were described by the SPC/E force field. [START_REF] Berendsen | The Missing Term in Effective Pair Potentials[END_REF] while. Two types of force fields were used for the ions and protein. First, we used conventional "full charge" force fields for Na + , [START_REF] Smith | Computer Simulations of NaCl Association in Polarizable Water Computer Simulations[END_REF] Ca 2+ 26 and Mg 2+ , [START_REF] Callahan | Solvation of Magnesium Dication: Molecular Dynamics Simulation and Vibrational Spectroscopic Study of Magnesium Chloride in Aqueous Solutions[END_REF] and the Amber99sb-ildn force field [START_REF] Lindorff-Larsen | Improved Side-Chain Torsion Potentials for the Amber ff99SB Protein Force Field[END_REF] for the protein. We also performed simulations with the recently developed "scaled charge" Empirical Continuum Correction (ECC) force fields. Within the ECC description, the fast electronic polarization is taken into account in a mean-field way via scaling the charges of the ions and charged protein residues by a factor 𝑜𝑓 1/√𝜀 !" = 0.75, where 𝜀 !" = 1.78 is the electronic part of water dielectric constant. In this work, we used the previously developed ECC force fields for the Na + , [START_REF] Kohagen | Accounting for Electronic Polarization Effects in Aqueous Sodium Chloride via Molecular Dynamics Aided by Neutron Scattering[END_REF] Ca 2+ 30 and Mg 2+ 31 cations, as well as for the Cl -counterions, [START_REF] Pluharova | Hydration of the Chloride Ion in Concentrated Aqueous Solutions Using Neutron Scattering and Molecular Dynamics[END_REF] and a new ECC version of the Amber99sb-ildn force field. The employed ionic parameters are summarized in Table 1 and details about the protein ECC force field are provided in the Supporting Information.

Results and discussion

Measurements at pH 3.0

The initial experiments involved measurements of the insulin mobility in acidic background electrolytes of pH 3.0. Under these conditions, insulin was strongly positively charged and migrated towards the cathode. The mobility was measured in buffers prepared from formic or phosphoric acids by titration with hydroxides of sodium, potassium, magnesium, or calcium, or with zinc oxide. The ionic strength of all the buffers was kept at 53 mmol L -1 . As can be seen in

Figure 1, the nature of the cation did not have a significant effect on the mobility of insulin.

Insulin was injected at two concentration levels of 4.3 or 57.3 µmol L -1 . Given the range of estimates for the dimerization constant depending on pH and ionic strength found in the literature, [START_REF] Dodson | The Role of Assembly in Insulin's Biosynthesis[END_REF][START_REF] Weiss | Two-Dimensional NMR and Photo-CIDNP Studies of the Insulin Monomer: Assignment of Aromatic Resonances with Application to Protein Folding, Structure, and Dynamics[END_REF][START_REF] Jeffrey | An Equilibrium Ultracentrifuge Study of the Self-Association of Bovine Insulin[END_REF][START_REF] Rupley | Concentration Difference Spectra in the Dimerization of Insulin[END_REF][START_REF] Goldman | Zinc Binding, Circular Dichroism, and Equilibrium Sedimentation Studies on Insulin (Bovine) and Several of Its Derivatives[END_REF][START_REF] Strazza | The Thermodynamics of Bovine and Porcine Insulin and Proinsulin Association Determined by Concentration Difference Spectroscopy[END_REF] we expect the insulin, in the absence of divalent ions, to be mainly monomeric at 4.3 µmol L -1 and to be mainly dimeric at 57.3 µmol L -1 . Mobility of insulin at the latter concentration tends to be slightly higher than that at the former one, which may be caused by a stronger tendency for dimerization at higher protein concentration. Compared to the monomer, the electric charge of the dimer increases more than its hydrodynamic radius, which leads to increased electrophoretic mobility. Changing the buffer composition from formate to phosphate affected the mobility only slightly without any systematic trend. Driving voltage was 30 kV, at a pressure of 5 kPa applied to the inlet end of the capillary at a temperature of 25 °C.

Measurements at pH 6.3

As a next step, we measured the electrophoretic mobilities of insulin at higher pH values. At high pH conditions, the carboxylate groups of the glutamate residues (four per insulin molecule) and of C-terminal residues are deprotonated and one may thus expect the effects of different cations to be more distinct. However, some of the tested divalent cations precipitated at higher pH with the phosphate buffer anions. We thus had to use different buffers at high pH than in the previous experiments at pH = 3.0. Since we wanted to avoid carboxylate-based buffers, which contain the same negatively charged group as the protein, we used the MES and MOPS buffers, both based on sulfonic acid groups. The range of accessible pH conditions was limited by the fact that insulin is only marginally soluble between pH = 5-6, and that the solubility of calcium and magnesium hydroxides strongly decreases with increasing pH. The highest accessible pH with calcium and magnesium cations was thus found to be pH = 6.3 for an ionic strength kept at 53 mmol L -1 . Buffers were prepared by titrating MES and MOPS acids with the appropriate hydroxides. This procedure, however, did not work with zinc oxide due to salt precipitation, where we had to use much lower zinc concentrations with a background sodium electrolyte and go for higher pH to increase the insulin solubility in the presence of Zn 2+ (see next Section).

The results obtained for sodium, potassium, magnesium, and calcium are shown in Figure 2.

At pH = 6.3, the net charge of the insulin molecule is negative and the protein thus tends to migrate toward the anode. Note that the absolute values of mobilities in MES buffers are slightly higher than those in MOPS buffers. This is probably due to a slight difference in the interactions of the two buffers with insulin. Nevertheless, in both buffers, the same ordering is seen among the individual cations. While mobilities in sodium vs potassium solutions do not significantly differ from each other, the mobilities decrease (i.e., become less negative) upon adding divalent cations (Figure 2). The presence of magnesium significantly slows down the migration of insulin with the effect of calcium being even stronger and the mobility of insulin approaching zero (Figure 2). With the MES buffer, the mobility of insulin with calcium cations even changes sign, indicating a weak overcharging of insulin by the interaction with calcium cations (Figure 2).

These results overall suggest stronger binding of divalent ions to insulin, and stronger binding of

Ca 2+ than Mg 2+ .
Mobilities in all buffers were measured at two insulin concentrations. As previously, the two concentrations were chosen so that the insulin is expected to be predominantly monomeric at the lower concentration, while being mainly dimeric at the higher one. However, note that in the MES buffer, the highest concentration level had to be lowered to 34.4 µmol l -1 due to the limited solubility of insulin. Consistently with the experimental data shown in Figure 1, higher insulin concentrations result in slightly higher (i.e., more negative) mobilities. Measurements at pH 7.5

In order to investigate the effect of zinc cations on the mobility of insulin, one has to choose a pH range where insulin is more soluble. We thus chose to work at pH = 7.5. As already discussed, such a high pH favors the precipitation of calcium and magnesium hydroxides, so that we had to change experimental conditions compared to the previously described experiments.

We used in all cases a 10 mM HEPES/NaOH buffer to adjust the pH, and added into it only small amounts of calcium, magnesium, or zinc nitrates. Figure 3 shows that the presence of zinc cations in concentrations as low as units of µmol L -1 causes a significant increase in the mobility of insulin, while concentrations of calcium cations up to 50 µmol L -1 have no significant effect on the mobility of insulin. This increased mobility of insulin in presence of zinc cations can be (Figure 3A) and 34.4 µmol L -1 (Figure 3B). The mobility obtained for the lower concentration of insulin steeply increases with increasing zinc concentrations and the effect of zinc rapidly saturates. For the higher concentration of insulin, the increase in mobility is more gradual, which is in agreement with the lower zinc-to-insulin ratio. In this case, the zinc concentration range for which insulin monomers, dimers, and hexamers co-exist is broader. Also, the starting mobility of insulin is more negative with higher concentration of insulin. As already discussed, this is probably due to partial dimerization of insulin even in the absence of zinc.

We see from Figure 3 and Figure 4 that the addition of calcium and magnesium nitrates only starts to influence significantly the mobility of insulin at almost two orders of magnitude higher concentrations than for zinc. The observed effect of calcium and magnesium ions on the mobility is opposite to that of zinc, with the absolute value of the insulin mobility decreasing. This can be attributed to the interaction of insulin with calcium and magnesium cations and consequent screening of its negative charge. Since calcium, magnesium, and zinc nitrates were added in varied amounts to the HEPES buffer of constant concentration, the increase in ionic strength caused by the addition of the nitrates may also have affected the insulin mobility (the concentration of the added divalent salt was up to 1 mmol L -1 and the ionic strength thus increased from 5.5 mmol L -1 to 8.5 mmol L -1 ). To investigate the effect of such a change of ionic strength, we measured the mobility of insulin in the HEPES buffer with addition of 3 mmol L -1 sodium nitrate. The mobility in this background electrolyte did not significantly differ from that measured in pure HEPES buffer. We can thus conclude that the observed change in mobility with addition of calcium or magnesium is predominantly caused by the interaction of insulin with these cations and does not originate from the increased ionic strength.

The differences in the sign of the effect and in the concentration range of cation leading to observable effects support the conclusion that the influence of zinc on insulin mobility is based on a different mechanism (i.e., insulin aggregation) than that of calcium and magnesium (i.e., charge screening).

Molecular Dynamics Simulations

Electrophoretic measurements allow us to characterize the specific effect of different cations on the insulin mobility, but probe only indirectly ion-protein interactions. We thus perform complementary force field molecular dynamics (MD) simulations of an insulin monomer in aqueous NaCl, CaCl 2 and MgCl 2 solutions in order to obtain further molecular details. Note that we do not attempt to model binding of zinc to insulin, since the present force fields do not capture the interaction of Zn 2+ with imidazole ligands, important also for the insulin hexamerization process. [START_REF] Dunn | Zinc-Ligand Interactions Modulate Assembly and Stability of the Insulin Hexamer -a Review[END_REF] Direct comparison of specific ion effects between all-atom MD simulations and electrophoretic experiments is difficult for several reasons. First, due to the limited size of the simulation box and need of proper sampling of ion binding within the simulation time scale, the salt concentration in the simulation of around 0.5 mol L -1 is about an order of magnitude higher than the highest salt concentration used in the experiments, because the low concentrations would be hard to simulate in a statistically converged way. Secondly, the simulated systems contain chloride salts instead of the experimental MOPS or MES buffers. Third, the electrophoretic mobility is a complex quantity, which is difficult to obtain from simulations. The electrophoretic mobility is connected to the electrostatic potential at the slipping plane (or slip plane), and thus only reflects the effective protein charge at the slipping plane. [START_REF] Delgado | Measurement and Interpretation of Electrokinetic Phenomena[END_REF][START_REF] Heikkilä | Atomistic Simulations of Functional Au 144 (SR) 60 Gold Nanoparticles in Aqueous Environment[END_REF] However, the specific value of the electrophoretic mobility is difficult to compare with simulation data, because the position of the slipping plane is not clearly defined and depends on other factors, such as particular effects in presence of an electric field. For small spherical solutes, the electrophoretic mobility can be approximately obtained from the weighted integral of the soluteions radial distribution functions. [START_REF] Wernersson | Counterion Condensation in Short Cationic Peptides: Limiting Mobilities beyond the Onsager-Fuoss Theory[END_REF] Yet, it is not clear to what extent such an approach is applicable to a non-spherical, non-homogeneously charged and flexible protein molecule. Taking into account the impact of protein flexibility has been attempted in the literature using coarsegrained peptide models, [START_REF] Allison | Conductivity and Electrophoretic Mobility of Dilute Ionic Solutions[END_REF][START_REF] Allison | Modeling the Electrophoresis of Oligolysines[END_REF] but such a coarse-grained approach is unlikely to be able to capture subtle specific ion effects.

We thus used the present simulations to characterize the locations and differences in ion bindings to insulin for the investigated salts. The ion density maps obtained with the standard Amber99sb-ildn force field (Figure 5A-C) exhibit extensive pairing of all the ions with any negatively charged protein group (i.e., glutamate and C-terminal residues). The ion pairing with carboxylate groups is particularly strong for divalent ions, with almost six Ca 2+ ion in direct contact with protein carboxylate groups (less than 4 Å away from the carbon atom), as shown by the radial distribution function pictured in Figure 6A. On average, almost all of the six carboxylate groups are thus ion paired with a calcium ion. The pairing is even stronger with Mg 2+ , with about ten Mg 2+ in direct contact with the protein's carboxylate groups. This means that more than half of the carboxylate groups are bound not only to one Mg 2+ ion but to two Mg 2+ ions simultaneously, which is an artifact of the simulations. A typical configuration of a doubly bound carboxylate group is shown in Figure 7. The effective protein charge (i.e., sum of the net charge of the protein and the charge of the surrounding ions within a given distance of the protein heavy atoms, which effectively define the protein surface) as obtained with the standard force field (Figure 8) is highly positive with both divalent cations and close to zero with the sodium ion, up to 8 Å from the protein surface. This is in clear contradiction with the experimental negative sign of the electrophoretic mobility for Na + and Mg 2+ and null mobility for Ca 2+ , and thus demonstrates the unphysical behavior of the standard force field, with overbinding cations.

The over-binding of cations is a known artifact of full charge force fields, [START_REF] Kohagen | Accounting for Electronic Polarization Effects in Aqueous Sodium Chloride via Molecular Dynamics Aided by Neutron Scattering[END_REF][START_REF] Martinek | Calcium Ions in Aqueous Solutions: Accurate Force Field Description Aided by Ab Initio Molecular Dynamics and Neutron Scattering[END_REF] even if attempts have been made to remedy this problem by careful re-parametrization and generalization of combination rules for van der Waals interactions. [START_REF] Mamatkulov | Force Fields for Divalent Cations Based on Single-Ion and Ion-Pair Properties[END_REF][START_REF] Kahlen | Interaction of Charged Amino-Acid Side Chains with Ions: An Optimization Strategy for Classical Force Fields[END_REF][START_REF] Satarifard | Is the Solution Activity Derivative Sufficient to Parametrize Ion-Ion Interactions? Ions for TIP5P Water[END_REF] Previous studies have shown that the lack of electronic polarization in non-polarizable force fields often leads to an overestimation of ion pairing in MD simulations. [START_REF] Kohagen | Calcium Binding to Calmodulin by Molecular Dynamics with Effective Polarization[END_REF][START_REF] Pluharova | Ion Pairing in Aqueous Lithium Salt Solutions with Monovalent and Divalent Counter-Anions[END_REF][START_REF] Vazdar | Aqueous Guanidinium -Carbonate Interactions by Molecular Dynamics and Neutron Scattering: Relevance to Ion -Protein Interactions[END_REF] Recently, we demonstrated in particular that calcium binding to carboxylate groups was overestimated in standard simulations compared to reference ab initio simulations. [START_REF] Martinek | Calcium Ions in Aqueous Solutions: Accurate Force Field Description Aided by Ab Initio Molecular Dynamics and Neutron Scattering[END_REF] Following previous theoretical developments, [START_REF] Leontyev | Electronic Continuum Model for Molecular Dynamics Simulations[END_REF][START_REF] Leontyev | Accounting for Electronic Polarization in Non-Polarizable Force Fields[END_REF] we showed that by effectively including electronic polarization effects via scaling the charges of ions by a factor of 0.75 allowed us to recover proper ion binding properties. This Electronic Continuum Correction (ECC) charge scaling procedure has been shown to be a physically well-justified mean-field way to incorporate electronic polarization in non-polarizable force fields. 49 in red) as obtained from the simulation performed with the A) full charges and B) scaled charge ECC force field. Integrated numbers of ions n(r) (dashes) are evaluated as 𝑛 𝑟 = 6𝑐 !"# 4𝜋𝑟 ! 𝑔 𝑟 𝑑𝑟. A multiplicative factor of six (number of carboxylate groups) has been used so that at short distances n(r) reflects the average total number of ions in direct contact with protein carboxylate groups. Effective protein charge or a standard full charge force field (dashes). In the case of NaCl, the ECC simulations were performed with either the same salt concentration as for the divalent ions (solid line), or, to be consistent with the experiments, with the same ionic strength (dot and dashes).

In order to obtain more realistic ion bindings, we thus repeated the simulations using an ECC force field for the ions and the protein, where the total charges of the side chains of charged amino acids are also scaled by a factor of 0.75. We checked for small pentapeptides composed of four glycine residues and a central charged amino acid (glutamate or arginine) that modifying the side chain charges did not significantly affect the geometric properties of the protein (see Supporting Information). In addition, we compared the insulin backbone flexibility in our 1.6 µs simulations employing either the standard full charge or the ECC force field. The root mean square fluctuations of backbone atoms averaged per residue did not evidence any significant change in protein flexibility or conformation, even though, as expected, detailed analysis of salt bridges showed reduced ion pairing between the involved residues (more details in the Supporting Information). Note that previous studies suggested [START_REF] Vazdar | Aqueous Guanidinium -Carbonate Interactions by Molecular Dynamics and Neutron Scattering: Relevance to Ion -Protein Interactions[END_REF] that the stability of salt bridges may be overestimated in standard force fields.

The ion density maps obtained with the scaled charges ECC force field show a dramatic reduction of ion pairing with protein residues compared to the previous full charges simulations (Figure 5D-F), the main ion-pairing sites remaining the protein carboxylate groups. Only about half a cation is found on average in direct contact (i.e., less than 4 Å away from the carbon atom)

with protein carboxylate groups (Figure 6B) at the investigated salt concentration. While every carboxylate group was in a contact pair with at least one divalent ion with the full charges force field, it is now only the case for less than one in six groups on average. We note that for sodium and calcium the peaks corresponding to the contact ion pair (between 2.5 and 4 Å) are split in two sub-peaks, which correspond to the bidentate (shorter distances) and monodentate (larger distance) binding geometries (Figure 7A-B). In contrast, Mg 2+ binds exclusively in monodentate geometries. If we include the ions forming solvent-shared ion pairs (a peak between 4 and 5 Å) (Figure 7C), we find that there are about two ions involved in an ion pair with the protein (for the six carboxylate groups). Also, the overall pairing seems to be slightly higher for Mg 2+ than for Na + and Ca 2+ , which exhibit very similar ion pairing properties with carboxylate groups in our ECC simulations. These ECC simulations lead to a negative apparent protein charge (Figure 8), in agreement with the sign of the experiment electrophoretic mobility, with little difference between the different electrolytes.

The experimental electrophoretic mobility data yield lower insulin mobilities (in absolute value) with divalent ions than with monovalent ions, which suggested a stronger binding of divalent ions compared to monovalent. This is because binding of the positively charged cations to the negatively charged protein would reduce its apparent charge and lead to lower mobility.

Similarly, experimental measurements suggest stronger binding for Ca 2+ than for Mg 2+ . Our simulation results focusing on the carboxylate groups, which are expected to be the strongest cation binding sites, are thus apparently inconsistent with the experimental findings. However, we note that the employed force field has been shown to be in very good agreement with ab initio calculations of the binding free energy between calcium ions and carboxylate groups. [START_REF] Martinek | Calcium Ions in Aqueous Solutions: Accurate Force Field Description Aided by Ab Initio Molecular Dynamics and Neutron Scattering[END_REF] In addition, a recent ab initio MD study found very similar binding free energies for Ca 2+ and Na + ions with carboxylate groups, [START_REF] Daily | Divalent Ion Parameterization Strongly Affects Conformation and Interactions of an Anionic Biomimetic Polymer[END_REF] which is consistent with the present simulations. Given the significant error bars (usually on the order of 0.5-1 kcal/mol) of AIMD simulations, our force field, while consistent with these simulations, may still somewhat underestimate the difference in binding between Ca 2+ and Na + . More importantly, even if carboxylate groups are expected to be the strongest binding sites, cations can also interact with backbone carbonyl groups.

We thus examined in more detail the ion binding to the backbone carbonyl groups constructing the radial distribution functions between any of the protein backbone carbonyl oxygen and the employed cations (Figure 9). Contrary to what was observed with carboxylate groups, ion binding to the carbonyl groups is not strongly affected by the use of the scaled charge ECC vs.

full charge force fields (except that the use of the ECC force field somewhat increases the amount of Mg 2+ at the backbone). With both full charges and ECC, sodium forms more direct contacts (distance lower than 3 Å) with carbonyl groups than calcium. This results is in disagreement with previous ab initio MD studies [START_REF] Daily | Divalent Ion Parameterization Strongly Affects Conformation and Interactions of an Anionic Biomimetic Polymer[END_REF][START_REF] Pluhařová | Aqueous Cation-Amide Binding: Free Energies and IR Spectral Signatures by Ab Initio Molecular Dynamics[END_REF] which consistently suggested, using model amide compounds, that Ca 2+ binds more strongly than Na + to backbone carbonyl groups. Even if the binding at carbonyl groups is weaker than at negatively charged carboxylate side chains, the number of exposed carbonyl groups is much higher than that of carboxylates. Hence, the fact that the employed force fields seem not to fully capture the specific ion bindings at the backbone may explain why our simulations do not reproduce the ordering in the specific ion binding deduced from electrophoresis measurements. Additional explanations may include the presence of the buffer in the experiment, particular effects in the presence of an explicit electric field, and/or the ability of divalent cations to bridge two carboxylic groups from different protein molecules (not considered in the present calculations). reflects the average total number of ions in direct contact with protein backbone carbonyl groups.

We performed additional simulations at pH = 3 (details in SI), where all the carboxylate groups are protonated, so that the cations are expected to interact only with the backbone carbonyl groups. As expected, the overall cation binding to insulin becomes weaker with the protein being now positively charged in overall. The specific bindings of the different ions follow the same trend as that observed at the carbonyl groups at high pH, with stronger binding of sodium than divalent ions, which is likely due to force field artifacts mentioned above.

Conclusions

Capillary electrophoresis reveals pH-dependent and cation-specific effects on the electrophoretic mobility of insulin in aqueous solutions. At pH = 3, at which all the carboxylate groups are protonated and the protein bears an overall positive charge, there is virtually no effect of sodium, potassium, magnesium, calcium, or zinc cations on the electrophoretic mobility of insulin. In contrast, at close to neutral conditions, when insulin in negatively charged and all carboxylates are deprotonated, divalent alkali earth cations (but not monovalent alkali cations) significantly reduce the absolute value of the electrophoretic mobility. For zinc, the effect on the mobility is by far the strongest, but it has the opposite sign due to induction of insulin hexamerization.

The electrophoretic measurements are complemented by molecular dynamics providing a detailed atomistic view on the ion-protein interactions. MD simulations employing standard nonpolarizable force fields grossly overestimate the affinities of the cations to the deprotonated carboxylate groups at the surface of insulin. This is largely remediated by inclusion of electronic polarization effects in a mean-field way via charge rescaling. Nevertheless, quantitative comparison to the experiment is difficult due to problems with accurate modeling of electrophoretic mobilities of large non-spherical solutes, as well as residual force field inaccuracies pertaining in particular to the cation-protein backbone interactions.

Figure 1 .

 1 Figure 1. Effective electrophoretic mobilities of insulin at I = 53 mmol L -1 in formate and

Figure 2 .

 2 Figure 2. Effective electrophoretic mobility of insulin in MES and MOPS buffers with sodium, potassium, calcium and magnesium cations, pH 6.3, I = 53 mmol L -1 . Driving voltage 30 kV,

  of insulin hexamers. Compared to the monomer, the electric charge of the hexamer increases more than the hydrodynamic radius resulting in increased (more negative) electrophoretic mobility. The electrophoretic mobility of insulin in the presence of varying concentrations of zinc and calcium was obtained for protein concentrations of 8.6 µmol L -1

Figure 3 .Figure 4 .

 34 Figure 3. Effective electrophoretic mobility of insulin in 10 mM HEPES/NaOH buffer, pH 7.5,

Figure 5 .

 5 Figure 5. Ion density maps around insulin for the three investigated cations (green: sodium,

Figure 6 .

 6 Figure 6. Radial distribution function g(r) (solid lines) between the carbon atom of the six protein carboxylate groups and the investigated cations (Na + in green, Ca 2+ in orange, and Mg 2+

Figure 7 .

 7 Figure 7. Snapshots of typical ion paired geometries with carboxylate groups: A) monodentate contact ion pair, B) bidentate contact ion pair, C) solvent-shared ion pair, and D) double Mg 2+

Figure 8 .

 8 Figure 8. Effective protein charge as a function of the distance from the protein heavy atoms (i.e.

Figure 9 .

 9 Figure 9. Radial distribution functions g(r) (solid lines) between the protein backbone carbonyl oxygen atoms and the ions (Na + in green, Ca 2+ in orange and Mg 2+ in red) as obtained from the

Table 1 .

 1 Details of the employed ion parameters (both full charges and ECC descriptions).

	ion	q (e)	Full charge force field σ (Å) ε (kJ/mol)	q (e)	ECC force field σ (Å) ε (kJ/mol)
	Na +	+1	2.3500	0.5439	+0.75	2.1150	0.5443
	Ca 2+	+2	2.9132	0.4184	+1.5	2.6656	0.5072
	Mg 2+	+2	1.8900	3.6610	+1.5	1.3600	3.6610
	Cl -	-1	4.4000	0.4184	-0.75	4.1000	0.4928

Acknowledgment This work has been supported by Charles University Research Centre program No. UNCE/SCI/014. P.J. acknowledges support from the Czech Science Foundation (grant no. 16-01074S). EDD acknowledges support from the EMBO and Marie Curie Actions (fellowship ALTF 952-2015). This work was performed using computer time allocated by the national supercomputing center IT4innovations in Ostrava (project OPEN-10-26) and MetaCentrum data storage facilities.

Supporting Information

Details about the development and testing of the ECC protein force field, additional results of simulations at low pH.

TOC graphics

Ca 2+

Insulin electrophoretic mobility

Na + Mg 2+ Ca 2+