Julie Puyo

Marie Juillé

Jérôme Hénin

Carine Clavaguéra

Elise Duboué-Dijon
email: duboue-dijon@ibpc.fr

A consistent picture of phosphate-divalent cation binding from models with implicit and explicit electronic polarization

The binding of divalent cations to the ubiquitous phosphate group is essential for a number of key biological processes, such as DNA compaction, RNA folding or interactions of some proteins with membranes. Yet, probing their binding sites, modes and associated binding free energy is a challenge for both experiments and simulations. In simulations, standard force fields strongly overestimate the interaction between phosphate groups and divalent cations. Here, we examine how different strategies to include electronic polarization effects in force fields-implicitly through the use of scaled charges or pair-specific Lennard-Jones parameters, or explicitly with the polarizable force fields Drude and AMOEBA-capture the interactions of a model phosphate compound, dimethylphosphate, with calcium and magnesium divalent cations. We show that both implicit and explicit approaches, when carefully parametrized, are successful in capturing the overall binding free energy, and that common trends emerge from the comparison of different simulation approaches. Overall, the binding is very moderate, slightly weaker for Ca 2+ than Mg 2+ , and the solvent-shared ion pair is slightly more stable than the contact monodentate ion pair. The bidentate ion pair is higher in energy (or even fully unstable for Mg 2+ ). Our results thus suggest practical ways to capture the divalent cations with biomolecular phosphate groups in complex biochemical systems. In particular, the computational efficiency of implicit models makes them ideally suited for large-scale simulations of biological assemblies, with improved accuracy compared to state-ofthe-art fixed-charge force fields.

Introduction

Phosphate groups are ubiquitous in biochemistry, [START_REF] Westheimer | Why nature chose phosphates[END_REF] as they not only form the backbone of nucleic acids and the head group of phospholipids, but are also present in metabolites such as NADPH, ADP, or ATP, and are added to protein side chains as post-translational modifications. Their interactions with divalent cations such as Mg 2+ or Ca 2+ play a key role in many fundamental biochemical processes. For instance, Mg 2+ ions have been shown to modulate nucleic acid structures and assist RNA folding. [START_REF] Draper | Ions and RNA Folding[END_REF][START_REF] Bowman | Cations in charge: Magnesium ions in RNA folding and catalysis[END_REF][START_REF] Woodson | Metal ions and RNA folding: A highly charged topic with a dynamic future[END_REF] They are also essential for the full catalytic activity of ribozymes. [START_REF] Hanna | Metal ions in ribozyme folding and catalysis[END_REF] In the cell, calcium ions can modulate the interactions of phospholipid membranes with proteins, and the formation of Ca 2+ -phosphate based structures is also essential for bone mineralization. [START_REF] Brini | Calcium in health and disease[END_REF] However, despite their high biological relevance, the binding modes and associated energetics of biological phosphate compounds with divalent cations are not fully characterized yet. For instance, Mg 2+ cations are usually found to interact with nucleic acid phosphate groups through their hydration water molecules (thus forming solvent-shared ion pairs or "outer-sphere" coordination), [START_REF] Egli | DNA-Cation Interactions: Quo Vadis?[END_REF][START_REF] Auffinger | Metal ion binding to RNA[END_REF] except at specific RNA binding sites, where direct "inner-sphere" interactions have been evidenced (Fig. 1). [START_REF] Bowman | Cations in charge: Magnesium ions in RNA folding and catalysis[END_REF][START_REF] Auffinger | Metal ion binding to RNA[END_REF][START_REF] Klein | The contribution of metal ions to the structural stability of the large ribosomal subunit[END_REF][START_REF] Freisinger | From nucleotides to ribozymes-A comparison of their metal ion binding properties[END_REF] The conditions allowing for direct magnesium-phosphate interactions are still not fully understood, also because the identification of magnesium cations in biological structures is sometimes controversial. [START_REF] Auffinger | Deflating the RNA Mg2+ bubble. Stereochemistry to the rescue![END_REF] There is thus a strong need to better understand phosphate-divalent cation interactions, to be able to predict their binding modes, sites, and energetics in biological systems.

A first step is to characterize the ion binding behavior to small model phosphate compounds. The thermodynamics of Mg 2+ binding to dihydrogen phosphate is well characterized experimentally, [START_REF] Verbeeck | Solubility of magnesium hydrogen phosphate trihydrate and ion-pair formation in the system Mg(OH)2-H3PO4-H2O at 25°C[END_REF] with a binding free energy of ∆G bind (Mg 2+ ) = -7.1 kJ/mol. The body of experimental data is thinner for Ca 2+ , but tends to suggest a comparable, slightly weaker binding free energy to dihydrogen phosphate (∆G bind (Ca 2+ ) = -5.9 kJ/mol [START_REF] Mcdowell | Solubility Study of Calcium Hydrogen Phosphate. Ion-Pair Formation[END_REF] ). The same trend between the two cations is observed for binding to more complex phosphate compounds (e.g. ATP, CTP, ADP). [START_REF] Smith | Critical Evaluation of stability constants for nucleotide complexes with protons and the accompanying enthalpy changes[END_REF] In contrast, the information regarding the structure of the corresponding ion pairs remains very partial. Very recently, 2D-IR experiments on dimethylphosphate (DMP) suggested a significant fraction of direct contact ion pairs with both Ca 2+ and Mg 2+ in 0.2 M concentrated solutions, [START_REF] Schauss | Contact Ion Pairs of Phosphate Groups in Water: Two-Dimensional Infrared Spectroscopy of Dimethyl Phosphate and ab Initio Simulations[END_REF] but a decomposition of the overall binding energetics into the contribution of the different ion pairs is still missing.

In this context, molecular dynamics simulations can complement experimental data The Pcation distance can be used to distinguish the different ion pairs, with, for Mg 2+ , typical values of 2.7-3.2 Å for the contact bidentate binding mode, 3.2-2.7 Å for the contact monodentate and 4.1-6.1 Å for the SShIP. Snapshots were prepared using the VMD software. [START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF] and directly investigate the relative binding free energies corresponding to different binding modes. [START_REF] Huston | Free energies of association for the sodium-dimethyl phosphate ion pair in aqueous solution[END_REF][START_REF] Kumar | Structure and thermodynamics of Mg:phosphate interactions in water: A simulation study[END_REF][START_REF] Villa | Classical Drude Polarizable Force Field Model for Methyl Phosphate and Its Interactions with Mg2+[END_REF] However, the interactions of divalent cations to phosphate groups have proved very challenging to capture in simulations, because of polarization and possible charge transfer effects. In addition, the very slow (µs) kinetics of water exchange in the magnesium cation hydration shell poses acute sampling problems. Standard fixed-charge force fields have been shown to exhibit strong overbinding artefacts, [START_REF] Grotz | Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation , Ion-Binding , and Water-Exchange Properties[END_REF][START_REF] Panteva | Force Field for Mg2+, Mn2+, Zn2+, and Cd2+ Ions That Have Balanced Interactions with Nucleic Acids[END_REF] in line with similar observations pertaining to the binding of divalent cations to other negatively charged biomolecular groups. [START_REF] Melcr | Accurate Biomolecular Simulations Account for Electronic Polarization[END_REF][START_REF] Duboué-Dijon | A practical guide to biologically relevant molecular simulations with charge scaling for electronic polarization[END_REF][START_REF] Oliveira | Binding of divalent cations to acetate: Molecular simulations guided by Raman spectroscopy[END_REF] Several approaches have been suggested to improve the simulated binding behavior, with either explicit or implicit inclusion of electronic polarization.

Polarizable force fields explicitly account for electronic polarization, which can be done in different ways. [START_REF] Halgren | Polarizable force fields[END_REF][START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF][START_REF] Bedrov | Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields[END_REF] The binding of Mg 2+ ions with phosphate groups has been investigated in the past using two different polarizable force fields on slightly different model compounds: AMOEBA, which uses permanent electrostatic multipole moments up to the quadrupole and includes explicit dipole polarization, on dihydrogenphosphate (DHP), [START_REF] Kumar | Structure and thermodynamics of Mg:phosphate interactions in water: A simulation study[END_REF] and the Drude polarizable model on methylphosphate. [START_REF] Villa | Classical Drude Polarizable Force Field Model for Methyl Phosphate and Its Interactions with Mg2+[END_REF] They were both successful in reproducing the overall binding free energy, and predicted solvent shared ion pairs (SShIP) to be predominant. These studies were however limited to Mg 2+ binding, and did not compare it with Ca 2+ .

Another strategy is to account for induced polarization implicitly, within the standard fixedcharge functional form of the force-field. Two such mean-field approaches have been suggested in the literature, either through the use of pair specific Lennard-Jones (LJ) parameters [START_REF] Grotz | Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation , Ion-Binding , and Water-Exchange Properties[END_REF][START_REF] Fyta | Ionic force field optimization based on single-ion and ion-pair solvation properties: Going beyond standard mixing rules[END_REF][START_REF] Yoo | New tricks for old dogs: Improving the accuracy of biomolecular force fields by pair-specific corrections to non-bonded interactions[END_REF][START_REF] Kashefolgheta | Developing force fields when experimental data is sparse: AMBER/GAFF-compatible parameters for inorganic and alkyl oxoanions[END_REF][START_REF] Hess | Cation specific binding with protein surface charges[END_REF] (the so-called NBFIX corrections in CHARMM), or through scaling the total charge of the ion (Electronic Continuum Correction, or ECC). [START_REF] Duboué-Dijon | A practical guide to biologically relevant molecular simulations with charge scaling for electronic polarization[END_REF][START_REF] Leontyev | Electronic continuum model for molecular dynamsimulations[END_REF][START_REF] Leontyev | Accounting for electronic polarization in non-polarizable force fields[END_REF][START_REF] Kirby | Charge Scaling Manifesto: A Way of Reconciling the Inherently Macroscopic and Microscopic Natures of Molecular Simulations[END_REF] The strength of both approaches is that they allow the use of the same software as standard non polarizable force fields, with no additional computational cost. Pairspecific parameters have been recently derived to describe the interactions of Mg 2+ and Ca 2+ with dimethylphosphate, enabling simulations to reproduce the overall binding free energy and hydration kinetics. [START_REF] Grotz | Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation , Ion-Binding , and Water-Exchange Properties[END_REF][START_REF] Cruz-León | Extended magnesium and calcium force field parameters for accurate ion -nucleic acid interactions in biomolecular simulations[END_REF] The scaled-charge ECC description has shown promising results for the interaction of polyvalent phosphate species with cations [START_REF] Bruce | Does an electronic continuum correction improve effective short-range ion-ion interactions in aqueous solution?[END_REF] but has yet to be tested on monovalent phosphate compounds.

In this work, our goal is to gain molecular insight into Mg 2+ and Ca 2+ binding to a model phosphate, DMP, using molecular dynamics simulations. We aim to compare how the relative binding of the two cations and the relative free energy of different binding modes are described by different force fields-non polarizable and with implicit or explicit electronic polarization. We examine how ion pairing free energies and binding modes are captured by a scaled-charge ECC force field, which we compare with a recently proposed force field based on pair-specific Lennard Jones parameters. [START_REF] Grotz | Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation , Ion-Binding , and Water-Exchange Properties[END_REF] For comparison, we have used two polarizable force fields-the Drude polarizable force field [START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF] and AMOEBA, [START_REF] Ponder | Current status of the AMOEBA polarizable force field[END_REF][START_REF] Zhang | AMOEBA Polarizable Atomic Multipole Force Field for Nucleic Acids[END_REF] -which have been recently used to study ion-phosphate interactions, [START_REF] Kumar | Structure and thermodynamics of Mg:phosphate interactions in water: A simulation study[END_REF][START_REF] Villa | Classical Drude Polarizable Force Field Model for Methyl Phosphate and Its Interactions with Mg2+[END_REF] and are both available for nucleic acids and proteins, so that they can be used for simulation of complex biomolecular systems. This makes it possible to identify robust trends in phosphatecation binding and compare the merits of different simulation strategies, which will pave the way for a better understanding of more complex, biologically-relevant systems.

Computational details

Molecular Dynamics simulations

Molecular dynamics simulations of DMP with a divalent cation-either Mg 2+ or Ca 2+ -were performed with different force fields to compare the ion pairing behavior of the two cations and assess its sensitivity to the employed force field. The simulation box was composed of one DMP anion, one cation (Mg 2+ or Ca 2+ ) and 1723 water molecules.

We performed simulations using force fields from different families, with either implicit or explicit descriptions of electronic polarization. The employed parameters and computational protocol are detailed below. All the parameter files and typical input files for the different softwares are provided in an online repository (https://zenodo.org/record/5971002#.YkrBDDw6-YA).

Non polarizable force fields. Several parametrizations of divalent cations using a non polarizable force field and standard integer charges (+2 for divalent cations) are available in the literature. In the main body of this paper, the simulations with standard "non polarizable" force fields were performed with one of the most commonly used parametrizations, namely the 12-6 Lennard-Jones parameters suggested by the Merz group for both Mg 2+ and Ca 2+ , choosing the SPC/E compatible versions (labelled "Mg_m" and "Ca_m" in Table 1). [START_REF] Li | Rational design of particle mesh ewald compatible Lennard-Jones parameters for +2 metal cations in explicit solvent[END_REF] In the supporting information, we provide additional comparison with simulations using the Mg 2+ force field suggested by D. Tobias and coworkers, [START_REF] Callahan | Solvation of magnesium dication: molecular dynamics simulation and vibrational spectroscopic study of magnesium chloride in aqueous solutions[END_REF] and the Ca 2+ force field developed in the Netz group. [START_REF] Mamatkulov | Force fields for divalent cations based on single-ion and ion-pair properties[END_REF] The dimethyl phosphate was described using the General Amber Force Field (GAFF), which presents similar charges as the phosphate backbone in current Amber nucleic acid force fields (see Supporting Information). These ion force fields were combined with the SPC/E model 42 for water.

All simulations with non polarizable force fields were performed with the Gromacs2019.6 software [START_REF] Berendsen | GROMACS: A message-passing parallel molecular dynamics implementation[END_REF][START_REF] Abraham | Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF] with a 2 fs integration step, using periodic boundary conditions, in the constant temperature/constant pressure (NpT) ensemble. We used a Parinello-Rahman barostat [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] with a 5 ps coupling time and the velocity rescaling thermostat [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] at 300 K with a 1 ps coupling time. Long-range electrostatic interactions were treated with Particle Mesh Ewald (PME), [START_REF] Darden | Particle mesh Ewald: An Nlog(N) method for Ewald sums in large systems[END_REF] with a 12 Å cutoff. Hydrogencontaining bonds were constrained using the LINCS algorithm, [START_REF] Hess | LINCS: A linear constraint solver for molecular simulations[END_REF] and the SETTLE algorithm for water molecules. [START_REF] Miyamoto | SETTLE: an analytical version of the SHAKE and RAT-TLE algorithm for rigid water models[END_REF] Electronic Continuum Correction. The Electronic Continuum Correction (ECC) is a mean-field theory that prescribes the scaling of ionic charges by a factor 1 √ el 0.75 (where el is the high frequency part of water dielectric constant) to implicitly account for electronic polarization. [START_REF] Leontyev | Accounting for electronic polarization in non-polarizable force fields[END_REF][START_REF] Kirby | Charge Scaling Manifesto: A Way of Reconciling the Inherently Macroscopic and Microscopic Natures of Molecular Simulations[END_REF] We used previously developed scaled-charge force fields with a 0.75 scaling factor for Mg 2+ 50 (labelled "Mg_s") and Ca 2+ 51 ("Ca_s"), which we combined with a scaled-charge version of the DMP GAFF ("GAFF_s" in Table 1) force field, where we simply scaled all atomic charges by the prescribed 0.75 factor (such a simple strategy was used successfully with acetate previously [START_REF] Martinek | Calcium ions in aqueous solutions: Accurate force field description aided by ab initio molecular dynamics and neutron scattering[END_REF] ). These ion force fields were used together with the SPC/E model 42 for water. Previous studies have shown that a 0.75 scaling factor often tends to underestimate the strength of the ion pairs. [START_REF] Oliveira | Binding of divalent cations to acetate: Molecular simulations guided by Raman spectroscopy[END_REF][START_REF] Kirby | Charge Scaling Manifesto: A Way of Reconciling the Inherently Macroscopic and Microscopic Natures of Molecular Simulations[END_REF] This probably originates from the fact that part of the electronic polarization is effectively wrapped into the employed empirical water model, whose dielectric constant is larger than what would correspond solely to nuclear polarization. [START_REF] Vega | Simulating water with rigid non-polarizable models: a general perspective[END_REF] Following earlier suggestions, we thus also designed ECC force fields with a milder 0.8 scaling factor, both for the cations 24 ("Mg_s2" and "Ca_s2" in Table 1) and for the DMP ("GAFF_s2").

Simulations with scaled-charge force fields were performed using Gromacs, following the same protocol as with standard force fields.

Pair-specific Lennard-Jones parameters. Prompted by the observation that none of the existing non polarizable force fields for Mg 2+ captures at the same time its hydration thermodynamics, its hydration water exchange kinetics, and its binding free energy to biomolecular phosphate groups, Schwierz and collaborators recently proposed 20 a new Mg 2+ force field (microMg) where they selected the best combination of 12-6 Lennard-Jones parameters to capture simultaneously the thermodynamics and kinetics of Mg 2+ solvation. In addition, to ensure proper binding properties with RNA, they went beyond the standard Lorentz-Berthelot combinaison rule and optimized specific Lennard-Jones parameters for the pairwise interaction between Mg 2+ and the phosphate non bridging oxygens, targeting the overall binding free energy. A similar strategy was later adopted for Ca 2+ and led to the so-called Ca2 force field. [START_REF] Cruz-León | Extended magnesium and calcium force field parameters for accurate ion -nucleic acid interactions in biomolecular simulations[END_REF] microMg and Ca2 parameters were initially optimized for TIP3P water, [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF] and a GAFF parametrization of dimethylphosphate, [START_REF] Grotz | Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation , Ion-Binding , and Water-Exchange Properties[END_REF] slightly different from ours, that we hence adopted in simulations with microMg. Note that this force field was very recently extended to be compatible with other water models. [START_REF] Grotz | Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation, Ion-Binding, and Water-Exchange Properties in SPC/E, TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D[END_REF] Simulations with these force fields followed the same protocol as with standard force fields. The Drude particles fluctuate around their parent atom during the course of the simulation, depending on the local electric field, thus modeling the dipole response. Together with the Drude polarizable water model SWM4-NDP, [START_REF] Lamoureux | A polarizable model of water for molecular dynamics simulations of biomolecules[END_REF] we used Drude parameters for DMP and Mg 2+ specifically refined [START_REF] Lemkul | Balancing the Interactions of Mg2+ in Aqueous Solution and with Nucleic Acid Moieties For a Polarizable Force Field Based on the Classical Drude Oscillator Model[END_REF] to reproduce Mg 2+ solvation properties and its binding to nucleic acid moieties.

In addition to Drude particles, this force field includes pair-specific Lennard-Jones interactions (NBFIX) and dipole screening (NBTHOLE) terms both for Mg 2+ -water and Mg 2+ -phosphate (non bridging oxygen) interactions. Simulations were performed with the NAMD2.14 program, 57 with a timestep of 1 fs, in the NpT ensemble at room temperature and pressure, using Langevin dynamics with a 1 ps -1 damping period and a Langevin Piston Nosé-Hoover barostat. [START_REF] Feller | Constant pressure molecular dynamics simulation: The Langevin piston method[END_REF] Periodic boundary conditions were used with a PME treatment. [START_REF] Darden | Particle mesh Ewald: An Nlog(N) method for Ewald sums in large systems[END_REF] The temperature of the Drude oscillators was maintained at 1 K with a separate Langevin thermostat with a 20 ps -1 damping coefficient.

AMOEBA. The AMOEBA [START_REF] Ponder | Current status of the AMOEBA polarizable force field[END_REF][START_REF] Ren | Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation[END_REF] force field models the electrostatic potential with atomic permanent electrostatic multipole moments up to the quadrupole, and includes explicit selfconsistent dipole polarization. It has been extended to describe phosphates (incl. DMP) [START_REF] Zhang | Polarizable multipole-based force field for dimethyl and trimethyl phosphate[END_REF] and nucleic acids. [START_REF] Zhang | AMOEBA Polarizable Atomic Multipole Force Field for Nucleic Acids[END_REF] Parametrization of the nucleic acid backbone in the most recent version of the force field (amoebanuc17.prm) is based on the initial DMP parametrization, [START_REF] Zhang | Polarizable multipole-based force field for dimethyl and trimethyl phosphate[END_REF] which we thus used here in combination with the water and Mg 2+ parameters currently recommended for nucleic acids. [START_REF] Zhang | AMOEBA Polarizable Atomic Multipole Force Field for Nucleic Acids[END_REF] Simulations with the AMOEBA force field were performed with Tinker 61 on GPU using the Tinker-OpenMM implementation, [START_REF] Harger | Tinker-OpenMM: Absolute and relative alchemical free energies using AMOEBA on GPUs[END_REF] in the NVT ensemble at 300 K. The RESPA integrator [START_REF] Tuckerman | Reversible multiple time scale molecular dynamics[END_REF] was used with a 2 fs time step. Long-range electrostatic interactions were treated with PME, using a realspace cutoff of 8 Å .

To examine the sensitivity of the binding free energy ∆G • bind to the set of multipoles, we refined new multipoles using clusters with explicit water molecules, that should thus be suited for simulation of DMP in solution. The geometry of a cluster of 26-water molecules solvating DMP was optimized at the DFT level using the ωB97X-D3 functional and the def2-SVP basis set, with the Orca program package. [START_REF] Neese | The ORCA program system[END_REF][START_REF] Neese | Software update: the ORCA program system, version 4.0[END_REF] Following the AMOEBA procedure, 66 the atomic multipolar distribution of the optimised DMPwater cluster was derived from the MP2/cc-pVDZ density using the Gaussian program [START_REF] Frisch | Gaussian09 -Revision D01[END_REF] and DMA analysis. [START_REF] Stone | Distributed multipole analysis Methods and applications[END_REF] The multipoles were then refined using an adjustment of atomic dipoles and quadrupoles on the electrostatic potential while keeping the charges fixed.

Binding free energy calculations

Binding free energy calculations were performed in each binding mode using the broadly adopted double decoupling theoretical framework. [START_REF] Jorgensen | Efficient computation of absolute free energies of binding by computer simulations. Application to the methane dimer in water[END_REF][START_REF] Hermans | The Free Energy of Xenon Binding to Myoglobin from Molecular Dynamics Simulation[END_REF] The corresponding thermodynamic cycle and notations are presented in details in Refs. 71 and 72. Two alchemical transformations are performed, where the cation is progressively decoupled from its environment, either in the bulk (∆G * bulk ) or in the ion pair (∆G * site ). The free energy contribution due to the change in net charge during the alchemical transformation (which induces a change in neutralizing "gellium" in our simulations with periodic boundary conditions and PME summation) is taken into account with an analytic correction ∆G PBC (eq (17) from Ref.

73). Ion pair geometries during the computation of ∆G *

site are maintained with flat-well harmonic restraints on the distance between the cation and the DMP phosphorus atom, which allows to distinguish the different binding modes (see Figure 1). The free energy of applying the restraints in the uncoupled state with respect to the standard state concentration, ∆G V • →rest decoupled , is estimated numerically as r) dr , where L is the box size and U rest the flat-well potential. In the coupled state, since the employed flat-well potential is very close to the infinite well which is taken as the definition of the bound state, the contribution of the restraints ∆G rest→site coupled (i.e. the free energy associated with removing the restraints in the coupled state, to have an unrestrained ion in the defined binding site) is almost exactly 0. Additional details on the employed thermodynamic site and meaning of the different free energy terms can be found in Ref. 72. The standard binding free energy for a given binding mode ∆G • bind is thus obtained by adding those contributions:

∆G V • →rest decoupled = -RT ln ∞ 0 4πr 2 e -βUrest(r) dr -RT ln L 0 4πr 2 e -βUrest(
∆G • bind = ∆G * bulk + ∆G PBC + ∆G V • →rest decoupled - ∆G * site + ∆G rest→site coupled
The overall binding free energy was obtained by combining the binding free energies obtained in each ion pair geometry, as follows:

∆G • bind,tot = -RT ln e -∆G • bind,mono /RT + e -∆G • bind,bi /RT + e - (1
) This alchemical approach for the calculation of binding free energies has been shown previously [START_REF] Deng | Computations of standard binding free energies with molecular dynamics simulations[END_REF][START_REF] Gumbart | Standard binding free energies from computer simulations: What is the best strategy?[END_REF] to be formerly equivalent to PMFbased determinations, if, as we do here, the free energy contributions of the restraints are reliably evaluated. In the present case, the high free energy barriers separating the different ion pairs could make obtaining a converged PMF challenging, especially with Mg 2+ . In contrast, the separate estimation of affinities for each of the binding modes reduces correlation times and accelerates convergence, hence our method choice here.

Prior to the alchemical transformations, short unrestrained equilibrations were started from different ion pairing geometries (Fig. 1) with each force field to estimate the range of Pcation distances corresponding to each ion pair. The detailed restraint potentials used in each case during the alchemical transformations are provided in the Supporting Information.

The computational details for the alchemical transformations depend on the simulation software, and thus the force field, as follows.

Non polarizable force fields. For all non explicitly polarizable force fields, alchemical transformations were carried on with Gromacs. Electrostatic interactions between the cation and its environment were first turned off (11 windows, λ elec = 1, 0.9, 0.8, ..., 0.1, 0) before switching off van der Waals interactions (4 windows with λ vdw = 0.7, 0.5, 0.1, 0). Starting in the coupled state, each window was equilibrated for 250 ps before starting the next one, and another 5 ns sampling was then added independently in each window. The Bennett Acceptance Ratio 76 (BAR) method, as implemented in Gromacs, was then used to obtain the free energy. The error bar provided on ∆G • bind is the sum of the error for each alchemical transformation, which are estimated from the average variance over 5 blocks as implemented by default in gmx bar. On a selected system, we replicated the alchemical transformations three times. The observed standard deviation of the computed free energies is consistent with the error estimate reported by "gmx bar". We also checked that the computed free energy is independent of the direction of the alchemical transformation (forward / backward), which is a good sign of convergence (see details in Supporting Information).

Drude force field. The alchemistry module in NAMD is not yet fully compatible with the Drude force field. Hence, we decomposed the decoupling in two stages, as suggested previously. [START_REF] Villa | Classical Drude Polarizable Force Field Model for Methyl Phosphate and Its Interactions with Mg2+[END_REF] First, the Mg 2+ charge and associated Drude particle charge were gradually reduced to zero, using scaling factors of 1, 0.95, 0.9, ..., 0.1, 0.05, 0. A 2 ns-simulation was performed with each of the 21 corresponding PSF files, printing out the structures and energies every 1 ps. In a post-processing step, the energy of each frame of the trajectories was recalculated with the potentials of the neighboring windows. These energies were then collected to yield a BAR 76 estimate of the free energy. In a second stage, the Lennard-Jones interactions of the uncharged Mg 2+ were gradually removed (λ LJ = 1, 0.9, 0.8, 0.6, 0.4, 0.2, 0) using NAMD FEP alchemistry module with Interleaved Double-Wide Sampling, which allows the use of the BAR estimator. For both stages, the BAR free energy estimation was performed in a post-processing step using a python jupyter notebook [START_REF] Kluyver | Jupyter Notebooks -a publishing format for reproducible computational workflows[END_REF] and the alchemlyb [START_REF] Shirts | Statistically optimal analysis of samples from multiple equilibrium states[END_REF][START_REF] Chodera | A Simple Method for Automated Equilibration Detection in Molecular Simulations[END_REF] library. Error estimates were obtained by subsampling the energy data to obtain a decorrelated set, then performing the MBAR error analysis as described by Shirts and Chodera [START_REF] Shirts | Statistically optimal analysis of samples from multiple equilibrium states[END_REF] and implemented in the pymbar library.

AMOEBA. Alchemical transformation with the AMOEBA force field were performed with the Tinker-openMM software. [START_REF] Harger | Tinker-OpenMM: Absolute and relative alchemical free energies using AMOEBA on GPUs[END_REF] The interactions of the cation with its environment were progressively decoupled using 24 4-ns long windows, first progressively reducing the electrostatic and polarization interactions (λ elec = 1, 0.9, ..., 0.1, 0) before decoupling the van der waals interactions (λ LJ = 1, 0.9, 0.8, 0.75, 0.7, 0.65, 0.6, 0.55, 0.5, 0.4, 0.2, 0.1, 0). The final free energy estimate and error was obtained with BAR as implemented in the bar_omm.x Tinker tool.

Computation of free energy landscapes with ABF

Free energy landscapes were sampled and computed using a variant of the Adaptive Biasing Force (ABF) method, 80 extended-system ABF (eABF). [START_REF] Lelièvre | Computation of free energy profiles with parallel adaptive dynamics[END_REF][START_REF] Lesage | Smoothed Biasing Forces Yield Unbiased Free Energies with the Extended-System Adaptive Biasing Force Method[END_REF] Consider a vector collective variable (CV) z = ξ(x) , where x is the vector of atomic Cartesian coordinates. We define the free energy landscape, up to an additive constant, as the quantity

A(z) ≡ -k B T ln ρ(z) ,
where ρ(z) is the marginal probability density of the collective variable. In ABF, the gradient ∇A of the free energy surface A(z) is estimated on-the-fly as the conditional average of a projected force, [START_REF] Darve | Calculating free energies using average force[END_REF][START_REF] Hénin | Exploring multidimensional free energy landscapes using time-dependent biases on collective variables[END_REF][START_REF] Comer | The adaptive biasing force method: everything you always wanted to know but were afraid to ask[END_REF] following a Thermodynamic Integration formalism. This running estimate is applied as a biasing force, which tends to level the local slope of the free energy landscape. At long times, as the estimate converges to the free energy gradient, the effective free energy landscape becomes flat and the biased distribution of z becomes uniform. In ABF, estimating the projected force imposes technical requirements on the CVs. [START_REF] Hénin | Exploring multidimensional free energy landscapes using time-dependent biases on collective variables[END_REF] These requirements are relaxed in eABF, by not applying the ABF algorithm directly to the CVs of interest, but to fictitious degrees of freedom λ that evolve according to Langevin dynamics and are harmonically coupled to the CVs through a potential V ext (x, λ) = k 2 |ξ(x) -λ| 2 with force constant k. [START_REF] Lelièvre | Computation of free energy profiles with parallel adaptive dynamics[END_REF][START_REF] Lesage | Smoothed Biasing Forces Yield Unbiased Free Energies with the Extended-System Adaptive Biasing Force Method[END_REF] The implementation used was the Collective Variables Module (Colvars), [START_REF] Fiorin | Using collective variables to drive molecular dynamics simulations[END_REF] which was recently interfaced with Gromacs. Gromacs patched with Colvars was downloaded from https://github.com/Colvars/gromacs/tags. Three collective variables were defined, corresponding to the three Cartesian coordinates of the divalent cation in a moving frame of reference tied to the phosphate group. An ideal, symmetric geometry of the PO 4 group was created and used to define the reference position to which those atoms were aligned transparently throughout the simulations, by invoking the fittingGroup option of Colvars. This allowed for the resolution of well-defined free energy landscapes reflecting the positional distribution of divalent cations with respect to the PO 4 group. Unbiased free energy gradients were obtained by the CZAR estimator, [START_REF] Lesage | Smoothed Biasing Forces Yield Unbiased Free Energies with the Extended-System Adaptive Biasing Force Method[END_REF] and integrated on the fly using a Poisson equation formalism, [START_REF] Hénin | Fast and accurate multidimensional free energy integration[END_REF] following the default settings of Colvars. For 2d visualization, free energy surfaces were read into Python notebooks using the colvars_grid Python module, [START_REF] Hénin | Fast and accurate multidimensional free energy integration[END_REF] then integrated over relevant slices of the collective variables using numpy and plotted using matplotlib. For 3d visualization, the DX files written by Colvars were read directly into VMD. [START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF] The collective variables were discretized on a grid of step size 0.2 Å. For eABF, the extended variables were coupled with an RMS deviation of 0.2 Å and a time constant of 200 fs. ABF forces were fully applied when 50 samples were collected in a given bin.

Quantum mechanical (QM) calculations

We selected from the molecular dynamics simulations 20 snapshots in the monodentate and SShIP geometries. We then extracted from each snapshot a cluster composed of one DMP anion, one Mg 2+ , and 20 water molecules, selected based on the closest distance d between their H atoms and the cation or DMP oxygen atoms. All water molecules with d < 3.5 Å were selected, complemented to 20 by selecting additional water molecules by increasing value of d. The number of 20 water molecules ensured that the entire first hydration shell of both Mg 2+ and the DMP oxygen atoms (the water molecules directly donating a hydrogen bond) were selected, including in the SShIP. The geometry of these clusters was optimized with Density Functional Theory (DFT) with the ORCA software, [START_REF] Neese | The ORCA program system[END_REF][START_REF] Neese | Software update: the ORCA program system, version 4.0[END_REF] in the gas phase, using the ωB97X-D3(BJ) functional with Becke-Johnson dispersion correction and the def2-SVP basis set. Single point en-ergies were then obtained with a larger basis set (as in Ref. 60), which proved key to obtain correct interaction energies, both for the full cluster and for the DMP and ion+water fragments. From there, the interaction energy was calculated as E int = E tot -E DMP -E cation+water . Single point energies and the corresponding interaction energies were then calculated on the same optimized fragments with different versions of the AMOEBA force field using the analyze program of Tinker. Detailed data are provided in Supporting Information.

Results

Phosphate-cation binding without explicitly polarizable force fields.

Cations can interact with phosphate groups in distinct binding modes. The interaction can take place either through a water molecule, thus forming a solvent-shared ion pair, or with direct contact between the cation and the DMP oxygen atoms. Such a contact ion pair can be either monodentate (interaction with only one of the two non-bridging phosphate oxygen atoms), or bidentate (Fig. 1). To investigate how different force fields capture the interactions of DMP with divalent cations, Mg 2+ and Ca 2+ , and how they predict the balance between these binding modes, we computed the standard binding free energy ∆G • bind between Mg 2+ or Ca 2+ and DMP in different ion pair geometries with different non polarizable force fields (Fig 2 and Table S2 in the Supporting Information). The overall binding free energy was compared with the experimental reference (available only for DHP, presumed very similar to DMP), -7.1 kJ/mol for Mg 2+ 12 and -5.9 kJ/mol for Ca 2+ . [START_REF] Mcdowell | Solubility Study of Calcium Hydrogen Phosphate. Ion-Pair Formation[END_REF] As discussed already in the literature [START_REF] Duboué-Dijon | Building intuition for binding free energy calculations: bound state definition, restraints, and symmetry[END_REF][START_REF] Luo | On the calculation of absolute macromolecular binding free energies[END_REF][START_REF] Mey | Best Practices for Alchemical Free Energy Calculations[END_REF] the definition of the "bound" state depends on the type of experiment performed, since different techniques are sensitive to different ion pairs. Here, following previous works, [START_REF] Kumar | Structure and thermodynamics of Mg:phosphate interactions in water: A simulation study[END_REF] we assume that the chosen reference free energies report on the formation of both contact and solventshared ion pairs. As we discuss later, this is not the case for spectroscopy-derived estimates.

We also note that some works [START_REF] Cruz-León | Extended magnesium and calcium force field parameters for accurate ion -nucleic acid interactions in biomolecular simulations[END_REF] chose different, weaker, reference values reported by Sigel and Sigel. [START_REF] Sigel | A stability concept for metal ion coordination to single-stranded nucleic acids and affinities of individual sites[END_REF] Those values result from an extrapolation of various data obtained on related, but different, systems. We, as others before, [START_REF] Kumar | Structure and thermodynamics of Mg:phosphate interactions in water: A simulation study[END_REF] thus preferred to use experimental estimates directly obtained on DHP. ∆G • bind for Mg 2+ is strongly overestimated, by more than 20 kJ/mol, with standard non polarizable force fields, which is consistent with the overbinding artefacts reported for Mg 2+ binding to biomolecules in the literature. [START_REF] Grotz | Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation , Ion-Binding , and Water-Exchange Properties[END_REF][START_REF] Duboué-Dijon | A practical guide to biologically relevant molecular simulations with charge scaling for electronic polarization[END_REF] This excessive binding is observed independently of the details of the chosen force field and is robust with respect to a change of water model (see Supporting Information). The overbinding is also observed for Ca 2+ , albeit to a much smaller extent, with the 12-6 force field by Merz and coworkers. [START_REF] Li | Rational design of particle mesh ewald compatible Lennard-Jones parameters for +2 metal cations in explicit solvent[END_REF] Other non polarizable force fields exhibit an overbinding for Ca 2+ almost as pronounced as for Mg 2+ (see Supporting Information). Additionally, with such force fields, both cations are predicted to interact preferentially in direct contact with the DMP non bridging oxygen atoms (rather than through their hydration shell) which, for Mg 2+ , goes against the usual picture of such interactions.

The artefacts exhibited by standard force fields are believed to originate from the lack of electronic polarization, which is particularly important to properly capture binding properties of highly charged cations. One possibility to overcome this limitation is to use a mean-field approach and take implicitly into account electronic polarization by scaling the ionic charges. [START_REF] Duboué-Dijon | A practical guide to biologically relevant molecular simulations with charge scaling for electronic polarization[END_REF][START_REF] Leontyev | Accounting for electronic polarization in non-polarizable force fields[END_REF][START_REF] Kirby | Charge Scaling Manifesto: A Way of Reconciling the Inherently Macroscopic and Microscopic Natures of Molecular Simulations[END_REF] This approach is expected to be transferable to all non-transition metal cations, interacting with any charged anion. It has already been successfully used for instance for Na + , Li + , K + , Ca 2+ , Mg 2+ , and Zn 2+ in interaction with simple anions (e.g. Cl -, Br -) and typical biomolecular charged groups such as carboxylates or phospholipids. [START_REF] Duboué-Dijon | A practical guide to biologically relevant molecular simulations with charge scaling for electronic polarization[END_REF] Following the prescriptions of the Electronic Continuum Correction (ECC) theory, we thus developed a scaled-charge force field both for DMP and the divalent cations, [START_REF] Duboué-Dijon | Hydration and Ion Pairing in Aqueous Mg2+ and Zn2+ Solutions: Force Field Description Aided by Neutron Scattering Experiments and Ab Initio Molecular Dynamics Simulations[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] S2). However, as observed previously with other solutes (e.g. acetate), the overall binding free energy is now underestimated. This led us to prefer a 0.8 scaling factor, which can be justified by the fact that some of the electronic polarization effects are already implicitly incorporated in the water model. [START_REF] Kirby | Charge Scaling Manifesto: A Way of Reconciling the Inherently Macroscopic and Microscopic Natures of Molecular Simulations[END_REF][START_REF] Vega | Simulating water with rigid non-polarizable models: a general perspective[END_REF] With this second version of the ECC force field, the binding free energy for Mg 2+ is now within less than 3 kJ/mol from the experimental value, probably within experimental (and computational) error bars, or slightly underestimated still. The overall binding free energy for Ca 2+ is slightly smaller than that of Mg 2+ , by about 2 kJ/mol, consistently with experimental data. Strikingly, compared to standard non polarizable force fields, in addition to a much weaker binding strength, the favoured binding mode is now the SShIP, for both cations. The contact binding mode is found slightly less stable, by 2-4 kJ/mol, the monodentate geometry being favored over bidentate for both cations. The bidentate binding mode is found unstable with Mg 2+ (it is not associated with a local free energy minimum and cannot be stabilized even with restraints), while it is metastable with Ca 2+ . The free energy difference between the monodentate and solvent-shared binding modes is very modest, so that the relative ordering of these two binding modes can be easily switched by changing, for instance, the water model (see discussion in Supporting Information). Another approach to eliminate overbinding artefacts is to parametrize pair-specific Lennard Jones parameters between the studied cations and the solute they interact with. [START_REF] Grotz | Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation , Ion-Binding , and Water-Exchange Properties[END_REF][START_REF] Fyta | Ionic force field optimization based on single-ion and ion-pair solvation properties: Going beyond standard mixing rules[END_REF][START_REF] Yoo | New tricks for old dogs: Improving the accuracy of biomolecular force fields by pair-specific corrections to non-bonded interactions[END_REF][START_REF] Kashefolgheta | Developing force fields when experimental data is sparse: AMBER/GAFF-compatible parameters for inorganic and alkyl oxoanions[END_REF][START_REF] Hess | Cation specific binding with protein surface charges[END_REF] Recently, parameters were developed specifically to capture Mg 2+ and Ca 2+ binding free energies with phosphate groups. [START_REF] Grotz | Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation , Ion-Binding , and Water-Exchange Properties[END_REF][START_REF] Cruz-León | Extended magnesium and calcium force field parameters for accurate ion -nucleic acid interactions in biomolecular simulations[END_REF] Since ∆G • bind was used as a target in the parametrization process-albeit with a different reference than ours for Ca bind (kJ mol -1 ) for different ion pair conformations (contact monodentate, contact bidentate and SShIP), together with the overall binding free energy ∆G • bind,tot (Eq. ??), computed either over only the contact binding modes, "total contact", or also including the solvent-shared ion pair, "total". The binding free energies were computed with different force fields, as described in the Methodology section. The specific combinaison of force fields used in each case is specified in Table 1. The reference experimental binding free energies ∆G • bind,expt are -7.1 kJ/mol for Mg 2+ 12 and -5.9 kJ/mol for Ca 2+ . [START_REF] Mcdowell | Solubility Study of Calcium Hydrogen Phosphate. Ion-Pair Formation[END_REF] The light red zone indicates the range of binding free energy for the contact ion pairs as estimated from Ref. 15. Unstable ion pair geometries, which are not associated with a local free energy minimum and cannot be stabilized even with restraints, are marked with a * symbol. recently derived for Mg 2+ , [START_REF] Cruz-León | Extended magnesium and calcium force field parameters for accurate ion -nucleic acid interactions in biomolecular simulations[END_REF] that binds more strongly than Ca 2+ (more details in Supporting Information). The relative binding of the two cations thus sensitively depends on details in the parametrization. Here, we decomposed the binding free energy into the contribution of the different ion pairs (Fig 2 and Table S2). As with the ECC force field, the most stable bind-ing mode is the SShIP, the contact monodentate ion pair being less stable by 1-2 kJ/mol. The contact bidentate ion pair, which, for Ca 2+ , was only a few kJ/mol less stable than the contact monodentate with ECC, is fully unstable (i.e. is not associated with a local free energy minimum) with this description for both cations. This probably comes from the repulsion caused by the pair-specific Lennard-Jones term.

Figure 3: Illustration of the monodentate and bidentate binding modes for the Mg 2+ :DMP ion pair, represented as an iso-free-energy surface computed using eABF for the full charge model. Prepared with VMD. [START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF] In order to get a better insight into the geometry and relative stability of the different binding modes with the different force fields, we computed, using the eABF enhanced sampling technique, [START_REF] Lelièvre | Computation of free energy profiles with parallel adaptive dynamics[END_REF][START_REF] Lesage | Smoothed Biasing Forces Yield Unbiased Free Energies with the Extended-System Adaptive Biasing Force Method[END_REF] the 3D free energy landscape around the DMP anion. 3D visualization of free energy isosurfaces provides a picture of the spatial extent of each of the binding modes (Figure 3), with symmetric lobes around each oxygen for the contact monodentate ion pair, and a thin lobe, orthogonal to the OPO plane, for the contact bidentate (when stable).

Projecting the 3D free energy data on the plane formed by the phosphorus central atom and the two non bridging oxygens, we obtained 2D free energy maps that enable easy comparison of the binding strength of each mode with the different force fields (Figure 4). While the overall shape of the binding modes is similar with all force fields, clear differences are visible in the binding strength. Full-charges non polarizable force fields yield much deeper free energy basins than the two implicit descriptions; and ECC, even with 0.8 scaling, yields the most shallow free energy minima, in line with the alchemical evaluation of ∆G • bind . This representation also confirms the absence of a metastable contact bindentate state (x = 0) with the use of pair-specific Lennard-Jones parameters, which contrasts with both ECC and full charge descriptions, where the contact bidentate state, albeit the least stable with ECC, indeed corresponds to a metastable state with Ca 2+ . In addition, even when the contact and SShIP basins are close in free energy, they can be separatede.g. with Mg 2+ and the ECC force field-by a high (> 15 kJ/mol) free energy barrier, which means that the kinetics of exchange remains slow between the states and won't be spontaneously sampled in standard non accelerated molecular dynamics simulations.

While scaled ECC force fields and the use of pair-specific LJ parameters are able to reproduce the overall binding free energy and both predict the SShIP to be the lowest free energy minimum, they exhibit differences regarding, for instance, the existence of a metastable bidentate binding mode. Since these force fields only implicitly account for electronic polarization, we now set out to compare their behavior with that of two different explicitly polarizable force fields, the Drude polarizable force field and AMOEBA.

Explicitly polarizable force fields

Drude polarizable force field. The Drude polarizable force field [START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF] takes into account electronic polarization through the Drude oscillator model, with negatively charged Drude particles fluctuating around their parent atom during the simulation. Lemkul and MacKerell 56 parametrized a Drude polarizable force field for DMP and Mg 2+ specifically designed to capture Mg 2+ hydration thermodynamics and structure, as well as its interactions with chloride anions and nucleic acid moieties. The inclusion of pair-specific LJ interactions (NBFIX) between Mg 2+ and non-bridging phosphate oxygen atoms, in addition to the Drude particles, was found necessary to reproduce the structure of Mg 2+ -phosphate (DMP) complexes. In addi- Figure 4: 2D free energy maps for the interaction of a Mg 2+ or Ca 2+ cation with dimethylphosphate, computed using eABF, integrated over a 2 Å-thick slice around the y = 0 OPO plane, and capped at ±20 kJ/mol. Rendered with matplotlib. [START_REF] Hunter | Matplotlib: A 2D graphics environment[END_REF] tion, Villa et al. [START_REF] Villa | Classical Drude Polarizable Force Field Model for Methyl Phosphate and Its Interactions with Mg2+[END_REF] refined in a similar way (using NBFIX terms), a Drude force field for the interactions of Mg 2+ cations with methylphosphate. They showed that it quantitatively reproduced the overall binding free energy and predicted the SShIP to be the most stable, the contact monodentate ion pair being significantly less stable by about 9.5 kJ/mol. Using the specifically parametrized force field for DMP and Mg 2+ , [START_REF] Lemkul | Balancing the Interactions of Mg2+ in Aqueous Solution and with Nucleic Acid Moieties For a Polarizable Force Field Based on the Classical Drude Oscillator Model[END_REF] we computed the DMP:Mg 2+ binding free energy for each binding mode (Fig 2 and Table S3). The overall binding free energy (-10.7 kJ/mol) is in good agreement with the experimental reference, and mainly stems from the solvent-shared binding mode, the contact monodentate mode being about 11 kJ/mol high in free energy. These results are fully in line with previous data obtained on methylphosphate. [START_REF] Villa | Classical Drude Polarizable Force Field Model for Methyl Phosphate and Its Interactions with Mg2+[END_REF] Note that the contact bidentate ion pair is not stable at all with this force field, presumably due to the NBFIX repulsion term, as is also the case with the non polarizable force field with pair-specific LJ terms (Fig 2 and Table S3).

The Drude polarizable force field for Ca 2+phosphate interactions has not yet been refined with NBFIX terms, as for Mg 2+ . For the sake of comparison, we still computed the binding free energy in the different modes using the default Ca 2+ force field. While the SShIP is associated with a reasonable binding free energy (-6.5 kJ/mol), slightly higher than that of Mg 2+ , which would be in line with experimental data, the contact ion pair is clearly misbehaving, both in terms of structure and energetics (see Supporting Information, Table S3). The monodentate binding mode is unstable, immediately falling in a pathological bidentate conformation, with an extremely small 2.5-2.9 Å Ca 2+ -phosphorus distance, to be compared with a reference 2.95-3.05 Å distance estimated on 10 clusters (see Methods section) optimized with DFT at the def2-SVP/wB97x-D3(BJ)/PCM level. This kind of deformation is in line with what was observed with Mg 2+ before the inclusion of the NBFIX term. [START_REF] Lemkul | Balancing the Interactions of Mg2+ in Aqueous Solution and with Nucleic Acid Moieties For a Polarizable Force Field Based on the Classical Drude Oscillator Model[END_REF] Unsurprisingly, the binding free energy corresponding to this bidentate mode is also totally unphysical, leading to an overall ∆G • bind = -73.2 kJ/mol. These results stress that without addition of pair-specific Lennard Jones repulsion terms, the Drude polarizable force field does not properly capture divalent ion binding with biologically relevant moieties, hence the need for careful testing before application to a new family of interaction groups. AMOEBA. AMOEBA is another family of polarizable force field which uses permanent electrostatic multipole moments up to the quadrupole and includes explicit dipole polarization. [START_REF] Ponder | Current status of the AMOEBA polarizable force field[END_REF] It has been parametrized for proteins [START_REF] Shi | Polarizable atomic multipole-based AMOEBA force field for proteins[END_REF] and, more recently, for nucleic acids, [START_REF] Zhang | AMOEBA Polarizable Atomic Multipole Force Field for Nucleic Acids[END_REF] which now enables its use for simulation of complex biochemical systems. Kumar et al. [START_REF] Kumar | Structure and thermodynamics of Mg:phosphate interactions in water: A simulation study[END_REF] studied the interactions of Mg 2+ ions with a model phosphate compound (dihydrogen phosphate, DHP), and showed that it was able to quantitatively reproduce the binding free energy. The most stable ion pair conformation was found to be the SShIP, with both contact ion pairs (monodentate and bidentate) very close in free energy, about 7 kJ/mol less stable than the solvent shared. More recently, Zhang et al. [START_REF] Zhang | Polarizable multipole-based force field for dimethyl and trimethyl phosphate[END_REF] derived AMOEBA parameters for dimethylphosphate, which later served as a basis for the nucleic acid backbone in the AMOEBA nucleic acid force field. [START_REF] Zhang | AMOEBA Polarizable Atomic Multipole Force Field for Nucleic Acids[END_REF] This DMP force field was validated for its interactions with Mg 2+ against the structure and interaction energy of cation:DMP dimer, in a bidentate geometry optimized at the QM (DFT) level.

With this recent AMOEBA force field for DMP, which should be a good proxy for the RNA/DNA interactions with divalent cations with the AMOEBA force field for nucleic acids, we examined the binding free energy associated with each type of ion pair (Fig 2 and Table S3). The overall Mg 2+ -DMP binding (∆G • bind = -20.8 kJ/mol) appears significantly too strong compared to experimental data, especially for the contact monodentate ion pair, which is the most stable. Meanwhile, the SShIP exhibiting a binding free energy much more similar to that obtained with other force fields (∆G • bind,sship = -7.8 kJ/mol). This contrasts with the results previously obtained on DHP, both quantitatively and qualitatively.

We thus endeavor to better understand the origin of this difference, analyzing how sensitive the binding free energy was to different force field parameters. First, we derived new multipoles for DMP (see Methods), taking into account explicit water molecules in the derivation. However, the binding free energy ∆G • bind proved very little sensitive to the change of multipoles (see Supporting Information). We then noticed that a significant difference between the early DHP force field [START_REF] Kumar | Structure and thermodynamics of Mg:phosphate interactions in water: A simulation study[END_REF] and the current DMP force field [START_REF] Zhang | Polarizable multipole-based force field for dimethyl and trimethyl phosphate[END_REF] was in the value of the atomic polarizability of the phosphorus and non bridging oxygen atoms. For DMP, they were tuned to reproduce the overall DMP polarisability. [START_REF] Zhang | Polarizable multipole-based force field for dimethyl and trimethyl phosphate[END_REF] Computing the binding free energy for different polarizability values (see Supporting Information) evidenced that the binding free energy of the monodentate binding mode was extremely sensitive to the non-bridging oxygen polarizability and very little to that of the phosphorus. In contrast, the SShIP binding mode was very lit-tle affected by changes in both polarizabilities.

We then went back to the spirit of the initial development of the AMOEBA force field with the idea of benchmarking some parameters against QM calculations on the structures and energies of small clusters. The polarizability values were thus adjusted on clusters in the monodentate and SShIP binding modes (details in the Supporting Information). We found out that the set of polarizabilities that minimize the error on relative interaction energies (α O = 0.837 Å 3 and α P = 3.3 Å 3 , close to those used in the early DHP force field [START_REF] Kumar | Structure and thermodynamics of Mg:phosphate interactions in water: A simulation study[END_REF] ) yield binding free energies for both Mg 2+ and Ca 2+ in very good agreement with experimental data (see Fig 2 and Table S3), slightly more negative with Mg 2+ than with Ca 2+ . In both cases, the SShIP is the most stable ion pair, and the contact monodentate ion pair about 6-10 kJ/mol higher in energy. The contact bidentate ion pair is not stable at all, which is different from what was earlier observed with DHP, but in agreement with the picture provided by most of the other tested force fields. Refinement of the DMP polarizabilities thus had a major impact on binding free energies and led to an ion binding behavior in line with experimental data as well as with simulations with other force fields. Whether a similar change in large scale simulations of nucleic acids would improve their cation binding behavior without deteriorating the properties of the nucleic acid remains to be investigated.

Discussion

The interactions of DMP with Mg 2+ and Ca 2+ were examined in details with different force fields that use distinct, implicit or explicit, approaches to take into account electronic polarization. We show that all approaches, after careful refinement, are able to capture the overall binding free energy. From a structural point of view, the obtained binding distances (first peak of the radial distribution function between the DMP non bridging oxygen and the cations, see Table S3 andS4) are all compatible with experimental estimates [START_REF] Kolev | Interaction of Na+, K+, Mg2+ and Ca2+ counter cations with RNA[END_REF] (1.92-2.14 Å for Mg 2+ and 2.2-2.5 Å for Ca 2+ ). From there, a common picture emerges for divalent cation-phosphate interactions: the binding is overall quite moderate, slightly weaker for Ca 2+ than Mg 2+ , and the (slightly) most stable binding mode, for both cations, seems to be the SShIP (outersphere). As for the relative stability of the two contact binding modes, monodentate and bidentate, all force fields favor the monodentate binding mode with both cations. However, whether the contact bidentate binding mode can still form a minority population is still an open question. For Mg 2+ , the contact bidentate binding mode is associated with a local free energy minimum with none of the tested force fields (except for the standard non polarizable full charge, where the strong contact binding is clearly artefactual). This tends to suggest that this binding mode may not exist, which contrasts with the conclusion of an earlier simulation study. [START_REF] Kumar | Structure and thermodynamics of Mg:phosphate interactions in water: A simulation study[END_REF] For Ca 2+ , our results are less conclusive, as a (meta)stable state is found with ECC force fields, but does not exist with the tested pair-specific LJ parameters, nor with the AMOEBA polarizable force field. Hence, we cannot fully rule out the existence of this binding mode with Ca 2+ , which is also difficult to probe experimentally and to distinguish from the monodentate contact ion pair.

In addition, the different tested force fields differ in the balance between SShIP and contact ion pairs. With both implicit approaches (ECC and pair-specific LJ) and for both cations, the contact monodentate ion pair is only slightly (2-4 kJ/mol) less stable than the SShIP. In contrast, with the Drude polarizable force field and AMOEBA, the contact monodentate ion pair is more destabilized with respect to the SShIP (by about 6-10 kJ/mol, more for Mg 2+ than Ca 2+ ), which means that contact ion pairs would be much less abundant.

Very recently, using 2D-infrared spectroscopy experiments on DMP Mg 2+ or Ca 2+ solutions complemented by QM simulations,Schauss et al. [START_REF] Schauss | Contact Ion Pairs of Phosphate Groups in Water: Two-Dimensional Infrared Spectroscopy of Dimethyl Phosphate and ab Initio Simulations[END_REF] suggested the existence of contact ion pairs between DMP and these divalent cations, making use of the fact that the blueshift observed in IR measurements are only caused by the formation of contact ion pairs, and not with solvent-shared. They estimated the con-centration of DMP:cation contact ion pairs in C 0 = 0.2 M solutions between 0.23C 0 and 0.37C 0 for Mg 2+ , and between 0.12C 0 and 0.3C 0 for Ca 2+ . These estimates can be compared with that provided by simulations with different force fields. Direct counting in simulations run at the same concentration as in the experiment is not practical because of the very long exchange times between the different ion pairs, which prevents proper convergence of the contact ion pair concentration. Instead, we can convert the estimated contact ion pair concentrations into the binding free energy associated with contact pairs, ∆G bind,contact . This thus provides an estimated ∆G bind,contact = -7.8 to -4.3 kJ/mol for Mg 2+ and ∆G bind,contact = -5.9 to -1.5 kJ/mol for Ca 2+ . ∆G bind,contact can be directly compared with that computed with different force fields, from the exponential average of the ∆G bind associated with the two contact pairs (see Fig. 2). All force fields (among those giving a correct overall binding free energy) seem to underestimate the fraction of contact ion pairs. The contact ion pair binding free energy is slightly underestimated by about 2 to 4 kJ/mol with both ECC and pair-specific LJ descriptions, except for Ca 2+ , where the pair specific LJ description falls within experimental error bars. Both polarizable force fields (Drude and our modified AMOEBA) underestimate even more the contact ion pair free energy, by at least 4.5 kJ/mol for Mg 2+ .

Overall, implicit and explicit descriptions of electronic polarization yield quite similar pictures of the ion pairing behavior of phosphate groups with Mg 2+ and Ca 2+ . However, with the Drude polarizable force field, correct ion binding behavior crucially depends on added pairspecific interaction; with AMOEBA, it is highly sensitive to the DMP atomic polarizabilities of the non-bridging oxygen atoms. Both implicit strategies (ECC and pair-specific LJ) can be directly used with standard simulation software at no additional cost, which makes them ideally suited for large-scale simulations of biological assemblies. However, they also come with their drawbacks. Pair-specific Lennard Jones parameters lack transferability, as they need to be parametrized specifically for each pair of interacting ions and the relative behavior the two cations proved quite sensitive to details in the parametrization. This approach compensates the short-range overbinding artefacts by an ad hoc repulsive term, thus recovering a correct binding behavior, but does not fix the missing electronic term in long-range electrostatics. In contrast, the scaled-charge ECC force fields introduce the missing electronic polarizability in a physically well-grounded mean-field approach. However, this strategy requires scaling the charges of both the cations and the anionic solute, which means that in biosimulations of nucleic acids, the charges of the RNA or DNA moiety should also in principle be scaled. This will affect the nucleic acid-water interactions and could alter the biomolecule conformation. Additional developments and tests are thus needed to be able to use such force fields for more complex systems. Also, for non polarizable force fields, we focused here on the cation force fields and did not explore how modifying the DMP force field (beyond a uniform charge scaling factor) would impact the binding free energy, which should be explored in future works.

Finally, even if we don't explicitly examine the molecular details of the conversion between inner-sphere (contact) and outer-sphere (SShIP) binding modes, as has been done elsewhere, [START_REF] Neumann | Artificial Intelligence Resolves Kinetic Pathways of Magnesium Binding to RNA[END_REF] we can see that even with the force fields that properly capture the binding free energy, exchange kinetics between the different ion pairs remains slow (compared to typical simulation timescales). This is evidenced, for instance, by the large barriers separating the contact and solvent-shared binding modes in Fig. 4. Hence, special care should always be taken to properly sample the distinct binding modes.

Conclusion

While standard non-polarizable force fields cannot properly capture the binding structure and free energy of divalent cations to phosphate moieties, we have shown that force fields that include electronic polarization in an implicit mean-field way-either through the use of scaled charges (ECC theory) or pair-specific Lennard-Jones parameters-successfully reproduce the modest binding free energy of Mg 2+ and Ca 2+ to a model isolated phosphate compound, DMP. Similar results are obtained with explicitly polarizable force fields, provided extra care is taken during the parametrization process, with the addition of pair-specific repulsion for the Drude polarizable model, and fine tuning of the atomic polarizabilities for AMOEBA. From the comparison of these different approaches, a common picture of cationphosphate binding emerges: with both Mg 2+ and Ca 2+ cations, the most stable binding mode is the SShIP, and among the two contact geometries, the monodentate binding mode seems the most stable, even though the various force fields differ in the relative stability of the monodentate and solvent shared binding modes. Comparison with recent experimental estimates of the fraction of contact ion pairs however suggests that all the force fields may underestimate the stability of the contact binding mode. Given that force fields which only implicitly take into account electronic polarization yield very similar descriptions of the cationphosphate binding modes and binding free energies as explicitly polarizable force fields, these approaches seem very promising for simulation of large scale biochemical assemblies with divalent cations. The two strategies however pose their own challenges, the lack of transferability for the pair-specific LJ approach, and the scaling of the biomolecular charges for ECC. It also remains to be examined whether those different approaches are able to describe the formation of strong binding sites, as found in RNA for instance, that involve multiple phosphate groups.

Figure 1 :

 1 Figure1: Molecular representation of the different cation-DMP ion pairing modes with the surrounding layer of water molecules. The Pcation distance can be used to distinguish the different ion pairs, with, for Mg 2+ , typical values of 2.7-3.2 Å for the contact bidentate binding mode, 3.2-2.7 Å for the contact monodentate and 4.1-6.1 Å for the SShIP. Snapshots were prepared using the VMD software.[START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF]

Figure 2 :

 2 Figure 2: Cation-DMP binding free energy ∆G •bind (kJ mol -1 ) for different ion pair conformations (contact monodentate, contact bidentate and SShIP), together with the overall binding free energy ∆G • bind,tot (Eq. ??), computed either over only the contact binding modes, "total contact", or also including the solvent-shared ion pair, "total". The binding free energies were computed with different force fields, as described in the Methodology section. The specific combinaison of force fields used in each case is specified in Table1. The reference experimental binding free energies ∆G • bind,expt are -7.1 kJ/mol for Mg 2+ 12 and -5.9 kJ/mol for Ca 2+ .[START_REF] Mcdowell | Solubility Study of Calcium Hydrogen Phosphate. Ion-Pair Formation[END_REF] The light red zone indicates the range of binding free energy for the contact ion pairs as estimated from Ref.15. Unstable ion pair geometries, which are not associated with a local free energy minimum and cannot be stabilized even with restraints, are marked with a * symbol.

Table 1

 1

	summarizes the combination of force
	fields used in each simulation with non explic-
	itly polarizable force fields.
	Table 1: Combination of force fields used
	in each simulation with non explicitly po-
	larizable force fields, using labels for each
	force fields as defined in the text, where
	we provide the corresponding references.
	Label in Fig. 2 Cation	water	DMP
	Mg 2+
	full charges	Mg_m SPC/E GAFF
	ECC 0.75	Mg_s	SPC/E GAFF
	ECC 0.8	Mg_s2 SPC/E GAFF
	pair LJ	microMg TIP3P GAFF
	Ca 2+
	full charges	Ca_m	SPC/E GAFF
	ECC 0.75	Ca_s	SPC/E GAFF
	ECC 0.8	Ca_s2 SPC/E GAFF
	pair LJ	Ca2	TIP3P GAFF
	Drude polarizable force field. The Drude
	polarizable force field 26 takes into account elec-
	tronic polarization through the Drude oscil-
	lator model, with negatively charged parti-
	cles attached to each heavy atom via har-
	monic springs.

Acknowledgement This work was supported by the "Initiative d'Excellence" program from the French State (Grants "DYNAMO", ANR-11-LABX-0011, and "CACSICE", ANR-11-EQPX-0008). Computational work was performed using local HPC resources managed by G. Letessier, as well as HPC resources from GENCI (project A0090711021). JP ac-knowledges a PhD grant from the "DYNAMO" Labex (ANR-11-LABX-0011). We thank Dr. J. Lemkul (Virginia Tech) and Prof. Thomas Simonson (Ecole Polytechnique) for their help to set up simulations and alchemical transformations with the Drude polarizable force field, and Hubert Santuz (CNRS, LBT) for help with interfacing Colvars with Gromacs.

Supporting Information Available

Supporting Information is available, with additional computational details, and complementary analyses. Input and parameter files are shared on a Zenodo public folder (https://zenodo.org/ record/5971002#.YkrBDDw6-YA).

TOC Graphic