Force Field for Monovalent, Divalent, and Trivalent Cations Developed

Oct 27, 2012 - (15) Only the last 17.5 ps of each window are used to compute the free energy change for the increment in the dispersion interaction. T...
1 downloads 10 Views 363KB Size
Download PDF
Article pubs.acs.org/JPCA

Force Field for Monovalent, Divalent, and Trivalent Cations Developed under the Solvent Boundary Potential Youngdo Won* Department of Chemistry, Hanyang University, Seoul 133-792, Korea S Supporting Information *

ABSTRACT: Presented are parameters for mono-, di-, and trivalent cations compatible with the CHARMM additive force field and the TIP3P water model. Thermodynamic perturbation molecular dynamics simulations were performed for the cations located at the center of a TIP3P water sphere under a solvent boundary potential. A series of perturbations generated free energies of hydration indexed by the two Lennard-Jones parameters, ε and Rmin. Interpolating the experimental free energies of hydration showed that multiple combinations of ε and Rmin values reproduced the free energies of hydration for each ion. To overcome this nonunique parameter problem, the hydration shell model in combination with an empirical scaling parameter was applied to assign values for each ion. Rmin values were then identified via interpolation of the calculated free energies of hydration. The presented parameters are anticipated to be of utility for simulations of ions, including ions complexed to proteins. rule that prescribes εij = (εiεj)1/2 and Rmin,ij = 1/2(Rmin,i + Rmin,j). In this model, only two parameters, εi and Rmin,i, determine the atomic characteristics of an ion. While most popular force fields include molecular building blocks for proteins, lipids, nucleic acids, and sugars, they provide parameters for only a limited set of metal ions.2,3 Previous parameter developments covered mostly alkali metal cations and halide anions.4−8 The additive CHARMM force field includes monovalent Li+, Na+, K+, Rb+, and Cs+ ions and divalent Mg2+, Ca2+, Ba2+, Zn2+ ,and Cd2+ ions.9 Ideally, the force field for ions is developed to reproduce experimentally observed properties of ionic solutions with the use of molecular dynamics (MD) simulations. The primary target property is the free energy of hydration, and other properties that may be considered include the hydration entropy, the radial distribution function of water oxygen atoms around the ion, the lattice constant and the lattice energy of ionic crystals, and osmotic coefficients. Data from ab initio quantum mechanical calculations of monohydrates to water clusters of ions are also considered in balancing the dispersion parameters. Although a multitude of properties have been utilized in the optimization of ion parameters, the choice of the dispersion parameters was not clearly understood.4−8 This work presents a force field for a series of mono-, di-, and trivalent cations for which experimental free energy of hydration was compiled by Marcus.10 Thermodynamic perturbation MD simulations11 are performed with the TIP3P12 explicit water model to generate free energies of hydration as a function of the force field parameters. An interpolation of the experimental

1. INTRODUCTION Molecular simulations allowed for detailed insights into microscopic events taking place in chemical and biological systems. The structural information required to model macromolecular assemblies is often available from public domain databases. The protein data bank (PDB)1 provides atomic coordinates of enzymes, structural proteins, and macromolecular complexes of biological interests. A great number of structures in the PDB contain metal ions, which are an integral part of the molecular structure and often play a key role in their biological activity. A partial list of metal ions found in the PDB entries includes monovalent Li, Na, K, Rb, Cs, Cu, and Au ions, divalent Mg, Ca, Sr, Ba, V, Mn, Fe, Co, Ni, Cu, Zn, Pd, Cs, Pt, Hg, Pb, Eu, and Yb ions, and trivalent Al, Y, La, Ce, Pr, Sm, Eu, Gd, Tb, Ho, Er, Yb, Lu, Cr, Mn, Fe, Co, and Au ions.1 Thus, in order to perform molecular simulations on many proteins, it is necessary to have empirical force field parameters for the ions they contain. In additive empirical force fields, a cation is represented as a charged van der Waals (vdW) atom that interacts with all atoms in the molecular system, with the Coulomb and Lennard-Jones (LJ) potentials normally used for the electrostatic and the dispersion interactions, respectively. The pairwise interaction energy of an atom i of charge qi and an atom j of charge qj is thus 12 ⎡⎛ ⎛ R min , ij ⎞6 ⎤ R min , ij ⎞ ⎢ ⎟⎟ − 2⎜⎜ ⎟⎟ ⎥ Eij = + εij⎢⎜⎜ ⎥ 4πε0 r r ⎝ ij ⎠ ⎦ ⎣⎝ ij ⎠

qiqj

(1)

where ε0 is the permittivity of vacuum. The potential well depth εij and the distance Rmin,ij at which the potential reaches the minimum are obtained by combining the atomic parameters. An example of the combining rule is the Lorentz−Berthelot © 2012 American Chemical Society

Received: September 14, 2012 Revised: October 11, 2012 Published: October 27, 2012 11763

dx.doi.org/10.1021/jp309150r | J. Phys. Chem. A 2012, 116, 11763−11767

The Journal of Physical Chemistry A

Article

Figure 1. Hydration free energy map is shown on the ε−Rmin plain for monovalent ions. Contour lines are drawn from −650 to −250 in every 50 kJ mol−1. The proposed force field values presented in Table 1 are located on the map.

Figure 2. Hydration free energy map is shown on the ε−Rmin plain for divalent ions. Contour lines are drawn from −2200 to −1200 in every 200 kJ mol−1. The proposed force field values presented in Table 1 are located on the map.

free energies of hydration gives a series of ε and Rmin pair. The solvation shell model of ions is adopted to devise a scaling formula for the well depth parameter as a function of the radius and the shell thickness. The final outcome is a force field compatible to the existing one for sodium and potassium ions.13

2. THEORETICAL METHODS A monatomic cation is an atomic species that has the charge q (+1 for monovalent, +2 for divalent, and +3 for trivalent cations considered in this work) for the electrostatic interaction and the Lennard-Jones potential with the parameters ε and Rmin for the exchange-repulsion and dispersion interaction. The simulation system contains the cation placed at the origin of a water sphere of the radius 1.6 nm. The sphere includes 572 TIP3P water molecules equilibrated at 298.15 K under the solvent boundary potential (SBP).13 The hydration free energy is defined as the reversible work required in the process of changing the state from one mole of ideal gas to the 1 m standard solution. Note that the SBP of the simulation includes the contribution due to the phase potential of the air−water boundary.7 The free energy of a cation is calculated in three steps. First, one mole of ideal gas of volume 0.024788 m3 is confined into the volume of 0.001 m3 at 298.15 K and 1 bar. The entropic contribution is ΔG1 = −RT ln(0.001/0.024788) = 7.958 kJ mol−1, where R is the gas constant and T is the absolute temperature. Second, the vdW particle of given Rmin and ε = 0.0 is inserted at the center of the TIP3P water sphere under SBP. Its mass is set to 1 amu, and its position is fixed at the origin. The hydration free energy, ΔG2, of the fictitious uncharged atom is calculated by thermodynamic perturbation with the soft core potential14 where ε is increased to the value of the particle in steps of 41.84 J mol−1. Langevin dynamics at constant temperature and pressure is performed for 20 ps per each window with the integration time step of 1 fs. The Langevin friction coefficient is set to 25 ps−1 for TIP3P oxygen atoms. The O−H and H−H distances of water are fixed using SHAKE.15 Only the last 17.5 ps of each window are used to compute the free energy change for the increment in the dispersion interaction. The recommended nonbonded cutoff options are applied to

Figure 3. Hydration free energy map is shown on the ε−Rmin plain for trivalent ions. Contour lines are drawn from −4400 to −2600 in every 200 kJ mol−1. The proposed force field values presented in Table 1 are located on the map. The points are so congested that the ion is denoted vertically above its position.

the system under SBP: electrostatic switching and LJ switching functions with cut-on, cut-off, cut-nonbond distances of 1.8, 2.0, and 2.2 nm, respectively. Third, the vdW particle of given Rmin and ε is charged in steps of +0.1e to yield cations with +1, +2, and +3 charges in sequence. Eleven charging steps (11 windows ending 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, and 1.0 positive charge, respectively) are employed for each unit charge. The first 11 charging thermodynamics perturbation steps reach at the monovalent cation with the free energy change ΔG31 (i.e., 0 to +1). The next 11 charging steps reach at the divalent cation with ΔG32 (i.e., +1 to +2), and the following 11 steps at the trivalent cation with ΔG33 (i.e., +2 to +3). Each charging step is performed in a thermodynamic perturbation window of 20 ps, of which the first 2.5 ps is used to equilibrate to the new charge state, and the last 17.5 ps is used to calculate the free energy change. The same Langevin dynamics protocol is used 11764

dx.doi.org/10.1021/jp309150r | J. Phys. Chem. A 2012, 116, 11763−11767

The Journal of Physical Chemistry A

Article

Table 1. Compiled Experimental Gibbs Free Energy of Hydration, Marcus Radii, and the Hydration Shell Thickness of Ions with Derived Force Field and Calculated Gibbs Free Energy of Hydrationa ion

ΔGhyda

+

−475 −365 −295 −275 −250 −525 −430 −575 −300 −2395 −1830 −1505 −1380 −1250 −1250 −1825 −1850 −1760 −1840 −1915 −1980 −2010 −1955 −1910 −1850 −1755 −1490 −1960 −1760 −1425 −1375 −1385 −1510

Li Na+ K+ Rb+ Cs+ Cu+ Ag+ Au+ Ti+ Be2+ Mg2+ Ca2+ Sr2+ Ba2+ Ra2+ V2+ Cr2+ Mn2+ Fe2+ Co2+ Ni2+ Cu2+ Zn2+ Pd2+ Ag2+ Cd2+ Sn2+ Pt2+ Hg2+ Pb2+ Sm2+ Eu2+ Yb2+

r 69 102 138 149 170 77 115 137 150 45 72 100 113 136 143 79 82 83 78 75 69 73 75 86 89 95 93 80 102 118 119 117 102

Δr 172 176 74 64 49 156 97 75 63 303 227 171 150 118 109 212 205 203 213 220 233 224 220 197 180 178 183 209 168 143 138 145 168

εb 0.1670 0.2007 0.3600 0.4022 0.4910 0.1852 0.2870 0.3563 0.4067 0.1016 0.1485 0.2016 0.2288 0.2830 0.3019 0.1609 0.1667 0.1685 0.1595 0.1539 0.1434 0.1505 0.1539 0.1742 0.1853 0.1926 0.1879 0.1631 0.2054 0.2394 0.2448 0.2368 0.2054

Rmin/2c

ΔGcalcdd

107.67 160.76 198.81 213.15 233.24 86.20 121.43 63.14 193.07 4.00 110.50 162.06 185.59 209.86 208.99 110.37 105.92 120.37 108.01 92.76 78.03 72.51 81.44 92.71 104.90 119.62 165.80 79.71 118.18 175.80 185.65 184.04 160.90

−469 −365 −291 −276 −250 −528 −432 −572 −300 −2392 −1825 −1503 −1381 −1248 −1247 −1831 −1846 −1758 −1840 −1913 −1973 −2009 −1960 −1895 −1851 −1760 −1490 −1958 −1757 −1423 −1375 −1380 −1516

errore 6 0 4 −1 0 −3 −2 3 0 3 5 2 −1 2 3 −6 4 2 0 2 7 1 −5 15 −1 −5 0 2 3 2 0 5 −6

ion

ΔGhyda

r

Δr

εb

Rmin/2c

ΔGcalcdd

errore

3+

−4525 −3795 −3450 −3145 −3200 −3245 −3280 −3250 −3325 −3360 −3375 −3400 −3425 −3470 −3495 −3515 −3570 −3515 −3205 −3235 −4015 −4220 −4010 −4410 −4265 −4495 −4515 −4350 −3980 −4420 −3970 −3480

53 75 90 105 101 100 98 97 96 95 94 92 91 90 89 88 87 86 104 101 67 64 62 65 65 61 62 67 79 70 89 101

324 262 228 203 207 209 212 214 216 218 220 223 226 228 231 233 235 237 201 207 282 291 296 288 288 299 296 282 253 275 233 205

0.1066 0.1410 0.1656 0.1896 0.1841 0.1823 0.1792 0.1775 0.1757 0.1740 0.1723 0.1693 0.1673 0.1656 0.1636 0.1620 0.1604 0.1588 0.1896 0.1841 0.1285 0.1236 0.1206 0.1252 0.1252 0.1191 0.1206 0.1285 0.1473 0.1330 0.1629 0.1850

22.10 109.69 141.34 177.27 171.15 165.91 161.99 165.72 156.67 152.30 150.50 147.50 144.44 139.07 136.73 134.86 129.52 135.12 170.12 166.97 82.89 61.42 84.31 38.44 56.63 25.81 22.69 46.13 86.63 36.54 86.93 136.82

−4507 −3805 −3451 −3130 −3194 −3243 −3275 −3237 −3327 −3363 −3366 −3402 −3421 −3462 −3495 −3505 −3554 −3498 −3205 −3235 −4007 −4218 −4016 −4403 −4272 −4484 −4498 −4352 −3975 −4412 −3969 −3473

18 −10 −1 15 6 2 5 13 −2 −3 9 −2 4 8 0 10 16 17 0 0 8 2 −6 7 −7 11 17 −2 5 8 1 7

Al Sc3+ Y3+ La3+ Ce3+ Pr3+ Nd3+ Pm3+ Sm3+ Eu3+ Gd3+ Tb3+ Dy3+ Ho3+ Er3+ Tm3+ Yb3+ Lu3+ U3+ Pu3+ Ti3+ V3+ Cr3+ Mn3+ Fe3+ Co3+ Ga3+ Rh3+ In3+ Au3+ Tl3+ Bi3+

a Compiled data are from ref 10, where the hydration model depicts the ion as a sphere with the radius r that is surrounded by the hydration shell with the thickness Δr. ΔGhyd, ΔGcalcd, and ε are in kJ mol−1 and r, Δr, and Rmin are in pm. bε is evaluated in eq 2 based on r and Δr. cRmin is interpolated from the values in Table S5−S7 of the Supporting Information for the given ε value except Be2+, for which Rmin is adjusted to reproduce the experimental free energy of hydration. dΔGcalcd is the free energy of hydration calculated by the thermodynamic perturbation simulation with the force field ε and Rmin for the ion. eError is ΔGcalcd − ΔGhyd in kJ mol−1.

CHARMM force field parameter file, these correspond to Rmin/2 = 0.2 to 2.4 Å and ε = 0.02 to 0.24 kcal mol−1, respectively. The results of 144 simulations for ΔG2, ΔG31(0 to +1), ΔG32(+1 to +2), and ΔG33(+2 to +3) are, respectively, collected in Tables S1−S4 of the Supporting Information. The thermodynamic perturbation windows are combined to yield the hydration free energy ΔGhydration(n) = ΔG1 + ΔG2 + ∑ni = 1ΔG3i for monovalent (n = 1), divalent (n = 2), and trivalent (n = 3) cations, where ΔG1 is the entropic contribution from one mole of ideal gas to the 1 m standard solution. Tables S5−S7 of the Supporting Information list the hydration free energies in a matrix form indexed by the ε and Rmin values. The range of free energies covers −700 to −230 kJ mol−1 for +1 ions, −2300 to −1050 kJ mol−1 for +2 ions, and −4450 to −2440 kJ mol−1 for +3 ions. The hydration free energy data are depicted as the contour map on the ε−Rmin plane in Figures 1−3, respectively. The ordinate spans the ε values, and the abscissa spans the Rmin values used in the free energy simulations. A contour line traces Rmin as a function of ε, which yields the same hydration free energy. As is evident, any

as in the second stage that inserts the vdW particle, and the weighted histogram analysis method is employed to obtain each ΔG3.16 The sum of three components, ΔGhydration(n) = ΔG1 + ΔG2 + ∑ni = 1ΔG3i (n = +1, +2, and +3), is the hydration free energy of the cation with the dispersion parameters, Rmin and ε. The series of free energy calculations scanning the Rmin and ε space yield the table of ΔGhydration(n), which is used to determine the dispersion parameters for a variety of cations with known free energy of hydration. All calculations were performed with the CHARMM program package version c35b1.3

3. RESULTS AND DISCUSSION Thermodynamic perturbation simulations were performed to create the vdW spheres of Rmin and ε and to then subsequently charge it into a mono-, di-, or trivalent cation. A total of 144 pairs of the dispersion parameters were selected to cover the Rmin and ε space: Rmin spanning from 40 to 480 pm with the interval of 40 pm and ε spanning from 0.08368 kJ mol−1 to 1.00416 kJ mol−1 with the interval of 0.08368 kJ mol−1. In the 11765

dx.doi.org/10.1021/jp309150r | J. Phys. Chem. A 2012, 116, 11763−11767

The Journal of Physical Chemistry A

Article

Rmin and ε pair on the line for the given free energy of hydration is a feasible candidate for the force field of the ion. The experimental hydration free energy of an ion is interpolated from the values of Tables S5−S7 of the Supporting Information to determine the Rmin and ε pairs for the ion. A cubic spline is used in the interpolation procedure. As there is no unique ε,Rmin pair that corresponds to a free energy hydration, a series of Rmin values are obtained for the range of ε values. Table S8 of the Supporting Information lists the interpolated Rmin values for monovalent ions, Table S9 for divalent ions, and Table S10 for trivalent ions. Thus, the present analysis yields nonunique sets of parameters for the free energies of hydration, consistent with the previously described parameter correlation problem.17 Accordigly, it is necessary to identify a unique set of parameters for the final ion models. Marcus compiled experimental values of the hydration free energies of a variety of ions. On the basis of this compilation, he suggested a simple model that accounted for the thermodynamics of hydration of ions. The ion is modeled as a charged sphere of radius r that is surrounded by a spherical hydration shell of thickness Δr. The Marcus model parameters are presented in Table 1 along with the hydration free energies. To overcome the parameter correlation problem leading to the nonunique parameters encountered in the present study, the hydration shell model of ions is adopted to relate the dispersion parameter ε to the size r and the shell thickness Δr. The hydration shell basically identifies the layer of water molecules more tightly bound to the ion than bulk ones. The size of the hydration shell is roughly proportional to the volume around the ion accessible by the tightly bound water molecules. In the present study, a scaling formula is empirically devised to relate ε to r and Δr. The coefficient is adjusted to yield ε values close to the well-optimized CHARMM force field values for sodium and potassium ions; 0.196 and 0.364 kJ mol−1 respectively.13

ε = 0.063 r /Δr

variety of cations. The resultant free energies of hydration were organized by the dispersion parameters. For the ion with the hydration energy known experimentally, the force field may be determined on the table by interpolation. Unfortunately, due to the parameter correlation problem, no unique pair of ε and Rmin is located on the free energy table for any of the ions. The interpolation of the table values gives Rmin as a function of ε. For each ion with its experimental hydration free energy, the ranges of ε and Rmin pairs are reported in Tables S8−S10 of the Supporting Information. To overcome the presence of the nonunique parameters, the hydration shell model of Marcus was adopted. By relating the hydration of an ion as a sphere with the hydration shell of a defined thickness, ε is scaled to a function of the ratio of the radius and the shell thickness as in eq 2. The formulation is shown to yield close agreement with the published ε values of both sodium and potassium ions in the CHARMM force field. The rather arbitrary scaling aims at compatibility to the existing force field and suggests a set of force fields for the wide spectrum of ions listed in Table 1. The presented force field should be useful for molecular models solvated with TIP3P water.

■

ASSOCIATED CONTENT

S Supporting Information *

Hydration free energy data (Tables S1 to S7) and force fields for mono-, di-, and trivalent cations (Tables S8 to S10). This material is available free of charge via the Internet at http:// pubs.acs.org.

■ ■

AUTHOR INFORMATION

Notes

The authors declare no competing financial interest.

ACKNOWLEDGMENTS I thank Alexander D. MacKerell Jr. for helpful discussions and comments. This work was supported by Research Institute for Natural Sciences, Hanyang University.

(2)

The formula based on the hydration shell model yields ε values of 0.2007 kJ mol−1 for the sodium ion and 0.3600 kJ mol−1 for the potassium ion. The hydration free energy and the ε values based on eq 2 are used to interpolate the data in Tables S5−S7 of the Supporting Information using a cubic spline to obtain the Rmin value for each ion that yields the experimental free energy of hydration. Table 1 presents the resulting LJ parameters that are compatible with the TIP3P explicit water solvent. Note that the hydration free energy of Be2+ is off the range and that its final Rmin value was adjusted to reproduce the experimental free energy of hydration. The same free energy simulation protocol was applied to each ion with the suggested force field. The calculated free energies of hydration are presented in Table 1 with errors from the experimental values. The difference between the experimental and the calculated free energy is averaged to −1 kJ mol−1 with the standard deviation of 3 kJ mol−1 for monovalent cations. The differences are averaged to −1 kJ mol−1 with the standard deviation of 5 kJ mol−1 for divalent cations and −5 kJ mol−1 with the standard deviation of 7 kJ mol−1 for trivalent cations. All simulations reproduce the experimental free energy with less than 1% error.

■

REFERENCES

(1) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. Nucleic Acids Res. 2000, 28, 235. (2) Case, D. A.; Cheatham, T. E., III; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M., Jr.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. Comput. Chem. 2005, 26, 1668. (3) Brooks, B. R.; Brooks, C. L., III; MacKerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. J. Comput. Chem. 2009, 30, 1545. (4) Joung, I. S.; Cheatham, T. E., III J. Phys. Chem. B 2008, 112, 9020. (5) Jensen, K. P.; Jorgensen, W. L. J. Chem. Theory Comput. 2006, 2, 1499. (6) Aqvist, J. J. Phys. Chem. 1990, 94, 8021. (7) Lamoureux, G.; Roux, B. J. Phys. Chem. B 2006, 110, 3308. (8) Horinek, D.; Mamatkulov, S. I.; Netz, R. R. J. Chem. Phys. 2009, 130, 124507. (9) Yu, H.; Whitfield, T. W.; Harder, E.; Lamoureux, G.; Vorobyov, I.; Anisimov, V. M.; MacKerell, A. D., Jr.; Roux, B. J. Chem. Theory Comput. 2010, 6, 774. (10) Marcus, Y. Biophys. Chem. 1994, 51, 111. (11) Mezei, M.; Beveridge, D. L. Ann. N. Y. Acad. Sci. 1986, 482.

4. CONCLUSIONS A minimal set of thermodynamic perturbation simulations was performed to span the LJ parameter space of ε and Rmin for a 11766

dx.doi.org/10.1021/jp309150r | J. Phys. Chem. A 2012, 116, 11763−11767

The Journal of Physical Chemistry A

Article

(12) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (13) Beglov, D.; Roux, B. J. Chem. Phys. 1994, 100, 9050. (14) Zacharias, M.; Straatsma, T. P.; McCammon, J. A. J. Chem. Phys. 1994, 100, 9025. (15) Ryckaert, W. E.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (16) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. J. Comput. Chem. 1992, 13, 1011. (17) MacKerell, A. D., Jr. J. Comput. Chem. 2004, 25, 1584.

11767

dx.doi.org/10.1021/jp309150r | J. Phys. Chem. A 2012, 116, 11763−11767