Kinetic pathways of water exchange in the first hydration shell of magnesium.

Water exchange between the coordination shells of metal cations in aqueous solutions is fundamental in understanding their role in biochemical processes. Despite the importance, the microscopic mechanism of water exchange in the first hydration shell of Mg2+ has not been resolved since the exchange dynamics is out of reach for conventional all-atom simulations. To overcome this challenge, transition path sampling is applied to resolve the kinetic pathways, to characterize the reaction mechanism and to provide an accurate estimate of the exchange rate. The results reveal that water exchange involves the concerted motion of two exchanging water molecules and the collective rearrangement of all water molecules in the first hydration shell. Using a recently developed atomistic model for Mg2+, water molecules remain in the first hydration shell for about 40 ms, a time considerably longer compared to the 0.1 ms predicted by transition state theory based on the coordinates of a single water molecule. The discrepancy between these timescales arises from the neglected degrees of freedom of the second exchanging water molecule that plays a decisive role in the reaction mechanism. The approach presented here contributes molecular insights into the dynamics of water around metal cations and provides the basis for developing accurate atomistic models or for understanding complex biological processes involving metal cations.


INTRODUCTION
Water exchange between the first and second hydration shell around metal ions is fundamental for a large variety of processes ranging from simple chemical reactions in aqueous solutions to the structure and function of biomolecules. [1][2][3][4] In particular, the mechanism governs any type of reaction involving the replacement of strongly bound hydration water from the first hydration shell, which is essential, for instance, in metalloenzyme catalyzed reactions, regulatory biochemical processes, and the transport of metal ions across cell membranes. 3,[5][6][7] The molecular nature of water gives rise to an intriguing interplay of hydrogen bonding, packing, and orientational effects leading to complex dynamics in which the concerted motion of solvent molecules becomes important. 2,[8][9][10] It is, therefore, not surprising that water relaxation times in the first hydration shell surrounding ions can be significantly altered compared to bulk. 11 A particularly striking example is the mean lifetime of water molecules in the first coordination shell of metal cations. For different cations, the lifetimes span more than 18 orders of magnitude, ranging from a few picoseconds for Cs + to hundreds of years for Ir 3+ . 7 Given this tremendously broad timescale, the question arises how mechanistic insights can be obtained at the molecular level.
Experimental techniques such as dielectric relaxation, femtosecond mid-infrared, far-infrared (terahertz), and x-ray adsorption spectroscopy provide insight into the structural and dynamical properties of water around ions, 2,4,8,12 while water exchange rates are typically determined by nuclear magnetic resonance. 7,13,14 Since the structural changes at the microscopic level during an exchange process are not directly accessible from these experiments, only indirect classifications of the reaction mechanism according to stoichiometry or activation volume exist. 7,15 Here, simulations can contribute important insights by characterizing the solvent behavior and by providing a unique atomistic description of the dynamics. In particular, ab initio quantum mechanical calculations can provide unbiased insights into the structure of the first hydration shell and the physical origin of water-ion interactions. 12,16 However, due to their high computational costs, these methods are limited to ion-water clusters and to exchange processes on the picosecond timescale. Even though classical all-atom simulations are less expensive, up to date, water exchange can only be simulated for weakly hydrated ions such as monovalent alkali metal ions for which the exchange process is sufficiently fast. 17,18 For strongly hydrated ions, water exchange becomes so rare that the process cannot be simulated with sufficient statistics or, in many cases, not at all.
The present study focuses on water exchange in the first hydration shell of Mg 2+ as an intriguing example for exchange dynamics on the micro-to millisecond timescale. 7,13,14,19,20 Mg 2+ ions are ubiquitous in nature and play a decisive role in a large variety of biological processes including enzyme catalysis or ribozyme chemistry and have, therefore, received considerable scientific attention. 19,[21][22][23][24] In particular, transition state theory (TST) has been used to compare different or develop new force field parameters for Mg 2+19,25 or to derive a universal relationship between the water exchange time and charge density of different metal cations. 20 The first hydration shell around Mg 2+ ions contains six water molecules in an octahedral arrangement. 7,16,26 However, the most elementary process of exchanging these strongly bound hydration waters is so rare that it remains out of reach for conventional molecular dynamics (MD) simulations. Consequently, the microscopic mechanism of water exchange has not been resolved so far, and the kinetic pathways of exchange remain elusive.
To fill this gap, transition path sampling is applied as a particularly powerful sampling strategy to provide unbiased microscopic insight into the dynamics and to characterize the mechanism of water exchange. The approach reveals that water exchange involves the concerted motion of two water molecules. The molecular void provoked by the leaving water molecule is immediately filled by an entering water molecule and collective rearrangement of the first hydration shell. This detailed mechanistic insight contradicts the simplified picture in which water exchange is described by the distance between the cation and one water molecule alone. As a consequence, the rate predicted by transition state theory (TST) with the distance of one water molecule as a reaction coordinate overestimates the rate by more than two orders of magnitude.
In order to make progress in these complex many body systems, the mechanism of water exchange is characterized by collecting a large number of exchange trajectories using path sampling. Subsequently, an improved reaction coordinate is defined, which includes the coordinates of the two exchanging water molecules. Finally, the rate constant of water exchange is calculated and a detailed comparison of the accuracy of different methods is given.

Atomistic model and simulation setup
The system consists of a single Mg 2+ ion a cubic simulation box (L = 25 Å) filled with 506 water molecules. For Mg 2+ , recently optimized force field parameters were used 17 in combination with the TIP3P water model. 27 It should be noted that other water models are superior in capturing the properties of pure water. 28 Yet, out of the large number of rigid, non-polarizable water models, the most popular ones are the simple point charge (SPC) 29 or the transferable intermolecular potential TIP3P. 27 In addition, the amber force fields for proteins and nucleic acids have been optimized in combination with TIP3P water explaining its widespread use in biomolecular simulations and the selection of TIP3P in our current work. Regarding the choice of the force field parameters for the Mg 2+ ion, our choice was motivated by the fact that the optimized parameters reproduce experimental solution properties including hydration free energies and activity coefficients. 17 All simulations were performed using GROMACS 30 with periodic boundary conditions. Particle mesh Ewald summation was used with tin foil boundary conditions and a Fourier spacing of 0.12 nm and a grid interpolation up to order four to handle long-range electrostatic forces and the electro-neutrality condition. Close Coulomb real space interactions were cut off at 1.2 nm and Lennard-Jones (LJ) interactions after 1.2 nm, respectively. Long-range dispersion corrections for energy and pressure were applied to account for errors stemming from truncated LJ interactions. An implicit background charge was used to neutralize the system via the tin foil boundary conditions. Since the system is homogeneous, a uniform background charge is an appropriate model for a converged distribution of counterions. 31 In addition, the uniform background charge allows us to avoid possible artifacts due to ion pairing and sampling problems in the counterion distribution. The initial energy minimization was performed with the steepest descent algorithm. For each simulation, a NVT equilibration was performed for 1 ns, controlling the temperature at 300 K with the Berendsen thermostat. 29 All production runs and the transition path sampling were carried out in the NVT ensemble at a temperature of 300 K using the velocity rescaling thermostat with a stochastic term 32 and a time step of 2 fs. Note that the velocity rescaling thermostat generates the canonical ensemble and ensures detailed balance in the canonical ensemble of transition paths. 33

Rate constant at high temperatures
Water exchange rates at high temperatures were determined from temperature replica exchange molecular dynamics (TREMD) simulations in combination with a kinetic rate model. 34 64 replicas of the system were used in a temperature range T = 300-895 K with a temperature series generated from the temperature generator. 35 The simulations were performed for 90 ns, and exchanges between replicas were attempted every 2 ps leading to an overall exchange probability of 0.3. The rate constants of water exchange were calculated from 34 where N is the total number of transitions between the two stable states and t b/ub is the total time a water molecule spends inside or outside of the first hydration shell. Further details on the transition counting can be found in the supplementary material. To validate the approach further, a straightforward 200 ns MD simulation was performed at T = 694 K. The standard errors were calculated from the transition counts. 34 The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp

Free energy profiles
The one-dimensional free energy profile as a function of the ion-water distance was taken from our previous work. 17 The twodimensional free energy profiles as a function of the distance between Mg 2+ and the two exchanging water molecules were calculated from umbrella sampling using PLUMED. 36 A force constant k b = 100 000 kJ/(mol nm 2 ), and a window spacing of 0.01 nm was used. Each umbrella simulation was performed for 5 ns discarding 500 ps for equilibration. The free energy profiles were calculated using the weighted histogram analysis method (WHAM). 37 Without further restraints, the system has four stable states. In order to consider only transitions between the two relevant states (corresponding to the state in which one of the two exchanging water molecules is bound), an additional biasing potential on the hydration number was applied (see the supplementary material for further details). In experiments, the measured rate constant k corresponds to the exchange of one water molecule in the first hydration shell by any other water molecule. To facilitate direct comparison, the rate constant is calculated from Here, kA→B is the rate in the restrained system, N is the number of water molecules in the simulation box, and 6 is the coordination number of Mg 2+ . Hence, (N − 6) corresponds to the number of choices to select the second exchanging water molecule. Herby, it is assumed that subsequent transitions are uncorrelated (as confirmed by the clear separation of timescales obtained from the simulations as discussed later).

Transition path sampling
Transition path sampling 33,38 was used to harvest an ensemble of rare trajectories that connect the two stable states. Starting from an initial reactive pathway generated at high temperature, new trial trajectories were created by randomly selecting a time slice, randomizing the velocities and integrating the equations of motion forward and backward in time. Using standard two way shooting moves with a fixed length of 1 ps, trial moves are accepted if they connected the two stable states and rejected otherwise. The trajectories generated in the transition path sampling are true dynamic trajectories free of any bias.

Committor analysis and transition states
The commitment probability p A was calculated from the fraction of trajectories initiated with randomized velocities that reach the product state A. As starting points, 160 conformations were selected from the one-dimensional umbrella sampling at the barrier top (r 1 = 0.29 nm). Similarly, 600 conformations obtained from the two-dimensional umbrella sampling were selected along the separatrix (r 1 = r 2 ) with weights according to the equilibrium distribution. For each conformation, 100 trajectories were initiated with velocities drawn from a Maxwell-Boltzmann velocity distribution and run forward and backward for 1 ps. From these 60 000 fleeting trajectories, the probability p A to relax back to state A and its probability distribution p(p A ) were calculated. In addition, 4800 conformations along more than 2500 independent pathways obtained from transition path sampling were used. For each conformation, 100 trajectories were initiated with velocities drawn from a Maxwell-Boltzmann velocity distribution and run forward and backward for 1 ps. A conformation was identified as a transition state if half of the trajectories relaxed into each stable state. Note that the fleeting trajectories initiated from the restrained ensemble generated by umbrella sampling were also used to calculate the correction factor in modified TST (see below). The transition state ensemble was characterized further by calculating the symmetry of the two exchanging and the five spectator water molecules in the first hydration shell. In particular, the mirror symmetry is quantified by calculating the root mean square deviation of the transition state and its mirror image, where ri are the positions of the oxygen atoms of the n = 7 water molecules closest to Mg 2+ in the transition state and rj are the corresponding closest image positions obtained by reflection at the mirror plane perpendicular to the plane containing the two exchanging water molecules and Mg 2+ .

Projection of the free energy landscape onto the reaction coordinate
In order to quantify the progress along a reaction pathway, the coordinates of the two exchanging water molecules are combined into a single reaction coordinate λ = arctan ((r 2 − r 0 ) (r 1 − r 0 )). The projection of the two-dimensional free energy landscape onto the one-dimensional reaction coordinate λ is obtained by integrating over all coordinates as follows: where β = 1/(k B T), F(r 1 , r 2 ) is the two-dimensional free energy landscape, and δ(x) is the Dirac delta function used to project the free energy onto λ. R corresponds to the radius of a sphere containing the same number of water molecules as in the cubic simulation box.

Rate calculation from conventional and modified transition state theory
In conventional transition state theory (TST), the rate constant follows from 39,40 where λ * is the position of the barrier top, θ(λ) is the Heaviside step function, and ⟨⋯⟩c denotes the average over the restrained ensemble of trajectories initiated from an equilibrium ensemble of phase points on the dividing surface. Note that the free energy rather than the potential of mean force enters into the equation. If potentials of mean force are used, additional modifications are required TST gives an accurate rate only if the reaction coordinate is exactly known. Otherwise, TST gives a rough upper limit for the rate constant and further corrections are required. Therefore, three different methods were applied to calculate the accurate correction factor that accounts for recrossings including non-reactive trajectories. First, the rate constant is calculated from the modified reactive flux (RF) method, 41,42 where cA denotes the probability of being in state A. peq is the normalized equilibrium distribution defined further below. Second, the rate constant is calculated from the positive flux (PF) method as follows: 43 Third, the rate constant is calculated from the positive weight (PW) method as follows: 44 where the sum is over all intersections of each path with the dividing surface andλi is the velocity normal to the dividing surface at intersection point i. θ TP is one if forward and backward trajectory form a transition path and zero otherwise. In Eqs. (6)-(8), the normalized equilibrium distribution peq is given by Quantifying the limits of transition state theory TST is the most popular theory to calculate reaction rates. However, in complex systems, TST could fail due to the violation of the non-recrossing hypothesis, which forms the cornerstone of the theory. By applying reactive flux methods as detailed earlier, we can quantify whether the non-crossing hypothesis is satisfied or not. This is achieved by calculating the transmission coefficient κ, which is defined as the ratio between the true rate constant (kA→B) and that obtained from TST (k TST A→B ), In particular, κ is unity if there are no recrossings and the rate obtained from TST is exact. Otherwise, κ is smaller than unity, and the corrections due to recrossings can become important.
In addition, κ can be used to test the quality of different reaction coordinates [namely, r 1 and λ = arctan ((r 2 − r 0 ) (r 1 − r 0 ))]. Rewriting Eq. (7) yields with Since the rate constant kA→B does not change under a coordinate transformation from λ to r 1 , a kinetic quasi-universality can be written as follows: Therefore, combining a dynamic prefactor α(λ) = κ(λ) ⟨λ 0 θ(λ 0 )⟩c with the free energy profile yields the rate constant for arbitrary reaction coordinates. The value of κ calculated for different reaction coordinates provides direct insight into the quality of the coordinate. In particular, the three different scenarios are differentiated in the following: (i) κ(λ) = 1 corresponds to an ideal reaction coordinate.
Here, the reactant and product basins are separated such that trajectories never recross the dividing surface at the barrier top. Therefore, the distribution p(p A ) of the probability to relax back to state A has a sharp peak at pA ≈ 1/2. In this scenario, TST is valid without further corrections. (ii) κ ≈ 0.1 corresponds to a good reaction coordinate.
Here, p(p A ) shows a peak at pA ≈ 1/2 but is broader compared to the ideal distribution. In this scenario, few recrossings occur at the dividing surface and TST yields a reasonable upper estimate of the rate constant (within one order of magnitude). Using the reactive flux methods allows one to take the recrossings into account and to provide a more accurate estimate of the rate constant. (iii) κ < 0.01 corresponds to a poor choice of the reaction coordinate. Here, p(p A ) has two peaks, one at pA = 0 and pA = 1. In this scenario, TST overrates the rate significantly (more than two orders of magnitude). Therefore, additional degrees of freedom should be included into the reaction coordinate. Still, it is possible to calculate the rate constant based on reactive flux methods. However, the relative error of κ from reactive flux scales with 45 where Ntr is the number of trajectories initiated from the constrained ensemble. Consequently, an enormously large number of trajectories would be required. In the following, kA→B is calculated from the three different reactive flux methods using Eqs. (6)-(8) and the reaction coordinate λ = arctan ((r 2 − r 0 ) (r 1 − r 0 )). Hereby, the equilibrium distribution peq is obtained from umbrella sampling, while the dynamic prefactors are obtained using 60 000 fleeting trajectories. The transmission coefficient κ(λ) was calculated directly using Eq. (12). Finally, the transmission coefficient κ(r 1 ) was calculated from Eq. (13).

ARTICLE scitation.org/journal/jcp
All values and the errors obtained from block averaging are listed in Table S1 in the supplementary material.

RESULTS AND DISCUSSION
The aim of this work is to contribute microscopic insights into the kinetic pathways of water exchange in the first hydration shell of Mg 2+ and to provide an accurate estimate of the exchange rate at room temperature. Initially, the free energy landscapes of water exchange are presented. Subsequently, transition path sampling is applied to gain insight into the reaction mechanism and to identify the solvent molecules, which play a decisive role in the exchange reaction. Based on the reaction mechanism obtained from transition path sampling, a reaction coordinate is defined, which incorporates the concerted motion of the two exchanging water molecules. Finally, the rate constant of water exchange is calculated, and a comparison of the accuracy of different theoretical methods is given.

Free energy landscape of water exchange
In water exchange reactions, one of the six water molecules in the first hydration shell of Mg 2+ is replaced by a water molecule from the second hydration shell. In order to understand the underlying mechanism, it is essential to first identify the solvent molecules, which are relevant for the exchange reaction.
In the simplest case, only the coordinates of the leaving water molecule are considered, while all other degrees of freedom of the many body system are integrated out. The resulting onedimensional free energy profile as a function of the distance r 1 between Mg 2+ and the leaving water molecule is shown in Fig. 1(a). The free energy profile has two stable states (denoted as state A and B), which are separated by a high free energy barrier. In order to provide insight whether r 1 yields an adequate description of water exchange, a committor analysis is applied by shooting off trajectories from conformations selected at the barrier top (r 1 = r * 1 ). If r 1 was an appropriate choice for the reaction coordinate, about half of the trajectories initiated from a certain conformation should relax back to each stable state. As a consequence, the distribution p(p A ) of the probability to relax back to state A is expected to have a sharp peak at pA ≈ 1/2 (where the width of the distribution depends on the number of trajectories). Comparison of the resulting distribution to the expected ideal distribution clearly reveals that this is not the case [ Fig. 1(b)]. The vast majority of conformations are committed either to stable state A or B while only a vanishingly small fraction corresponds to pA ≈ 1/2. Hence, p(p A ) has two peaks, one at pA = 0 and the other at pA = 1 as the basins of attraction are not well separated along r 1 [ Fig. 1(c)]. Hence, the coordinates of the leaving water molecule alone serve as a very poor choice for a reaction coordinate and are, therefore, insufficient to capture the mechanism of water exchange.
In order to truly capture the mechanism, additional solvent degrees of freedom are required. In particular, water exchange involves the concerted motion of two water molecules in which the leaving water molecules is immediately replaced by an entering water molecule. The free energy landscape as a function of the distances of leaving (r 1 ) and entering (r 2 ) water molecules is shown in Fig. 1(c). From the two-dimensional representation, the failure of  has to involve at least the coordinates of the two exchanging water molecules.

Kinetic pathways from transition path sampling
To gain a clear mechanistic interpretation of water exchange, transition path sampling is applied to sample a large number of independent reactive pathways. Four representative transition paths are shown in Fig. 2(a). The distribution of the exchange angles [ Fig. 2(b)] along reactive pathways indicates that two alternative exchange pathways exist: In the indirect exchange mechanism, the entering and leaving water molecules occupy different positions on the water octahedron [ Fig. 2(e)]. During activation, one water molecule from the second hydration shell enters the molecular void leading to the concerted motion of another water molecule out of the first hydration shell. The distances of the leaving and entering water molecules in the transition states are elongated, while the distances of the remaining spectator water molecules remain relatively unchanged (Table I). At the same time, the spectator water molecules rearrange such that the transition state has mirror symmetry. Entering and leaving water molecules are symmetrically equivalent and indistinguishable. In the direct exchange mechanism, the exchange takes place via the attack of the incoming water molecule onto the edge of the water octahedron. Here, entering and leaving water molecules occupy the same positions on the water octahedron [ Fig. 2(f)]. In equilibrium, the indirect exchange mechanism is observed much more frequently (92%) compared to the direct mechanism (8%). Note that, in a few cases (<1% of all trajectories), three water molecules are involved in the exchange typically by a short-time rebound of the third water molecule.
The sum of all Mg 2+ -oxygen distances ΔΣ during the activation process can be used to classify the exchange mechanism. 46,47 Based on ΔΣ shown in Fig. 2(c), both exchange mechanisms correspond to interchange dissociative (I d ) processes in agreement with experimental results. 7,14 Figure 2(d) shows the distribution of transition times. For both reaction pathways, the exchanging water molecules spend less than 0.4 ps in transition. Therefore, a clear separation of timescales between the duration of a transition event (ps) and the dwell time TABLE I. Properties of the transition state ensemble: Mg 2+ -oxygen distances of the exchanging water molecules r 1,2 , Mg 2+ -oxygen distances of non-exchanging water molecules in first hydration shell rs, change of the sum of all Mg 2+ -oxygen distances ΔΣ, angle αex between the exchanging water molecules and Mg 2+ , average transition time τ, and probability p for indirect and direct exchange mechanism in equilibrium. Standard errors are indicated using a sample size of 327 or 137 transition states for the indirect and direct exchange mechanisms, respectively. Mechanism

Reaction coordinate
In the following, the mechanistic insight obtained from transition path sampling is used to provide a simple scalar reaction coordinate. Initially, the commitment probability as a function of the coordinates of the two exchanging water molecules is calculated [ Fig. 3(a)]. Transition states along reactive pathways are distinguished as conformations in which the exchanging water molecules have elongated but equal distances from the Mg 2+ ion and are, hence, distributed along the diagonal section with r 1 = r 2 [ Fig. 3(a)]. A simple scalar reaction coordinate λ that allows us to monitor the progress during the reaction is given by an angular function that combines the coordinates of the two exchanging waters, where r 0 = 0.2 nm. Note that based on likelihood maximization, 48 the angular function provides a better choice than any linear combination of r 1 and r 2 (data not shown). Subsequently, the onedimensional free energy profile along λ is obtained by projecting the two-dimensional free energy landscape onto λ according to Eq. (4).
In the resulting free energy profile along λ [ Fig. 3(b)], reactant and product states are identical, and the free energy difference is zero as expected since forward and backward reaction must be symmetrical. The putative reaction coordinate takes into account the coordinates of the two exchanging water molecules and allows us to provide a simple mechanistic interpretation of water exchange. In the concerted motion of the two exchanging water molecules, the molecular void provoked by the leaving water molecule is immediately filled by the entering water molecule.

Water exchange rate
Finally, the question how the rate of water exchange can be calculated accurately is addressed. To provide insight into the accuracy of different methods, the rate is calculated based on conventional and modified TST and compared to the results from TREMD (see the section titled Methods). Figure 4 shows the rate constant as a function of temperature. At high temperatures, TST based on r 1 as reaction coordinate  Table II. At room temperature, water exchange becomes so rare that the acceleration by TREMD is insufficient to provide a reliable rate estimate since no exchange events are observed on the 100 ns timescale. In order to determine the rate, the concept of TST is combined with the information obtained from dynamical trajectories using the reactive flux (RF), positive flux (PF), and positive weight approach (PW) to correctly account for recrossings (see the section titled Methods). In addition, the improved reaction coordinate λ based on the coordinates of the two exchanging water molecules is used [Eq. (15)]. To validate the approach, the rate constant is calculated at high temperatures (T = 649 K) and compared to the results from straightforward MD simulations and TREMD. Using this improved approach, the results from all methods (RF, PF, and PW) agree, within errors, with the benchmark results from MD and TREMD (Table II and Fig. S2 in the supplementary material).
At 300 K, the rate constant calculated from the PF approach is 24.0 ± 8.8 s −1 . The rate is more than two orders of magnitude lower compared to the value predicted by TST based on r 1 ( Fig. 4; Table II). This corresponds to a very low transmission coefficient of κ(r 1 ) = 0.003, which reflects the reaction coordinate problem discussed above. TST based on λ provides significant improvement (κ(λ) = 0.12). Still, moving the system reversibly from the state with λ = λ * to the true transition state (κ = 1) requires an amount of work of about 2 k B T. This additional barrier is connected to the rearrangement of the waters in the first hydration shell. Therefore, the coordinate λ fails to provide a complete description of the transition state since it does not include the collective rearrangement of all water molecules in the first hydration shell.  The quality of the reaction coordinate can be assessed further by the committor analysis shown in Fig. 5: Conformations with λ = λ * lead to a broad distribution. However, if the rearrangement of the water molecules is considered by taking the mirror symmetry of the transition state into account, the committor distribution is peaked at pA ≈ 1/2. In summary, providing an accurate estimate of the rate constant of water exchange in the first hydration shell of Mg 2+ from reaction rate theory is challenging due to the concerted motion of the exchanging water molecules and the additional collective rearrangements in the first hydration shell. Neglecting these solvent degrees of freedom leads to an overestimate of the true rate by more than two orders of magnitude.
On the other hand, a reliable estimate of the rate can be obtained by incorporating the coordinates of the two exchanging water molecules into the reaction coordinate and accounting for further collective rearrangement via the transmission coefficient obtained from dynamic trajectories used in reactive flux approaches.
The results clearly reveal that the exchange rate using recently developed force field parameters for Mg 2+ is much smaller compared to experimental results (Table II). In progressing toward improved force fields for Mg 2+ , it is important to reproduce experimental properties including solvation free energies, the activity derivatives, and the rate of water exchange. 17,19 For the latter, the approach presented here will greatly enhance the accuracy of the theoretical rate calculations and therewith the agreement between simulation and experiments.

CONCLUSION
Water exchange between the coordination shells of Mg 2+ is fundamental in understanding its reactivity and binding kinetics in aqueous solutions. Yet, simulating the elementary process of exchanging one strongly bound hydration water molecule from the The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp first hydration shell is tremendously challenging due to the long timescales involved. To overcome this challenge, transition path sampling is applied to gain detailed insights into the exchange dynamics and to resolve the molecular mechanism. The results reveal that the majority of pathways follow an indirect exchange mechanism, in which the entering and leaving water molecules occupy different positions on the water octahedron. A minority of pathways follows a direct exchange mechanism, in which the exchange takes place via the attack of the incoming water onto the edge of the water octahedron. In both pathways, one water molecule from the second hydration shell enters the molecular void leading to the concerted motion of another water molecule out of the first hydration shell. At the same time, the spectator molecules in the first hydration shell rearrange collectively such that the transition state has mirror symmetry.
In order to capture this reaction mechanism, a good reaction coordinate must include the degrees of freedom of all water molecules that play a decisive role in the exchange dynamics. Therefore, finding a perfect reaction coordinate as judged by the commitment analysis (and which yields the reaction rate without any further corrections i.e., κ = 1) is challenging due to the entangled contributions from the first hydration shell or even beyond. In the future, modern machine learning methods may lead to further progress in this complex many body problem. 49 Given the complexity of the exchange dynamics, it is not surprising that transition state theory overestimates the rate of water exchange by more than two orders of magnitude. On the other hand, the improved reaction coordinate presented here in combination with reactive flux methods allows us to provide an accurate value (as verified by high temperature simulations). Providing a more accurate estimate of the rate constant can contribute substantially to the development of improved Mg 2+ force fields for biomolecular simulations. 17,19,25 At the present time, the exchange rate using nonpolarizable force fields is significantly smaller (24 s −1 ) compared to experimental results (5.3 ⋅ 10 5 s −113 ). The discrepancy between exchange rates from simulations and experiments may result from the insufficient accuracy of the atomistic models for the water molecules and the Mg 2+ ions. For instance, the diffusion coefficient of the TIP3P water model differs by a factor of three from experiments. 28 Flexible water models might capture the dynamics of the hydrogen bonding network better. In addition, using polarizable models for the water molecules and the Mg 2+ ions could help to improve the agreement between simulations and experiments further. 24 The too slow exchange kinetics obtained with current rigid, non-polarizable atomistic force fields has significant consequences on metal cation binding to active sites on proteins or nucleic acids since the binding kinetics is limited by water exchange. However, with improved force fields and powerful path sampling approaches, direct sampling of the slow transitions including the binding of metal cations to RNA on the hundred millisecond timescale 50 should become possible in the future. 51

SUPPLEMENTARY MATERIAL
See the supplementary material for further discussion of the two-dimensional free energy landscape, the Jacobian correction, and the rate calculation at 694 K.