Adaptive resolution simulation of polarizable supramolecular coarse-grained water models

Multiscale simulations methods, such as adaptive resolution scheme, are becoming increasingly popular due to their significant computational advantages with respect to conventional atomistic simulations. For these kind of simulations, it is essential to develop accurate multiscale water models that can be used to solvate biophysical systems of interest. Recently, a 4-to-1 mapping was used to couple the bundled-simple point charge water with the MARTINI model. Here, we extend the supramolecular mapping to coarse-grained models with explicit charges. In particular, the two tested models are the polarizable water and big multiple water models associated with the MARTINI force field. As corresponding coarse-grained representations consist of several interaction sites, we couple orientational degrees of freedom of the atomistic and coarse-grained representations via a harmonic energy penalty term. This additional energy term aligns the dipole moments of both representations. We test this coupling by studying the system under applied static external electric field. We show that our approach leads to the correct reproduction of the relevant structural and dynamical properties. C 2015 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http: // dx.doi.org / 10.1063 / 1.4923008]


I. INTRODUCTION
Biological systems are profoundly affected by the surrounding solvent.Water, as the most important solvent, plays a crucial role in phenomena such as protein folding and stabilization of nucleic acid structure.Simulations of such systems therefore require the use of an accurate water model.Commonly used atomistic (AT) water models, such as simple point charge (SPC) model, 1 perform well and can reproduce many experimental properties, but at the same time they are computationally expensive.
7][8] However, there is even greater gain when multiple (typically 3-5) water molecules are modeled as a single site model.The most used CG water model of this kind is the one developed by Marrink et al. in the MARTINI force field. 9,10In this model, the water is described as a neutral bead that represents four water molecules and has been parametrized using partitioning free energies.Due to the lack of explicit charges, the electrostatic interactions are unfortunately not properly screened.As a consequence, artifacts arise at the interfaces between water and other phases and in the vicinity of charged particles. 11,12o address these issues, two CG water models with explicit charges have been recently introduced in combination a) s.j.marrink@rug.nlb) praprot@cmm.ki.si with the MARTINI force field, i.e., the polarizable water (PW) 13 and the big multipole water (BMW) 11 model.Both models use three sites per bead that, like the standard MAR-TINI, represent four water molecules.Although they are computationally slower, their performance is much better in representing a number of physical phenomena.For example, the PW model shows improvement in reproduction of the air/water surface tension, the freezing temperature, 13 and the free energy estimate of the pore formation in a lipid membrane, 14 while the BMW model gives the correct dipole potential at the membrane-water interface 11 and can be regarded as orientationally polarizable.
However, there are cases where AT resolution is required in some part of the system, e.g., for an accurate description of solvation shells around biomolecules.The emerging multiscale simulation methods 12,[15][16][17][18][19][20] are the optimal way to tackle such situations as they bypass the trade-off between accuracy and efficiency of a model.2][23][24][25] In the AdResS approach, different resolution regions coexist within one simulation and molecules can freely move from one region to another without any barrier.Recently, a multiscale simulation of water was performed where four semi-harmonically connected SPC waters were mapped to a MARTINI water site. 26Simulation of a protein embedded in such solvent 27 was shown to give no notable differences compared to the fully atomistic simulation.
In this work, we broaden the applicability of the AdResS simulation by coupling the SPC water model to supramolecular CG water models that incorporate orientational polarizability.The multiscale solvent is compatible with the popular MAR-TINI force field and can therefore be readily used to solvate biomolecular systems.Due to the polarizable nature of the multiscale solvent, we expect particular useful applications to systems that feature strong electrostatic interactions, for example, multiscale simulations of lipid bilayers.

II. MULTISCALE MODEL
In our multiscale model, four atomistic water molecules are mapped to one polarizable coarse-grained bead, as illustrated in Figure 1.The multiscale character is achieved by splitting the simulation box along the x-coordinate into regions of different resolutions.High resolution is employed for the central AT region.On both ends, the AT region is flanked by the intermediate hybrid (HY) and subsequently low resolution CG regions.In the AT region, the water clusters, each with four water molecules, are modeled by the atomistic bundled-SPC 26,28 model.In this model, an attractive semi-harmonic potential is added between all oxygen atoms within a cluster.The potential has the following form: where r and r 0 = 0.3 nm are the current and the equilibrium distance between two given oxygen atoms, respectively.Thus, the attractive bonds between water molecules as defined by Eq.
(1) start acting smoothly from the onset of r 0 .The k s represents the force-constant, which is set to 1000 kJ/(mol nm 2 ).In addition, the oxygen-oxygen Lennard-Jones (LJ) parameters of the standard SPC water model are modified to match the density of SPC water.In particular, the C 12 -parameter is equal to 3.25 × 10 −6 kJ mol −1 nm 12 .The use of attractive potential between molecules within a cluster is a convenient solution to problems that arise with multi-molecule coupling. 28,29If one wishes to map more than one AT molecule always to the same CG molecule, the motion of those AT molecules with respect to the CG site must be FIG. 1. Schematic representation of the simulated multiscale water system.From the simulation box's center to the edges, the regions of different resolution change in the horizontal direction from the AT to CG. Above each region, the corresponding water representation is shown.In the AT region, we use the bundled-SPC model, 26,28 where water cluster consists of four water molecules that are connected by semi-harmonic springs.Either PW 13 or BMW 11  restricted.Otherwise, the AT molecules will diffuse in time too far away from each other and the correspondence between AT and CG coordinates will be lost.Note that the use of bundled-SPC water model brings about some limitations.In a recent paper by Gopal et al., 30 the performance of the bundled-SPC model with respect to standard SPC was critically assessed.In general, correspondence between the two models was very good, in particular concerning the thermodynamic behavior of amino acid residues and conformational analysis of short peptides.Differences between the water models were found for the structure of a coiled-coil dimer, and hydration of the active site of a serine protease, however.
In the low resolution region, we employ two CG water models with explicit charges, the PW 13 and BMW 11 models.Both are three-site models consisting of a van der Waals center that is modeled with a Lennard-Jones potential in the PW model and with a Born-Mayer-Huggins function in the BMW model.The central particle is bonded with two interaction sites that only interact via electrostatic interactions.The BMW model is rigid with all three particles charged, thus forming a multipole.On the other hand, in the PW model the angle between the two bonds that connect the charged site with the center particle is not constrained, making the model polarizable.

III. ADRESS AND COMPUTATIONAL DETAILS
To allow the water clusters to freely move across different regions and change their resolution, and consequently degrees-of-freedom (DOFs), on-the-fly we employ the AdResS method. 21,22In this method, the smooth transition between the two resolutions is achieved through HY region, where the intercluster forces between the two resolutions are interpolated with a weighting function w(x).We use a smooth weighting function, where w = 1 corresponds to the AT region and w = 0 to the CG region. 26The values of w in the HY region are between 0 and 1.The total force on a water cluster α is given by 21,22 where F at α β and F cg α β are the forces between clusters α and β, obtained from the AT and CG potentials, respectively.Note that the bonded forces within a cluster are left unweighted.To correct for unwanted density artifacts in the HY region, we apply an external thermodynamic (TD) force 31,32 which acts in the HY region.The weighting function w and the TD force depend on the X position of the cluster's center of mass (COM).
The TD force F T D is derived in detail in Refs.23, 31, and 32.Here, for completeness, we briefly summarize its derivation.To enforce a uniform density profile, we have to compensate for the chemical potential differences between the high and low resolution molecular models.Therefore, F T D is defined as a negative gradient of the effective excess chemical potential due to the intermolecular interactions. 31Numerically, this translates into an iterative formula, 23,32 where C is an appropriately chosen numerical prefactor.Note that the TD force plays a similar role for eliminating unwanted density fluctuations as effective boundary forces in hybrid atomistic-continuum methods. 33,34hen a cluster crosses a boundary between the different regimes it will gain or lose the DOFs.A local thermostat is needed to supply or remove the latent heat caused by the switch of resolution 22 in the AdResS method.
We prepare the initial system of clustered water molecules as in our previous work, 26 i.e., we first equilibrate a system of 2796 SPC water molecules.Afterwards, we run series of short atomistic simulations.At each step, the distances between oxygen atoms are measured and, if they are small enough and the oxygen atoms are connected to less than three other oxygen atoms, semiharmonic springs are added between them.This process is repeated until all the four-water bundles are formed.We then map the AT clusters to CG particles in such a way that the COM and dipole moments of both representations match.We simulate 2796 water molecules in an orthorhombic simulation box of size 11.2 × 2.8 × 2.8 nm 3 as depicted in Figure 1.The AdResS regions are set so that the width of AT and HY regions is 4.2 and 1.4 nm, respectively.
The simulations were performed using the ESPResSo++ software package. 35Newton's equations of motion are solved using the standard velocity-Verlet integrator with a time step of 1 fs.Water geometry is constrained with the SETTLE 36 algorithm.The electrostatic interactions are treated with the reaction-field correction. 37,38The dielectric constants are equal to ϵ 1 = 1, 2.5, and 1.3 and ϵ 2 = 54, 28 75.6, 13 and 74 11 for the SPC, PW, and BMW, respectively.The dielectric constants' indices 1 and 2 correspond to 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. 37,38The cutoff distance for all non-bonded interactions is set to r c = 1.4 nm.Furthermore, the atomistic LJ interactions were capped for distances smaller than 0.2 nm.The constant temperature of 300 K is achieved with the Langevin thermostat with a coupling constant 4.0/ps.All simulations were run for 10 ns under NVT conditions with periodic boundary conditions and minimum image convention.For reference, we performed also fully atomistic and fully coarse-grained simulations.
The technical implementation of this kind of multiscaling is not straightforward and deserves a few notes.In particular, we used the following integration procedure in the HY region.First, we compute a given cluster's CG COM force as , where summation runs over all CG particles in a cluster.This force is distributed to the corresponding AT particles.The total force acting on ith AT particle of a cluster is thus , where F at i is the force obtained from the AT potentials.Next, we integrate both the AT and CG particles' positions and velocities according to the velocity-Verlet algorithm.Then, the COM position and velocity of the CG beads are matched to the AT ones.This implementation scheme ensures that the COMs of the AT and CG representations are always superimposed.In addition, whenever a new representation comes into play, i.e., upon entering the HY region, the orientation is mapped to the existing one.Finally, in the HY region, thermostatting is applied to both representations.

IV. COUPLING OF ROTATIONAL DEGREES OF FREEDOM
So far, the AdResS has only been applied to systems where the CG molecular representation consists of a single particle with one interaction site and is thus rotationally invariant.In the present case, this is no longer true.Therefore, we have to consider in our resolution coupling also the orientational degrees of freedom.The question then arises: how can one introduce the rotational coupling between the AT and CG representations into the existing scheme and is such coupling really necessary?To address this issue, we first discuss some associated difficulties.
First, the cluster's orientation of either AT or CG representation can be defined through the dipole moment.It therefore seems reasonable to map the CG representation to the AT one in such a way that the dipole moments are superimposed.However, since the dipole moment is invariant under any rotation about its axis, such mapping is not unique.In other words, many orientations of the CG configuration correspond to the same AT configuration and vice versa.Additionally, the dipole magnitudes of the AT and CG models are by construction not equal.
Second, while the molecular masses of both representations are the same, this is not the case with the moments of inertia.Consequently, it is not clear how one should transfer the torques in the HY region in order to couple the clusters' orientations.Furthermore, inclusion of torques is not convenient from the algorithmic point of view anyway.
Third, as the TD force is an external central force acting on clusters' COMs, it cannot change orientations of clusters.Hence, it cannot be used to eliminate deviations in dipole orientations of the AT and CG models in the HY region.
Because of the above, we opted to introduce the rotational coupling in the HY region through an additional energy penalty term of the form where θ α is the angle between the AT and CG dipole moments of a given cluster α.The potential depends on all the AT and CG particles' coordinates of the cluster and is therefore able to orient both models so that the two dipoles point in the same direction.At the same time, it leaves the cluster's centerof-mass intact.In this way, only the direction of the dipole moment is coupled, while the magnitudes still differ.The corresponding forces acting on AT and CG particles are and FIG. 2. The average cosine value of the angle formed by the AT and CG representations dipole vectors.Labels 1 to 5 correspond to increased strength of the rotational coupling.From no rotational coupling and a random initial orientation (1) followed by no rotational coupling with a matched initial orientation (2) to rotational coupling with the force constant equal to 50, 100, and 500 kJ/(mol rad 2 ) (labels 3, 4, and 5, respectively).
respectively.The dipole moments are calculated as d at =  N at i=0 e i r i and d cg = N c g j=0 e j r j , where e i and e j are partial electric charges of the ith AT and jth CG particle, respectively, and N at = 12 and N cg = 3 are numbers of the AT and CG particles within the cluster.
To find the appropriate strength of rotational coupling, we test several different force constants K.In Figure 2, we plot the average cosine value of the angle formed by the dipole orientations of the AT and CG representations.In the case of no rotational coupling (labels 1 and 2), we observe no correlation FIG. 3. The properties of bundled-SPC/PW (top) and bundled-SPC/BMW (bottom) systems subject to external electric field.The plots show the average cosine value of the angle formed by the dipole vector of water cluster and vector normal to boundary and pointing toward the CG region.Left plots show the results with no rotational coupling, while the right plots are for the case with rotational coupling and the force constant K = 100 kJ/(mol rad 2 ).To allow greater insight, we plot separately the AT (dotted black) and CG (blue and green for the PW and BMW model, respectively) dipole components and the interpolated dipole moment with weighting function w(x) (red).The results from reference full blown AT and CG simulations are shown as circles.
between the orientation of the AT and CG representations of clusters.This is also true even if orientations are initially aligned when clusters enter the HY region.On the other hand, drastic improvement can be seen with incorporation of rotational coupling.For the value of the force constant K, we have chosen K = 100 kJ/(mol rad 2 ) (label 4) as it is moderate, 39 yet sufficient.
The significance of the rotational coupling can be further examined by the response of the system to an external electric field (Figure 3).We plot the average cosine values of angle α formed by the dipole and normal vectors pointing toward the CG region as a function of the x-coordinate.The static electric field is applied in the y-direction and has the intensity of 1.0 V/nm.The results are presented for two cases, i.e., without (left) and with (right) rotational coupling.In both cases, we find that under the external electric field the clusters' dipole moments align with the direction of the field (see Figure 7 for comparison for the case without the electric field).The extent of alignment is not equal in all AdResS domains; for both PW and BMW couplings, it is higher in the AT than in the CG region due to a higher dielectric constant of the bundled-SPC water. 11,13,26However, while in both cases, the orientation is correctly reproduced in the AT and CG regions, in the HY region, fewer artifacts are observed with the rotational coupling.In particular, the magnitude of the clusters' alignment is noticeably reduced.From results in Figure 3, it can be seen that the deviations of orientations in the HY region are the most pronounced near the AT/HY boundaries.This indicates that this effect is related to the atomistic water in the presence of the boundary.As discussed in Ref. 6, this effect could be further reduced by an increase of HY region's thickness although we want to keep the HY region as thin as possible.Note that if the electric field is applied, the TD force has to be further iterated in the presence of the field due to the difference in the electrochemical potential of both models.Similarly, further iteration is also required if one changes the force constant K, but the final TD force changes by a small amount as it is center-of-mass based and does not have any influence on the orientation.

V. RESULTS AND DISCUSSION
We coupled the bundled-SPC model to two polarizable CG water models and will now show that the structural and dynamical properties of the simulated system are well reproduced by the AdResS multiscale simulation.The results hereafter include the rotational coupling described in Section IV.The radial distribution functions (RDFs) of oxygen atoms (O-O), central CG particles (W-W), and cluster's center-ofmass (COM-COM) are shown in Figure 4 for the bundled-SPC model coupled to the PW (left) and BMW (right) models.All AdResS RDFs very closely match the reference results from fully atomistic and coarse-grained simulations.In the case of AdResS simulations, the RDFs are computed locally, i.e., we only consider the molecules either in the all-atom (AdResS AT) or the coarse-grained (AdResS CG) domain.Interestingly, considering that neither of the CG models is parametrized on the basis of structural properties, the COM-COM RDF for the PW model more closely resembles the AT results (see Figure 4, bottom plots).
Next, we examine the density across different AdResS domains.To this end, we compute the normalized density profile (NDP), i.e., the local density divided by the bulk density, along the direction of resolution change.Figure 5 shows the uncorrected (top plot) and corrected (middle plot) NDP for both polarizable water models.
Previous studies have shown that coupling different representations with different chemical potentials induces density artifacts in the HY region.Interestingly, apart from oscillations in the HY region, pronounced density oscillations are also observed in the CG region in the case of PW model without TD force.Normally, such density oscillations occur if the system is confined.The obtained results suggest that the HY/CG interface induces the structural order (in a similar manner as a FIG. 4. RDFs of oxygen atoms (top), central CG particles (middle), and cluster's COM (bottom) for bundled-SPC model coupled to PW (left) and BMW (right).The reference (overlaid circles) fully atomistic and coarse-grained RDFs are well reproduced by the AdResS simulation.Labels "AdResS AT" and "AdResS CG" denote the AdResS domain, i.e., the atomistic and coarse-grained regions, respectively.FIG. 5. Normalized density profile along the xcoordinate (the direction of resolution change).The boundaries between AT, HY, and CG regions are denoted by vertical lines.The inhomogeneities (top) that arise due to differences in the chemical potential at different levels of resolution are eliminated when the appropriate TD force is used (middle).Dotted blue and solid green lines indicate the PW and BMW models, respectively.The bottom plot shows the applied TD force.wall would). 28Similar wall-induced layering problems occur also in other multiscale methods where the unwanted effect has been studied, see, for example, Refs.40-42.However, the TD force removes this wall effect of HY/CG interface and flattens out the density profile.The bottom plot shows the converged TD force that acts on COMs of all clusters that are located in the HY region.From Eq. ( 2), it follows that the hybrid forces are nonzero in the part of the AT region that is within a cutoff distance away from the HY region.For this reason, a smoother density profile is achieved when TD force is extended into AT and CG regions.In particular, the TD force is applied in the range d AT − d sk in < r < d AT + d HY + d sk in , where d sk in = 0.3 nm.Similar extension of the TD force was used in previous AdResS simulations. 43As can be observed from Figure the TD force can have a complicated form and therefore cannot be described with simple analytical functions.
In Figure 6, we plot the distribution of the dipole moment's magnitude for water clusters of the bundled-SPC/PW system.The distribution from the AT simulation is compared to the AdResS simulation, where we average only over the clusters within the AT region.Similarly, the distribution obtained from the CG simulation is compared to the distribution averaged over clusters within the CG region of the AdResS simula-tion.Our results show that the distributions are very well preserved in the AdResS simulation.The fixed dipole moment (µ = 5.8 D) of the BMW model is shown with a vertical line.
Figure 7 shows the cluster's average magnitude (top) and average orientation (middle) of the dipole moment across the resolution regions.The average orientation is quantified by the average cosine values of angle α formed by the dipole vector and normal vector pointing toward CG region.In the HY region, we calculate the clusters' dipole moment in both AT and CG representations.The obtained values are weighted by the function w(X), defined earlier.The obtained results indicate that the orientation is random in both AT and CG regions as in the bulk water.The HY region therefore sufficiently removes all orientational artifacts that might take place at the AT/HY or CG/HY interface.
The difference between normal and tangential components of the pressure across the simulation box is shown in Figure 7 (bottom).The local pressure is calculated using the Irving and Kirkwood method, with the pressure tensor given by 44 where A is the surface area normal to the x direction and ρ(x k ) is the number density in slab k.The step function θ(x) is defined by θ(x) = 0 when x < 0 and θ(x) = 1 when x ≥ 0 so that molecules i and j contribute to a given slab if the line joining the molecules crosses, starts, or finishes in the slab.The normal component p N of the pressure is equal to p x x , while the tangential component p T is calculated as (p y y + p z z )/2.The term p N − p T is directly related to the surface tension, 44 which is non-zero in the HY region (see Figure 7).In fact, the surface tension at the AT/HY interface is higher than at the CG/HY transition boundary.The normal component p N variations are compensated by the work done by the TD force. 32owever, the tangential component p T oscillations cannot be controlled by the TD force.They reflect the fact that we are coupling different water models with intrinsically different surface tensions. 44Note that the compressibilities of the AT and CG water models are not identical either as their RDFs do not match.When water clusters move from one AdResS region to another they could in principle, due to resolution change, experience a potential barrier (see Figure 5, top panel).As a consequence, the clusters would rebound off the transition boundary.However, results in Figure 8 show that this does not occur in our case as the clusters can freely move across the boundaries. 45e show three sequential diffusion profiles for clusters' COM, i.e., at times: t = 0, 50, and 100 ps.The clusters are initially in a slab within the AT, CG, or HY region but diffuse in time throughout the whole simulation box.The obtained diffusion profiles are slightly asymmetrical because the lower friction in the CG region allows a slightly faster diffusion.The estimated diffusion constants, which are extracted by fitting the Green's function for the diffusion equation to the diffusion profiles, are 1.9, 2.0, and 2.2 × 10 −9 m 2 s −1 (error bars are roughly 10%) for the bundled-SPC, PW, and BMW models, respectively.They are in good agreement with our previous published value of 1.8 × 10 −9 m 2 s −1 for the bundled-SPC water. 26

VI. CONCLUSIONS
In this paper, our aim was to develop a polarizable multiresolution water model that can be used in multiscale simulations of biological systems.We employed the AdResS method, which concurrently couples the atomistic and coarse-grained representations and is based on a force interpolation scheme.The 4-to-1 mapping applied in this study was made possible by using the bundled-SPC model in the high resolution region.We successfully coupled the atomistic system to two different supramolecular, explicit-charge, CG water models-the PW and BMW models.Our approach, however, is not limited to 4-to-1 mapping but can be applied to any multi-molecule mapping.In this regard, one interesting future application would be to use the two-site water model developed by Riniker and van Gunsteren 46 that represents five atomistic waters.
FIG.1.Schematic representation of the simulated multiscale water system.From the simulation box's center to the edges, the regions of different resolution change in the horizontal direction from the AT to CG. Above each region, the corresponding water representation is shown.In the AT region, we use the bundled-SPC model,26,28 where water cluster consists of four water molecules that are connected by semi-harmonic springs.Either PW13 or BMW11 model (both with three charged particles) is used in the CG region.The water clusters (shown as light blue beads) can change their resolution on-the-fly as they move through the HY region.
Top plot depicts the average absolute value of dipole moment along the x-coordinate of the simulation box.The average cosine value of the angle formed by the dipole vector of water cluster and vector normal to boundary and pointing toward the CG region is shown in the middle plot, while the difference between the normal and tangential component of the pressure is shown in the bottom plot.The dashed-dotted blue and solid green lines denote the PW and BMW models, respectively.