# Implicit Electronic Polarization via Charge Scaling in Amber OLx Nucleic Acid Force Fields Improves Ion Pairing without Backbone Reparametrization

Julie Puyo-Fourtine, Marie Zgarbová, Laetitia Kantin, Petr Jurečka, Elise Duboué-Dijon

► To cite this version:

Julie Puyo-Fourtine, Marie Zgarbová, Laetitia Kantin, Petr Jurečka, Elise Duboué-Dijon. Implicit Electronic Polarization via Charge Scaling in Amber OLx Nucleic Acid Force Fields Improves Ion Pairing without Backbone Reparametrization. Journal of Chemical Theory and Computation, 2025, 10.1021/acs.jctc.5c01262. [DOI: 10.1021/acs.jctc.5c01262](https://doi.org/10.1021/acs.jctc.5c01262)

HAL Id: hal-05227535

https://hal.science/hal-05227535v1

Submitted on 28 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Implicit electronic polarization via charge
scaling in Amber OLx nucleic acid force fields
improves ion pairing without backbone
reparametrization

Julie Puyo-Fourtine,† Marie Zgarbová,‡ Laetitia Kantin,† Petr Jurečka,‡ and Elise

Duboué-Dijon*,†

†Université Paris Cité, CNRS, Laboratoire de Biochimie Théorique, 13 rue Pierre et Marie
Curie, 75005, Paris, France

‡Department of Physical Chemistry, Faculty of Science, Palacky University, 77146
Olomouc, Czech Republic

E-mail: elise.duboue-dijon@cnrs.fr
Abstract

While divalent ions are known to be involved in key biological processes such as RNA
folding or DNA-histone interactions, these interactions are poorly captured in molecular
dynamics simulations with empirical force fields, which suffer from strong overbinding
artefacts. Hence, there is a high need for improved descriptions of (divalent) ions in
nucleic acid simulations. In this work, we explore the possibility to improve ion-binding
properties of the popular Amber-OL15 force field using the Electronic Continuum Cor-
rection (ECC) approach, which includes electronic polarization through charge scaling,
limited here to the phosphate backbone. This strategy yields very promising results,
with essentially no degradation of the conformational properties of selected DNA (and
a ds-RNA) sequences, and a strong improvement of both monovalent ion retention in
G-quadruplexes and divalent ion pairing. As the ECC modification appears mostly or-
thogonal to force field refinements focused on backbone dihedral parameters, this work
suggests a systematic way to improve ion pairing properties of nucleic acids in all-atom
MD simulations.

**Introduction**

Nucleic acids—DNA (deoxyribonucleic acid) and RNA (ribonucleic acid)—are long nega-
tively charged biopolymers. In a biological context, their phosphate backbone charge is
screened by countercations,¹ which affect their structure and dynamics,²,³ a prominent ex-
ample being the shift of DNA duplexes from a B-DNA helix at physiological salt conditions to
an A-form at high salt concentration.⁴ Cations interact loosely with nucleic acids—forming
the so-called ion atmosphere¹,³,⁵—but specific ion binding sites have also been characterized,
with both monovalent²,⁶⁻¹⁰ and divalent cations such as Ca²⁺ and Mg²⁺.¹¹⁻¹⁵ Such divalent
ions play a key role in a number of biological processes, e.g., in the formation of DNA–protein
or DNA–drug complexes,¹⁶ or RNA folding.¹⁷⁻²⁰ Mg²⁺ ions also modulate DNA compaction
in the chromatin, mediate histone–DNA interactions,²¹ and are associated with the catalytic
activity of DNA- and RNA-processing enzymes²²,²³ and ribozymes.²⁴

Molecular dynamics investigations of these fundamental biological processes thus de-
mands a proper description of ion-nucleic acid interactions in molecular simulations. Com-
monly employed nucleic acid force fields (typically of the Amber or CHARMM families)
have been extensively refined over the years to properly capture base-base or base-water
interactions as well as conformational properties of nucleic acids (i.e. equilibrium between
backbone substates),²⁵⁻²⁸ but do not guarantee an accurate description of the cation-nucleic
acid interactions, with a typical overbinding of cations next around the negatively charged
phosphate groups.²⁹ Yet, with monovalent ions, combining state-of-the-art DNA force fields
with various ion force fields has been shown not only to give a good description of DNA
duplex conformations but also to capture fine details of Na⁺ or K⁺ locations within the
grooves, as identified in high-resolution crystal structures,³⁰⁻³³ even if problems may arise at
higher ion concentrations.

In contrast, it is well acknowledged in the literature³⁴,³⁵ that divalent ions, such as Ca²⁺ and Mg²⁺, and their interactions with nucleic acids remain beyond the reach of usual fixed-charge force fields. This is due to large polarization and charge transfer effects that are not included in standard force fields, where ions are described as van-der-Waals spheres with a fixed point charge. The long exchange times for water and ligands in the cations' hydration sphere—> µs for Mg²⁺³⁶,³⁷—pose an additional sampling challenge for simulations.³¹ Plain MD simulations will often not sample the most thermodynamically stable arrangements of ions around nucleic acids but instead be trapped in local minima with high kinetic barriers to overcome.

There is thus a strong need for improved descriptions of (divalent) cations in nucleic
acid simulations, which is reflected in the large body of recent work in this direction.²⁹,³⁸⁻⁴⁴
Divalent ions van der Waals parameters—either the traditional 12-6 Lennard-Jones potential,
or more recently, a 12-6-4 Lennard-Jones variant⁴¹,⁴⁴—are often tuned on thermodynamic
and structural bulk solution properties such as the solvation free energy, position and height
of the first peak of cation-water radial distribution functions. However, several works have shown that even divalent ions with good bulk properties do not properly describe binding to nucleic acid sites<sup>45</sup>—and, more generally, charged biomolecular groups—even if 12-6-4 potentials could be refined to better capture cation binding to phosphate groups.<sup>44</sup> These artifacts of fixed charge force fields are thought to originate from the lack of electronic polarization—which should screen ion-ion interactions—and charge transfer in common MD simulations. The development of polarizable force fields for nucleic acids<sup>46-48</sup> offers a first solution for ion-rich nucleic acid systems. For instance, the Drude polarizable force field offers a Mg<sup>2+</sup> force field that can reproduce both the thermodynamic and kinetic properties of Mg<sup>2+</sup> hydration,<sup>48</sup> and yields an improved description of G-quadruplexes.<sup>49</sup> However, despite these significant advances and recent developments, including key software developments to limit their computational cost,<sup>50</sup> polarizable force fields are not yet broadly used by the biosimulation community. Among possible reasons, we can mention a higher number of parameters, a more tedious parametrization procedure, the requirement of specialized software with smaller user bases, and their slower computational performance which hampers their application to large biomolecular systems.

Several alternative approaches have thus been developed to take into account electronic polarization "implicitly", incorporating the effects of electronic polarization within the standard functional form of non-polarizable force fields, through a modification of either charges or van der Waals parameters. These approaches allow to effectively correct ion-ion interactions without any additional computational cost with respect to standard force fields.

In particular, a strategy adopted by several groups to improve the ion pairing behavior with nucleic acids is the use, instead of the usual combination rules, of pair-specific Lennard-Jones parameters for key interactions, and in particular for the ion-phosphate oxygen pair.<sup>29,40,43</sup> These pair-specific terms—usually referred to as Non-Bonded FIX (NBFIX) corrections as in the CHARMM package—can be adjusted on model moieties to reproduce ion-pair sensitive experimental data, such as osmotic coefficients<sup>29,51</sup> or ion binding free en-
ergies.40,43 For Mg2+, the force field by Yoo and Aksimentiev29 assumes no exchanges in
the Mg2+ hydration shell, and thus is not designed to describe the equilibrium between solvent-
mediated (outer-sphere) and direct contact (inner sphere) ion pairing with phosphates, which
limits its use to understand the formation of contact ion pairs, common in RNA systems.
In contrast, the parametrizations by Panteva et al.40 and Grotz et al.43 are designed to
accurately capture the formation of contact ion pairs. Overall, the NBFIX strategy has been
shown to successfully remove the artificially strong pairing of divalent ions to nucleic acid
phosphate groups and to yield a much improved description of the ion atmosphere around
nucleic acids, as validated against, for instance, ion counting experiments52 or the location
of inner-sphere Mg2+ binding sites in the add A-riboswitch.43 Despite these important suc-
cesses, the NBFIX approach relies on pair-specific interactions, so that specific parameters
need to be defined for the interaction of each ion with each possible interaction site, with
no a priori transferability of parameters between sites. Parametrization could thus proove
tedious in large systems containing different biomolecules (e.g. proteins, sugars, phospho-
lipids), as a Mg2+ ion optimized with NBFIX to interact with DNA phosphates would not be
expected to have good interaction properties with, for instance, protein carboxylate groups.

In this work, we thus proceed to investigate an alternative approach, the so-called Elec-
tronic Continuum Correction (ECC),53,54 that relies on charge scaling instead of tuning van
der Waals interactions. The idea underlying this physically well-grounded mean-field ap-
proach is to treat the fast electronic part of the polarization as a dielectric continuum. This
is equivalent to scaling all ionic charges, with a physics-derived prescription for the scaling
factor of 1/√ε∞ ≃ 0.75, where ε∞ is the fast electronic part of water dielectric response (the
slow component is already present and well-described in MD simulations). ε∞ (in contrast
to εr) presents only small variations in condensed phases, so that pragmatically a single
scaling factor can be used for all biomolecular systems. Recent theoretical work has however
demonstrated that a small variation of the optimal scaling factor (from 0.75 to 0.78-0.83)
can be expected in the case of contact ion pairs,55 which suggests that an optimal choice
of scaling factor would be somewhere in the range 0.75-0.85 and may depend in part on
the employed water model, as observed and discussed in previous works.<sup>54,56–58</sup> The ECC
approach has been successfully tested in various salt solutions<sup>56,59,60</sup> and biological systems,
including proteins and phospholipids.<sup>61–64</sup> Since the ECC strategy has been previously shown
to strongly improve the description of ion pairing in electrolytes with singly<sup>58</sup> or multiply
charged phosphates,<sup>65</sup> we now want to extend it to nucleic acids, DNA and RNA. As the
charged moiety in nucleic acids is directly located on the phosphate backbone, a key question
is whether charge modification of the phosphate group affects conformational properties and
thus, whether it implies the need for extensive force field parametrization. Note that charge
scaling is already used—with a 0.9 charge scaling factor—in the latest DES-Amber force
fields for proteins and nucleic acids,<sup>66,67</sup> where it was found essential to reach a good fit of
target osmotic coefficients used for parametrization.

Our ambition here is not to develop from scratch a new ECC force field for nucleic
acids, but rather to investigate whether we can improve ion pairing properties of the state-
of-the-art nucleic acid force fields by applying the ECC charge-scaling prescription to the
nucleic acid phosphate backbone—the main interaction site for cations—and counterions,
without degrading their overall performance for non ion-specific structural properties. We
consider this an exploratory evaluation of whether the ECC strategy can be applied on top
of existing force fields without the need for complete reparameterization. We will use as
a starting point the popular Amber-OL15 nucleic acid force field,<sup>68</sup> and design an ECC
variant with local charge modifications on the phosphate backbone to ensure more accurate
ion binding with both monovalent and divalent ions. We then demonstrate that this ECC
force field variant does not yield unwanted degradation of conformational behavior of various
DNA and RNA sequences, including challenging non-canonical Z-DNA and G-quadruplex
structures in addition to canonical double helical duplexes. We further characterize how
much our ECC variant affects ion pairing with monovalent ions and show preliminary results
with divalent ions, where sampling poses an additional challenge.
# Methods

## Molecular simulations of a model compound for phosphate backbone: DMP

Following common practice in the literature,⁴⁶,⁶⁷ the dimethylphosphate anion (DMP) was used as a simple model for the nucleic acid backbone to evaluate the impact of different charge distributions on ion pairing properties.

## Force fields

Our strategy was to use DMP as a proxy for nucleic acid backbones, so we started with building a full charge DMP force field that would have the same phosphate moiety as the Amber-OL15 DNA backbone. To this end, we designed a DMP force field ("DMPamb") that mimicks the atom types and charges from the Amber-OL15⁶⁸ force field for DNA (see Table 1). The OMe oxygens were assigned the average charge of the DNA O3' and O5' atoms, and the methyl carbon atoms were assigned the averaged charge of the DNA C3' and C5' atoms. The overall charge was finally adjusted to -1 equally dividing the remaining charge between the 6 methyl hydrogen atoms.

A first "naive" ECC version of the Amber-OL15 force field for DNA, "OL15-ECC-1", was derived by applying a 0.8 charge scaling factor—that seemed to lead to the best ion pairing properties in previous work— to all 5 atoms of the backbone phosphate group (P, O1P, O2P, O3' and O5'). As the phosphate group does not have a charge of exactly -1 in the original force field, this strategy does not yield a -0.8 charge for the nucleotide, as desired. The total charge was then adjusted to -0.8 by distributing the excess charge equally between O3' and O5' atoms, which are presumably less essential for ion binding properties than the non bridging oxygen atoms. The resulting charges are summarized in Table 1. A DMP proxy for this DNA force field, "DMPecc1", was then built following the same strategy as for the standard OL15 force field (see Table 1).
A second variant of the ECC force field for DNA ("OL15-ECC-2" in Table 1 and Table S3 for terminal residues) was then designed to increase its binding to divalent ions, which seemed underestimated in the first version. This was achieved through an increase of the polarity of the PO bonds, setting the charge on the non-bridging oxygen atoms to -0.7 and increasing the charge of the phosphorus accordingly. The corresponding DMP proxy ("DMPecc2") was designed as before.

Since this charge distribution was found to provide the best agreement with experimental data used for parametrization (see Results section), we then simply call this DNA force field "OL15-ECC". As backbone charges slightly differ between DNA OL15 and RNA OL3 force fields, we then derived an OL3-ECC RNA variant, imposing the same charges on the phosphorus atom and non-bridging phosphate oxygen atoms, while the charges of the O3' and O5' atoms were adjusted to yield an overall charge of -0.8 for each nucleotide (Supporting Information, Tables S4 and S5).

<table><caption>Table 1: Detailed charges (e) of the different DNA force fields and their respective DMP proxies.</caption><thead><tr><th></th><th>OL15</th><th>OL15-ECC-1</th><th>OL15-ECC-2</th><th></th><th>DMPamb</th><th>DMPecc1</th><th>DMPecc2</th></tr></thead><tbody><tr><td>P</td><td>1.1659</td><td>0.93272</td><td>1.09096</td><td>P</td><td>1.1659</td><td>0.93272</td><td>1.09096</td></tr><tr><td>O1P/O2P</td><td>-0.7761</td><td>-0.62088</td><td>-0.7</td><td>OP</td><td>-0.7761</td><td>-0.62088</td><td>-0.7</td></tr><tr><td>O3'</td><td>-0.5232</td><td>-0.45905</td><td>-0.45905</td><td>O(Me)</td><td>-0.5093</td><td>-0.44793</td><td>-0.44793</td></tr><tr><td>O5'</td><td>-0.4954</td><td>-0.43681</td><td>-0.43681</td><td></td><td></td><td></td><td></td></tr><tr><td>C3'</td><td>0.0713</td><td>0.0713</td><td>0.0713</td><td>C</td><td>0.0322</td><td>0.0322</td><td>0.0322</td></tr><tr><td>C5'</td><td>-0.0069</td><td>-0.0069</td><td>-0.0069</td><td></td><td></td><td></td><td></td></tr><tr><td>total</td><td>-1</td><td>-0.8</td><td>-0.8</td><td></td><td>-1</td><td>-0.8</td><td>-0.8</td></tr></tbody></table>

These variants of the DMP force field (with full and scaled charges) were used in combination with SPC/E water,⁶⁹ which is a good 3-point water model that has been shown to yield good conformational ensembles for DNA when used in combination with the Amber-OL15 force field.⁷⁰ Further, this choice allowed us to take advantage of previously derived ECC force fields for ions, compatible with SPC/E.⁵⁸

Full charge DMP or DNA force fields were combined with the Joung and Cheatham force
field for monovalent ions (K⁺, Na⁺, Cl⁻),⁷¹ while ECC force fields were used together with 0.8-scaled versions of these ions, derived herein, that include—following the usual ECC ions parametrization strategy which has been found to improve the description of electrolytes in previous works⁷²⁻⁷⁴—fine tuning of the van der Waals radius in addition to charge scaling to recover proper ion-water structural properties (Supporting Information Table S1). While these ECC0.8 force field variants for ions are not as optimized as state-of-the-art force fields such as the Madrid2019⁵⁶ ions (which use a 0.85 scaling and thus could not be employed here), we checked against recent neutron scattering data with isotopic substitution on 4 m KCl solutions⁷⁵ that this ECC0.8 force field variant for ions yields a good structure for bulk KCl electrolyte, very similar to that of the original full charge force field (Supporting Information Fig S1). We then more specifically assessed the cation-phosphate interaction, which is the core of our interest here.

## Osmotic coefficient calculations

We estimated with each force field variant the osmotic coefficient of KDMP and NaDMP solutions based on the procedure initially described by Luo and Roux,⁷⁶ and adapted for a spherical constraint by Piana et al.⁶⁷ The simulated system consisted of 32 DMP:K or DMP:Na ion pairs in a water box containing 23692 SPC/E water molecules (size of about ≈ 90 Å).

A 1 ns equilibration phase under NPT conditions (using the stochastic velocity rescaling thermostat⁷⁷ and Parinello-Rahman barostat⁷⁸) was followed by 500 ps equilibration under NVT conditions. A flat-bottom harmonic potential was then added to restrain the position of the ions relative to the center of the box, with a force constant of 10 kcal/mol/Å². This harmonic restraint was applied gradually over 200 ps to reach target concentrations of 0.3 M, 0.6 M for K⁺ (or 0.5 M for Na⁺) and 1 M inside the restraint volume, corresponding to restraint radii of 34.84 Å 27.65 Å (29.38 Å for Na⁺), and 23.32 Å respectively. 250 ns-long NVT simulations were then performed to evaluate the osmotic coefficient at each concen-
tration, using 10 replicas with different initial velocities to reach proper convergence and estimate error bars.

The osmotic pressure was calculated from the average force applied per unit area on the harmonic wall. The osmotic coefficient $\phi$ was then determined by the ratio $\phi = P/P_0$, where $P_0$ represents the reference pressure for a 1:1 salt solution at 1 mol/L. The error on this coefficient is estimated from the standard deviation calculated on the 10 replicas and varies between 0.006 and 0.01, depending on the simulations.

## Binding free energy calculations

The standard binding free energy between a Mg$^{2+}$ cation and the DMP anion, $\Delta G^{\circ}_{bind}(Mg:DMP)$, was computed for the various force fields using the well-established alchemical simulation strategy,⁷⁹⁻⁸³ within the double decoupling theoretical framework.⁸³⁻⁸⁷ Binding free energy calculations were performed independently in each binding mode (monodentate, bidentate, solvent-shared). They rely primarily on two alchemical transformations, that gradually decouple the cation from its environment, either in bulk ($\Delta G_{bulk}$) or within the ion pair ($\Delta G_{site}$), where flat-bottom harmonic restraints on the cation-phosphorus distance maintain the ion pair in the desired binding mode throughout the transformation. A correction term $\Delta G_{PBC}$ accounts for the free energy contribution due to net charge changes during the alchemical transformation,⁸⁸ and the free energy contribution due to applying restraints in the decoupled state, relative to the standard state concentration, $\Delta G^{V\circ \rightarrow rest}_{decoupled}$, was estimated numerically, as previously described.⁵⁸ Detailed discussion of this thermodynamic framework and free energy terms can be found in⁸⁹ and.⁹⁰

The total binding free energy $\Delta G^{\circ}_{bind}(Mg:DMP)$ was finally derived by combining the free energies for each ion-pair geometry as follows:

$$ \Delta G^{\circ}_{bind,tot} = -RT \ln \left( e^{-\Delta G^{\circ}_{bind,mono}/RT} + e^{-\Delta G^{\circ}_{bind,bi}/RT} + e^{-\Delta G^{\circ}_{bind,sship}/RT} \right) \quad (1) $$
Alchemical transformations were performed with the different force fields using Gro-
macs2023.2 with a simulation box that contained one DMP anion, one Mg²⁺ cation and
1723 water molecules. Electrostatic interactions between the cation and its surroundings
were first removed using 11 windows with λ<sub>elec</sub> values ranging from 1 to 0 in increments of
0.1, followed by the removal of van der Waals interactions across 4 windows with λ<sub>vdw</sub> values
of 0.7, 0.5, 0.1, and 0.

Starting from the coupled state, each window was equilibrated for 250 ps before proceed-
ing to the next, with an additional 5 ns of sampling per window. The Bennett Acceptance
Ratio (BAR) method,⁹¹ as implemented in Gromacs, was then used to compute the free
energy. The reported error on ΔG°bind is the sum of the errors for each alchemical transfor-
mation, estimated from the variance over 5 blocks. For a selected system, the alchemical
transformations were replicated three times. The standard deviation of the calculated free
energies was consistent with the error estimate reported by gmx bar. We also verified that
the computed free energy was independent of the direction of the alchemical transformation
(forward/reverse), which is a good indicator of convergence.

Molecular simulations of nucleic acids

Molecular dynamics simulations

To assess the performance of our AmberOL15-ECC and AmberOL21-ECC force fields to
capture the conformations of nucleic acids, simulations were carried on 3 different DNA se-
quences: the well characterized double-stranded Dickerson-Drew dodecamer,⁹² d(CGCGAATTCGCG)₂
(PDB ID 1NAJ⁹³); the antiparallel guanine quadruplex (G-DNA) formed by the sequence
d(G₄T₄G₄) (PDB ID 2GWQ, chain A); a left-handed Z-DNA hexamer duplex d(CGCGCG)₂
(PDB ID 1ICK⁹⁴; not including spermine and Mg²⁺ ions); a Jun-Fos oligomer duplex
d(GCATTCTGAGTCAG)⁹⁵ (starting structure prepared using the nucgen⁹⁶ module of AM-
BER); a short A-tract d(CGCA₃T₃GCG)₂ (PDB ID 1S2R⁹⁷); a sequence with a purine-
rich strand d(GAAGAGAAGC) (PDB ID 1AXP;⁹⁸ first NMR model) a G₄C₄ sequence
d(GGGGCCCC)$_2^{99}$ (starting structure prepared by the nucleic acid builder (NAB) module of AMBER;$^{100}$ and a DNA/RNA hybrid r(GAAGAGAAGC)/d(GCTTCTCTTC) (PDB ID 1RRD;$^{101}$ first NMR model). We also simulated, as a preliminary test, a double stranded A-RNA decamer r(GCACCGUUGG)$_2$ (1QC0', excised from the longer sequence PDB ID 1QC0$^{102}$). A list of all MD simulations is provided in Supporting Information Table S6.

Using Gromacs tools, all systems were explicitly solvated in a cubic box of SPC/E water molecules (size of about ≈ 70 Å for G-DNA, Dickerson-Drew dodecamer and A-RNA, and of about ≈ 60 Å for Z-DNA, ≈ 78 Å for Jun-Fos oligomer, ≈ 71 Å for A-tract, ≈ 61 Å for 1AXP, ≈ 54 Å for G$_4$C$_4$ sequence, ≈ 65 Å for DNA/RNA hybrid) and neutralized with potassium ions with addition of an extra 150 mM KCl background, except for Z-DNA, experimentally observed at high salt concentrations, which was simulated in a 2 M KCl solution. The B-DNA dodecamer simulations were repeated with a 150 mM NaCl background for comparison.

For each system and each tested force field, an initial minimization step was followed by a 20 ns equilibration in the NPT ensemble, using the Parinello-Rahman barostat$^{78}$ with a 5 ps coupling time and a stochastic velocity rescaling thermostat$^{77}$ at 300 K with a 1 ps coupling time. Three 1 µs-long replicas were then propagated in the NPT ensemble using different initial velocities. All simulations were carried out using Gromacs2019.6$^{103,104}$ with a 2 fs integration time step, using periodic boundary conditions in the isothermal-isobaric ensemble (NPT). We employed a Parinello-Rahman barostat$^{78}$ with a 5 ps coupling time and a stochastic velocity rescaling thermostat$^{77}$ at 300 K with a 1 ps coupling time. Long-range electrostatic interactions were treated with the particle-mesh Ewald (PME) method,$^{105}$ with a 12 Å cutoff. Bonds containing hydrogens were constrained using the LINCS algorithm,$^{106}$ and the SETTLE algorithm was used for water molecules.$^{107}$

Each system was simulated using at least two force fields: the original broadly used Amber-OL15 force field for DNA—i.e. amber ff99 with bsc0,$^{108}$ χ$_{OL4}$,$^{109}$ ϵζ$_{OL1}$ $^{110}$and β$_{OL1}$ $^{68}$ corrections—and the OL3 force field for RNA—Amber99ff-bsc0 with χ$_{OL3}$ $^{26}$ corrections—
together with the corresponding ECC variants proposed in this work. The original full
charge force field was used in combination with the Joung and Cheatham force field (SPCE-
compatible version) for potassium (or sodium) and chloride ions.⁷¹ In simulations with the
ECC-force field for nucleic acids, we also used ECC (0.8-scaling) force fields for the ions,
which were derived herein based on the full charge parametrization, targeting the K-OW
(or Na-OW) and Cl-OW correlation functions to tune the van der waals radius after charge
scaling, typically decreasing it by about 5% (see Supporting Information). Z-DNA and G-
DNA simulations were also repeated in triplicate with the more recent OL21 force field—that
improves the description of α/γ conformational substates in non canonical DNA forms—and
its OL21-ECC variant. Force field parameters are shared on a public Github folder¹¹¹ and
trajectories are shared on a public Zenodo folder.¹¹²

## Simulations with divalent ions

Simulations of the Dickerson Dodecamer⁹² d(CGCGAATTCGCG)₂ (PDB ID 1NAJ⁹³) were prepared as previously described, neutralizing its charge with either Ca²⁺ or Mg²⁺ ions and adding a background 150 mM concentration of KCl. The system was equilibrated for 1 µs, before a 1 µs production run, performed in 3 replicas with different initial velocities. The simulations were performed both with the original Amber-OL15 force field and full charge divalent ions with 12-6 Lennard-Jones parameters parametrized by Merz and coworkers,¹¹³ and with our Amber-OL15-ECC 0.8-scaled charge force field variant for nucleic acids associated with a consistent ECC0.8 parametrization for Ca²⁺ and Mg²⁺ ions.⁵⁷,⁵⁸

# Analyses of nucleic acid conformation and ion binding

## Structural analyses.

All molecular dynamics simulations were analyzed using CPPTRAJ¹¹⁴ to extract helical pa-
rameters and dihedral angles (see definitions on Fig. 1). Typical analysis scripts are shared
on a public Zenodo folder.¹¹²
<figure><img src="image_2.png" /><figcaption>Figure 1: Scheme of dihedral angles on the DNA backbone (in orange) and the glycosidic torsion angle (in purple) within a nucleotide unit</figcaption></figure>

**B-DNA** The BI/BII conformational states were determined using the same definition as Hartmann et al.,¹¹⁵ based on the difference between the ϵ and ζ backbone dihedral angles: BI state when ϵ − ζ < 0°, and BII when ϵ − ζ > 0°. Sugar puckers were assigned to the North (N) conformation if the pseudorotation angle P was between -90° and 90°; otherwise, they were classified as South (S).

**G-DNA** Guanines in *anti* and *syn* conformations are treated separately.⁶⁸ Guanines in the *anti* conformation are not discussed, as they adopt canonical angles. For *syn*-guanines, the analysis focuses on the *α*/*γ* states, which are usually non-canonical. To calculate the population of *α*/*γ* states, the *α* and *γ* dihedral angles are calculated for all *syn* guanines (G3, G9, G11, G15, G21, G23). The different backbone substates are classified using the following criteria on *α* and *γ* dihedrals: gauche+ (g+): 0-150°, trans (t): 150-250°, and gauche− (g−): 250-360°. Reference experimental backbone angles and substate populations of G-DNA were extracted from the X-ray structures with a resolution of 2 Å or better without ligands (1JPQ, 2GWQ, 2HBN, 4R45, 4R47), as done in previous work.⁶⁸
**Z-DNA** Dinucleotide conformations (GpC steps: ZI, ZII; CpG steps: native) were assigned based on the nomenclature by Černý et al.,¹¹⁶ using nucleotide conformational classes (NtC) from the DNATCO database (ZI = NtC ZZS1; ZII = NtC ZZS2; CpG native = NtC ZZ1S). Each dinucleotide step was assigned the class with smallest Euclidean distance in the nine-dimensional space of backbone and glycosidic dihedral angles (δ₁, χ₁, ϵ₁, ζ₁, α₂, β₂, γ₂, δ₂, χ₂) to the DNATCO-defined classes. The analysis focuses on ZI/ZII; a minor fraction (3-7 %) of other NtCs appeared in GpC steps but is not considered further.

**Ion density maps** Ion density maps were computed using a custom Python script, leveraging the MDAnalysis package.¹¹⁷ Each system was first aligned to a reference structure ensuring structural consistency across replicas. The K⁺ ion density was calculated with a spatial resolution of 1 Å using MDAnalysis.density. The density values (ions/Å³) were then averaged over three independent replicas and converted to molar concentrations (mol/L). The final density maps, exported as .dx files are shared on a public Zenodo folder¹¹² and can be visualized using VMD.¹¹⁸

# Results

## Assessment of ion-pairing to the phosphate backbone on a model system: DMP

The negatively charged phosphate backbone is the primary interaction site for cations with nucleic acids, even if high cation densities can be found at specific locations, for example in DNA grooves.⁶,¹⁰,¹¹⁹,¹²⁰ To assess the performance of various force fields to capture the strength of ion-nucleic acid interactions, we thus chose to focus on the phosphate backbone and evaluate ion binding with a model phosphate compound, the dimethylphosphate (DMP, Figure 2), often considered as the simplest proxy for nucleic acid backbone and sometimes used as a building block for nucleic acid force fields.⁴⁶
<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic representation of ions bound to DNA, with a zoom on the phosphate backbone and a 3D chemical representation of the employed proxy, dimethylphosphate (DMP). Prepared with VMD.¹¹⁸</figcaption></figure>

In this work, we assessed DMP ion pairing behavior with both $K^+$—as the most abundant monovalent cation inside a cell— and with the divalent $Mg^{2+}$ ion, known to play a key role in numerous nucleic acid related processes. More specifically and following earlier work on phosphate-ion pairing,⁶⁷,¹²¹,¹²² we used as a reference the experimental osmotic coefficients for KDMP solutions at concentrations between 0.3 and 1 M.¹²³ Osmotic coefficients quantify the deviation of a solution from an ideal behavior; for electrolytes, they reflect the tendency of ions to associate in solution¹²⁴,¹²⁵ and are thus sensitive probes of ion pairing. They have been often used to assess the quality of force fields to describe simple electrolytes, or even as targets for force field development.⁵¹,⁶⁵,⁶⁷,¹²⁵⁻¹²⁷ Since charge scaling is particularly important for binding to divalent ions, we also target the experimental estimates for the binding free energy between dihydrogen phosphate (DHP) and $Mg^{2+}$ ions, $\Delta G°_{bind}(Mg^{2+}) = -7.1 \text{ kJ mol}^{-1}$.¹²⁸ Unfortunately, there is no experimental estimate of the binding free energy with DMP, which we assumed similar to that with DHP. We note however, that a slightly weaker estimate was obtained by theoretical extrapolation of various experimental values.¹²⁹ The limits and uncertainties associated with this experimental comparison are discussed in details in the Discussion. Overall, the experimental data pictures a weak binding of
Mg²⁺ ions to isolated phosphate groups, which—combined with the relatively low solubility of magnesium-phosphate salts—probably explains why experimental estimates of such an important ion pairing process are so scarce.

As we used DMP as a proxy for the nucleic acid backbone phosphate group, we derived a DMP force field ("DMPamb") directly from the existing AmberOL15 and AmberOL21 force fields for nucleic acids—that differ only in refined parameters for the backbone torsional angles α and γ but use identical partial charges (see Methods section). With this model, we computed both the osmotic coefficient for KDMP and NaDMP solutions, experimentally available up to 1 mol/L, (Figure 3 and Supporting Information Table S2) and the standard binding free energy of a DMP:Mg²⁺ or DMP:Ca²⁺ ion pair (Table 2), using for the ions standard full charge force fields often used in combination with Amber force fields. The osmotic coefficients of our model KDMP and NaDMP solutions are not only too low with the standard force field—suggesting some overbinding of potassium ions to phosphate groups—but also diminishes (more markedly for Na⁺ than K⁺) with increasing concentrations, which is at odds with the experimental trend. As expected from earlier studies,⁵⁸ the binding of magnesium ions is dramatically overestimated, with a binding free energy more than 20 kJ/mol more negative than the experimental reference. Together, these results highlight the poor description of ion pairing to nucleic acids by state-of-the-art full charge force fields, and stress the need for improvement.

Following the Electronic Continuum Correction (ECC) strategy⁵³,⁵⁴ and earlier works on its application to various electrolytes,⁵⁷,⁵⁸,⁶⁰ we chose to design a scaled-charge version of DMP, using a 0.8 scaling factor, building it so that the charges of the phosphate group would be directly transferable to the DNA force field. The ion force fields are combined with the SPC/E water model,⁶⁹ which we used in previous work on related electrolytes and which has been shown to be compatible with the Amber force fields for nucleic acids.⁷⁰ A "naive" rescaling of all the 5 partial charges of the PO₄ group by this 0.8 factor (DMPecc1) brings the KDMP osmotic coefficient (0.93) significantly closer to the experimental reference
(0.98) at 0.3 M, even if it does not significantly increases with concentration, as it should, but rather stays constant (see Figure 3). Charge scaling has a very striking impact on the binding free energy with Mg²⁺, reducing it by almost an order of magnitude, to -2.8 kJ/mol (see Table 2). While the full charge model was dramatically overbinding divalent ions, this first scaled charge model thus seems to slightly underestimate it.

The ECC approach has a clear prescription for the total charge of the ions, but there is no consensus yet as how to derive the best charge distributions for molecular ions. We thus tested, keeping an overall charge of -0.8 different charge distributions on the phosphate group, varying in particular the polarity of the PO bonds—the charge of the non-bridging oxygen atoms being intuitively a key parameter for interaction with ions. The osmotic coefficient of this DMPecc2 model with K⁺ is similar to that of DMPecc1 and significantly improved with respect to the initial full charge model, even if there is still room for improvement, with our osmotic coefficient still smaller than the experimental estimates. The improvement is similar with Na⁺, where ECC also avoids the marked unphysical decrease with increasing concentration, which was observed with the full charge model. While the osmotic coefficient and the stability of the solvent-shared ion pairs (SShIP) are almost independent of (small) changes in the charge distribution, the binding free energy to Mg²⁺ ions is greatly modulated by the choice of partial charges. ΔG°bind(Mg²⁺) computed with our best model, DMPecc2 (see Methods section, Table 1) falls within the error bars of the experimental estimate for DHP (-8.0 ± 2.0 vs -7.1 kJ/mol, see Table 2). The same DMPecc2 model for DMP also yields a binding free energy with Ca²⁺ in good agreement with experimental estimates for DHP (-3.4 ± 2.0 vs -5.9 kJ/mol, see Table 2). We note that a theoretical estimate¹²⁹ of these binding constants for DMP (from data on various related compounds) yields weaker binding constants with both cations, which would then be closer to that of our DMPecc1 model for Mg²⁺ and in perfect agreement with our DMPecc2 model for Ca²⁺. We selected DMPecc2 as our best model for the rest of the study.
<figure><img src="image_4.png" /><figcaption>Figure 3: Osmotic coefficients $\phi$ for KDMP (circles and solid lines) and NaDMP (diamonds and dashes) solutions at increasing concentrations computed with different force fields— full charge (blue), DMPecc1 (orange) and DMPecc2 (red)— compared with experiment by Tamaki et al.¹²³ (black). The estimated error of 0.009 with all force fields. Lines are guides for the eye.</figcaption></figure>

<table><caption>Table 2: Standard binding free energy $\Delta G^{\circ}_{\text{bind}}$ (kJ/mol) for a DMP:Mg²⁺ and DMP:Ca²⁺ ion pair computed with different force fields, with full or scaled charges. $\Delta G^{\circ}_{\text{bind}}$ is calculated independently for each binding mode (contact monodentate or bidentate and solvent-shared). The estimated error bar for all binding free energies is around 2 kJ/mol. These calculated values are compared with both the experimental references pertains to cation binding to a closely related compound, dihydrogen phosphate (DHP), and theoretical extrapolations of experimental data on neighbor compounds by Sigel and Sigel¹²⁹ (the reported binding constant at 0.1 M ionic strength was extrapolated to $I = 0$ using a Davies-like equation).</caption><thead><tr><th></th><th></th><th>bidentate</th><th>monodentate</th><th>SSHIP</th><th>Total</th></tr></thead><tbody><tr><td rowspan="4">Mg<sup>2+</sup></td><td>DMPamb</td><td>-22.8</td><td>-29.6</td><td>-5.9</td><td>-29.8</td></tr><tr><td>DMPecc1</td><td>7.72</td><td>-0.3</td><td>-1.6</td><td>-2.8</td></tr><tr><td>DMPecc2</td><td>-2.5</td><td>-7.4</td><td>-2.2</td><td>-8.0</td></tr><tr><td><strong>Expt. <sup>128</sup> (DHP)</strong></td><td></td><td></td><td></td><td><strong>-7.1</strong></td></tr><tr><td></td><td><strong>Theoretical extrapol.</strong><sup>129</sup></td><td></td><td></td><td></td><td><strong>-3.2</strong></td></tr><tr><td rowspan="3">Ca<sup>2+</sup></td><td>DMPamb</td><td>-9.3</td><td>-8.4</td><td>-4.5</td><td>-10.8</td></tr><tr><td>DMPecc2</td><td>+1.3</td><td>-0.5</td><td>-1.9</td><td>-3.4</td></tr><tr><td><strong>Expt. <sup>130</sup> (DHP)</strong></td><td></td><td></td><td></td><td><strong>-5.9</strong></td></tr><tr><td></td><td><strong>Theoretical extrapol.</strong><sup>129</sup></td><td></td><td></td><td></td><td><strong>-2.9</strong></td></tr></tbody></table>
ECC force field variant for DNA and RNA

We directly transferred the charge distribution that yielded the best ion pairing properties for DMP (DMPecc2) on the nucleic acid phosphate backbone to build a first scaled-charge ECC variant of the Amber-OL15 force field, Amber-OL15-ECC-2 (see Table 1 Methods), which we will from now on simply call Amber-OL15-ECC. This force field differs from its full charge parent only in the partial charges assigned to the 5 atoms (PO₄) of each phosphate group, each nucleotide now carrying a −0.8 charge. This approach limits the changes strictly to the phosphate groups, in the hope to improve ion pairing properties without disturbing all the other qualities of the parent force field, optimized to properly capture the conformational equilibrium of various nucleic acid sequences.

The first step to validate the performance of this force field variant is to check the impact
of the charge modification on DNA conformational properties, that are well optimized in
the Amber-OL15 force field. To this aim, we selected 3 model DNA sequences with distinct
helical conformations (Figure 4): the famous B-DNA Dickerson dodecamer, a left-handed
Z-DNA hexamer duplex d(CGCGCG)₂ and the antiparallel guanine quadruplex (G-DNA)
formed by the sequence d(G₄T₄G₄). These structures are well characterized experimentally
and have been previously used in the literature to assess the performance of nucleic acid
force fields.⁶⁸,¹¹⁰ Designing force fields for RNA is even more challenging than DNA due to
the broad range of possible conformations. Preliminary tests of the applicability of our ECC
corrected force field to RNA have been performed using a model double-stranded sequence
in the typical A-form.
<figure><img src="image_5.png" /><figcaption>Figure 4: Representations of selected simulated nucleic acid sequences: B-DNA Dickerson-Drew dodecamer, Z-DNA and a DNA antiparallel G-quadruplex (G-DNA) and a ds-RNA system. Prepared with Chimera.¹³¹</figcaption></figure>

## Dickerson-Drew DNA Dodecamer

The Dickerson-Drew dodecamer $d(CGCGAATTCGCG)_2$ is probably the best characterized B-DNA duplex, with multiple NMR and high-resolution X-ray structures available in the literature.⁷,⁹²,¹³²⁻¹³⁴ It often serves as a benchmark for force field development and has been extensively studied with multi-μs-long molecular dynamics simulations.²⁷,³³,¹¹⁰ In particular, modern DNA force fields have been refined to properly capture B-DNA conformational properties, including the helical twist and the equilibrium between B1 and B11 backbone
substates (see definition in Methods section)—with early force fields underestimating the BII population.²⁷,¹¹⁰

We performed µs-long simulations of the dodecamer in both 150 mM KCl and 150 mM NaCl with the original OL15 force field and our OL15-ECC variant, and systematically compared the distributions of backbone dihedral angles—that define several conformational substates—and helical parameters (Table 3, Figure 5 and Supporting Information Figure S5 for KCl and Table S7 and Figs S7 to S9 for NaCl). The impact of the ECC variant being very similar with both electrolytes, we focus in the main text on KCl, while the analysis of simulations in NaCl is provided in the Supporting Information. All the dihedral and helical parameters are strikingly similar between the two force field variants, which shows that our charge scaling modifications on the phosphate backbone do not deteriorate the capability of the original force field to describe DNA substates. The only noticeable change in backbone dihedral angle distributions (Fig. 5) pertains to the ϵ angle, and reflects a small change in the BI/BII conformational equilibrium, even if the estimated BII population with OL15-ECC remains similar to that predicted by OL15. The overall BII population is within the experimental range (see Table 3) and the magnitude of the change between OL15 and OL15-ECC is much smaller than the variation between two experimental determinations. A closer examination of BII populations per each base step and comparison with experimental estimates (Supporting Information Figure S6) however suggests that OL15 and OL15-ECC both underestimate the BII population in the central part of the sequence.

As expected, our ECC force field variant leads to a decrease in the interaction between monovalent ions and the DNA phosphate backbone—in line with the increase in osmotic coefficient of the model KDMP system. The number of solvent-shared ion pairs—mediated by a single water molecule— is significantly reduced, from on average 5.4 to 3.9 K⁺ ions in solvent-shared ion pairs with phosphate groups (Table 4). In contrast, the number of contact ion pairs with phosphates does not change, while the total number of contacts slightly increases, thus reflecting a very limited increase in ion pairing to nucleobases upon
<table><caption>Table 3: Average structural parameters and standard deviation over the concatenated trajectory of the 3 1µs-replicas of Dickerson-Drew DNA dodecamer simulations, with the OL15 and OL15-ECC force fields. The two terminal base pairs at each end were excluded from the analysis.</caption><thead><tr><th></th><th>NMR<sup>a</sup></th><th>X-ray<sup>b</sup></th><th>OL15</th><th>OL15 ECC</th></tr></thead><tbody><tr><td>minor groove/Å</td><td>10.5</td><td>9.9 ± 1.1</td><td>11.1 ± 1.7</td><td>10.7 ± 1.6</td></tr><tr><td>major groove/Å</td><td>17.6</td><td>17.6 ± 1.0</td><td>17.7 ± 1.6</td><td>17.9 ± 1.6</td></tr><tr><td>shift/Å</td><td>0.0</td><td>0.0</td><td>0.0 ± 0.8</td><td>0.0 ± 0.8</td></tr><tr><td>slide/Å</td><td>-0.2</td><td>0.0 ± 0.5</td><td>0.0 ± 0.7</td><td>0.1 ± 0.7</td></tr><tr><td>rise/Å</td><td>3.2</td><td>3.3 ± 0.2</td><td>3.3 ± 0.3</td><td>3.3 ± 0.3</td></tr><tr><td>tilt/deg</td><td>0.0</td><td>0.2 ± 1.7</td><td>0.0 ± 5.0</td><td>0.0 ± 5.0</td></tr><tr><td>roll/deg</td><td>3.0</td><td>2.0 ± 4.0</td><td>2.2 ± 5.8</td><td>1.4 ± 6.0</td></tr><tr><td>twist/deg</td><td>35.7</td><td>33.6 ± 5.9</td><td>35.3 ± 6.6</td><td>36.0 ± 6.1</td></tr><tr><td>shear/Å</td><td>0.0</td><td>0.0</td><td>0.0 ± 0.3</td><td>0.0 ± 0.3</td></tr><tr><td>buckle/deg</td><td>0.0</td><td>0.9 ± 6.8</td><td>0.1 ± 10.5</td><td>0.0 ± 13.4</td></tr><tr><td>stretch/Å</td><td>-0.3</td><td>-0.1 ± 0.1</td><td>0.0 ± 0.1</td><td>0.0 ± 0.2</td></tr><tr><td>propeller/deg</td><td>-17.5</td><td>-11.8 ± 3.7</td><td>-11.7 ± 9.8</td><td>-11.9 ± 9.9</td></tr><tr><td>stagger/Å</td><td>-0.1</td><td>0.1 ± 0.1</td><td>0.0 ± 0.4</td><td>0.0 ± 0.4</td></tr><tr><td>opening/deg</td><td>-1.1</td><td>2.1 ± 2.4</td><td>-0.3 ± 4.2</td><td>-0.3 ± 8.8</td></tr><tr><td>X-displacement/Å</td><td>-0.8</td><td>-0.5 ± 0.6</td><td>-0.3 ± 1.5</td><td>-0.2 ± 1.4</td></tr><tr><td>Y-displacement/Å</td><td>0.0</td><td>0.1 ± 1.0</td><td>0.0 ± 1.2</td><td>0.0 ± 1.2</td></tr><tr><td>inclination/deg</td><td>5.0</td><td>4.3 ± 8.4</td><td>3.9 ± 10.2</td><td>2.6 ± 10.1</td></tr><tr><td>tip/deg</td><td>0.0</td><td>0.5 ± 3.4</td><td>0.0 ± 8.2</td><td>0.0 ± 8.3</td></tr><tr><td>helical twist/deg</td><td>36.1</td><td>34.0 ± 5.4</td><td>36.2 ± 6.3</td><td>36.9 ± 5.8</td></tr><tr><td>% BII</td><td>22.3,<sup>133</sup><br>31.1<sup>135</sup></td><td></td><td>24.0 ±0.5<sup>d</sup></td><td>23.8 ±<br>0.4<sup>d</sup></td></tr><tr><td>RMSD (all)<sup>c</sup></td><td></td><td></td><td>1.15 ±<br>0.20</td><td>1.16 ±<br>0.22</td></tr><tr><td>RMSD<br>(backbone)<sup>c</sup></td><td></td><td></td><td>1.22 ±<br>0.24</td><td>1.22 ±<br>0.26</td></tr></tbody></table>

<sup>a</sup> Average values from the 1NAJ structure<sup>93</sup> <sup>b</sup> Average values from well resolved crystal structures (resolution better than 1.5 Å, PDB IDs 1FQ2, 1JGR, 3U2N, 4C64, 355D) as in previous work.<sup>68</sup> <sup>c</sup> RMSD (mass-weighted) for all atoms and backbone calculated against the 1NAJ structure. <sup>d</sup> The error is estimated as the standard deviation between the averages calculated independently on each replica.

force field modification (Table 4). In addition, ion density maps computed with both force fields around the dodecamer (Supporting Information Figure S2) show that the ions occupy very similar locations with similar density in the grooves, so that the use of scaled charges
<figure><img src="image_6.png" /><figcaption>Figure 5: Distribution of backbone dihedral angles averaged over 3 1-μs replicas of simulations performed with the OL15 (blue) or OL15-ECC (red) force field.</figcaption></figure>

for the cations has little effect on their binding to nucleobases.

Antiparallel guanine DNA quadruplex

DNA G-quadruplexes (G4s) are non canonical structures formed by guanine-rich single-
stranded DNA sequences. Hoogsteen hydrogen bonding between four guanines yields planar
G-quartets, that further stack into G4s, with additional stabilization by monovalent ions,
primarily K+. G4s are implicated in key biological functions, including regulation of gene
<table><caption>Table 4: Average number (standard deviation between the 3 replicas in parenthesis) of K⁺ ions in various interaction modes with ds-DNA, ds-RNA and Z-DNA: direct contact (< 3.2 Å to any nucleic acid heavy atom), direct contact with a phosphate group and solvent-shared ion pair with a phosphate group (4 < d(P-K) < 6 Å)</caption><thead><tr><th rowspan="2"></th><th colspan="2">total direct contacts</th><th colspan="2">phosphate direct contacts</th><th colspan="2">SShIP w. phosphates</th></tr><tr><th>OL15</th><th>OL15-ECC</th><th>OL15</th><th>OL15-ECC</th><th>OL15</th><th>OL15-ECC</th></tr></thead><tbody><tr><td>B-DNA dodecamer</td><td>5.8 (0.0)</td><td>6.6 (0.1)</td><td>2.5 (0.0)</td><td>2.3 (0.0)</td><td>5.4 (0.0)</td><td>3.9 (0.0)</td></tr><tr><td>Z-DNA</td><td>10.2 (0.0)</td><td>10.3 (0.1)</td><td>7.2 (0.0)</td><td>6.7 (0.0)</td><td>8.3 (0.0)</td><td>7.2 (0.0)</td></tr><tr><td>G-DNA</td><td>9.3 (0.1)</td><td>8.3 (0.0)</td><td>5.1 (0.2)</td><td>2.6 (0.0)</td><td>5.9 (0.1)</td><td>3.1 (0.0)</td></tr><tr><td>ds-RNA</td><td>8.4 (0.0)</td><td>7.2 (0.0)</td><td>3.4 (0.0)</td><td>2.2 (0.0)</td><td>6.1 (0.0)</td><td>3.9 (0.0)</td></tr></tbody></table>

expression, maintenance of telomere integrity, and genome stability. They are suggested to play crucial roles in cancer, aging, and are potential therapeutic targets.¹³⁶,¹³⁷

G4s often present unusual backbone topologies and are thus challenging test cases for force field refinement,¹⁰⁹,¹³⁸,¹³⁹ with, for instance, a complex equilibrium of α/γ conformational substates. The conformation of the single-stranded loops of the G4s is also very sensitive to force field accuracy. In addition, non polarizable force fields all fail to properly stabilize the closely-spaced monovalent ions inside G-quadruplexes,⁶⁸,¹⁴⁰ so that the use of explicitly polarizable force fields is thought required to recover proper ion binding and G4 structures.⁴⁹,¹⁴⁰

Here, we performed long simulations (3 replicas, 1 µs) of the d(G₄T₄G₄) (PDB ID 2GWQ) antiparallel DNA G-quadruplex (Fig. 6a), that has been previously simulated with different force fields.⁶⁸,¹⁰⁹ It is composed of two G₄T₄G₄ strands that form 4 stacked G-quartets with two four-thymidine diagonal loops. A total of 5 K⁺ ions are resolved in the structure: 3 in the G-stem, located between tetrads, and one at each stem-loop junction. This arrangement of five cations in the quadruplex core region was previously found unstable in all simulations performed with non polarizable force fields, presumably because of the overestimated repulsion between cations.⁶⁸,¹⁰⁹,¹³⁸

Our OL15-ECC force field variant remarkably stabilizes the G4 structure, in particular in the loops, with much smaller RMSD deviations with respect to the PDB structure (Fig. 6b). With Amber-OL15, significant fluctuations are observed in the loops in all replicas, with an
a

<figure><img src="image_7.png" /></figure>

b

<figure><img src="image_8.png" /><figcaption>Figure 6: a) Schematic representation of the modeled antiparallel DNA quadruplex d(G<sub>4</sub>T<sub>4</sub>G<sub>4</sub>). The 5 monovalent ions observed in the crystal structures are shown in blue and numbered. Guanines in syn glycosidic conformation are shown in pink, while those in anti are shown in gray. b) RMSD evolution for one replica using OL15 (blue) and OL15-ECC (red) force fields, for the entire structure (ALL), guanines and the loop regions (LOOPS). RMSDs were computed by aligning the snapshots to the heavy atoms of PDB 2GWQ, chain A.</figcaption></figure>
increase of the RMSD that indicates a departure from the initial crystallographic structure.
In contrast, with the OL15-ECC variant, the loops maintain their initial structure throughout
the simulations.

This increased stability of the G4 structure is probably due to the much improved ion
retention observed with OL15-ECC: while at least 2 ions quickly escape the guanine core in
simulations with the parent full charge OL15 force field (as previously reported), the location
of the 5 potassium ions is stably conserved during most of our simulations with OL15-ECC (
Fig. 7) with only transient escapes from their reference positions. This strong improvement
in ion retention probably originates from the reduced charge on K⁺ ions within the ECC
framework. These results thus demonstrate that G4 structures can be stabilized without
explicitly polarizable force field, but do require electronic polarization to be properly taken
into account—even if only implicitly.

<figure><img src="image_9.png" /><figcaption>Figure 7: Distribution of conformations with 3, 4, or 5 potassium ions in the G-DNA structure simulated for the three replicas of each force field: the standard force field (OL15, in blue) and its ECC variant (OL15-ECC, in red).</figcaption></figure>

The unusual backbone conformational substates adopted by G-quadruplexes constitute
an additional stringent test for force field accuracy. *Syn* guanines can be divided in two
groups in this quadruplex: those that belong to the outer tetrads (G9 and G21), that adopt
an *α*/*γ* =g+/t backbone conformation in the X-ray structure; and those of the inner tetrads
(G3, G11, G15 and G23, see Fig. 6a), found in two distinct conformations: g+/g− (most populated) and g−/g−. The population of these conformational substates converges slowly in simulations, so they were quantified on very long (10 µs) simulations. Table 5 shows that the quality of description of the G4 backbone conformational substates slightly decreases when using the OL15-ECC variant, compared to its parent Amber-OL15 force field. For outer guanines, the population of the probably spurious g−/g+ substate (not observed in X-ray structures) increases; for inner guanines, the population of g+/g− increases—getting closer to the experimental value— but the population of the substate g−/g−, already too small in OL15, gets even smaller.

We have noticed that these conformational substates have been substantially improved by a reparametrization of the α/γ backbone dihedrals,¹⁴¹ which is implemented in the so-called Amber-OL21 force field. We thus designed an OL21-ECC variant of this force field, using the same charge modifications as previously described for OL15-ECC. An improved description of backbone substate is observed between OL21-ECC and OL15-ECC, similar to that between OL21 and OL15: the population of the g+/t dominant substate is increased for outer guanines; for inner guanines, the OL21 reparametrization increases the population of the minority g−/g− substate to bring it closer to experimental estimates, even if the gain seems smaller for our ECC variant than for the original force field, with a still too populated g−/g+ substate. Interestingly, the OL21 force field significantly increases the loop stability (Supporting Information Figure S12), even if it does not improve potassium ion retention within the quadruplex (Supporting Information Figure S10).

Overall, the OL15/OL21-ECC variants of the OL15/OL21 force fields provide a very good description of this DNA G-quadruplex, with a most significant improvement in K⁺ retention within the structure and loop stabilization achieved by the ECC variant, with only a small degradation in backbone conformational substates.
<table><caption>Table 5: Population (%) of various α/γ substates of syn-χ guanines in the stem of the antiparallel d(G₄T₄G₄) G-DNA quadruplex, computed on a 10 µs-long simulation with each force field.ᵃ</caption><thead><tr><th rowspan="2"></th><th rowspan="2">syn guanines</th><th colspan="4">α/γ</th></tr><tr><th>g-/g+</th><th>g+/g-</th><th>g+/t</th><th>g-/g-</th></tr></thead><tbody><tr><td rowspan="2">OL15</td><td>G9, G21</td><td>12</td><td>0</td><td>88</td><td>0</td></tr><tr><td>G3, G11, G15, G23</td><td>9</td><td>63</td><td>0</td><td>27</td></tr><tr><td rowspan="2">OL15-ECC</td><td>G9, G21</td><td>26</td><td>0</td><td>74</td><td>0</td></tr><tr><td>G3, G11, G15, G23</td><td>18</td><td>73</td><td>0</td><td>9</td></tr><tr><td rowspan="2">OL21</td><td>G9, G21</td><td>3</td><td>0</td><td>97</td><td>0</td></tr><tr><td>G3, G11, G15, G23</td><td>7</td><td>54</td><td>0</td><td>39</td></tr><tr><td rowspan="2">OL21-ECC</td><td>G9, G21</td><td>6</td><td>0</td><td>94</td><td>0</td></tr><tr><td>G3, G11, G15, G23</td><td>28</td><td>57</td><td>1</td><td>14</td></tr><tr><td rowspan="2">X-ray</td><td>G9, G21</td><td>0</td><td>0</td><td>100</td><td>0</td></tr><tr><td>G3, G11, G15, G23</td><td>0</td><td>78</td><td>0</td><td>22</td></tr></tbody></table>

ᵃ Gauche+, gauche–, and trans states of the α and γ dihedral angles are denoted g+, g–, and t, respectively. Green backgrounds indicate native states, white backgrounds indicate non-native states, and light green indicates states that are less populated in the X-ray structures. Populations do not sum to 100% because some rare substates are not shown in the table.

## Z-DNA

Z-DNA is a challenging molecule for force fields due to its backbone containing multiple non-canonical conformational substates. Since these substates involve torsions around the phosphate group, on which ECC was applied, we will focus on them in the following. In the simulated Z-DNA hexamer d(CGCGCG)₂ there are only two types of base pair steps: CpG and GpC. Experiments show that two conformational substates coexist in the GpC steps, both non-canonical, generally denoted as ZI and ZII.¹⁴² In the CpG step, experiments show only one native conformation, also non-canonical, characterized by the α/γ = g+/t conformation. Modelling these Z-DNA substates in the Amber ff family has been challenging for a long time and several consecutive ff modifications were necessary to establish a picture consistent with experimental data,⁶⁸,¹¹⁰ including the latest OL21 refinement.¹⁴¹ Thus, in the case of Z-DNA, in addition to OL15 simulations, we also performed simulations using the newer OL21 force field.

The Z-DNA molecule was generally stable in all our simulations, both with and without
ECC charges. RMSD of the stem is comparable in all simulations (1.00, 1.01, 1.02, and 1.01 Å for OL15, OL21, OL15-ECC, and OL21-ECC, respectively) (Supporting Information Figure S13). The more pronounced effect of ECC was on the populations of the non-canonical backbone conformations. In combination with both OL15 and OL21, ECC increased the population of the ZII substate. Since ZII is significantly populated in the X-ray structures (though ZI is dominant), and has been suggested to prevail in hydrated films,¹⁴³ a non-negligible ZII population is expected in solution. Unfortunately, accurate experimental data is unavailable—for more discussion, see Ref. ⁶⁸ However, populations of 11 % and 12 % found in OL15 and OL21 simulations (Table 6), respectively, seem too low. The ECC-variant increases these populations to 17 % (OL15-ECC) and 30 % (OL21-ECC), which is thus likely a beneficial effect of the ECC charge scaling.

<table><caption>Table 6: Populations (%) of the native backbone substates in the GpC and CpG steps of Z-DNA, 3' terminal base pairs excluded</caption><thead><tr><th>Step</th><th>Conformation</th><th>OL15</th><th>OL15-ECC</th><th>OL21</th><th>OL21-ECC</th></tr></thead><tbody><tr><td rowspan="3">GpC</td><td>ZI</td><td>85</td><td>77</td><td>84</td><td>63</td></tr><tr><td>ZII</td><td>11</td><td>17</td><td>12</td><td>30</td></tr><tr><td>ZI + ZII</td><td>97</td><td>94</td><td>96</td><td>93</td></tr><tr><td>CpG</td><td>Native</td><td>64</td><td>48</td><td>80</td><td>77</td></tr></tbody></table>

In the case of CpG steps, Z-DNA exhibits an unusual (native) $\alpha/\gamma = g+/t$ conformation. However, simulations with force fields prior to OL21 significantly populated a spurious conformation with $\alpha/\gamma = t/g+$, which is not experimentally supported.¹⁴¹ As shown in Table 6, ECC somewhat decreases the percentage of the native $\alpha/\gamma = g+/t$ conformation, from 64 % with OL15 to 48 % with OL15-ECC, primarily due to the spurious $\alpha/\gamma = t/g+$ substate. With the OL21 refinement, which stabilizes the native conformation, its decrease with ECC is less pronounced, from 80 to 77 %. Thus, in OL21 simulations, the native conformation dominates, and ECC only slightly reduces its population. Importantly, this indicates that the use of ECC does not change the trends observed with the underlying force fields. For example, if OL21 improves the description of the CpG step compared to OL15, a similar improvement is seen in OL21-ECC compared to OL15-ECC. This further supports our belief
that ECC can be combined with other force field variants without introducing undesirable behaviors.

In our simulations, Z-DNA exhibits some tendency for terminal bases to fray (Supporting Information Figure S13). Fraying, or the transient opening of the terminal bases, is a natural phenomenon in B-DNA.<sup>144</sup> However, there are no experimental data available on frequency of fraying in Z-DNA, so it is unclear whether fraying is described better with or without ECC charges.

In contrast with the effect of ECC on K<sup>+</sup> ion pairing in the DNA dodecamer—where it mainly reduced solvent-shared ion pairs with phosphates and had little effect on the formation of contact ion pairs—for Z-DNA, we observe a small decrease also in the number of direct contacts, from 7.2 on average with OL15 to 6.7 with OL15-ECC (see Table 4 and density maps Supporting Information Figure S3).

Overall, ECC did not introduce any new (non-native) substates in Z-DNA backbone, as seen from the histograms of the backbone angles (Supporting Information Figure S14), and it only slightly altered the existing equilibria, mainly increasing the population of the ZII backbone substate. Therefore, ECC can be used in combination with both OL15 and OL21 force fields, although the latter is recommended.

## Other B-DNA molecules and RNA/DNA hybrid

In further testing, we included a short B-DNA A-tract, A<sub>3</sub>T<sub>3</sub>, with a known X-ray structure<sup>97</sup> (PDB ID 1S2R), a Jun-Fos sequence in which its BI/BII equilibrium was characterized by NMR spectroscopy,<sup>95</sup> a GC-rich sequence, G<sub>4</sub>C<sub>4</sub>, for which NMR data on S/N sugar pucker equilibrium is available,<sup>98</sup> a DNA/RNA hybrid<sup>101</sup> (PDB ID 1RRD), in which higher A-form content is induced in the DNA strand by the A-form RNA complement, and a B-DNA duplex sharing the same purine-rich DNA strand as the hybrid<sup>98</sup> (PDB ID 1AXP).

Overall, the ECC correction had only a minor impact on the behavior of these systems, with simulations using OL15 and OL15-ECC producing very similar results. RMSDs are
compared in Fig. S15, where the ECC variants exhibited slightly lower RMSDs than OL15. Terminal base pairs remained mostly stable, with only occasional fraying events observed in both OL15 and OL15-ECC simulations (Fig. S16). As with the Dickerson-Drew dodecamer, the BII population and its sequence dependence in the Jun-Fos sequence were largely unaffected by the ECC charges (Fig. S17).

The structural parameters of the short $A_3T_3$ A-tract are compared with the X-ray structure in Table S8. Similar to the Dickerson-Drew dodecamer, simulations with OL15-ECC yielded helical parameters very similar to those of OL15, with a slightly narrower minor groove, marginally wider major groove and a slightly larger helical twist. Other helical parameters and backbone dihedral angle values were almost unchanged.

In B-DNA duplexes, the S/N sugar puckering equilibrium, which—together with the *anti* conformation of the glycosidic torsion $\chi$ defines A- and B-form DNA—is usually shifted toward the B-form, with only minor A-DNA content. It has been shown that current AMBER force fields, including OL15, substantially underestimate A-form population,¹⁴⁵ an issue corrected only in the most recent OL24¹⁴⁶ force field. As shown in Fig. S18, for two test sequences with available NMR data on the N sugar pucker population, the underestimation of the N pucker was similar for OL15 and OL15-ECC, possibly somewhat more pronounced with ECC. Sequence dependence of the pucker equilibrium remained almost unchanged.

In conclusion, the ECC correction has only a modest effect on the geometry, dynamics and conformational equilibria of B-DNA molecules. While this may appear surprising given the substantial reduction in the net charge of the DNA duplexes, it is important to consider that the screening effects of the dissolved ions are also reduced. The relatively small structural and dynamical changes observed are therefore likely the result of a compensatory balance between reduced molecular charge and reduced electrostatic screening.
Preliminary RNA simulations

While DNA adopts predominantly a double-stranded structure, RNA is more typically single-
stranded and adopts a much richer variety of secondary and tertiary structures than DNA,
with frequent non canonical motifs and complex folds. ¹⁴⁷,¹⁴⁸ Validation of RNA force fields
is thus an even more challenging task than for DNA, and full validation of our ECC variant
is beyond the scope of the present work. However, we performed preliminary µs-long simu-
lations of a short model double-stranded RNA r(GCACCGUUGG)₂ (1QC0', see Methods),
in its typical A-form helix. In this A-form helix, the phosphate groups on both side of the
RNA deep groove are significantly closer than in the DNA B-form helix,¹⁴⁹,¹⁵⁰ so that one
could expect there a larger impact of scaling the phosphate backbone on the helical geome-
try. However, no significant difference between OL3 and OL3-ECC force fields was found on
backbone dihedral angles nor helical parameters for this model ds-RNA (see Table 7). In this
simple double-stranded RNA form, our OL3-ECC force field yields stable ds-RNA structures,
with very similar structural properties as the original force field, which is currently regarded
as standard—though not perfect—in the community. The OL3-ECC variant thus seems to
preserve the structural qualities of the original force field, while, by design, guaranteeing a
much improved description of the RNA interactions with ions, in particular divalent ones.
With monovalent ions, our ECC force field variant reduced both the number of potassium
ions in direct contact with phosphate groups (see Table 4) and that of solvent-shared ion
pairs —interactions mediated by water molecules. As for ds-DNA, ion locations and occu-
pancies within the grooves are not significantly affected by ECC (Supporting Information
Figure S4).
<table><caption>Table 7: Average of helical parameters and dihedral angles for the 3 replicas of the two force fields studied, compared to the X-ray structure for the studied RNA model system. The RMSD is calculated with respect to the 1QC0 structure.</caption><thead><tr><th></th><th>X-RAY<sup>102</sup></th><th>OL3</th><th>OL3 ECC</th></tr></thead><tbody><tr><td>Shallow minor groove/Å</td><td>15.4</td><td>17.2 ± 0.8</td><td>17.0 ± 0.7</td></tr><tr><td>Deep Major groove/Å</td><td>14.7</td><td>17.2 ± 1.6</td><td>18.3 ± 1.7</td></tr><tr><td>χ/deg</td><td>-162.9</td><td>-153.0 ± 39.9</td><td>-153.4 ± 44.9</td></tr><tr><td>α/deg</td><td>-65.7</td><td>-75.7 ± 17.9</td><td>-75.3 ± 16.3</td></tr><tr><td>β/deg</td><td>170.6</td><td>174.0 ± 9.9</td><td>173.8 ± 9.4</td></tr><tr><td>γ/deg</td><td>56.7</td><td>61.9 ± 18.9</td><td>62.4 ± 18.8</td></tr><tr><td>δ/deg</td><td>79.4</td><td>79.3 ± 9.2</td><td>78.2 ± 8.5</td></tr><tr><td>ε/deg</td><td>-152.4</td><td>-153.6 ± 35.5</td><td>-155.1 ± 28.8</td></tr><tr><td>ζ/deg</td><td>-71.5</td><td>-68.2 ± 9.2</td><td>-67.7 ± 8.8</td></tr><tr><td>shift/Å</td><td>-0.1</td><td>0.0 ± 0.6</td><td>0.0 ± 0.6</td></tr><tr><td>slide/Å</td><td>-1.7</td><td>-1.5 ± 0.5</td><td>-1.7 ± 0.5</td></tr><tr><td>rise/Å</td><td>3.2</td><td>3.3 ± 0.4</td><td>3.3 ± 0.4</td></tr><tr><td>tilt/deg</td><td>-1.2</td><td>0.3 ± 4.6</td><td>0.3 ± 4.5</td></tr><tr><td>roll/deg</td><td>8.1</td><td>10.1 ± 6.9</td><td>8.6 ± 6.7</td></tr><tr><td>twist/deg</td><td>30.8</td><td>30.8 ± 3.9</td><td>30.2 ± 3.9</td></tr><tr><td>shear/Å</td><td>0.1</td><td>0.0 ± 0.4</td><td>0.0 ± 0.4</td></tr><tr><td>buckle/deg</td><td>1.2</td><td>0.2 ± 12.1</td><td>0.1 ± 11.2</td></tr><tr><td>stretch/Å</td><td>-0.3</td><td>-0.1 ± 0.1</td><td>-0.1 ± 0.1</td></tr><tr><td>propeller/deg</td><td>-12.5</td><td>-12.9 ± 9.0</td><td>-11.8 ± 8.9</td></tr><tr><td>stagger/Å</td><td>0.1</td><td>0.0 ± 0.4</td><td>0.0 ± 0.4</td></tr><tr><td>opening/deg</td><td>1.4</td><td>0.0 ± 4.6</td><td>0.0 ± 4.3</td></tr><tr><td>X-displacement/Å</td><td>-4.4</td><td>-4.2 ± 1.3</td><td>-4.4 ± 1.4</td></tr><tr><td>Y-displacement/Å</td><td>-0.1</td><td>0.0 ± 1.3</td><td>0.0 ± 1.2</td></tr><tr><td>inclination/deg</td><td>15.2</td><td>17.8 ± 11.6</td><td>15.5 ± 11.6</td></tr><tr><td>tip/deg</td><td>2.1</td><td>-0.6 ± 8.3</td><td>-0.5 ± 8.4</td></tr><tr><td>helical twist/deg</td><td>32.3</td><td>33.3 ± 4.5</td><td>32.4 ± 4.3</td></tr><tr><td>RMSD (all) /Å</td><td></td><td>1.42 ± 0.31</td><td>1.36 ± 0.35</td></tr><tr><td>RMSD (backbone) /Å</td><td></td><td>1.49 ± 0.36</td><td>1.43 ± 0.40</td></tr></tbody></table>

## ECC impact on divalent ion pairing

Overall, except in the case of G-DNA with a marked improvement in ion retention inside the quadruplex, the OL15-ECC variant does not cause major changes in K⁺ ion pairing with the studied nucleic acid sequences. Scaling charges only slightly decreases the amount of ions forming water-mediated (solvent-shared) ion pairs with the phosphate backbone, while it
has no significant effect on the number of monovalent ions in direct contact with phosphate groups and nucleobases—except in Z-DNA, as previously discussed. As the overbinding artefacts of standard full charge force fields are much more obvious with divalent ions, we expect a much stronger impact of the ECC variant on their pairing with nucleic acids.

Proper characterization of ion pairing to nucleic acid is however much more challenging with divalent ions than with monovalent ions, due to the (very) slow timescales for ligand exchange in their hydration shells, in particular for Mg²⁺ with timescales slower than µs.³⁶,³⁷ This precludes proper sampling of various ion binding arrangements within the typical timescale of our molecular dynamics simulations, even with our ECC force field variant. Supporting Information Fig S21 shows that starting from different ion pairing configurations, very few ion exchanges are observed in the course of 5 µs-long simulations, which demonstrate that plain MD simulations with Mg²⁺ do not properly sample ion pairing arrangements. Extensive characterization of divalent ion binding and their impact on nucleic acid structures would thus require first the development of specific enhanced sampling strategies, which is beyond the scope of the present work.

However, we present here a preliminary analysis of the impact of our OL15-ECC variant on Ca²⁺ ion binding to our model double-stranded B-DNA, the Dickerson-Drew dodecamer. Ligand exchanges are much faster in the Ca²⁺ hydration shell, which allows sampling of multiple ion binding and unbinding events within the course of our µs-long simulations (Supporting information Figures S22 and S23). In sharp contrast with the behavior with monovalent ions, the OL15-ECC scaled charge force field variant now causes a pronounced reduction in the number of direct contact ion pairs with phosphate groups—from 3.2 on average with the standard OL15 full charge force field (out of a total of 11 divalent ions) to only 1.1 with the ECC variant (Table 8). Note that in these simulations we do not see significant formation of contact ion pairs with other groups than phosphates (bases, sugars). In parallel there is also a decrease in the number of solvent-shared ion pairs with phosphate groups, from 2.8 to 1.7. While an extensive analysis of the impact of these interactions with
divalent ions is beyond the scope of the present work, these preliminary results demonstrate the pronounced impact of properly talking into account electronic polarization on divalent ion pairing properties with nucleic acids. The effect is expected to be even much stronger with Mg²⁺—its binding free energy with phosphate groups being even more overestimated by standard force fields (see Table 2).

<table><caption>Table 8: Average number (standard deviation between the 3 replicas in parenthesis) of Ca²⁺ ions in various interaction modes with ds-DNA: direct contact (< 3.2 Å to any nucleic acid heavy atom), direct contact with a phosphate group and solvent-shared ion pair with a phosphate group (3.9 < d(P-Ca²⁺) < 6 Å). The simulated system contains a total of 11 Ca²⁺ ions with an additional 150 mM KCl background.</caption><thead><tr><th></th><th>total direct contacts</th><th>phosphate direct contacts</th><th>SShIP w. phosphates</th></tr></thead><tbody><tr><td>OL15</td><td>2.8 (0.1)</td><td>2.8 (0.1)</td><td>2.8 (0.0)</td></tr><tr><td>OL15-ECC</td><td>1.1 (0.1)</td><td>1.1 (0.1)</td><td>1.7 (0.1)</td></tr></tbody></table>

## Discussion

This work investigates the possibility to improve ion pairing properties of widely used nucleic acid force fields by adopting the ECC charge-scaling strategy, limited to the negatively charged phosphate backbone (5 atoms per nucleotide). Starting from the popular Amber-OL15 force field, we use a scaling factor of 0.8, as suggested in previous work,⁵⁷,⁵⁸ and a charge distribution on the phosphate group optimized to reproduce ion pairing properties of the model DMP compound with both K⁺ and Mg²⁺. We show that the resulting Amber-OL15-ECC force field variant yields a very similar performance to that of the original Amber-OL15 force field and is able to capture conformational substates in various DNA sequences, including non canonical forms such as the challenging Z-DNA. In addition, for DNA G-quadruplexes, Amber-OL15-ECC strongly improves ion retention—and thus, loop stability—leading to a much better agreement between simulations and X-ray structures, at the same cost as usual MD simulations and without the need for more complex explicitly polarizable force fields. Results are also promising for RNA, with preliminary tests on a model double-
stranded A-RNA. How this would affect single-stranded RNA remains to be tested in future
work.

Interestingly, our simulations showed that the impact of scaling the charge of the phosphate backbone seems mostly orthogonal to that of refinements focused on backbone dihedral angles: the same correction can be used for the OL15 force field for DNA, and the OL3 force field for RNA, which differ only by some backbone dihedral angles. Besides, the OL21 force field refined α and γ dihedral angles to better capture Z-DNA substates with respect to OL15, and a similar improvement is observed between their ECC variants. This surprising insensitivity of backbone conformation to phosphate charge may be explained by a compensation of two effects: both the repulsion between phosphate groups and the screening by the background electrolyte are reduced. This result suggests that a similar charge scaling strategy could be used on top of various nucleic acid force fields, for instance Amber-bsc1, the recent OL24 variant for DNA or CHARMM nucleic acid force field, even if the performance of ECC variants would need to be carefully validated before use.

Overall, the use of scaled charges (for both the nucleic acid and counterions) has only a small effect on monovalent ion binding, except inside the DNA G-quadruplex, where the inter-ion repulsion is strongly overestimated in common force fields. The ECC correction leads to a reduced number of both contact and solvent-shared ion pairs with phosphate groups, as consistent with the increase in osmotic coefficient for the model DMP:K solution, getting closer to experimental data. In B-DNA and Z-DNA, we observe a small increase (smaller than 1 excess ion on average in the studied sequences) in the number of K⁺ in direct contact with nucleobases, while the preferred ion locations within grooves remain unchanged. In the absence of quantitative measure of the population of monovalent ions within the grooves, we cannot assess which behavior is more realistic.

In contrast, the ECC charge scaling has a much more pronounced effect on the binding of divalent ions, as we exemplified with preliminary simulations with Ca²⁺. There, ECC caused a significant decrease in the number of both contact and solvent-shared ions pairs, providing
a probably much more realistic (given the computed binding free energies, see Table 2) description of ion pairing in these systems. An even stronger improvement is expected for Mg²⁺ pairing, which is dramatically overestimated by usual force fields. However, the very slow exchange of ligands in the Mg²⁺ hydration shell often hides these overbinding artefacts, as the formation of contact ion pairs (thermodynamically favored by standard full charge force fields) is rarely observed in the timescale of typical MD simulations due to high kinetic barriers. For realistic simulations of divalent ions, there is thus a strong need for efficient enhanced sampling strategies, even with our ECC force field variant.

In addition, binding properties of divalent ions sensitively depend on the chosen charge distribution on the phosphate group, even with a given ECC overall charge scaling factor. Here, we chose a charge distribution that yields an accurate binding free energy for Mg²⁺ ions. We should keep in mind however that the employed experimental reference¹²⁸—which was also used by other groups before us¹²²—pertains to DHP and not DMP, which would be a closer proxy for the DNA backbone, but where no direct experimental measurements are available. We however cannot exclude that DMP may have a slightly different binding free energy, even if we expect the order of magnitude to be similar. Notably, a theoretical model for ion binding at different nucleic acid sites—based on extrapolation of various experimental binding free energies and pKas of related compounds—predicts log stability constants for the formation of a DMP⁻:Mg²⁺ [respectively DMP⁻:Ca²⁺] ion pair of log(K<sub>bind</sub>) = 0.45 [resp. 0.4] for an ionic strength *I* = 0.1 M,¹²⁹ which has been used in previous work⁴⁰,⁴³ as a target for tuning divalent ions binding to phosphates. Once extrapolated to *I* = 0 with a Davies equation,¹⁵¹ this estimate would yield ΔG°<sub>bind</sub> = -3.2 kJ mol⁻¹ for Mg²⁺ and -2.9 kJ/mol for Ca²⁺, weaker than the DHP binding free energies. Hence, in the absence of experimental measurements on DMP directly, the best we can do is probably assess that the real binding free energy is somewhere between those two numbers. This would suggest our charge distribution to yield a very good binding free energy for Ca²⁺ and maybe a bit too strong for Mg²⁺. A yet additional source of uncertainty relates to the mode of ion pairing
these various experimental estimates report on: contact pairing exclusively, as suggested in Ref.,¹²⁹ or the overall binding averaged over different binding modes, as we assumed here? There is thus an acute need for clean and detailed characterization of the relative binding free energies of different binding modes on model compounds to use as targets for force field refinement, from either experiments or high level theoretical work.

Overall, the present work should be viewed as a first promising proof-of-concept of the possibility to improve ion pairing behavior of state-of-the-art nucleic force fields with the ECC0.8 charge scaling strategy applied to the phosphate backbone. While we so far only tested it with a single water model (SPC/E) and provide only a limited set of ECC compatible 0.8-scaled ions to pair it with—K⁺, Na⁺, Ca²⁺ and Mg²⁺, not as thoroughly optimized and tested as the most popular ion force fields, such as those of Merz and coworkers⁴¹ or the Madrid2019 ions by Vega and coworkers⁵⁶—our longer-term goal is to provide the community with a full set of ECC0.8 biomolecular force fields, including in a consistent and compatible way nucleic acids, proteins and phospholipids, validated with a modern water force field and together with a full set of compatible ECC optimized ions. This is a huge already ongoing task,⁶³,¹⁵² and we think that this work demonstrates the potential of this approach for nucleic acids, showing that the ECC strategy can really be used on top of existing force fields, while keeping the benefit of a lot of previous refinements.

## Conclusion

This work is a promising proof-of-concept of the possibility to improve ion pairing properties of non polarizable nucleic acid force fields using the ECC strategy. The physics-based mean-field ECC approach takes implicitly into account electronic polarization—neglected in non-polarizable force fields—via charge scaling. Here, we design an ECC variant of the popular Amber-OL15 force field for nucleic acids, Amber-OL15-ECC, that uses a local charge scaling of the negatively charged phosphate backbone—targeted against experimental data for model
moieties. We show that this variant does not deteriorate the overall performance of the
starting force field on a broad range of structures (ds-DNA, ds-RNA, Z-DNA, G-DNA)
while yielding a much improved description of ion retention in G-quadruplexes and ion-
pairing with divalent ions. Further, the effects of the ECC scaling seem mostly orthogonal to
force field refinements of the backbone dihedral terms, which suggests that a similar strategy
could be used to improve ion-pairing properties of various force fields. While this work was
limited to few counterions with a specific water model (SPC/E), on-going work will aim to
provide the community with a full set of ECC force fields for various ions and biomolecules,
combined with the most popular water models. With divalent ions, the upcoming challenge
will be to derive specific strategies to allow extensive sampling of ion pairing and proper
characterization of their impact on nucleic acid structure and related processes.

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). JPF acknowledges a PhD grant from the DYNAMO Labex. EDD and LK acknowledge support from the ANR JCJC MUSIRICAT (ANR-22-CE29-0004-01). This project was provided with computing HPC and storage resources by GENCI at CINES, TGCC and IDRIS (projects A0150711021 and A0170711021). We thank Philip E. Mason (IOCB, Prague) for his help with comparing electrolyte structure with neutron scattering signals, and the whole Jungwirth group (IOCB, Prague) for stimulating exchanges.

**Supporting Information Available**

Supporting Information is available, with additional computational details and complemen-
tary analyses. In addition, input and trajectory files for all simulations are shared on a
Zenodo public folder,¹¹² and parameter files for our Amber-OL15-ECC force field are shared
on a public GitHub folder.¹¹¹
References

(1) Lipfert, J.; Doniach, S.; Das, R.; Herschlag, D. *Understanding Nucleic Acid-Ion Interactions*.; 2014; Vol. 83.

(2) Hud, N. V.; Polak, M. DNA-cation Interactions: The Major and Minor Grooves Are Flexible Ionophores. *Curr. Opin. Struct. Biol.* **2001**, *11*, 293-301.

(3) McFail-Isom, L.; Sines, C. C.; Williams, L. D. DNA Structure: Cations in Charge? *Curr. Opin. Struct. Biol.* **1999**, *9*, 298-304.

(4) Nishimura, Y.; Torigoe, C.; Tsuboi, M. Salt Induced B-A Transition of Poly(dG).Poly(dC) and the Stabilization of A Form by Its Methylation. *Nucleic Acids Res.* **1986**, *14*, 2737-2748.

(5) Jacobson, D. R.; Saleh, O. A. Counting the Ions Surrounding Nucleic Acids. *Nucleic Acids Res.* **2017**, *45*, 1596-1605.

(6) Shui, X.; Sines, C. C.; McFail-Isom, L.; VanDerveer, D.; Williams, L. D. Structure of the Potassium Form of CGCGAATTCGCG: DNA Deformation by Electrostatic Collapse around Inorganic Cations. *Biochemistry* **1998**, *37*, 16877-16887.

(7) Shui, X.; Mcfail-isom, L.; Hu, G. G.; Williams, L. D. The B-DNA Dodecamer at High Resolution Reveals a Spine of Water on Sodium The B-DNA Dodecamer at High Resolution Reveals a Spine of Water on Sodium. *Structure* **1998**, *37*, 8341-8355.

(8) Denisov, V. P.; Halle, B. Sequence-Specific Binding of Counterions to B-DNA. *Proc. Natl. Acad. Sci. U.S.A.* **2000**, *97*, 629-633.

(9) Tereshko, V.; Wilds, C. J.; Minasov, G.; Prakash, T. P.; Maier, M. A.; Howard, A.; Wawrzak, Z.; Manoharan, M.; Egli, M. Detection of Alkali Metal Ions in DNA Crystals Using State-of-the-Art X-ray Diffraction Experiments. *Nucleic Acids Res.* **2001**, *29*, 1208-1215.

(10) Auffinger, P.; D'Ascenzo, L.; Ennifar, E. *The Alkali Metal Ions: Their Role for Life*; Springer International Publishing: Cham, 2016; Vol. 16; pp 167-201.
(11) Minasov, G.; Tereshko, V.; Egli, M. Atomic-Resolution Crystal Structures of B-DNA Reveal Specific Influences of Divalent Metal Ions on Conformation and Packing. *J. Mol. Biol.* **1999**, *291*, 83–99.

(12) Leonarski, F.; Henning-Knechtel, A.; Kirmizialtin, S.; Ennifar, E.; Auffinger, P. Principles of Ion Binding to RNA Inferred from the Analysis of a 1.55 Å Resolution Bacterial Ribosome Structure – Part I: Mg²⁺. *Nucleic Acids Res.* **2025**, 53, gkae1148.

(13) Kolev, S. K.; Petkov, P. S.; Rangelov, M. A.; Trifonov, D. V.; Milenov, T. I.; Vayssilov, G. N. Interaction of Na⁺, K⁺, Mg²⁺ and Ca²⁺ Counter Cations with RNA. *Metallomics* **2018**, *10*, 659–678.

(14) Guéroult, M.; Boittin, O.; Mauffret, O.; Etchebest, C.; Hartmann, B. Mg²⁺ in the Major Groove Modulates B-DNA Structure and Dynamics. *PLoS ONE* **2012**, *7*, e41704.

(15) Auffinger, P.; Ennifar, E.; D'Ascenzo, L. Deflating the RNA Mg2⁺ Bubble. Stereochemistry to the Rescue! *RNA* **2021**, *27*, 243–252.

(16) Palù, G.; Valisena, S.; Ciarrocchi, G.; Gatto, B.; Palumbo, M. Quinolone Binding to DNA Is Mediated by Magnesium Ions. *Proc. Natl. Acad. Sci. U.S.A.* **1992**, *89*, 9671–9675.

(17) Draper, D. E. RNA Folding: Thermodynamic and Molecular Descriptions of the Roles of Ions. *Biophys. J.* **2008**, *95*, 5489–95.

(18) Bowman, J. C.; Lenz, T. K.; Hud, N. V.; Williams, L. D. Cations in Charge: Magnesium Ions in RNA Folding and Catalysis. *Curr. Opin. Struc. Biol.* **2012**, *22*, 262–272.

(19) Woodson, S. A. Metal Ions and RNA Folding: A Highly Charged Topic with a Dynamic Future. *Curr. Opin. Chem. Biol.* **2005**, *9*, 104–109.

(20) Draper, D. E.; Grilley, D.; Soto, A. M. Ions and RNA Folding. *Annu. Rev. Biophys. Biochem. Struct.* **2005**, *34*, 221–243.
(21) Yang, Z.; Hayes, J. J. The Divalent Cations Ca²⁺ and Mg²⁺ Play Specific Roles in Stabilizing Histone-DNA Interactions within Nucleosomes That Are Partially Redundant with the Core Histone Tail Domains. *Biochemistry* **2011**, *50*, 9973–9981.

(22) Yang, W. Nucleases: Diversity of Structure, Function and Mechanism. *Q. Rev. Biophys.* **2011**, *44*, 1–93.

(23) Palermo, G.; Cavalli, A.; Klein, M. L.; Alfonso-Prieto, M.; Dal Peraro, M.; De Vivo, M. Catalytic Metal Ions and Enzymatic Processing of DNA and RNA. *Acc. Chem. Res.* **2015**, *48*, 220–228.

(24) Hanna, R.; Doudna, J. A. Metal Ions in Ribozyme Folding and Catalysis. *Curr. Opin. Chem. Biol.* **2000**, *4*, 166–170.

(25) MacKerell, A. D.; Banavali, N.; Foloppe, N. Development and Current Status of the CHARMM Force Field for Nucleic Acids. *Biopolymers* **2000**, *56*, 257–265.

(26) Zgarbová, M.; Otyepka, M.; Šponer, J.; Mládek, A.; Banáš, P.; Cheatham, T. E.; Jurečka, P. Refinement of the Cornell et al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles. *J. Chem. Theory Comput.* **2011**, *7*, 2886–2902.

(27) Ivani, I. et al. Parmbsc1: A Refined Force Field for DNA Simulations. *Nat. Methods* **2016**, *13*, 55–58.

(28) Liebl, K.; Zacharias, M. The Development of Nucleic Acids Force Fields: From an Unchallenged Past to a Competitive Future. *Biophys. J.* **2023**, *122*, 2841–2851.

(29) Yoo, J.; Aksimentiev, A. Improved Parametrization of Li⁺, Na⁺, K⁺, and Mg²⁺ Ions for All-Atom Molecular Dynamics Simulations of Nucleic Acid Systems. *J. Phys. Chem. Lett.* **2012**, *3*, 45–50.

(30) Jayaram, B.; Beveridge, D. L. Modeling DNA in Aqueous Solutions: Theoretical and Computer Simulation Studies on the Ion Atmosphere of DNA. *Annu. Rev. Biophys. Biomol. Struct.* **1996**, *25*, 367–394.
(31) Rueda, M.; Cubero, E.; Laughton, C. A.; Orozco, M. Exploring the Counterion Atmosphere around DNA: What Can Be Learned from Molecular Dynamics Simulations? *Biophys. J.* **2004**, *87*, 800–811.

(32) Noy, A.; Soteras, I.; Luqueb, F. J.; Orozco, M. The Impact of Monovalent Ion Force Field Model in Nucleic Acids Simulations. *Phys. Chem. Chem. Phys.* **2009**, *11*, 10541–10542.

(33) Dans, P. D.; Danilāne, L.; Ivani, I.; Dršata, T.; Lankaš, F.; Hospital, A.; Walther, J.; Pujagut, R. I.; Battistini, F.; Gelpí, J. L.; Lavery, R.; Orozco, M. Long-Timescale Dynamics of the Drew–Dickerson Dodecamer. *Nucleic Acids Res.* **2016**, *44*, 4052–4066.

(34) Ditzler, M. A.; Otyepka, M.; Šponer, J.; Walter, N. G. Molecular Dynamics and Quantum Mechanics of RNA: Conformational and Chemical Change We Can Believe In. *Acc. Chem. Res.* **2010**, *43*, 40–47.

(35) Šponer, J.; Banáš, P.; Jurečka, P.; Zgarbová, M.; Kührová, P.; Havrila, M.; Krepl, M.; Stadlbauer, P.; Otyepka, M. Molecular Dynamics Simulations of Nucleic Acids. from Tetranucleotides to the Ribosome. *J. Phys. Chem. Lett.* **2014**, *5*, 1771–1782.

(36) Bleuzen, A.; Helm, L.; Merbach, E. Water Exchange on Magnesium (II) in Aqueous Solution: A Variable Temperature and Pressure. *Magn. Reson. Chem.* **1997**, *35*, 765–773.

(37) Neely, J.; Connick, R. E. The Rate of Water Exchange from Hydrated Magnesium. *J. Am. Chem. Soc.* **1970**, *92*, 3476–3478.

(38) Babu, C. S.; Lim, C. Empirical Force Fields for Biologically Active Divalent Metal Cations in Water. *J. Phys. Chem. A* **2006**, *110*, 691–699.

(39) Allner, O.; Nilsson, L.; Villa, A. Magnesium Ion - Water Coordination and Exchange in Biomolecular Simulations. *J. Chem. Theory Comput.* **2012**, *8*, 1493–1502.

(40) Panteva, M. T.; Giambasu, G. M.; York, D. M. Force Field for Mg²⁺, Mn²⁺, Zn²⁺, and Cd²⁺ Ions That Have Balanced Interactions with Nucleic Acids. *J. Phys. Chem. B* **2015**, *119*, 15460–15470.
(41) Li, P.; Merz, K. M. Taking into Account the Ion-Induced Dipole Interaction in the Nonbonded Model of Ions. *J. Chem. Theory Comput.* **2014**, *10*, 289–297.

(42) Mamatkulov, S.; Schwierz, N. Force Fields for Monovalent and Divalent Metal Cations in TIP3P Water Based on Thermodynamic and Kinetic Properties. *J. Chem. Phys.* **2018**, *148*, 074504.

(43) Grotz, K. K.; Cruz-León, S.; Schwierz, N. Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation, Ion-Binding, and Water-Exchange Properties. *J. Chem. Theory Comput.* **2021**, *17*, 2530–2540.

(44) Koca Fındık, B.; Jafari, M.; Song, L. F.; Li, Z.; Aviyente, V.; Merz, K. M. Binding of Phosphate Species to Ca²⁺ and Mg²⁺ in Aqueous Solution. *J. Chem. Theory Comput.* **2024**, *20*, 4298–4307.

(45) Panteva, M. T.; Giambaşu, G. M.; York, D. M. Comparison of Structural, Thermodynamic, Kinetic and Mass Transport Properties of Mg²⁺ Ion Models Commonly Used in Biomolecular Simulations. *J. Comput. Chem.* **2015**, *36*, 970–982.

(46) Zhang, C.; Lu, C.; Jing, Z.; Wu, C.; Piquemal, J.-P.; Ponder, J. W.; Ren, P. AMOEBA Polarizable Atomic Multipole Force Field for Nucleic Acids. *J. Chem. Theory Comput.* **2018**, *14*, 2084–2108.

(47) Lemkul, J. A.; MacKerell, A. D. Polarizable Force Field for RNA Based on the Classical Drude Oscillator. *J. Comput. Chem.* **2018**, *39*, 2624–2646.

(48) Lemkul, J. A.; MacKerell, A. D. 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. *J. Phys. Chem. B* **2016**, *120*, 11436–11448.

(49) Ratnasinghe, B. D.; Salsbury, A. M.; Lemkul, J. A. Ion Binding Properties and Dynamics of the Bcl-2 G-Quadruplex Using a Polarizable Force Field. *J. Chem. Inf. Model.* **2020**, *60*, 6476–6488.
(50) Lagardère, L.; Jolly, L. H.; Lipparini, F.; Aviat, F.; Stamm, B.; Jing, Z. F.; Harger, M.; Torabifard, H.; Cisneros, G. A.; Schnieders, M. J.; Gresh, N.; Maday, Y.; Ren, P. Y.; Ponder, J. W.; Piquemal, J. P. Tinker-HP: A Massively Parallel Molecular Dynamics Package for Multiscale Simulations of Large Complex Systems with Advanced Point Dipole Polarizable Force Fields. *Chem. Sci.* **2018**, *9*, 956–972.

(51) Yoo, J.; Aksimentiev, A. New Tricks for Old Dogs: Improving the Accuracy of Biomolecular Force Fields by Pair-Specific Corrections to Non-Bonded Interactions. *Phys. Chem. Chem. Phys.* **2018**, *20*, 8432–8449.

(52) Yoo, J.; Aksimentiev, A. Competitive Binding of Cations to Duplex DNA Revealed through Molecular Dynamics Simulations. *J. Phys. Chem. B* **2012**, *116*, 12946–12954.

(53) Leontyev, I.; Stuchebrukhov, A. Accounting for Electronic Polarization in Non-Polarizable Force Fields. *Phys. Chem. Chem. Phys.* **2011**, *13*, 2613–2626.

(54) Kirby, B. J.; Jungwirth, P. Charge Scaling Manifesto: A Way of Reconciling the Inherently Macroscopic and Microscopic Natures of Molecular Simulations. *J. Phys. Chem. Lett.* **2019**, *10*, 7531–7536.

(55) Kostal, V.; Jungwirth, P.; Martinez-Seara, H. Nonaqueous Ion Pairing Exemplifies the Case for Including Electronic Polarization in Molecular Dynamics Simulations. *J. Phys. Chem. Lett.* **2023**, *14*, 8691–8696.

(56) Zeron, I.; Abascal, J. L. F.; Vega, C. A Force Field of Li, Na, K, Mg, Ca, in Aqueous Solution Based on the TIP4P / 2005 Water Model and Scaled Charges for the Ions in Aqueous Solution Based on the TIP4P / 2005 Water Model and Scaled Charges for the Ions. *J. Chem. Phys.* **2019**, *151*, 134504.

(57) De Oliveira, D. M.; Zukowski, S. R.; Palavec, V.; Hénin, J.; Martinez-Seara, H.; Ben-Amotz, D.; Jungwirth, P.; Duboué-Dijon, E. Binding of Divalent Cations to Acetate: Molecular Simulations Guided by Raman Spectroscopy. *Phys. Chem. Chem. Phys.* **2020**, *22*, 24014–24027.
(58) Puyo-Fourtine, J.; Juillé, M.; Hénin, J.; Clavaguéra, C.; Duboué-Dijon, E. A Consistent Picture of Phosphate – Divalent Cation Binding from Models with Implicit and Explicit Electronic Polarization. *J. Phys. Chem. B* **2022**, *126*, 4022–4034.

(59) Kohagen, M.; Pluharová, E.; Mason, P. E.; Jungwirth, P. Exploring Ion-Ion Interactions in Aqueous Solutions by a Combination of Molecular Dynamics and Neutron Scattering. *J. Phys. Chem. Lett.* **2015**, *6*, 1563–1567.

(60) Duboué-Dijon, E.; Javanainen, M.; Delcroix, P.; Jungwirth, P.; Martinez-Seara, H. A Practical Guide to Biologically Relevant Molecular Simulations with Charge Scaling for Electronic Polarization. *J. Chem. Phys.* **2020**, *153*, 050901.

(61) Melcr, J.; Martinez-Seara, H.; Nencini, R.; Kolafa, J.; Jungwirth, P.; Ollila, O. H. Accurate Binding of Sodium and Calcium to a POPC Bilayer by Effective Inclusion of Electronic Polarization. *J. Phys. Chem. B* **2018**, *122*, 4546–4557.

(62) Duboué-Dijon, E.; Delcroix, P.; Martinez-Seara, H.; Hladilkova, J.; Coufal, P.; Krizek, T.; Jungwirth, P. Binding of Divalent Cations to Insulin: Capillary Electrophoresis and Molecular Simulations. *J. Phys. Chem. B* **2018**, *122*, 5640–5648.

(63) Nencini, R.; Tempra, C.; Biriukov, D.; Riopedre-Fernandez, M.; Cruces Chamorro, V.; Polák, J.; Mason, P. E.; Ondo, D.; Heyda, J.; Ollila, O. H. S.; Jungwirth, P.; Javanainen, M.; Martinez-Seara, H. Effective Inclusion of Electronic Polarization Improves the Description of Electrostatic Interactions: The prosECCo75 Biomolecular Force Field. *J. Chem. Theory Comput.* **2024**, *20*, 7546–7559.

(64) Hui, C.; De Vries, R.; Kopec, W.; De Groot, B. L. Effective Polarization in Potassium Channel Simulations: Ion Conductance, Occupancy, Voltage Response, and Selectivity. *Proc. Natl. Acad. Sci. USA* **2025**, *122*, e2423866122.

(65) Bruce, E. E.; van der Vegt, N. F. A. Does an Electronic Continuum Correction Improve Effective Short-Range Ion-Ion Interactions in Aqueous Solution? *J. Chem. Phys.* **2018**, *148*, 222816.
(66) Tan, D.; Piana, S.; Dirks, R. M.; Shaw, D. E. RNA Force Field with Accuracy Comparable to State-of-the-Art Protein Force Fields. *Proc. Natl. Acad. Sci. USA* **2018**, *201713027*.

(67) Piana, S.; Robustelli, P.; Tan, D.; Chen, S.; Shaw, D. E. Development of a Force Field for the Simulation of Single-Chain Proteins and Protein-Protein Complexes. *J. Chem. Theory Comput.* **2020**, *16*, 2494–2507.

(68) Zgarbová, M.; Šponer, J.; Otyepka, M.; Cheatham, T. E.; Galindo-Murillo, R.; Jurečka, P. Refinement of the Sugar–Phosphate Backbone Torsion Beta for AMBER Force Fields Improves the Description of Z- and B-DNA. *J. Chem. Theory Comput.* **2015**, *11*, 5723–5736.

(69) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. *J. Phys. Chem.* **1987**, *91*, 6269–6271.

(70) Kuhrova, P.; Otyepka, M.; Sponer, J.; Banáš, P. Are Waters around RNA More than Just a Solvent?–An Insight from Molecular Dynamics Simulations. *J. Chem. Theory Comput.* **2014**, *10*, 401–411.

(71) Joung, I. S.; Cheatham, T. E. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. *J. Phys. Chem. B* **2008**, *112*, 9020–9041.

(72) Kohagen, M.; Mason, P. E.; Jungwirth, P. Accounting for Electronic Polarization Effects in Aqueous Sodium Chloride via Molecular Dynamics Aided by Neutron Scattering. *J. Phys. Chem. B* **2016**, *120*, 1454–1460.

(73) Duboué-Dijon, E.; Mason, P. E.; Fischer, H. E.; Jungwirth, P. Hydration and Ion Pairing in Aqueous Mg²⁺ and Zn²⁺ Solutions: Force Field Description Aided by Neutron Scattering Experiments and Ab Initio Molecular Dynamics Simulations. *J. Phys. Chem. B* **2017**, *122*, 3296–3306.

(74) Laage, D.; Stirnemann, G. Effect of Ions on Water Dynamics in Dilute and Concentrated Aqueous Salt Solutions. *J. Phys. Chem. B* **2019**, *123*, 3312–3324.
(75) Mason, P. E.; Tavagnacco, L.; Saboungi, M.-L.; Hansen, T.; Fischer, H. E.; Neilson, G. W.; Ichiye, T.; Brady, J. W. Molecular Dynamics and Neutron Scattering Studies of Potassium Chloride in Aqueous Solution. *J. Phys. Chem. B* **2019**, 123, 10807–10813.

(76) Luo, Y.; Roux, B. Simulation of Osmotic Pressure in Concentrated Aqueous Salt Solutions. *J. Phys. Chem. Lett.* **2010**, *1*, 183–189.

(77) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. *J. Chem. Phys.* **2007**, *126*, 014101.

(78) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. *J. Appl. Phys.* **1981**, *52*, 7182–7190.

(79) Mobley, D. L.; Klimovich, P. V. Perspective: Alchemical Free Energy Calculations for Drug Discovery. *J. Chem. Phys.* **2012**, *137*, 230901.

(80) Mey, A. S. J. S.; Allen, B.; Macdonald, H. E. B.; Chodera, J. D.; Kuhn, M.; Michel, J.; Mobley, D. L.; Naden, L. N.; Prasad, S.; Rizzi, A.; Scheen, J.; Shirts, M. R.; Tresadern, G.; Xu, H. Best Practices for Alchemical Free Energy Calculations. *Living J. Comp. Mol. Sci.* **2020**, 1–51.

(81) Mobley, D. L.; Gilson, M. K. Predicting Binding Free Energies: Frontiers and Benchmarks. *Annu. Rev. Biophys.* **2017**, *46*, 531–558.

(82) Song, L. F.; Merz, K. M. Evolution of Alchemical Free Energy Methods in Drug Discovery. *J. Chem. Inf. Model.* **2020**, *60*, 5308–5318.

(83) Deng, Y.; Roux, B. Computations of Standard Binding Free Energies with Molecular Dynamics Simulations. *J. Phys. Chem. B* **2009**, *113*, 2234–2246.

(84) Gilson, M. K.; Given, J. A.; Bush, B. L.; McCammon, J. A. The Statistical-Thermodynamic Basis for Computation of Binding Affinities: A Critical Review. *Biophys. J.* **1997**, *72*, 1047–1069.
(85) Mobley, D. L.; Chodera, J. D.; Dill, K. A. On the Use of Orientational Restraints and Symmetry Corrections in Alchemical Free Energy Calculations. *J. Chem. Phys.* **2006**, *125*, 084902.

(86) Pranata, J.; Jorgensen, W. L. Monte Carlo Simulations Yield Absolute Free Energies of Binding Theory and the OPLS Potential Functions Lar Dynamics. *Tetrahedron* **1991**, *41*, 2491-2501.

(87) Hermans, J.; Shankar, S. The Free Energy of Xenon Binding to Myoglobin from Molecular Dynamics Simulation. *Israel J. Chem.* **1986**, *27*, 225-227.

(88) Simonson, T.; Roux, B. Concepts and Protocols for Electrostatic Free Energies. *Mol. Sim.* **2016**, *42*, 1090-1101.

(89) Salari, R.; Joseph, T.; Lohia, R.; Hénin, J.; Brannigan, G. A Streamlined, General Approach for Computing Ligand Binding Free Energies and Its Application to GPCR-Bound Cholesterol. *J. Chem. Theory Comput.* **2018**, *14*, 6560-6573.

(90) Duboué-Dijon, E.; Hénin, J. Building Intuition for Binding Free Energy Calculations: Bound State Definition, Restraints, and Symmetry. *J. Chem. Phys.* **2021**, *154*, 204101.

(91) Bennett, C. H. Efficient Estimation of Free Energy Differences from Monte Carlo Data. *J. Comput. Phys.* **1976**, *22*, 245-268.

(92) Drew, H. R.; Dickerson, R. E. Structure of a B-DNA Dodecamer. *J. Mol. Biol.* **1981**, *151*, 535-556.

(93) Wu, Z.; Delaglio, F.; Tjandra, N.; Zhurkin, V. B. Overall Structure and Sugar Dynamics of a DNA Dodecamer from Homo- and Heteronuclear Dipolar Couplings and ³¹P Chemical Shift Anisotropy. *J. Biomol. NMR* **26**, 297-315.

(94) Dauter, Z.; Adamiak, D. A. Anomalous Signal of Phosphorus Used for Phasing DNA Oligomer: Importance of Data Redundancy. *Acta Crystallogr. D Biol. Crystallogr.* **2001**, *57*, 990-995.
(95) Heddi, B.; Foloppe, N.; Oguey, C.; Hartmann, B. Importance of Accurate DNA Structures in Solution: The Jun-Fos Model. *J. Mol. Biol.* **2008**, *382*, 956-970.

(96) Bansal, M.; Bhattacharyya, D.; Ravi, B. NUPARM and NUCGEN: Software for Analysis and Generation of Sequence Dependent Nucleic Acid Structures. *Bioinformatics* **1995**, *11*, 281-287.

(97) Woods, K. K.; Maehigashi, T.; Howerton, S. B.; Sines, C. C.; Tannenbaum, S.; Williams, L. D. High-Resolution Structure of an Extended A-Tract: [*D(CGCAAATTTGCG)]<sub>2</sub>. *J. Am. Chem. Soc.* **2004**, *126*, 15330-15331.

(98) Gyi, J. I.; Lane, A. N.; Conn, G. L.; Brown, T. Solution Structures of DNA·RNA Hybrids with Purine-Rich and Pyrimidine-Rich Strands: Comparison with the Homologous DNA and RNA Duplexes. *Biochemistry* **1998**, *37*, 73-80.

(99) Štefl, R.; Trantírek, L.; Vorlíčková, M.; Koča, J.; Sklenář, V.; Kypr, J. A-like Guanine-Guanine Stacking in the Aqueous DNA Duplex of d(GGGGCCCC). *J. Mol. Biol.* **2001**, *307*, 513-524.

(100) *Modeling Unusual Nucleic Acid Structures*; ACS Symposium Series; American Chemical Society: Washington, DC, 1997; Vol. 24; pp 379-393.

(101) Gyi, J. I.; Lane, A. N.; Conn, G. L.; Brown, T. Solution Structures of DNA·RNA Hybrids with Purine-Rich and Pyrimidine-Rich Strands: Comparison with the Homologous DNA and RNA Duplexes. *Biochemistry* **1998**, *37*, 73-80.

(102) Klosterman, P. S.; Shah, S. A.; Steitz, T. A. Crystal Structures of Two Plasmid Copy Control Related RNA Duplexes: An 18 Base Pair Duplex at 1.20 Å Resolution and a 19 Base Pair Duplex at 1.55 Å Resolution. *Biochemistry* **1999**, *38*, 14784-14792.

(103) Berendsen, H. J.; van der Spoel, D.; van Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. *Comput. Phys. Comm.* **1995**, *91*, 43-56.
(104) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. Gromacs: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. *SoftwareX* **2015**, *1-2*, 19-25.

(105) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An Nlog(N) Method for Ewald Sums in Large Systems. *J. Chem. Phys.* **1993**, *98*, 10089-10092.

(106) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. *J. Comput. Chem.* **1997**, *18*, 1463-1472.

(107) Miyamoto, S.; Kollman, P. A. SETTLE: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. *J. Comput. Chem.* **1992**, *13*, 952-962.

(108) Pérez, A.; Marchán, I.; Svozil, D.; Sponer, J.; Cheatham, T. E.; Laughton, C. A.; Orozco, M. Refinement of the AMBER Force Field for Nucleic Acids: Improving the Description of α/γ Conformers. *Biophys. J.* **2007**, *92*, 3817-3829.

(109) Krepl, M.; Zgarbová, M.; Stadlbauer, P.; Otyepka, M.; Banáš, P.; Koča, J.; Cheatham, T. E.; Jurečka, P.; Šponer, J. Reference Simulations of Noncanonical Nucleic Acids with Different χ Variants of the AMBER Force Field: Quadruplex DNA, Quadruplex RNA, and Z-DNA. *J. Chem. Theory Comput.* **2012**, *8*, 2506-2520.

(110) Zgarbová, M.; Luque, F. J.; Šponer, J.; Cheatham, T. E.; Otyepka, M.; Jurečka, P. Toward Improved Description of DNA Backbone: Revisiting Epsilon and Zeta Torsion Force Field Parameters. *J. Chem. Theory Comput.* **2013**, *9*, 2339-2354.

(111) https://github.com/elisedubouedijon/Amber-ECC-ffs.

(112) Zenodo Dataset DOI:10.5281/Zenodo.15230600.

(113) Li, P.; Roberts, B. P.; Chakravorty, D. K.; Merz, K. M. Rational Design of Particle Mesh Ewald Compatible Lennard-Jones Parameters for +2 Metal Cations in Explicit Solvent. *J. Chem. Theory Comput.* **2013**, *9*, 2733-2748.
(114) Roe, D. R.; Cheatham, T. E. PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. *J. Chem. Theory Comput.* **2013**, *9*, 3084–3095.

(115) Hartmann, B.; Piazzola, D.; Lavery, R.; Thorique, L. D. B.; Ura, C.; Physico-chimique, I. D. B. BI - BII Transitions in B-DNA. *Nucl. Acids Res.* **1993**, *21*, 561–568.

(116) Černý, J.; Božíková, P.; Svoboda, J.; Schneider, B. A Unified Dinucleotide Alphabet Describing Both RNA and DNA Structures. *Nucleic Acids Res.* **2020**, *48*, 6367–6381.

(117) Michaud-Agrawal, N.; Denning, E. J.; Woolf, T. B.; Beckstein, O. MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. *J. Comput. Chem.* **2011**, *32*, 2319–2327.

(118) Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. *J. Molec. Graph-ics* **1996**, *14*, 33–38.

(119) Marincola, F. C.; Denisov, V. P.; Halle, B. Competitive Na⁺ and Rb⁺ Binding in the Minor Groove of DNA. *J. Am. Chem. Soc.* **2004**, *126*, 6739–6750.

(120) Lavery, R.; Maddocks, J. H.; Pasi, M.; Zakrzewska, K. Analyzing Ion Distributions around DNA. *Nucl. Acids Res.* **2014**, *42*, 8138–8149.

(121) Walker, B.; Jing, Z.; Ren, P. Molecular Dynamics Free Energy Simulations of ATP : Mg and ADP : Mg Using the Polarisable Force Field AMOEBA. *Mol. Sim.* **2020**, *0*, 1–10.

(122) Kumar, M.; Simonson, T.; Ohanessian, G.; Clavaguéra, C. Structure and Thermodynamics of Mg:Phosphate Interactions in Water: A Simulation Study. *ChemPhysChem* **2015**, *16*, 658–665.

(123) Tamaki, K.; Suga, K.; Tanihara, E. Solution Properties of Dialkyl Phosphate Salts. Apparent Molar Volumes, Viscosity *B* Coefficients, Heats of Solution, and Osmotic Coefficients. *BCSJ* **1987**, *60*, 1225–1229.

(124) Rasaiah, J. C.; Friedman, H. L. Integral Equation Methods in the Computation of Equilibrium Properties of Ionic Solutions. *J. Chem. Phys.* **1968**, *48*, 2742–2752.
(125) Kalcher, I.; Dzubiella, J. Structure-Thermodynamics Relation of Electrolyte Solutions. *J. Chem. Phys.* **2009**, *130*, 134507.

(126) Hess, B.; Holm, C.; Van Der Vegt, N. Osmotic Coefficients of Atomistic NaCl (Aq) Force Fields. *J. Chem. Phys.* **2006**, *124*, 164509.

(127) Fyta, M.; Kalcher, I.; Dzubiella, J.; Vrbka, L.; Netz, R. R. Ionic Force Field Optimization Based on Single-Ion and Ion-Pair Solvation Properties. *J. Chem. Phys.* **2010**, *132*, 024911.

(128) Verbeeck, R. M.; De Bruyne, P. A.; Driessens, F. C.; Verbeek, F. Solubility of Magnesium Hydrogen Phosphate Trihydrate and Ion-Pair Formation in the System Mg(OH)₂-H₃PO₄-H₂O at 25°C. *Inorg. Chem.* **1984**, *23*, 1922–1926.

(129) Sigel, R. K.; Sigel, H. A Stability Concept for Metal Ion Coordination to Single-Stranded Nucleic Acids and Affinities of Individual Sites. *Acc. Chem. Res.* **2010**, *43*, 974–984.

(130) Mcdowell, H.; Brown, W. E.; Sutter, J. R. Solubility Study of Calcium Hydrogen Phosphate. Ion-Pair Formation. *Inorg. Chem.* **1971**, *10*, 1638–1643.

(131) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera—A Visualization System for Exploratory Research and Analysis. *J. Comput. Chem.* **2004**, *25*, 1605–1612.

(132) Tereshko, V.; Minasov, G.; Egli, M. The Dickerson-Drew Dodecamer Revisited at Atomic Resolution. *J. Am. Chem. Soc.* **1999**, *121*, 470–471.

(133) Schwieters, C. D.; Clore, G. M. A Physical Picture of Atomic Motions within the Dickerson DNA Dodecamer in Solution Derived from Joint Ensemble Refinement against NMR and Large-Angle X-ray Scattering Data. *Biochemistry* **2007**, *46*, 1152–1166.

(134) Heddi, B.; Oguey, C.; Lavelle, C.; Foloppe, N.; Hartmann, B. Intrinsic Flexibility of B-DNA: The Experimental TRX Scale. *Nucl. Acids Res.* **2010**, *38*, 1034–1047.
(135) Tian, Y.; Kayatta, M.; Shultis, K.; Gonzalez, A.; Mueller, L. J.; Hatcher, M. E. ³¹P NMR Investigation of Backbone Dynamics in DNA Binding Sites. *J. Phys. Chem. B* **2009**, *113*, 2596-2603.

(136) Spiegel, J.; Adhikari, S.; Balasubramanian, S. The Structure and Function of DNA G-Quadruplexes. *Trends in Chemistry* **2020**, *2*, 123-136.

(137) Monsen, R. C.; Trent, J. O.; Chaires, J. B. G-Quadruplex DNA: A Longer Story. *Acc. Chem. Res.* **2022**, *55*, 3242-3252.

(138) Fadrná, E.; Špačková, N.; Sarzyňska, J.; Koča, J.; Orozco, M.; Cheatham, T. E.; Kulinski, T.; Šponer, J. Single Stranded Loops of Quadruplex DNA As Key Benchmark for Testing Nucleic Acids Force Fields. *J. Chem. Theory Comput.* **2009**, *5*, 2514-2530.

(139) Havrila, M.; Stadlbauer, P.; Islam, B.; Otyepka, M.; Šponer, J. Effect of Monovalent Ion Parameters on Molecular Dynamics Simulations of G-Quadruplexes. *J. Chem. Theory Comput.* **2017**, *13*, 3911-3926.

(140) Li, N.; Gao, Y.; Qiu, F.; Zhu, T. Benchmark Force Fields for the Molecular Dynamic Simulation of G-Quadruplexes. *Molecules* **2021**, *26*, 5379.

(141) Zgarbová, M.; Šponer, J.; Jurečka, P. Z-DNA as a Touchstone for Additive Empirical Force Fields and a Refinement of the Alpha/Gamma DNA Torsions for AMBER. *J. Chem. Theory Comput.* **2021**, *17*, 6292-6301.

(142) Wang, A. H.-J.; Quigley, G. J.; Kolpak, F. J.; van der Marel, G.; van Boom, J. H.; Rich, A. Left-Handed Double Helical DNA: Variations in the Backbone Conformation. *Science* **1981**, *211*, 171-176.

(143) Rauch, C.; Pichler, A.; Trieb, M.; Wellenzohn, B.; Liedl, K. R.; Mayer, E. Z-DNA's Conformer Substates Revealed by FT-IR Difference Spectroscopy of Nonoriented Left-Handed Double Helical Poly(dG-dC). *J. Biomol. Struct. Dyn.* **2005**, *22*, 595-614.

(144) Zgarbová, M.; Otyepka, M.; Sponer, J.; Lankas, F.; Jurečka, P. Base Pair Fraying in Molecular Dynamics Simulations of DNA and RNA. *J. Chem. Theory Comput.* **2014**, *10*, 3177-3189.
(145) Zgarbová, M.; Jurečka, P.; Šponer, J.; Otyepka, M. A- to B-DNA Transition in AMBER Force Fields and Its Coupling to Sugar Pucker. *J. Chem. Theory Comput.* **2018**, *14*, 319–328.

(146) Zgarbová, M.; Šponer, J.; Jurečka, P. Refinement of the Sugar Puckering Torsion Potential in the AMBER DNA Force Field. *J. Chem. Theory Comput.* **2025**, *21*, 833–846.

(147) Conn, G. L.; Draper, D. E. RNA Structure. *Curr. Opin. Struc. Biol.* **1998**, *8*, 278–285.

(148) Westhof, E.; Auffinger, P. In *Encyclopedia of Analytical Chemistry*, 1st ed.; Meyers, R. A., Ed.; Wiley, 2000.

(149) Auffinger, P.; Westhof, E. Water and Ion Binding around r(UpA)12 and d(TpA)12 Oligomers - Comparison with RNA and DNA (CpG)12 Duplexes. *J. Mol. Biol.* **2001**, *305*, 1057–1072.

(150) Auffinger, P.; Westhof, E. Water and Ion Binding around RNA and DNA (C,G) Oligomers1. *J. Mol. Biol.* **2000**, *300*, 1113–1131.

(151) Davies, C. W. The Extent of Dissociation of Salts in Water. Part VIII. An Equation for the Mean Ionic Activity Coefficient of an Electrolyte in Water, and a Revision of the Dissociation Constants of Some Sulphates. *J. Chem. Soc. 1938*, 2093–2098.

(152) Fan, S.; Mason, P. E.; Chamorro, V. C.; Shanks, B.; Martinez-Seara, H.; Jungwirth, P. Charge Scaling Force Field for Biologically Relevant Ions Utilizing a Global Optimization Method. *ChemRXiv* **2025**, 10.26434/chemrxiv-2025-k63tq.
TOC Graphic

<figure><img src="image_10.png" /></figure>