{'fraction_non_alphanumeric': 0.029309559757471182, 'fraction_numerical': 0.023270304598940732, 'mean_word_length': 4.616906716920178, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Contact pairs of RNA with magnesium ions-electrostatics beyond the Poisson-Boltzmann equation

The electrostatic interaction of RNA with its aqueous environment is most relevant for defining macromolecular structure and biological function. The attractive interaction of phosphate groups in the RNA backbone with ions in the water environment leads to the accumulation of positively charged ions in the first few hydration layers around RNA. Electrostatics of this ion atmosphere and the resulting ion concentration profiles have been described by solutions of the nonlinear Poisson-Boltzmann equation and atomistic molecular dynamics (MD) simulations. Much less is known on contact pairs of RNA phosphate groups with ions at the RNA surface, regarding their abundance, molecular geometry, and role in defining RNA structure. Here, we present a combined theoretical and experimental study of interactions of a short RNA duplex with magnesium (Mg2+) ions. MD simulations covering a microsecond time range give detailed hydration geometries as well as electrostatics and spatial arrangements of phosphate-Mg2+ pairs, including both pairs in direct contact and separated by a single water layer. The theoretical predictions are benchmarked by linear infrared absorption and nonlinear two-dimensional infrared spectra of the asymmetric phosphate stretch vibration which probes both local interaction geometries and electric fields. Contact pairs of phosphate groups and Mg2+ ions are identified via their impact on the vibrational frequency position and line shape. A quantitative analysis of infrared spectra for a range of Mg2+-excess concentrations and comparison with fluorescence titration measurements shows that on average 20–30% of the Mg2+ ions interacting with the RNA duplex form contact pairs. The experimental and MD results are in good agreement. In contrast, calculations based on the nonlinear Poisson-Boltzmann equation fail in describing the ion arrangement, molecular electrostatic potential, and local electric field strengths correctly. Our results underline the importance of local electric field mapping and molecular-level simulations to correctly account for the electrostatics at the RNA-water interface.


INTRODUCTION
RNA structures display a high negative electric charge density, originating from the closely spaced phosphate groups in the sugar-phosphate backbone. Electrostatic interactions of the phosphates with mono-and divalent ions and water mol-ecules in the native aqueous environment have a direct impact on RNA structure and on the multifaceted function of RNAs in cells (1)(2)(3). At the molecular level, hydrated RNA represents a complex many-body ensemble of coupled charged and polar entities with a structure undergoing thermally activated fluctuations on a femto-to picosecond timescale (4). The variety of phosphate-ion-water arrangements and the strength and frequency spectrum of local electric fields in such geometries are understood only in part, with a particular lack of quantitative information. Basic open questions have been summarized in (3). Studies on DNA suggest electric field strengths on the order of 90 megavolt per centimeter (MV/cm) with fluctuation amplitudes of some 25 MV/cm (5) and partly ordered water structures were suggested to make up the hydration shell of doublestranded RNA (6,7).
The interaction of the charged RNA backbone with ions in the aqueous environment is frequently modeled with the help of Poisson-Boltzmann (PB) theory (8,9). Such work includes the application of size-modified PB approaches based on a lattice gas model (10,11) and generalizations to fluctuating counterion distributions (12). In most cases, the nonlinear Poisson-Boltzmann equation (PBE) has been solved for a given charge density rðrÞ of RNA and different (counter-)ion species i with charge z i and concentration c i . The ions are embedded in a water dielectric continuum, characterized by the static dielectric constant εðrÞz80. Inside the RNA structure, the dielectric constant is reduced to a much smaller value (ε ¼ 4 À 10) (13) and the interface between the inner and outer region is defined with the help of local geometries, as described by, e.g., atomic van der Waals radii (14). PB theories predict an accumulation of positively charged (counter-)ions around RNA, resulting in a pronounced gradient of ion concentration within 20-30 Å along the radial distance from the RNA surface (3). This condensation of (counter-)ions in a so-called diffuse ion atmosphere reduces the overall free energy of the ensemble and stabilizes hydrated RNA. However, the neglect of hydrogen bonding and electric polarizabilities throughout the ensemble together with the restriction to static molecular geometries limit the accuracy of PB pictures substantially.
Molecular dynamics (MD) simulations include the different constituents of hydrated RNA at the atomic and/ or molecular level and, thus, represent an intrinsically microscopic approach that allows to grasp dynamical properties (15)(16)(17). Electrostatic interactions are typically treated with the help of force fields and combined with molecular water models. The large numerical effort limits the size of the molecular simulation box while long-time trajectories into the microsecond time range are required for an equilibration of the ion atmosphere (18,19). MD simulations have confirmed the occurrence of ion condensation around RNA structures with the calculated radial profiles of ion concentration depending on the particular force field (16,17,19,20). Challenges arise in particular for doubly charged ions (e.g., Mg 2þ and Ca 2þ ) where, because of the neglect of electric polarizabilities and charge transfer among molecular constituents, fixed charge force fields are expected to overestimate the ion-phosphate interaction at the RNA/water interface. Continuous effort is devoted to the reparametrization of ion (and water) force field parameters, e.g., via additional charge-dipole interaction terms or a rescaling of atomic charges used in the simulations (21)(22)(23)(24). MD simulations allow for distinguishing different local interaction geometries of phosphate groups, ions, and water molecules and account for directed hydrogen bonds.
Testing the validity of the different, partly contradicting theoretical predictions (25) requires benchmark experimental data. Both x-ray small-angle scattering (26,27) and chemical methods such as fluorescence titration (28,29) and ion counting (30,31) have provided insight in the timeaveraged ion distribution around RNA. Recent results from ion counting experiments on short DNA and RNA duplexes suggest that PB simulations of RNA electrostatics fail to explain part of the experimental observations. Ultrafast fluctuation dynamics and specific phosphate-ion interaction geometries have been mapped experimentally by noninvasive vibrational probes, in particular phosphate stretching vibrations at the RNA-water interface and their linear and nonlinear two-dimensional infrared (2D-IR) spectra (7,32,33). An analysis of the time-dependent spectroscopic line shapes together with MD simulations has allowed for quantifying the strength and fluctuation timescale of electric fields acting on the polarizable phosphate groups (5). Characteristic spectral blue shifts of the asymmetric phosphate stretching vibrations have been established as a quantitative reporter of contact ion interactions (34)(35)(36).
Here, we combine MD simulations and PB calculations with stationary and ultrafast time-resolved infrared spectroscopy of RNA phosphate stretching vibrations for developing a broader quantitative understanding of local ion geometries and electric interactions close to the RNA-water interface. The interaction of short, hydrated RNA duplex structures with an excess of magnesium (Mg 2þ ) ions is studied in a wide range of ion concentrations. We establish the existence of phosphate-Mg 2þ contact pairs in short RNA duplexes that are spectroscopically distinguished from waterseparated geometries. This behavior is accounted for by the microscopic MD simulations but not by the PB treatment. The origin of such discrepancies is identified and underlines the necessity for a theoretical treatment at the molecular level.

MATERIALS AND METHODS
MD simulations of duplex RNA containing 23 adenine-uracil (AU) basepairs were performed with AMBER 18 (37) employing the ff99bsc0 force field with c OL3 corrections (38) and used the TIP4P-FB (39) water model and the 12-6-4 Lennard-Jones-type nonbonded model for Na þ and Mg 2þ (21,40) (TIP4P-FB/12-6-4 LJ). MD simulations of a single RNA double strand were performed using two different box sizes in a truncated octahedral simulation cell with total bounding box size of 171.4 Å (40 Å buffer region, 78,733 water molecules) and 98 Å (10 Å buffer region, 13,792 water molecules), respectively. The large simulation cell accommodates a cylinder around the RNA helical center of maximal radius r max z 55 Å which assures the presence of an extensive bulk region in the simulation. We find that the Mg 2þ ion atmosphere converges on a length scale rz20 Å from the helical axis (see Supporting materials and methods). Thus, the small simulation cell with a maximal cylinder radius r max z 32 Å contains a >10 Å bulk region. We have assured that the presented results are independent of the size of the simulation box. Multiple simulations were performed with a different number of added Mg 2þ ions R ¼ c(Mg 2þ )/ c(RNA) ¼ [0, 5,10,20,25,30], with R defining the ratio of the total Mg 2þ and RNA concentration c(Mg 2þ ) and c(RNA), and for different force-field-ion parameter combinations (see Table 1). For R ¼ 20, the total Mg 2þ ion concentration corresponds to $77 mM for the standard simulation box and $13 mM for the large simulation box, closely resembling the experimental conditions. Additional simulations were performed with a background concentration c 32,65,97,130], to explore the impact of singly charged ions on contact ion pair (CIP) formation with Mg 2þ (see Supporting materials and methods). Results from MD simulations are compared to PB simulations performed with the Adaptive Poisson-Boltzmann Solver (version 1.4.2) (41) (see Table 2).
Vibrational frequencies of the asymmetric phosphate stretching vibration n AS ðPO À 2 Þ were simulated on the mixed quantum mechanical/molecular mechanical (QM/MM) level of theory, following the procedure described in detail in (32). QM/MM simulations were performed with the NWchem program (43) and take into account the first solvation shell around PO 2 À groups, contact ions and the first solvation shell around contact ions in the QM region. A detailed description of the theoretical methods and a comparison of different water model and ion force field combinations is given in the Supporting materials and methods.
The short RNA duplexes (supplier Integrated DNA Technologies (IDT), Coralville, Iowa) studied in the experiments contain 23 alternating AU basepairs and form an A-helix. They were dissolved in ultrapure H 2 O (Rotipuran Ultra, Carl Roth, Karlsruhe, Germany). Samples with RNA duplex concentrations of c(RNA) ¼ 0.8-7 mM were used for linear infrared and nonlinear 2D-IR spectroscopy. To such solutions, MgCl 2 was added with different excess concentrations. In the aqueous environment, MgCl 2 dissociates in separately solvated Mg 2þ and Clions (44). In the following, the Mg 2þ concentration c(Mg 2þ ) is given in units of c(RNA) using the ratio R ¼ c(Mg 2þ )/ c(RNA). The IDT-supplied lyophilized duplex RNA samples contain residual acetate (CH 3 CO 2 À ) buffer with singly charged ions. We have characterized the rehydrated duplex RNA samples by inductively coupled plasma optical emission spectrometry (ICP-OES) and an Na þ -selective electrode (Supporting materials and methods). Per duplex RNA, there is a number of 50 5 5 Na þ and 5 5 2 K þ ions and, thus, an excess of $10 singly charged ions per duplex RNA (see Supporting materials and methods). Additionally, ICP-OES confirms negligible magnesium concentrations (R < 0.5) for the duplex RNA samples before the addition of Mg 2þ ions. For determining the fraction of Mg 2þ ions interacting with RNA duplexes, a fluorescence titration method was applied (28,29). Solutions of MgCl 2 in a buffer containing the chromophore 8-hydroxyquinoline sulfonic acid (8-HQS; c ¼ 0.5 mM), sodium cacodylate (c ¼ 10 mM), and NaCl (c ¼ 22 mM) were prepared in a range of c(Mg 2þ ) ¼ 0.1-150 mM. The measurement of fluorescence intensity from 8-HQS as a function of Mg 2þ content provides a reference to which the emission intensity from a sample additionally containing RNA (c(RNA) ¼ 0.26-1.1 mM) is compared. In the latter, the formation of RNA/Mg 2þ complexes reduces the number of 8-HQS/Mg 2þ complexes and, consequently, the emission intensity. A comparison with the reference measurement gives the concentration of Mg 2þ ions interacting with RNA (32).
In the different infrared experiments, liquid samples of a 12-mm (2D-IR) or 25-mm thickness (linear infrared absorption) were cast in a cell with two 1-mm thick CaF 2 or BaF 2 windows. Infrared absorption spectra were recorded with a Fourier transform infrared spectrometer (Bruker Vertex v80, spectral resolution 2 cm -1 ). The 2D-IR spectra were measured with a 3-pulse photon-echo method in which three femtosecond pulses interact with the sample and a fourth synchronized pulse serves for heterodyning the nonlinear signal diffracted from the sample. The first two pulses are separated in time by the coherence time t and the third pulse generates the nonlinear signal after a waiting or population time T relative to the second pulse. All 2D-IR spectra reported in the following were recorded at T ¼ 300 fs. The nonlinear signal is heterodyned with a fourth synchronized pulse, spectrally dispersed, and detected by an infrared detector array, defining the detection frequency n 3 (spectral resolution 2 cm -1 ). From signals measured for different coherence times t, frequency domain signals are derived by a Fourier transform, thus defining the excitation frequency n 1 . The femtosecond pulses were centered at 1230 cm -1 with a spectral full width at half maximum of 150 cm -1 and a duration of 80 fs. Further details of the method and the experimental setup have been reported in (45).

Simulation results
We first present results of atomistic MD simulations of duplex RNA in a water environment with R ¼ 0-30 up to the microsecond (ms) timescale. Fig. 1 a shows the instantaneous molecular RNA structure for R ¼ 20 after 0.6 ms simulation time. The RNA duplex forms an A-helix with a narrow and comparably deep major groove and a wide shallow minor groove. Populations of ions interacting with TABLE 1 Parameters and results of MD simulations: the different MD simulations are numbered in the first column, all simulations were performed with a truncated octahedral simulation cell, the simulation time is 600 ns and is preceded by a 400-ns equilibration phase the A-helix are predominantly located in the narrow major groove and, to much lesser extent, in the minor groove. Of the 20 Mg 2þ ions present per RNA duplex, $16 interact with the double-helix and $4 are solvated in the water environment. Among the interacting ions, 10 form CIPs with phosphate groups in the backbone via the PO 2 À subunit which contains the two nonbridging oxygen atoms. MD simulations were performed in presence of Mg 2þ and Clions, as well as 44 charge-neutralizing Na þ ions, whereas in experiments additionally $10 singly charged Na þ and K þ ions are contained in the samples because of the presence of CH 3 CO 2 À buffer. To closely resemble the experimental conditions, we have verified in extensive simulations for background concentrations of K þ ions in a range R K þ ¼ 0-130 that the ratio R is a key parameter for CIP formation between Mg 2þ ions and PO 2 À groups. In the experimentally relevant range R K þ ¼ 0-50 of K þ excess ions, we observe a minor influence of K þ ions on the number of CIPs (see Fig. S3). The simulations with a large simulation cell find a comparable number of CIPs (10 5 1) (see Table 1 and Supporting materials and methods), which suggests that within the experimentally investigated concentration range, formation of CIPs is largely independent of the absolute Mg 2þ and RNA concentration but mostly determined by the concentration ratio R, in agreement with experimental observations (see below, Fig. 4).
In a CIP, the Mg 2þ ion is surrounded by 5 H 2 O molecules and the vacant position of the octahedral first hydration shell is occupied by one of the nonbridging oxygen atoms of the PO 2 À group ( Fig. 1 a, middle inlay). Vice versa, the  (42) together with partial charges of the amber force field, see Supporting materials and methods for details. Data are shown for simulations J and G, simulation parameters are summarized in Tables 1 and 2.
Contact pairs of RNA with magnesium ions hydration of the PO 2 À group is distorted. Only one of the nonbridging oxygen atoms of the PO 2 À group is accessible for H 2 O molecules that form a prototypical tetrahedral hydration structure. The number of CIPs among the interacting Mg 2þ ions depends critically on the force field used in the MD simulations as will be discussed below. The remaining six interacting Mg 2þ ions are accommodated as solventshared ion pairs (SSIPs) with a water-mediated interaction between the Mg 2þ ion and the PO 2 À group. In an SSIP, the sixfold water coordination of the Mg 2þ ion is preserved in the octahedral hydration shell (Fig. 1 a, bottom inlay). The H 2 O molecules in the first solvation shell around the ion are simultaneously part of the tetrahedral solvation shell around the PO 2 À group. The narrow major-groove width of the RNA A-helix (see Supporting materials and methods) facilitates the stabilization of CIP and SSIP hydration geometries by interstrand interaction with more than one PO 2 À group. The Mg 2þ ions in the CIPs are immobile on the timescale of the simulations, their average residence time is longer than 1 ms. In contrast, SSIPs display a limited mobility with major-groove residence times extending into the tens of nanoseconds timescale. The major-groove residence time distribution of SSIPs is best represented by a power-law distribution (see Supporting materials and methods). The submicrosecond timescale mobility is evident from the delocalized character of the ion density averaged over a 0.6-ms period (gray contour in Fig. 1 a).
The relative stability of the CIP and SSIP hydration structures is determined by the complex electrostatic many-body scenario at the RNA phosphate-water interface. Fig. 1 b (upper panel) presents the potential of mean force (PMF) of Mg 2þ -PO 2 À interactions as a function of the coordinates d and a (upper inlay of Fig. 1 a), averaged over the different PO 2 À positions of the duplex RNA. For a P.Mg 2þ distance d < 4 Å , the 2D-PMF shows two distinct minima at polar angles a $ 5 and a $ 111 (a ¼ : O1¼P.Mg 2þ ) that reflect an approximately linear P¼O.Mg 2þ arrangement in CIPs with the O1 and O2 atoms of the PO 2 À group (inlay Fig. 1 a, O1¼P¼O2 angle $103 in PO 2 À group). For d > 4 Å , a pronounced SSIP minimum is identified at a $ 62 . In the range 4 < d < 6 Å , the 2D-PMF is particularly corrugated. Here, displacements among different SSIP structures occur due to thermally excited fluctuations (DE$1.5 k B T). The 2D-PMF is not symmetric in the angular coordinate a because the steric constraints of the RNA helical surface allow for increased accessibility of the O1¼P positions compared with O2¼P positions of the PO 2 À group. O1 oxygen atoms point toward solvent accessible volume, whereas the O2 oxygen atoms point into the groove direction.
The 2D-PMF derived from the MD simulations is compared with the electrostatic potential f evaluated on the PB level that determines the free energy of interaction between Mg 2þ and RNA (Fig. 1 b, bottom panel). The potential f shows no distinct minima that distinguish SSIP or CIP states. Instead, f exhibits a smooth surface without major corrugations. Details of f are determined by the molecular surface employed in the simulations. Specifically, the magnitude of f at the solute-solvent boundary is crucially affected by the molecular radii employed. Nevertheless, the asymmetry of O1 and O2 oxygen atoms is reflected in f due to an increased electrostatic potential within the groove, originating from the high PO 2 À charge density.
Slices of the 2D-PMF along the angular coordinate a for different d allow for a quantitative comparison of the stability of MD-derived CIP and SSIP geometries to the electrostatic potential of angular ion translation predicted by PB theory (Fig. 1 c). We find that a CIP involving the O2 oxygen atom of the PO 2 À group (a $ 111 ) is by $0.6 k B T more stable than a CIP formed with the O1 oxygen atom (a $ 5 ) despite the increased accessibility of the latter. The energetics of the CIP and SSIP structures are found to be within <1 k B T, resulting in a similar population in the MD simulations for R ¼ 20. In the CIP geometries, the angular PMF is strongly confined which prohibits thermal angular translation of ions around the PO 2 À group. As a result, MD simulations predict quasilinear CIP structures. A transformation among CIP and SSIP species is suppressed with an estimated barrier of 6.5 5 0.5 k B T (see Supporting materials and methods). The corrugated angular PMF from the MD simulations contrasts markedly with findings from PBE simulations in which almost free angular translation of ions is possible in the plane of the PO 2 À group. The $1 k B T/e free energy difference along a reflects the different electrostatic potentials inside and outside the major groove (a > 80 ). Fig. 2 compares the ion population along the angular coordinate a obtained from MD and PBE simulations. There are distinct maxima of ion population in MD simulations where the maxima at a $ 5 and a $ 111 arise from CIP structures and maxima at a $ 62 and a $ 140 correspond to SSIP structures. The asymmetry of the CIP maxima reflects the relative population of CIP at O1 and O2 positions of the PO 2 À group. The asymmetry of the SSIP population maxima is inverted compared to CIPs. SSIPs at the entrance of the major groove (a $ 120 ) can be stabilized via coordination by a second PO 2 À group on the opposite side of the groove (cf. Fig. 1 a, bottom inlay). The PBE calculations give a significantly different ion population. For a>80 , the larger ion population is induced by the differences in electrostatic potential inside and outside the major groove (a<80 ).
Phosphate stretching vibrations are sensitive quantitative reporters of their local hydration environment (5,33,34). Fig. 3 a compares the linear infrared absorption spectrum of the asymmetric PO 2 À stretching vibration n AS ðPO À 2 Þ of the RNA duplex simulated on the QM/MM level of theory to the experimental spectrum. The QM/MM simulations reveal a broad distribution of n AS ðPO À 2 Þ frequency positions due to different local hydration geometries of PO 2 À groups and specific ion interactions at different positions of the RNA duplex. We find excellent agreement in frequency position of n AS ðPO À 2 Þ covering a range from $1180 to 1280 cm -1 , whereas some deviation in lineshape is recognized. The density of states (Fig. 3 b) reveals the contributions of CIP and SSIP geometries of Mg 2þ ions with the PO 2 À group to the linear absorption spectrum. The predominant contribution of CIPs arises in the frequency range n AS ðPO À 2 Þ ¼ 1243-1282 cm -1 , which closely mimics experimental observations (see below). SSIPs in general appear red shifted compared with CIPs with n AS ðPO À 2 Þ covering a range $1183-1244 cm -1 . Very few CIP and SSIP show absorption at a much lower or higher frequency n AS ðPO À 2 Þ ¼ 1223 cm -1 and n AS ðPO À 2 Þ ¼ 1260-1273 cm -1 , due to particular geometric constraints in the groove of the RNA duplex.
The local electric fields imposed on the PO 2 À groups affect the frequency position of n AS ðPO À 2 Þ via a fieldinduced red-shift of the transition frequencies (33). In our analysis, we consider the electric field E projected on the bisector of the O1¼P¼O2 angle (Fig. 3 c). The MD simulations allow for isolating PO 2 À groups embedded in a neat water shell (red bars in Fig. 3 c) with electric fields up to some 100 MV/cm. The strong phosphate-water interaction is due to the particularly strong PO 2 À -H 2 O hydrogen bond and the rigid hydration structure around the PO 2 À group. The ordering of the first few water shells leads to an anisotropic water distribution in which the intermolecular interactions arising from different solvent molecules do not compensate, even upon the temporal averaging over configuration space. The electric field distribution including all PO 2 À groups extends up to 200 MV/cm (Fig. 3 c, blue bars). Such higher fields occur at sites where the Mg 2þ ions of CIPs and adjacent PO 2 À groups contribute to the total electric field.
The electric field in the PBE calculations contains the intermolecular contributions from dielectric and ionic boundary forces as well as reaction field forces (46,47) and reveals appreciable forces experienced by the interfacial PO 2 À group. The field distribution shows a maximum at E ¼ 20-25 MV/cm (Fig. 3 d), i.e., significantly smaller   $ 62 and a $ 140 ). Shown are results from MD simulations using different simulation box sizes (10 Å buffer size: 13,792 water molecules, Mg 2þ /RNA concentration ratio R ¼ 20; 40 Å buffer size: 78,733 water molecules, R ¼ 30). The number density was averaged over the radial interval of CIP and SSIP (d ¼ 3-9 Å ; cf. Fig. 1 b). The number density evaluated with the PBE method is modulated because of the higher electrostatic potential within the major groove (a > 80 ) compared with outside the groove region (a <80 ). Data are shown for simulations J, N, and G, simulation parameters are summarized in Tables 1 and 2. than predicted by the MD simulations. Contributions from dielectric boundary and reaction field forces are comparable, whereas ionic boundary contributions are negligible at the used ionic strength (0.405 M).

EXPERIMENTAL RESULTS
The asymmetric phosphate stretching vibrations n AS ðPO À 2 Þ of duplex (AU) 23 RNA serve as probes of local interaction geometries between RNA, ions, and water molecules, and of the related electric field strengths (33). This approach allows to benchmark the contradicting theoretical predictions of MD and PB simulations via the spectral positions and shapes of vibrational bands.
Linear infrared absorption spectra of n AS ðPO À 2 Þ for different ratios R ¼ c(Mg 2þ )/c(RNA) are presented in Fig. 4 a, showing the absorbance A(R) as a function of wavenumber. The spectra exhibit two strong components with maxima at 1220 and 1245 cm -1 and a weaker component around 1278 cm -1 . Upon addition of Mg 2þ ions, the absorption component at 1245 cm -1 shows a decrease, whereas absorption increases between 1260 and 1290 cm -1 . From such spectra, the differential absorbance DAðRÞ ¼ AðRÞ À cAðR ¼ 0Þ was derived (Fig. 4 b) in which c ¼ AðR; 1245Þ=AðR ¼ 0; 1245Þ is a correction factor accounting for the decrease in concentration of phosphate groups absorbing at 1245 cm -1 (AðR; 1245Þ, AðR ¼ 0; 1245Þ: absorbance at 1245 cm 1 for R and R ¼ 0). The difference spectra show a prominent increase in the frequency range 1260-1290 cm -1 with R.
In Fig. 4 c, the differential absorbance DA at 1278 cm -1 , normalized to the absorbance A at 1245 cm -1 , is plotted as a function of R (solid squares, right ordinate scale, c(RNA) ¼ 6.9 mM). There is a nearly linear increase up to R ¼ 20 and a pronounced saturation in the range from R 2D-IR experiments give insight in ultrafast frequency fluctuations and resulting spectral lineshapes of the asymmetric PO 2 À stretching vibration n AS ðPO À 2 Þ. Fig. 5 summarizes 2D-IR spectra of the RNA duplexes for concentration ratios R between 0 and 32. The absorptive 2D signal, i.e., the real part of the sum of the rephasing and nonrephasing signal is shown as a function of excitation frequency n 1 (ordinate) and detection frequency n 3 (abscissa). Yellow-red contours represent 2D-IR signals which arise from the v ¼ 0 / 1 transitions of the n AS ðPO À 2 Þ vibrations and are due to bleaching of the v ¼ 0 ground state and stimulated emission from the v ¼ 1 state. Blue contours characterize the transient v ¼ 1 / 2 absorption band from excited oscillators.
Three components at 1220, 1245, and 1278 cm -1 are resolved as individual peaks in the 2D-IR spectra and correspond to the three contributions to the linear infrared absorption spectra (Fig. 4). Crosspeaks are absent in the 2D-IR spectra, pointing to a negligible coupling between the three components. The elliptic lineshapes of the diagonal peaks are elongated along the diagonal, an indication of inhomogeneous broadening because of different local environments of the contributing phosphate groups. For increasing R, we observe, similar to the linear infrared  absorption spectra, an increase of the component at 1278 cm -1 . Cuts along a diagonal line parallel to n 1 ¼ n 3 (Fig. 5 e) show the increase of the 2D-amplitudes of the 1278 cm -1 component relative to the 1245 cm -1 component. The Mg 2þ -induced change in amplitude ratio of the 1278 relative to the 1245 cm -1 component of up to $0.1 is substantially higher than in the infrared absorption spectra.

DISCUSSION
The results of MD and PB simulations differ fundamentally in the description of electrostatics at the phosphate-water interface of RNA. The MD simulations predict the formation of distinct molecular species, i.e., hydrated CIP and SSIP geometries of Mg 2þ ions and charged PO 2 À groups, and a highly corrugated PMF of PO 2 À -Mg 2þ interactions. Moreover, there are strong local electric fields originating from water molecules and ions in different solvation arrangements. In contrast, the PB treatment suggests an elec-trostatic stabilization of RNA mediated via the interaction with ionic species fully solvated in water (48). The PBE electrostatic potential f in the vicinity of the PO 2 À group exhibits a smooth surface without formation of minima or maxima. Thus, a single ionic species is predicted to stabilize the RNA structure via long-range electrostatic interactions. The interfacial electric fields calculated in the PB approach (cf. Fig. 3) originate from the discontinuity of the dielectric constant at the RNA surface, i.e., the model-dependent dielectric boundary conditions, and the reaction field due to dielectric polarization in response to the high-charge density of RNA. The high solvent dielectric constant efficiently screens ionic field contributions. We now discuss the predictions of the MD and PB simulations by relating them to key results of the spectroscopic experiments.
Both the infrared absorption spectra and the 2D-IR spectra of the n AS ðPO À 2 Þ vibrations display three different components corresponding to distinct configurations of PO 2 À groups, ions, and water molecules at the RNA surface. This conclusion is supported by the absence of crosspeaks in 2D-IR spectra, i.e., the absence of vibrational couplings between the different sites. Our experimental and simulation work has allowed for the following assignments (32,49): the component with maximum at 1220 cm -1 is due to PO 2 À groups in prototypical tetrahedral hydration geometries that are fully exposed to water and are surrounded by up to six water molecules in the first solvation shell. The contribution around 1245 cm -1 is due to PO 2 À groups that are, due to the steric constraints of A-helical RNA, undercoordinated in the number of water molecules. Exemplary hydration geometries are chain-like arrangements of water molecules that bridge neighboring PO 2 À groups. The high polarizability of phosphate groups makes the frequency position of the n AS ðPO À 2 Þ vibration particularly sensitive to electric fields originating from the environment. Except for CIPs, this electrostatic interaction induces a redshift of n AS ðPO À 2 Þ compared with dehydrated PO 2 À groups without a water shell (50,51). High level simulations of RNA and DNA suggest a linear dependence of the n AS ðPO À 2 Þ frequency on the field strength E with Stark tuning rates of $0.5-0.6 cm -1 /(MV/cm) (34,52,53). Thus, the 25-cm -1 frequency separation of the 1220 and 1245 cm -1 components in the infrared spectra of the RNA duplex translates into a difference of the relevant electric fields on the order of 50 MV/cm. The MD simulations and previous work (32,53) reproduce this range of electric fields by explicitly including the water dipoles. In contrast, the PB treatment gives field strengths that are too small to explain the relative spectral positions (Fig. 3). This failure originates from the use of a dielectric continuum rather than a molecular approach.
The component around 1278 cm -1 is a hallmark of CIP geometries in which one of the PO 2 À oxygen atoms is integrated in the octahedral solvation shell of a Mg 2þ ion, thus replacing a water molecule (35,36,49). The frequency a b d c e FIGURE 5 2D-IR spectra in the frequency range of the n AS ðPO À 2 Þ vibration. The absorptive 2D signal is plotted as a function of excitation frequency n 1 and detection frequency n 3 . Yellow-red contours represent signals on v ¼ 0 / 1, blue contours on v ¼ 1 / 2 transitions. The signal is normalized to the respective maximum of the v ¼ 0 / 1 component, the signal change between neighboring contour lines is 7.5%. 2D spectra are shown for ratios R of Mg 2þ to RNA concentration of (a) R ¼ 0, (b) R ¼ up-shift of the n AS ðPO À 2 Þ bands of CIPs compared to the other two species is due the strong interaction of one of the PO 2 À oxygens with a Mg 2þ ion at a short distance < 2:2 Å . Because of this asymmetric distortion of the PO 2 À binding potential, elongations of the n AS ðPO À 2 Þ vibration probe the repulsive part of the intermolecular interaction potential, leading to a blue-shift of the vibration. This mechanism and the underlying molecular geometry are grasped by the QM/MM level of theory but not by the PB approach. It should be noted that the n AS ðPO À 2 Þ frequencies of SSIPs cover a broad range from 1200 to 1260 cm -1 , which overlaps with the two main components in the spectra (32) and, thus, prevents a spectroscopic distinction of SSIPs.
Our experimental results give quantitative insight in the abundance of CIPs. As shown in Fig. 4 c and Fig. S9, the total number of Mg 2þ ions interacting with an RNA duplex saturates at $15 for high ion concentrations, corresponding to roughly one-third of the number of phosphate groups. A subset of interacting Mg 2þ ions is accommodated as CIPs the number of which can be estimated from the spectroscopic data. The differential absorption DA 1278 shown in Fig. 4 c saturates at a value of 5.4% of the peak absorption at 1245 cm -1 . Assuming a similar infrared absorption strength of the different species, one estimates a lower limit of CIPs per RNA duplex of 2-3 in the range of R ¼ 20 to 30 (some 5% of all 46 phosphate groups). On the other hand, the analysis of the relative amplitude ratio at 1245 and 1278 cm -1 in the 2D-IR spectra for varying R (cf. Supporting materials and methods) suggests the formation of 55 1 CIPs per RNA duplex as an upper limit.
The number of CIPs per RNA duplex derived from the experimental data is about a factor of two smaller than the predictions of the MD simulations, giving a number of CIPs per RNA duplex of 10 for R ¼ 20. As discussed in more detail in the Supporting materials and methods, the calculated number of CIPs depends critically on the force field employed in the MD simulations. The present force fields neglect the impact of electric polarizability on the Coulomb interaction strength between phosphate groups and Mg 2þ ions. As a result, the attractive electrostatic forces and, thus, the CIP abundance may be overestimated. To test this hypothesis, we performed MD simulations with a force field in which the electric charge of the Mg 2þ ions was rescaled (z Mg 2þ ¼ 1.5), i.e., electric polarizability is accounted for in a mean field approach (22,24). This model shows a substantially reduced tendency to form CIP (see Supporting materials and methods), which indicates the importance of electronic polarization for the equilibrium between CIP and SSIP at the phosphate-water interface of RNA. Recent progress in the development of polarizable force fields is promising (54-56) but challenges persist because of extensive computational demands arising from large system size and long propagation times on the microsecond timescale.
We now address the impact of CIPs and SSIPs on the helix structure of the RNA duplex. The MD simulations sug-gest that the narrow major-groove geometry (groove width $10 Å ) mediates the stabilization of CIPs and SSIPs via the outersphere coordination by more than one PO 2 À group. Similar bridging hydration structures have been reported in the x-ray structures of mixed A-form DNA/RNA duplexes in which hexahydrated Mg 2þ ions in SSIP structures bridge the outer mouth of the duplexes (57). The intergroove electrostatic stabilization of helical RNA via the interaction of Mg 2þ ions with two different PO 2 À groups is expected to be self-enforcing and reveals new ion correlations. The accumulation of ion population in the groove region leads to a contraction of the groove (see Supporting materials and methods) that in turn increases the strength of the PO 2 À -ion interaction. The corrugated nature of the 2D-PMF in the range 4 < d < 8 Å reflects the coordination of SSIP across the narrow major groove. Such coordination suppresses a free diffusive movement of ions but imposes distinct transitions between solvation states and structures, activated by excitations via thermal fluctuation ($1.5 k B T). Major-groove residence times of SSIPs are found to extend well into the tens of nanoseconds timescale. The mutual interactions of both RNA strands with hydrated Mg 2þ ions over the narrow major groove (cf. inlay in Fig. 1  a) impose particularly complex dynamics and long residence times of SSIP in the major-groove region. Along the angular coordinate a, distinct ion population maxima are identified because of CIP and SSIP structures with Mg 2þ ions.
In comparison, PB simulations predict a qualitatively different picture of the ion population in the groove region. Here, the ion population is predicted to decrease in the opening region of the groove (a>90 ) as the electrostatic potential is reduced toward the solvent accessible bulk region. Structure-ion correlations are absent and, because of the neglect of structural relaxation, the impact of negative charge density from adjacent and intergroove PO 2 À groups on the electrostatic potential is additive. Consequently, macromolecular SSIP stabilization via correlated ion-backbone motion is absent. Details of the electrostatic potential f, like the ion-phosphate interaction strength at close contact, depend on the employed molecular surface. The magnitude of ion accumulation at the solute-solvent boundary is, thus, crucially affected by the molecular radii employed in simulations. The identified limitations of the PB approach are expected to be relevant for extensions of the PB method, like the size-modified PB approach (10,11) where deficiencies in the description of Mg 2þ ions have been identified (25).

CONCLUSIONS
We have identified clear spectroscopic signatures of CIPs of Mg 2þ and PO 2 À groups for short double-helical RNA via the high sensitivity of noninvasive infrared spectroscopy of phosphate backbone reporter modes. CIP formation is quantified via the systematic increase of infrared absorption around 1278 cm -1 for increasing Mg 2þ concentration which allows to derive a number of 2-5 CIPs per RNA duplex for a Mg 2þ /RNA concentration ratio of R ¼ 20-30. MD simulations provide a most detailed view of hydration geometries and show the coexistence of CIP and SSIP hydration structures of Mg 2þ ions and PO 2 À groups. Electrostatic stabilization of helical RNA via the intergroove coordination of Mg 2þ ions by two PO 2 À groups was identified to lead to structural correlations of the sugar-phosphate backbone and ions which are absent in PB approaches. Our results reveal the limited reliability of PB approaches in the description of interfacial properties. Deficiencies were identified for local ion arrangements, the molecular electrostatic potential and local electric field strengths determining peak positions in infrared absorption experiments. Our work established the relevant interactions and solvation geometries for modeling of the ionic properties at the biological relevant phosphate-water interface. Our findings underscore the importance of balanced molecular level descriptions of ionic interactions, the desolvation of ions (58) and phosphate groups, coupled to structural relaxation of the RNA helix.

DATA AVAILABILITY
The experimental data sets and calculation results generated and/or analyzed in the current study are archived at the Max-Born-Institute and are available from the corresponding author upon reasonable request.

AUTHOR CONTRIBUTIONS
T.E. and B.P.F. designed the research. J.S. and A.K. performed the experiments and analyzed the data together with T.E. B.P.F. carried out all theoretical calculations and simulations. All authors discussed the data and wrote the article.