Efficient Osmotic Pressure Calculations Using Coarse-Grained

Letter Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

Efficient Osmotic Pressure Calculations Using Coarse-Grained Molecular Simulations Jörg Sauter and Andrea Grafmüller* Department of Theory and Biosystems, Max Planck Institute of Colloids and Interfaces, Potsdam 14424, Germany ABSTRACT: Osmotic pressure data is increasingly used to parametrize all-atom simulation Force Fields (FFs), leading to large computational cost for larger molecules. Here, we show that the osmotic pressure can be calculated precisely using transferable coarse-grained FFs obtained from short atomistic simulations using an inhomogeneously regularized coarsegraining procedure. This is demonstrated for carbohydrates, where compared to the equivalent atomistic system, an increase of the computational efficiency by a factor of ≈500 is achieved.

T

namic quantities, such as the activity coefficients, has been applied to parametrize a number of small molecules.1−5 Similarly, the McMillan and Mayer theory6 has been used recently, to optimize the MARTINI model7 to reproduce the second virial coefficient of the osmotic pressure. An alternative is the direct use of osmotic pressure data as a target property. The osmotic pressure method (OPM) of Lou and Roux8 that measures the osmotic pressure directly is increasingly being used to optimize FF parameters. Osmotic pressure data has been used to improve FF parameters for ion solutions8−10 and amine-phosphate interactions,11 multisite,12 and implicit solvent13 ion models. Application to the larger molecules GLY1, GLY2, and GLY3 has highlighted that agreement with experimental data becomes poorer with increasing molecules size, underlining the limited transferability of parameters from small model molecules to larger ones.2 Recently, the OPM has been also applied to optimize amino acid−amino acid interactions14 and to improve the performance of different carbohydrate FFs with respect to their ability to reproduce solute−solute aggregation15,16 and amino acid− sugar interactions.17 To correct the trend of the osmotic pressure as a function of the degree of polymerization for saccharides, osmotic data for oligomers up to trimers has been matched for various concentrations in the optimization of the 16,18 GLYCAM06TIP5P OSMOr14 FF. The above examples illustrate that standard mixing rules often perform poorly at capturing the balance of solute−solute and solvent interactions correctly. As the osmotic pressure data is used to optimize only the solute−solute interactions, the agreement with the previously matched target data utilized to obtain the solute−solvent interactions will remain largely

he predictive power of molecular simulations depends primarily on how well the force field used can describe the system of interest. Therefore, force field parametrization is a key aspect in the molecular modeling process. There are different strategies for force field development, and balancing the contributions of different interactions remains an ongoing challenge. Many macromolecular and solution properties depend critically on the description of the nonbonded interactions. These are most often described by a combination of electrostatic and Lennard-Jones interactions. To represent the electrostatic properties, atoms carry partial charges, which are typically derived based on fits to quantum mechanical calculation data. The Lennard-Jones interaction potential takes the form ⎡⎛ σ ⎞12 ⎛ σ ⎞6 ⎤ U (r ) = 4ε⎢⎜ ⎟ − ⎜ ⎟ ⎥ ⎝r⎠ ⎦ ⎣⎝ r ⎠

where the free parameters ε and σ are typically adjusted to match properties of single compounds such as the enthalpy of vaporization or sublimation (target parameter ε) and of binary mixtures, for example the mass density (target parameter σ) of small model compounds. However, these target observables contain little information on the balance of solute and solvent interactions. In addition, the interactions between different atom types, and especially between solutes, are often obtained based on mixing rules for combining these parameters. Thus, they do not necessarily reflect the correct effective aggregation behavior. Furthermore, parameters cannot always be readily transferred to larger molecules, as cooperative or screening effects due to the specific molecular structure may arise. Different strategies to include experimentally accessible solution properties in the parametrization of the solute−solute interactions have been developed. The Kirkwood-Buff theory, which relates the radial distribution functions to thermody© XXXX American Chemical Society

Received: December 5, 2017 Published: February 5, 2018 A

DOI: 10.1021/acs.jctc.7b01220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Letter

Journal of Chemical Theory and Computation unaffected. This has been shown for instance for the mass density and molecular structure of small sugar molecules.16 While such computations are still feasible using atomistic resolution for small polysaccharides, the computational cost soon becomes too high for larger molecule sizes, as the necessary sampling time to reach converged results increases rapidly due to the increasingly slow dynamics of the polymer systems. For the glucose trimer studied in ref 16, for instance, a simulation time of 1.5 μs was required in order to obtain converged results. The high computational cost of measuring osmotic pressure for larger molecules renders the process of scanning the parameter space for force field parametrization to match osmotic pressure trends very time-consuming. In addition, it prevents the application of the method to systems involving larger and more complex molecules. These limitations may be overcome by the use of coarse-grained (CG) models with reduced number of degrees of freedom. However, the effective CG interactions are typically set up to represent one state point, and reproducing thermodynamic data, such as the osmotic pressure, is often a challenge. Recently, we have demonstrated that CG-FFs based on the Multiscale CoarseGraining Method (MS-CG)19 can correctly reproduce systems of aqueous polysaccharides with respect to radial distribution functions.20 Since this CG procedure is based on reproducing the forces in the atomistic system with the CG interactions, and since the OPM calculates the osmotic pressure based on the force exerted on the solutes, this procedure appears to be a promising approach for CG osmotic pressure measurements. Here, we show that CG-FFs obtained with a similar, refined coarse-graining method using inhomogeneous regularization and employing separations of solute−solute, solute−solvent, and solvent−solvent interactions can be used to precisely calculate the osmotic pressure using an adapted OPM, optimized to work at the lower resolution. Finally, we demonstrate that this coarse-grained OPM (CG-OPM) can be used with other CG FFs and apply it to optimize the carbohydrate MARTINI21 force field.

pressure results. For CG simulations using the MARTINI FF,21,32 the standard input options for MARTINI 2.x were used, with no pressure coupling except during equilibration. Using a 4 fs time-step, 200 ns were sampled for the osmotic pressure measurements. 1.2. Coarse-Graining Method. We employ the coarsegraining procedure based on the MS-CG method19 described previously20 to derive the CG interactions. However, the following three modification to the procedure were required to obtain reliable results in the OPM setup: (i) We now treat solvent−solvent interactions separately in the MS-CG procedure, whereas previously, only solute−solute interactions were separated from interactions involving solvent. This is required, as in the osmotic pressure setup a separate pure water compartment is simulated, whereas in ref 20 the MS-CG method was applied to simulate mixed solute−solvent systems. Note that this separation only refers to the coarsegraining process; the sampling of the atomistic forces is performed in a mixed solution as before. (ii) The interaction range in the MS-CG method was extended from 1.0 to 1.8 nm, because of the significant cutoff sensitivity of the osmotic pressure.16 (iii) To be able to use a 1.8 nm cutoff despite the profound sensitivity to perturbations at long interaction ranges (from ≈1.0 nm) described previously for the CG solute− solvent systems, we use an inhomogeneous regularization in the MS-CG procedure: i.e. instead of solving

1. METHODS 1.1. Molecular Dynamics. Simulations were performed with a modified version of GROMACS 5.1.2,22 in which the wall force behavior was modified such that all solute particles experience the same wall force irrespective of the particle type. TIP5P For all atom simulations, the GLYCAM06OSMOr14 FF16 parameters were used with a 1.4 nm cutoff using the group cutoff-scheme. All covalent bonds involving hydrogen atoms were constrained with LINCS,23 and water molecules are kept rigid using Settle.24 Electrostatic interactions were treated with the Particle Mesh Ewald method25 using a 1.4 nm cutoff. For the osmotic pressure setup, the EW3DC26 method was used for long-range electrostatic interactions. Simulations were run using the Leap-Frog integrator27 with a 2 fs time-step. The temperature of 298.15 K was controlled using the Nosé−Hoover thermostat.28,29 Energy minimization and equilibration were performed using standard protocols. The pressure in the equilibration simulations was set to 1 bar with the Parinello-Rahmann barostat.30,31 All other simulations were performed in the NVT ensemble. Dispersion correction was used for energy and pressure, where applicable. For the CG simulation, using the MS-CG potentials, the tabulated interactions were cut off at 1.8 nm. 20 ns of sampling time was more than sufficient to obtain converged osmotic

α(r ) = 100

argmin FΦ − f

2

Φ

with the Euclidean norm ∥·∥, the basis matrix F, the atomistic force vector f, and the set of CG-FF parameters Φ, we now solve argmin( FΦ − f

2

+ diag(α)Φ 2 )

Φ

where diag(α) denotes the diagonal matrix with the diagonal entries given by a regularization vector α. The entries of the vector α are calculated according to the regularization function 2

( 1.8r )

+ 0.5 which is plotted in Figure 1. This coarse-graining procedure was performed with the MSCGFM 1.54 code, which has been modified such that the regularization

Figure 1. Regularization function α(r) used for the calculation of α. B


Letter

Journal of Chemical Theory and Computation

16 Figure 2. Calculated osmotic pressure using the GLYCAM06TIP5P and the CG-FF obtained at the highest OSMOr14 FF using atomistic resolution 35−37 is shown as well. concentration. The experimental data

is applied after preconditioning of F.33,34 The CG-FF was obtained at the highest concentration for all molecules and directly transferred to all lower concentrations. We chose the highest concentration, as this has previously been shown as the best choice to gain transferability.20 The results described below show that the application at lower concentrations still performs very well. 1.3. Coarse-Grained Osmotic Pressure Setup. The method we use to measure the osmotic pressure in the atomistic system is similar to the OPM by Luo and Roux8 as described in detail in ref 16. In this method, virtual walls exert a constant wall force on the atoms of the solute molecules if they move outside the central region of the box. The osmotic pressure is calculated from the average wall force. To apply this method to the CG system, the method parameters have to be adapted to take into account the reduced number of degrees of freedom. In the z-direction, the water molecules are kept in the box by a wall, exerting a force of 500(1−d)12 kJ mol−1nm−1 on any water molecule (1-site) within a distance d ≤ 1 nm from the nearest wall. This increase in the force constant by a factor 50 compared to the atomistic setup is necessary to keep the water confined. In addition, the water bath size was extended from 2 to 4 nm, to prevent boundary effects. A further increase produces virtually no deviation in the calculated osmotic pressure. Finally, the constant wall force experienced by solute beads was increased from 32 kJ mol−1 nm−1 to 64 kJ mol−1 nm−1 to ensure that the solutes are confined in the central region.

The simulation box for the CG osmotic pressure system using the MS-CG interactions was set up at the atomistic scale, using NPT simulations to equilibrate the mixed solute solvent system in a cubic box for 50 ns. Afterward the average box size obtained from the NPT simulation was extended in the z direction by 4 nm at either side, filled with water, and allowed to equilibrate in NVT. Water molecules were added until the correct water density was reached. Then the system was mapped onto the CG-coordinates. For the MARTINI Model the same setup procedure was used directly in the CG system.

2. RESULTS AND DISCUSSION 2.1. Coarse-Grained Osmotic Pressure. Osmotic pressure results using the MS-CG procedure and the CG-OPM are shown in Figure 2 for four test molecules. The CG model reproduces the atomistic results precisely. The plot of the calculated osmotic pressure as a function of simulation time, shown for the highest glucose concentration of 6.17m glucose in Figure 3, shows that convergence is achieved already within 2 ns. This represents a vast speed up compared to 400 ns required for convergence in the fully atomistic simulation.16 In addition to the ≈200× faster convergence, the CG simulation is faster by a factor of 2.4 per time-step using 12 cores, both for the atomistic and the coarse-grained system. Thus, an overall increase in computational efficiency by a factor of ≈480 can be achieved. The simulation time required to set up the CG interaction potentials, i.e. 100 ns of sampling, is much shorter than what is required to measure the osmotic pressure. In addition, the C


Letter


significantly with this model, which is equivalent to an overaggregation of glucose. This can be corrected by reducing the ε parameter of the solute−solute interactions by a factor of 0.95. The corresponding C6 and C12 parameters are listed in Table 1. The sugar− Table 1. New Nonbonded C6 and C12 Parameters for the Saccharide−Saccharide Interactions of the Reparametrized MARTINI FF

Figure 3. Average osmotic pressure versus simulation time for a 6.17m glucose solution.

atom type

atom type

C6

C12

P4 P1 P4

P4 P1 P1

0.204801 0.184319 0.184319

0.00220761 0.00198683 0.00198683

water interactions remain unchanged to keep the solute partitioning, for which MARTINI was parametrized, unchanged. The osmotic pressure data for these parameters, also shown in Figure 4a), is now in excellent agreement with the experimental results. Thus, a smaller change than suggested from the second virial coefficient7 is sufficient to correct the osmotic pressure. The solute−solute radial distribution functions for a 1m glucose solution for the two sets of parameters are shown in Figure 4b) together with the one obtained from the atomistic simulations, as well as the MS-CG simulations. The radial distribution functions clearly reflect the weaker tendency to aggregate with the new MARTINI parameter set, though the local solution structure reflected by the peak positions remains unchanged. On the other hand, despite leading to the same osmotic pressure, the radial distribution functions obtained both from the atomistic simulations and with the MS-CG interaction potentials for the GLYCAM06TIP5P OSMOr14 FF show much less local structuring, illustrating that even though the osmotic pressure is sensitive to details in the radial distribution function, inversely, different solution structures can produce the same osmotic pressure. As expected, many other solution properties such as the mass density and the diffusion coefficient remain virtually unchanged with the adjusted solute−solute interactions for glucose concentrations up to 4m.

simulated system size is smaller, since no water bath is simulated, and only one simulation at the highest concentration is required to calculate the osmotic pressure versus concentration curves. Thus, the combination of the MS-CG procedure and the OPM offers a very efficient approach to calculate concentration dependent osmotic pressure data. This opens up the possibility to explore the deviation from the ideal behavior for larger or more complex sugar molecules, than can be easily studied with atomistic resolution. Furthermore, such an approach could also be used to include concentration dependent osmotic pressure data in the force field parametrization process more efficiently, by scanning the parameter space at CG resolution. Even if the final parameter set is then evaluated at the atomistic scale, such an approach would save a significant amount of computational time. For example, for the solute−solute interaction parameters of the GLYCAM06TIP5P OSMOr14 FF, the effect of scaling several interdependent parameters was scanned over a range of concentrations and for several molecules. This computationally demanding procedure is exactly the type of application, which would greatly benefit from the CG procedure to calculate the osmotic pressure described here. 2.2. Optimized MARTINI Model. Similarly, the optimized CG-OPM can be used with other CG-FFs. For example, here we apply this method to calculate the osmotic pressure of glucose with the MARTINI model. The results plotted in Figure 4a) show that the osmotic pressure is underestimated

Figure 4. a) Osmotic pressure with both the original and reparametrized MARTINI CG-FF and b) glucose−glucose radial distribution functions g(r), for the GLYCAM06TIP5P OSMOr14 all atom FF, both MARTINI parameter sets, and the MS-CG interaction potentials. D


Letter


(12) Saxena, A.; García, A. E. Multisite ion model in concentrated solutions of divalent cations (MgCl2 and CaCl2): osmotic pressure calculations. J. Phys. Chem. B 2015, 119, 219−227. (13) Shen, J.-W.; Li, C.; van der Vegt, N. F.; Peter, C. Transferability of coarse grained potentials: Implicit solvent models for hydrated ions. J. Chem. Theory Comput. 2011, 7, 1916−1927. (14) Miller, M. S.; Lay, W. K.; Li, S.; Hacker, W. C.; An, J.; Ren, J.; Elcock, A. H. Reparametrization of Protein Force Field Nonbonded Interactions Guided by Osmotic Coefficient Measurements from Molecular Dynamics Simulations. J. Chem. Theory Comput. 2017, 13, 1812−1826. (15) Lay, W. K.; Miller, M. S.; Elcock, A. H. Optimizing Solute− Solute Interactions in the GLYCAM06 and CHARMM36 Carbohydrate Force Fields Using Osmotic Pressure Measurements. J. Chem. Theory Comput. 2016, 12, 1401−1407. (16) Sauter, J.; Grafmüller, A. Predicting the Chemical Potential and Osmotic Pressure of Polysaccharide Solutions by Molecular Simulations. J. Chem. Theory Comput. 2016, 12, 4375−4384. (17) Lay, W. K.; Miller, M. S.; Elcock, A. H. Reparameterization of Solute? Solute Interactions for Amino Acid−Sugar Systems Using Isopiestic Osmotic Pressure Molecular Dynamics Simulations. J. Chem. Theory Comput. 2017, 13, 1874−1882. (18) Sauter, J.; Grafmüller, A. Solution Properties of Hemicellulose Polysaccharides with Four Common Carbohydrate Force Fields. J. Chem. Theory Comput. 2015, 11, 1765−1774. (19) Izvekov, S.; Voth, G. A. A multiscale coarse-graining method for biomolecular systems. J. Phys. Chem. B 2005, 109, 2469−2473. (20) Sauter, J.; Grafmüller, A. Procedure for Transferable CoarseGrained Models of Aqueous Polysaccharides. J. Chem. Theory Comput. 2017, 13, 223−236. PMID: 27997210. (21) López, C. A.; Rzepiela, A. J.; De Vries, A. H.; Dijkhuizen, L.; Hunenberger, P. H.; Marrink, S. J. Martini coarse-grained force field: extension to carbohydrates. J. Chem. Theory Comput. 2009, 5, 3195− 3210. (22) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845−854. (23) Hess, B.; Bekker, H.; Berendsen, H. J.; Fraaije, J. G. LINCS: a linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463−1472. (24) Miyamoto, S.; Kollman, P. A. SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem. 1992, 13, 952−962. (25) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An NâŃ Ě log (N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089−10092. (26) Yeh, I.-C.; Berkowitz, M. L. Ewald summation for systems with slab geometry. J. Chem. Phys. 1999, 111, 3155−3162. (27) Birdsall, C. K.; Langdon, A. B. Plasma Physics Via Computer Simulation; McGraw-Hill, Inc.: New York City, NY, 1986. (28) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511−519. (29) Hoover, W. G. Canonical dynamics: equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695. (30) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182−7190. (31) Nose, S.; Klein, M. Constant pressure molecular dynamics for molecular systems. Mol. Phys. 1983, 50, 1055−1076. (32) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; De Vries, A. H. The MARTINI force field: coarse grained model for biomolecular simulations. J. Phys. Chem. B 2007, 111, 7812−7824. (33) MSCG-release/MSCGFM 1.7 Manual.txt. https://github.com/ uchicago-voth/MSCG-release (accessed 2017-09-04). (34) Lu, L.; Izvekov, S.; Das, A.; Andersen, H. C.; Voth, G. A. Efficient, regularized, and scalable algorithms for multiscale coarsegraining. J. Chem. Theory Comput. 2010, 6, 954−965.

3. CONCLUSIONS We have demonstrated that the osmotic pressure can be calculated precisely using transferable CG-FFs obtained by using an inhomogeneously regularized coarse-graining procedure. The osmotic pressure data from atomistic simulations for four different carbohydrates can be reproduced almost exactly with an increase in the computational efficiency by a factor of ≈500. Therefore, the method presented here may be used to significantly speed up osmotic pressure calculations and even the FF parametrization process. In addition, the osmotic pressure method can be used directly to tune existing coarsegrained FFs like the carbohydrate MARTINI Model to reproduce the osmotic pressure.

■

AUTHOR INFORMATION

Corresponding Author

*Phone: +49(0)331 567 9619. Fax: +49(0)331 567 9612. Email: [email protected]. ORCID

Andrea Grafmüller: 0000-0002-1671-3158 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The project was partly funded by the Deutsche Forschungsgemeinschaft (GR 3661/2-1). The authors thank Reinhard Lipowsky for helpful discussions.

■

REFERENCES

(1) Ploetz, E. A.; Bentenitis, N.; Smith, P. E. Developing force fields from the microscopic structure of solutions. Fluid Phase Equilib. 2010, 290, 43−47 Proceedings of the 17th Symposium on Thermophysical Properties.. (2) Karunaweera, S.; Gee, M. B.; Weerasinghe, S.; Smith, P. E. Theory and simulation of multicomponent osmotic systems. J. Chem. Theory Comput. 2012, 8, 3493−3503. (3) Satarifard, V.; Kashefolgheta, S.; Vila Verde, A.; Grafmüller, A. Is the Solution Activity Derivative Sufficient to Parametrize Ion−Ion Interactions? Ions for TIP5P Water. J. Chem. Theory Comput. 2017, 13, 2112−2122. (4) Kang, M.; Smith, P. E. A Kirkwood-Buff derived force field for amides. J. Comput. Chem. 2006, 27, 1477−1485. (5) Bentenitis, N.; Cox, N. R.; Smith, P. E. A Kirkwood- Buff Derived Force Field for Thiols, Sulfides, and Disulfides. J. Phys. Chem. B 2009, 113, 12306−12315. (6) McMillan, W. G., Jr; Mayer, J. E. The statistical thermodynamics of multicomponent systems. J. Chem. Phys. 1945, 13, 276−305. (7) Schmalhorst, P. S.; Deluweit, F.; Scherrers, R.; Heisenberg, C.-P.; Sikora, M. Overcoming the Limitations of the MARTINI Force Field in Simulations of Polysaccharides. J. Chem. Theory Comput. 2017, 13, 5039−5053. PMID: 28787166. (8) Luo, Y.; Roux, B. Simulation of osmotic pressure in concentrated aqueous salt solutions. J. Phys. Chem. Lett. 2010, 1, 183−189. (9) Yoo, J.; Aksimentiev, A. Improved parametrization of Li+, Na+, K +, and Mg2+ ions for all-atom molecular dynamics simulations of nucleic acid systems. J. Phys. Chem. Lett. 2012, 3, 45−50. (10) Yoo, J.; Wilson, J.; Aksimentiev, A. Improved model of hydrated calcium ion for molecular dynamics simulations using classical biomolecular force fields. Biopolymers 2016, 105, 752−763. (11) Yoo, J.; Aksimentiev, A. Improved Parameterization of Amine− Carboxylate and Amine−Phosphate Interactions for Molecular Dynamics Simulations Using the CHARMM and AMBER Force Fields. J. Chem. Theory Comput. 2016, 12, 430−443. E


Letter

Journal of Chemical Theory and Computation (35) Uedaira, H.; Uedaira, H. Activity coefficients of aqueous xylose and maltose solutions. Bull. Chem. Soc. Jpn. 1969, 42, 2137−2140. (36) Miyajima, K.; Sawada, M.; Nakagaki, M. Studies on aqueous solutions of saccharides. I. Activity coefficients of monosaccharides in aqueous solutions at 25 C. Bull. Chem. Soc. Jpn. 1983, 56, 1620−1623. (37) Miyajima, K.; Sawada, M.; Nakagaki, M. Studies on aqueous solutions of saccharides. II. Viscosity B-coefficients, apparent molar volumes, and activity coefficients of D-glucose, maltose, and maltotriose in aqueous solutions. Bull. Chem. Soc. Jpn. 1983, 56, 1954−1957.

F


Efficient Osmotic Pressure Calculations Using Coarse-Grained

Recommend Documents