Adaptive resolution simulation of salt solutions

We present an adaptive resolution simulation of aqueous salt (NaCl) solutions at ambient conditions using the adaptive resolution scheme. Our multiscale approach concurrently couples the atomistic and coarse-grained models of the aqueous NaCl, where water molecules and ions change their resolution while moving from one resolution domain to the other. We employ standard extended simple point charge (SPC/E) and simple point charge (SPC) water models in combination with AMBER and GROMOS force fields for ion interactions in the atomistic domain. Electrostatics in our model are described by the generalized reaction field method. The effective interactions for water–water and water–ion interactions in the coarse-grained model are derived using structure-based coarse-graining approach while the Coulomb interactions between ions are appropriately screened. To ensure an even distribution of water molecules and ions across the simulation box we employ thermodynamic forces. We demonstrate that the equilibrium structural, e.g. radial distribution functions and density distributions of all the species, and dynamical properties are correctly reproduced by our adaptive resolution method. Our multiscale approach, which is general and can be used for any classical non-polarizable force-field and/or types of ions, will significantly speed up biomolecular simulation involving aqueous salt.


Introduction
Solvents have a major effect on the functional properties of macromolecules, e.g. DNA molecules and proteins. Typically at least a monolayer of water is required for full protein functionality. Another important concern is also how a solvated macromolecule influences the local structure and dynamics of the surrounding solvent. Hence, a detailed study of interactions between a macromolecule and a solvent is required for an understanding of the macromolecule's structure, dynamics and function. To this end, chemistry specific interactions on the atomic level of detail have to be considered in biomolecular simulations. However, due to a large number of atoms and the corresponding degrees of freedom (DOFs) such systems are difficult to study using all-atom computer simulations. Moreover, the majority (up to 90%) of the simulation time is spent on the solvent, i.e. usually water or aqueous salt in biomolecular simulations, and not on the DNA molecule or the protein itself that are of our primary interest. To bridge the gap between the time and length scales accessible to simulations at an atomistic level of detail simplified models of the solvent have been developed.
Specifically, there have been presented many different simplified models for ionic solutions such as implicit solvent methods [1,2] where the solvent effects are described by a concentration-dependent dielectric constant. In another type of implicit solvent approaches, solvent free models for aqueous ionic solutions are derived using many-body interaction potentials and the generalized Langevin equations [3]. Simplified explicit solvent methods for ionic solutions, on the other hand, involve the development of coarse-grained models where only the relevant DOFs are retained. Typically coarse models are derived by a structure based ansatz [4][5][6][7][8][9][10][11][12][13][14][15][16] focusing on reproducing structural aspects or based on free energy mappings [17][18][19]. However, there are many situations in which resulting solvating phenomena manifest themselves simultaneously on atomistic, mesoscopic and macroscopic scales [20]. Multiscale modeling approaches that combine various simulation models represent the most efficient way to bridge many orders of magnitude in the spatial and temporal scales involved in such systems.
Interfacing descriptions on various length scales of molecular liquids and soft matter in a multiscale framework has attracted significant interest in recent years [21][22][23][24][25]. One type of these approaches focuses on coupling molecular dynamics with the continuum description of liquids using the Navier-Stokes equation [26][27][28][29][30][31][32][33][34][35][36][37]. Another type of method is concerned with concurrent coupling of particle-based models on different length scales, e.g. all-atom and coarse-grained molecular models [38][39][40][41][42][43]. In recent years, we have been developing the adaptive resolution scheme (AdResS) [24,44]. AdResS concurrently couples different particle descriptions of the molecular liquids from quantum to all the way to continuum scale [33,[44][45][46]. Very recently, a Hamiltonian version of the method H-AdResS was proposed, which also opens doors for adaptive resolution Monte Carlo simulations [47]. Efficient implementations of AdResS in molecular dynamics software packages, i.e. ESPResSo [48,49], GROMACS [50,51], ESPResSo + + [52], enable us to perform fast adaptive resolution simulations of soft matter and molecular liquids. This is especially important for simulations of ionic solutions where we need large systems to obtain acceptable statistics for the ions.
In this paper, we present the adaptive resolution simulations of aqueous salt (NaCl) using AdResS. For this purpose, we develop coarse-grained models of salt to be used with standard force fields and derive thermodynamic forces [53,54] to ensure the thermodynamic equilibrium distribution of all molecular species across the simulation box. Electrostatic interactions are modeled using the generalized reaction field method [55]. The developed multiscale model of salt solution will significantly speed up biomolecular simulations where we need atomistic descriptions of salt only in certain localized subdomains of the system, e.g. in the vicinity of a DNA molecule.
The paper is organized as follows. In section 2 our multiscale model of aqueous salt is presented. The simulation details are given in section 3. The results and discussion are reported in section 4, followed by conclusions in section 5.

Multiscale model of aqueous sodium chloride
We derive a hybrid all-atom/coarse-grained model for a salt solution, i.e. Na and Cl ions solvated in liquid water at ambient conditions, to be used in adaptive resolution simulations [24,44,56]. We first run atomistic simulations to obtain a reference system. We employ standard SPC/E [57] and SPC [58] water molecules in combination with AMBER [59] and GROMOS [60] force fields for ions (see section 3). The generalized reaction field method [55] is used for electrostatic interactions, where the force between point charges q i and q j of atoms i and j is calculated as 1 and 2 are dielectric permittivities of the inner (where the electrostatic interactions between point charges are explicitly considered up to the cutoff radius R c ) and outer (where water is considered as a dielectric continuum) regions, respectively, as explained in the literature [55]. The cutoff radius R c is set to 0.9 nm, 1 = 1, 2 = 80 and the Debye screening length κ = 3.25 nm −1 corresponding to a 1 M salt solution [55]. This concentration is higher than the physiological concentration (approx. 0.15 M), because we need an adequate number of ions for acceptable statistics. Moreover, many experiments are commonly performed at 1 M concentration. However, lower concentrations, e.g. 0.15 M, would not require a different method, just larger simulation systems. We have chosen the generalized reaction field method for electrostatics because it is very convenient for implementation in AdResS, i.e. it is pairwise and short-ranged. Instead, we could use the Ewald summation [61,62] for computing electrostatic interactions between ions. That approach, though not followed here, will be explored in future works.
After conducting the reference atomistic simulations, we set up a coarse-grained model by mapping one water molecule to a one-bead representation. For the coarse-grained model, we derive an effective pair potential for the water-water and water-ion interactions, such that structural properties of the underlying atomistic model are reproduced. We obtain the effective interaction potential following the structure based coarse-graining strategy using the iterative Boltzmann inversion (IBI) method (without the pressure correction) [5,7]. The corresponding coarse-grain force is computed in the standard way, as the negative gradient of the potential. For ion-ion interactions we also set 1 = 80 for the generalized reaction field, to properly screen the electrostatic interactions between ions in the coarse-grained model. Ion-ion Lennard-Jones interactions parameters are the same as in the atomistic simulation and charges are also the same as in the atomistic model. Ions thus interact via the same potentials in the atomistic and coarse-grained models aside from the changed dielectric permittivity ( 1 ).
With the atomistic and coarse-grained systems set up, we then employ the AdResS method. This multiscale method allows molecules to move freely between regions, changing their resolution and DOFs on-the-fly. The total force on a molecule α reads as [53,54] where X α and X β are center-of-mass positions of molecules α and β respectively, F atom αβ and F cg αβ are the forces between molecules α and β obtained from the potentials of atomistic and coarsegrained representations, respectively. A weighting function w(X ) ∈ [0, 1] ensures a smooth transition between the two regions. F th α is a thermodynamic force that compensates the effective chemical potential differences (since we do not have a Hamiltonian) across different levels of resolution and thus ensures a uniform density profile throughout the whole simulation box [25,53,54]. There are three thermodynamic forces corresponding to water molecules, Na and Cl ions. The molecular index α determines the corresponding thermodynamic force to be used. A local thermostat is used to supply or remove the latent heat associated with the switch of resolution [25,63].
Using the AdResS scheme given by equation (2) in combination with the generalized reaction field given by equation (1) we treat the electrostatic interactions in our adaptive resolution simulation similarly as this is done in fixed partitioning hybrid atomistic/coarsegrained methods [64,65]. However, since we use a one-to-one molecule mapping for the coarse-grained water model we use 1 = 80 for ion-ion interactions in the coarse-grained region instead of much lower values in case of four-to-one [65] and five-to-one [64] molecule mappings. Moreover, we do not adjust the value of neither atomistic ( 1 = 1) nor coarse-grained ( 1 = 80) relative dielectric permittivities for ion-ion electrostatic interactions between different resolution domains as it is done in [64,65]. Instead, we simply employ our AdResS scheme.

Simulation details
Simulations are performed with ESPResSo + + [52] which is an open source multiscale simulation package designed to perform parallel adaptive resolution scheme simulations of Figure 1. A schematic representation of the simulation box. The resolution changes along the horizontal direction, such that an explicit all-atom region is at the center of the simulation box surrounded by a coarse-grained region filled by blue one-cited coarse-grained water molecules and ions. Although ions (i.e. Na(yellow) and Cl(green)) are one sited in both representations, they are characterized by different interactions with water molecules and the screened interactions between themselves. soft matter. We simulate two systems, one with box size 16.08 × 4.02 × 4.02 nm with 8377 SPC/E water molecules, 157 Na and 157 Cl ions, which corresponds to 1 M NaCl solution. For the second system we use SPC waters instead (see the supplementary data, available from stacks.iop.org/NJP/15/105007/mmedia). For the integration a standard velocity Verlet with time step 1 fs is used [66,67]. The temperature is regulated at 300 K with a local Langevin thermostat, with the value of the friction constant 15 ps −1 . Standard periodic boundary conditions and the minimum image convention are employed. All simulations are run for 10 ns.
The all-atom systems are set up in the following way: for the system with SPC/E waters we use Lennard-Jones and electrostatic interactions with AMBER99 force field. For the system with SPC water molecules we use Lennard-Jones and electrostatic interactions with standard GROMOS force field parameters between atoms of water molecules and ions. The geometry of water molecules is constrained with the SETTLE [68] algorithm. Note that our multiscale model of aqueous salt solution could be used with any non-polarizable classical water model [69,70].
The interaction potentials for the coarse-grained water-water and water-ions are obtained using the VOTCA package [71,72]. As the initial approximation for water-water effective interactions the coarse-grained water model for the pure water case from the literature is chosen [54].
The AdResS simulation is set up with one-dimensional splitting along the x-coordinate, such that an explicit all-atom region is at the center of the simulation box sandwiched in between two coarse-grained domains as depicted in figure 1. The size of the atomistic region is set to 4 nm, and the width of the hybrid region is 3 nm on each side of the atomistic region. The thermodynamic forces are obtained with an iterative procedure as described in the literature [25,53,54]. Without thermodynamic forces, both the water molecules and the ions start to gather in the atomistic region. First we iterate the thermodynamic force for the water molecules until a satisfying normalized density profile (NDP) is reached. Here we assume that the ion distributions will have a small effect on the water density profile. Hence the thermodynamic force for water should be decoupled from the thermodynamic forces for ions. Our results (see below) indeed confirm this assumption. The obtained force is then used as the initial RDFs of water-water centers of mass from coarse-grained, AdResS and all-atom simulations (bottom). Effective pair potential for water-water center of mass interactions between molecules (top). Atomistic RDF cm is well reproduced by the fully coarse-grained and AdResS simulation (both in the explicit atomistic (AdResS ex) and coarse-grained (AdResS cg) regions). The results for the pure water case [54] are also presented for comparison.
thermodynamic force for both species of ions and we iterate further only the thermodynamic forces for ions [53,54]. Our multiscale model can also be used for adaptive resolution simulations with other packages like ESPResSo [48,49] and GROMACS [50,51].

Results and discussion
First, we calculate the radial distribution functions (RDFs) for atomistic, coarse-grained and AdResS simulations. In all simulations we consider centers of mass of water molecules and ions for RDF calculation. In the all-atom and coarse-grained simulations we average over all molecules in the box. To check the local structure in the AdResS simulations we average only over the molecules either in the all-atom (ex) or the coarse-grained (cg) domain. Figure 2 shows water-water center of mass RDFs cm (bottom) and effective pair potentials (top) between waters' centers of mass for the SPC/E system. The effective potentials obtained with the IBI method used for the coarse-grained molecules very well reproduce the atomistic RDFs cm . In particular, the water-water effective potential for the SPC/E system is similar to the water-water effective interaction potential for the pure water case [54] (see figure 2 (top)). However, despite the similarity all our salt solution simulations reproduce the known effect of adding salt, i.e. the first peak in the water-water center of mass RDFs cm weakens while the second peak moves inwards [73,74]. The water structure is modified due to the presence of the ions similarly as due to increased pressure [73]. Figure 3 shows water-Cl RDFs cm . Figure 4 shows water-Na SPC/E water -Na Cl Na SPC/E water pure SPC/E water Figure 5. Thermodynamic forces used for the corresponding molecules. The x-axis shows the distance from the center of the simulation box in x-direction. Vertical lines denote the boundaries of the hybrid region. The thermodynamic forces are symmetric about d = 0. The thermodynamic force for the pure water AdResS simulation [54] is also shown for comparison.
Next, to verify that the density is homogeneous across the different regions of the simulation box, we calculate the NDP along the x-coordinate for centers of mass water molecules and for the ions. Computing the thermodynamic force to ensure a homogeneous density of the water molecules is done in the same way as for the pure water system [54]. However, for the ions it requires slightly more effort due to electrostatics and poor statistics. Despite having a relatively large box, the number of ions in the system is small, so statistics for the NDP is worse than for water. Moreover, ions interact through electrostatics. Thus, one can expect correlations between the thermodynamic forces of different ion species due to correlated density fluctuations of the ion pairs. Indeed, applying a specific thermodynamic force to one type of ion causes changes in the NDP of the other ion type [53]. On the other hand, we observed in the thermodynamic force iteration procedure that ions do not have an impact on water molecules' density distribution. Figure 5 shows the thermodynamic forces used for the corresponding molecule types for the SPC/E system (see the supplementary data (available from stacks.iop.org/NJP/ 15/105007/mmedia) for the SPC thermodynamic forces). Water and Cl thermodynamic forces are approximately of the same magnitude, while the thermodynamic force for Na ions is larger indicating the larger effective chemical potential difference between the corresponding all-atom and coarse-grained models. Thermodynamic forces act mostly in the hybrid region. The thermodynamic force for water is almost identical to the pure bulk water case [54]. This facilitates the derivation of thermodynamic forces for ions as the thermodynamic force for water can be kept almost unchanged during the iterative procedure to determine the thermodynamic forces on ions. Figure 6 depicts NDPs cm for all three molecule types for the SPC/E system (see the supplementary material for SPC system NDPs). NPDs are shown along the x-coordinate, as this is the direction in which the resolution of the molecules changes. Shown are NDPs from three simulations: all-atom (bottom), AdResS with thermodynamic forces acting on all molecule types (middle) and AdResS with thermodynamic forces acting on water molecules only (top). One can clearly see that thermodynamic forces on ions are also needed for a flat density profile.  [75]. They are larger for the ions than water due to a lower number of the ions.
As can be seen in the middle of figure 6 the NDPs cm of water molecules have only very small deviations (up to 5%) from the ideally flat profile, while the NDPs of ions is not as smooth (deviations up to 20%). However, ion density fluctuations of the same order of magnitude are also present in the all-atom simulations shown at the bottom of figure 6. There these larger fluctuations are due to rather poor statistics (see the error bars). However, this does not prove that these fluctuations are due to poor statistics, only in the AdResS case. Nevertheless, we can conclude that within our error bars the density profiles are correctly reproduced by AdResS.
To check the dynamical properties we have monitored the mean square displacements (MSDs) of different molecular species. Figure 7 shows MSDs for all-atom, AdResS and coarsegrained simulations for the SPC/E system (see the supplementary material for SPC system MSDs). The corresponding diffusion constants are in agreement with those reported in the literature [76]. As expected, the diffusion constants of the coarse-grained models are higher than the all-atom counterparts due to softer interaction potentials [25]. Since in the AdResS simulations particles move in both resolution domains, the associated diffusion constant values are in between the atomistic and coarse-grained values. If we wanted to match the AdResS diffusion constants to all-atom ones we would have to resort to position-dependent local thermostats [63,70]. Finally, we estimate the speedup of AdResS compared to the full-blown atomistic simulation. The speedup is due to the coarse-grained model, which allows in our case for three times faster simulations using the same time step as in the atomistic case. The actual speedup of AdResS simulations thus depends on the ratio between the atomistic and coarse-grained domain sizes, i.e. up to three times.
To summarize, our results show that equilibrium structural and dynamical properties of aqueous salt solution at 1 M are well reproduced with adaptive resolution simulation. The derived effective coarse-grain interactions and thermodynamic forces ensure the thermodynamic equilibrium of all molecular species.

Conclusions
In this paper, we have introduced a hybrid atomistic/coarse-grained model of 1 M aqueous NaCl at ambient conditions for adaptive resolution simulations. We have derived the effective coarse-grained interactions and thermodynamic forces and shown that equilibrium structural and dynamical properties, e.g. radial and density distributions as well as diffusion constants, are well reproduced using AdResS. The presented multiscale approach is general and can be used for any type of ions and/or force field. The same method could be also used for lower concentrations than 1 M. However, larger systems would have to be simulated to obtain the acceptable statistics for ions. Our multiscale model allows us to perform efficient biomolecular simulations involving salt solutions, e.g. solvated DNA molecules. Work along these lines is underway.