Refined Dummy Atom Model of Mg2+

Article pubs.acs.org/jcim

Refined Dummy Atom Model of Mg2+ by Simple Parameter Screening Strategy with Revised Experimental Solvation Free Energy Yang Jiang, Haiyang Zhang,† Wei Feng, and Tianwei Tan* Beijing Key Lab of Bioprocess, College of Life Science and Technology, Beijing University of Chemical Technology, Beijing 100029, China S Supporting Information *

ABSTRACT: Metal ions play an important role in the catalysis of metalloenzymes. To investigate metalloenzymes via molecular modeling, a set of accurate force field parameters for metal ions is highly imperative. To extend its application range and improve the performance, the dummy atom model of metal ions was refined through a simple parameter screening strategy using the Mg2+ ion as an example. Using the AMBER ff03 force field with the TIP3P model, the refined model accurately reproduced the experimental geometric and thermodynamic properties of Mg2+. Compared with point charge models and previous dummy atom models, the refined dummy atom model yields an enhanced performance for producing reliable ATP/GTP-Mg2+-protein conformations in three metalloenzyme systems with single or double metal centers. Similar to other unbounded models, the refined model failed to reproduce the Mg−Mg distance and favored a monodentate binding of carboxylate groups, and these drawbacks needed to be considered with care. The outperformance of the refined model is mainly attributed to the use of a revised (more accurate) experimental solvation free energy and a suitable free energy correction protocol. This work provides a parameter screening strategy that can be readily applied to refine the dummy atom models for metal ions.

1. INTRODUCTION Metalloenzyme is a kind of enzyme that consists of one or more metal ions as its cofactor. Over one-third of the known enzymes are metalloenzymes.1 The metal ions in metalloenzymes have multiple functions and play an important role in the catalytic process.2 They can be a part of the active center and stabilize the conformation of the protein or act as a bridge to combine the enzyme with its substrate. They can also transport protons, electrons, or chemical groups in the catalytic process.2 In recent years, molecular dynamics (MD) simulation has been a useful tool in studying the conformation and activation mechanism of metalloenzymes.3−6 To simulate the metalloenzyme, one of the most important preparations is to get a set of suitable force field parameters of metal ions. Currently, the commonly used biomolecular simulation force field, like the AMBER force field,7 the CHARMM force field,8,9 and others, are of static charge and not able to account for multipolar electrostatics.10 In order to simulate metal ions using these classical force fields, three strategies for metal ion modeling have been developed: the covalent bonded model,3,4 the point charge model,6 and the nonbonded dummy atom model.11−14 The covalent bonded model has predefined covalent bonds or distance restrictions between metal and its ligands. The force field parameters of the metal−ligand complex are usually generated and optimized by quantum mechanics (QM) calculations. It can produce highly accurate geometries of the target metal centers. However, when the composition of the © XXXX American Chemical Society

metal site is changed, the covalent bonded model parameters have to be reoptimized to match the new metal site. Moreover, it does not allow the ligands exchange during the simulation time. The point charge model, also named “nonbonded softsphere model”, is the simplest one. The metal ion is described by a point charge surrounded by van der Waals (vdW) potentials. Only the parameters of vdW potentials need to be optimized. The point charge model is convenient to use and allows the ligands exchange. But when the simulation system has a complex metal center (two or more metal ions in one active site) or when transition metals are present, this model will be inadequate.12,15 The nonbonded dummy model makes a good balance between both previous models. It was originally presented by Åqvist and Warshel15 and improved by Pang,11 Oelschlaeger et al.,12 Saxena et al.,13 and Duarte et al.14 In this model, the metal ion has a metal core (MC) surrounded by six dummy atoms to form an octahedron (Figure 1). Each dummy atom possesses a charge of +δ, a small part of the total mass and a negligible vdW potential. The metal core possesses a charge of n-6δ (to ensure the total charge is +n), a major part of the total mass and a vdW potential. There are covalent bonds between metal core and dummy atoms, but no covalent bond between dummy atoms and ligands. Thus, the nonbonded dummy model allows the Received: May 15, 2015

A

DOI: 10.1021/acs.jcim.5b00286 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

the advantages of our refined model and then discuss the effect of different force fields and metal models on the simulation accuracy.

2. MATERIALS AND METHODS 2.1. Mg2+ Force Field Parameters Used in This Work. To refine the dummy model parameters of Mg2+, we used Duarte’s parameters as the initial values (shown in Table S1). The values of force constants (Kb for bonds and Kθ for angles), equilibrium bond length (b0) and angle (θ0), and mass and charge of Duarte’s model were still used in our model. The vdW parameters (Ai and Bi) were replaced and refined to be the vdW radius (Ri) and the 12−6 potential well depth (εi). The nonbonding interaction potential function of two atoms i and j in the AMBER ff03 force field has the following form:7

Figure 1. Geometry of nonbonded dummy model The metal center (MC) with a charge of n-6δ (−1 for Mg2+) is covalently bonded by six dummy atoms (D) with a charge of +δ (+0.5 for Mg2+) to form an octahedron. The total charge is +n.

ligands exchange and has a low computational cost. The dummy model of Zn2+ developed by Pang showed an excellent performance in simulating a phosphotriesterase with two Zn2+.16 Recently, Duarte et al. optimized dummy model parameters of seven metal ions (Mg2+, Mn2+, Fe2+, Co2+, Ni2+, Zn2+, Ca2+) using the OPLS-AA force field, and all the parameters were checked to be adaptable in simulating the metal substitution of metalloenzyme.14 The magnesium ion is of particular importance for nucleic acid systems and can stabilize the activated conformation of active sites in some enzymes and play a direct role in the catalysis.10 The dummy model of Mg2+ was first developed by Oelschlaeger et al.12 But Lu et al.6 reported that this parameter was inadaptable for simulating the multimetal center of kinase GSK3β and led to a serious conformational change. This parameter was then further refined by Duarte et al.14 in the OPLS-AA force field. Another Mg2+ dummy model was provided by Saxena et al.13,17 with zero charge of the metal core and small force constants of the bonds and angles. Although Saxena’s model was optimized under the AMBER force field, the simulation of Mg2+ in water showed a serious error in our test (discussed in detail in Section 1 of the Supporting Information). Moreover, using Duarte’s Mg2+ parameter in the AMBER force field may also become problematic (see the Discussion section of the main text). Above all, all the previous parameters of Mg2+ were optimized to reproduce the experimental solvation free energies that have been considered to be inaccurate (discussed in section 2.2). These existing drawbacks motivated us to refine the parameter of Duarte’s dummy model of Mg2+ to reproduce a revised experimental solvation free energy using a suitable free energy correction protocol. A simple parameter refining strategy is provided in the present work. Our refined model, optimized for the AMBER ff03 force field and the TIP3P water model, can accurately reproduce the experimental value of the Mg2+−O distance (d(Mg−O)), coordination number (CN), and the solvation free energy (ΔGexp sol ) in water. Finally, we demonstrate

Uijnonbond(rij) = UijvdW(rij) + Uijele(rij) ⎡⎛ ⎞12 ⎛ R ij ⎞6 ⎤ qiqj R ij ⎢ = εij⎢⎜⎜ ⎟⎟ − 2⎜⎜ ⎟⎟ ⎥⎥ + k r rij ⎝ rij ⎠ ⎦ ⎣⎝ ij ⎠

(1)

where rij is the distance between atom i and j in angstroms, qi is the charge of atom i in e, and k is a constant equal to 332. Rij and εij satisfy the Lorentz/Berthelot mixing rules: ⎧ 1 ⎪ R ij = (R ii + R jj) = R i + R j 2 ⎨ ⎪ε = ε ε = ε ε ii jj i j ⎩ ij

(2)

Two sets of point charge model parameters (from Åqvist et al.18 and Allnér et al.,19 respectively) and the Duarte’s dummy model parameters (shown in Tables 1 and S1) were used to compare with our refined dummy model parameters. To use Duarte’s vdW parameters in the AMBER force field, one method is modifying the frcmod file of metal ions with the corresponding R and ε. The nonbonding interaction potential function of OPLS-AA force field used by Duarte et al. has the following form:14,20 Uijnonbond(rij) = VijvdW(rij) + Uijele(rij) =

Aij rij12

−

Bij rij

6

+k

qiqj rij (3)

where the parameters Aij and Bij satisfy the geometric mixing rules

Table 1. Mg2+ Parameters Examined in This Work vdW parameters models Åqvist Allnér Duartea this workb

MC D MC D

Ri (Å)

εi (kcal/mol)

0.7926 1.5545

0.8947 0.00295

1.0360 1.3882

0.0660 1 × 10−8

Ai (kcal1/2 Å6/mol1/2)

Bi (kcal1/2 Å3/mol1/2)

63 0.05 20.3287 0.0458

9 0 3.2319 0.0030

a

Details of Duarte’s parameters are shown in Table S1 (SI). MC stands for metal center and D for dummy atoms. bRefined parameters in this work are in the form of Ri and εi; the values for Ai and Bi are computed using eq 5 and given here for comparison with Duarte’s model. B


Article

Journal of Chemical Information and Modeling Table 2. Calculated Geometry and Solvation Free Energy (kcal/mol) for Mg2+ Dummy Modelsa d(Mg−O) (Å) RMSD (Å)f RI (Å)g ΔGT1 M0→M2+ ΔGsoft‑core appearing ΔGPBC corr ΔGϕcorr ΔGcal sol

Duarte’s model

first eligible model of Iteration Ib

second eligible model of Iteration Ic

eligible model of Iteration IId

exp

2.126 0.036 1.97 −421.52 0.18 1.10 −35.04 −455.28

2.097 0.007 1.92 −430.70 0.99 1.06 −34.85 −463.50

2.100 0.010 1.92 −430.52 0.99 1.06 −34.85 −463.32

2.090 5 × 10−4 1.92 −433.81 1.24 1.06 −34.85 −466.36

2.0921

e

−468.2

The eligible models were obtained from the parameter screening strategy proposed in this work. bThe iteration gives RMC = 1.000 and εMC = 0.126. c RMC = 1.000 and εMC = 0.144. dIteration II gives only one model: RMC = 1.036 and εMC = 0.066. eThe average of six Mg−O distances was calculated from the last 10 ns trajectory; fRoot mean square deviation (RMSD) of d(Mg−O) from the experimental value. gIonic radius was estimated using the radial ion-oxygen distribution function.40 a

⎧A = A A = A A ii jj i j ⎪ ij ⎨ ⎪ Bij = Bii Bjj = Bi Bj ⎩

have a wrong direction in his work22 (ΔGo→* was added to ΔGosol). Thus, for the absolute solvation free energy provided by Marcus (ΔGMarcus * ), the conventional ionic solvation free energy ΔGconv sol can be derived by

(4)

In fact, the forms of eqs 1 and 3 can be converted to each other by the following formulas: ⎧ ⎛ 2Aij ⎞1/6 ⎪ ⎜ ⎟ R = ⎪ ij ⎜ B ⎟ ⎪ ⎝ ij ⎠ ⎨ ⎪ B2 ⎪ εij = ij . ⎪ 4Aij ⎩

conv * ΔGsol = ΔGMarcus − 1.9 + 252.3z

(6)

where ΔGMarcus * is under the 1 M solution standard state and z is the net charge of the ion. Then, our revised Mg2+ solvation free energy can be derived by using the accurate proton hydration energy: exp conv ΔGsol = ΔGsol − 265.9z − 1.9

(7)

The revised experimental solvation free energy for Mg2+ (−468.2 kcal/mol) was given in Table 2. Besides the Marcus’s solvation free energy, there are another set of experimental solvation free energies reported, such as those provided by Rosseinsky27 and Noyes.28 The latter two sets have very similar values (−455.5 kcal/mol for Rosseinsky and −454.2 kcal/mol for Noyes), which have some deviations from Marcus’s (−437.4 kcal/mol). Considering that the absolute solvation free energy from Noyes was derived by using a gas phase H2 hydration process, which is different from the proton hydration process used by Rosseinsky and Marcus, we only compared the sets from Rosseinsky and Marcus. By eliminating the effect of different proton hydration energies, i.e. by obtaining the conventional solvation free energies ΔGconv sol , we found a remarkable consistency between these two conventional solvation free energies of Mg2+ (both are 65.5 kcal/mol). Therefore, it is reliable that using Marcus’s free energy set to derive our revised experimental solvation free energy for Mg2+. 2.3. Parameter Refining Strategy. The first step was optimizing the vdW parameters of dummy atoms under the AMBER ff03 force field. Because of the negligible vdW potential of dummy atoms, the 12−6 potential well depth of dummy atoms εD was set to be close to zero (1 × 10−8 kcal/ mol). The vdW radius RD was optimized according to the following optimization problem:

(5)

Note that, eq 5 can be derived only if the values of Aij and Bij do not equal zero. Unfortunately, the Bi for dummy atoms in Duarte’s model are all zero (see Table 1), thus preventing the direct use of Duarte’s vdW parameters in the AMBER force field by modifying its frcmod file. Another method to use Duarte’s vdW parameters in AMBER force field is modifying the topology file (according to Saxena’s work13): Use any Ri and εi of dummy atoms to generate the Amber topology file; Then use Duarte’s vdW parameters to calculate the right Aij and Bij of each vdW interaction with metal ion (using the mixture rule in eq 4); Finally modify the wrong parameters in the topology file with the right values. The topology file needs to be modified every time a new simulation system is built. Therefore, this method is not as convenient as modifying the frcmod file. But it is the only feasible way to use Duarte’s vdW parameters in the AMBER force field. 2.2. Experimental Values Used in This Work. The experimental value of d(Mg−O) was obtained from Marcus’s work,21 in which d(Mg−O) was derived by averaging the distances measured by X-ray diffraction on different Mg2+ aqueous solutions, and equals 2.09 Å. ΔGexp sol was also obtained from Marcus’s work (−437.4 kcal/mol)22 and revised by the more accurate experimental value of the proton hydration free energy provided by many reports.23−25 The accurate proton hydration energy (−265.9 kcal/mol) was obtained by using the cluster pair approximation, which is considered to be more reliable and accurate than the extra-thermodynamic assumptions “TATB hypothesis” provided by Marcus (about −252.3 kcal/mol)26 and has been corrected by converting the 1-atm gas phase standard state ΔGosol to the 1 M solution standard state ΔG*sol (by subtracting ΔGo→* = 1.9 kcal/mol).25 Note that, Marcus’s solvation free energies also have been corrected to the 1 M solution standard state, but the correction seems to

min ∑ [VvdW(r , AD , BD) − UvdW(r , RD , εD)]2 RD

r

s.t. RD > 0

(8)

where VvdW(r, AD, BD) is the vdW potential of the dummy atom interacting with the oxygen atom of TIP3P water calculated using eq 3; UvdW(r, RD, εD) is the same vdW potential calculated using eq 1; and r is the distance between the dummy C


Article


Figure 2. Iterative scheme used to refine the vdW parameters of the metal core. The procedure starts from defining an initial hunting zone according to Duarte’s parameters. Step 1 produces six Mg−O distances all close to the experimental value (2.09 Å) that are retained after energy minimizations. Step 2 then produces an averaged Mg−O distance close to 2.09 Å after 10 ns MD simulations. Step 3 reproduces the experimental solvation free energy (ΔGsol) of Mg2+. If the relative error of ΔGsol is less than 1%, the screening process is terminated; if not, a narrower hunting zone according to the eligible parameters in step 2 is defined and a new iteration starts. The threshold values (T1 and T2) decrease during the iteration process according to an increasing accuracy requirement. Details about the parameter screening process are given in section 2.3.

atom and the oxygen atom (here, r ∈ [0.5, 3.0]). This optimization problem was solved by MATLAB.29 The optimized parameters of the dummy atom were then used to refine the metal core parameters. A simple parameter screening strategy was used and shown in Figure 2. At Step 1, a set of parameters was used to define the hunting zone of parameter refinement. The Duarte’s vdW parameters (AMC and BMC in Table S1) were transformed to RMC and εMC (here, RMC = 1.0736 Å and εMC = 0.4132 kcal/mol) by eq 5. Thus, the initial hunting zone was defined to be RMC ∈ (0, 2.0] and εMC ∈ (0, 0.9]. The hunting zone was then divided into 2500 discrete points (each parameter had 50 discrete points) to represent different parameter combinations. For each combination, a Mg2+ dummy model was parametrized, and the [Mg(H2O)6]2+ surrounded by a TIP3P water box (15 Å left to the ion) was then built. A 1000-step energy minimization was performed simultaneously for each system using Amber11 suite30−32 to get the distances between Mg2+ (the metal core MC of the dummy atom model) and coordinated oxygen atoms (d(Mg−O)). If all the six values of d(Mg−O) were close to the experimental value 2.09 Å21 with the threshold T1, i.e. |d(Mg− O) − 2.09| < T1, this set of parameters was retained and the strategy proceeded to Step 2. At Step 2, 10 ns MD simulations were performed simultaneously for the retained models to check the consistency of the calculated and the experimental d(Mg−O) values, i.e. the average value of the six distances (d(Mg−O)avg) was close to 2.09 Å with the threshold T2. Note that, the threshold values above (T1 and T2 in Figure 2) were not constant. They would decrease during the iteration process because of the increasing accuracy requirement. All the eligible

sets of parameters went to Step 3. At Step 3, the solvation free energy of each Mg 2+ model was calculated by the Thermodynamic Integration (TI) method. If the ΔGcal sol was cal exp close to the experimental value ΔGexp sol , i.e. |(ΔGsol − ΔGsol )/ | < 1%, this set of parameters was finally retained. ΔGexp sol Otherwise, a new hunting zone would be defined around the final sets of parameters with a narrower range than the previous one, and a new iteration would then start. The set of exp parameters that can produce the closest ΔGcal sol to ΔGsol was selected to be the final optimized parameter combination. 2.4. General Simulation Protocol. All simulations were performed using the AMBER11 suite with the all-atom ff03 force field.33 The solutes were embedded in a periodic TIP3P water box. Several ions were added to neutralize the simulation systems, except the systems in the parameter refining and the free energy calculation. An energy minimization was performed to eliminate improper contacts in each system. After energy minimization, each system was heated gradually from 0 to 300 K. During the heating steps, position restraints were imposed on the solutes with a force constant of 10.0 kcal·mol−1·Å−2. Each system was then equilibrated at constant temperature (300 K) and pressure (1 bar) conditions via the Langevin dynamics (the collision frequency is 2.0 ps−1), with a coupling constant of 0.2 ps for both parameters. All the production simulations were performed at 300 K and 1 bar using a time step of 2 fs. Electrostatic interactions were calculated using the particlemesh-Ewald (PME) algorithm. The cutoff distance for van der Waals interactions was 10.0 Å. The SHAKE algorithm34 was applied to the bonds involving hydrogen. D


Article


Figure 3. Metal centers in the crystal structures of the three proteins studied in this work. (a) Kinase GSK3β (1PYX); (b) yeast glutathione synthase (1M0W); and (c) GTPases OsRac1 (4U5X). Proteins are shown in transparent white and Mg2+ ions with pink spheres. Cofactor analogues (ANP/ GNP), water molecules, and residues contacting with Mg2+ are shown in “licorice” model. The interactions between Mg2+ and its ligand are represented by red dashed lines. The residues that form these hydrogen bonds (HBs) are marked. For 1PYX, HBs between Lys85 and ANP are also marked by red dashed lines.

The ⟨dU/dλ⟩ plot for our final refined model is shown in Figure S1b and the other systems show similar landscapes. The correction for the raw ionic charging free energies has been discussed and provided by many reports.37−40 There are two different correction terms should be added to the calculated solvation free energy. The first term ΔGPBC corr is the correction for using the periodic boundary condition (PBC) with the PME approach. PME is one of the most popular lattice-sum (LS) methods. For the PME approach, the ion−ion self-potential is already included in the instantaneous electrostatic potential of the ion.38 This is different from the general lattice-sum method. So the ion−ion self-potential term should be eliminated from the general finite-size correction. Thus, the correction for using PBC with PME approach is given by

The root-mean-square deviations (RMSDs) were calculated using a least-squares fit and starting from the crystal structure. The convergence of protein simulation systems were checked by the RMSDs of protein backbone atoms (C, CA, and N). The fluctuation of the metal ions was checked by the root-meansquare fluctuations (RMSFs) of Mg2+ (or the metal core in dummy model). The geometry of the metal site was checked using the distance between the two Mg2+ ions (d(Mg−Mg)) in the double-metal center, Mg-ligands distances and coordination numbers (CNs). The CN of Mg2+ was calculated by counting the oxygen atoms around Mg2+ within the distance of 2.5 Å. The conformation of the cofactor-protein system was checked using the hydrogen bond occupancy and the RMSF of the ATP and GTP. The hydrogen bond occupancy was calculated by counting the frames in which the distance between the donor and the acceptor was less than 3.5 Å and the angle of donor− H−acceptor was bigger than 120° and then was divided by the total number of frames. All the above analyses were calculated by Ptraj and Cpptraj module of AmberTools1.5. The visualization of conformations was generated by VMD software.35 2.5. Free Energy Calculations. The solvation free energy ΔGsol of Mg2+ was calculated using TI method in two steps. The first step was perturbing nothing to Mg atom (with no charge) in water solvent using the soft-core potential (“appearing” part) in Amber11 suite.36 ⟨dVsoft‑core/dλ⟩ was calculated from 14 independent simulations at λ = 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, and 0.95. The parameters of soft-core potential were set to default. Each simulation included 1 ns energy equilibration and 2 ns production. ⟨dVsoft‑core/dλ⟩ was obtained from the production simulation. The ⟨dVsoft‑core/dλ⟩ values at λ = 0 and 1 were extrapolated by using the cubic spline function. The free energy change ΔGsoft‑core appearing was then calculated using the trapezoidal integral rule. Both the extrapolation and integral were calculated using MATLAB. The ⟨dVsoft‑core/dλ⟩ plot for our final refined model is shown in Figure S1a and the other systems show similar landscapes. The second step was charging the atom Mg from 0 e to +2 e by a parameter λ. ⟨dU/dλ⟩ was calculated from 23 independent simulations at λ = 0, 0.01, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 0.99, and 1.00. Each simulation included 1000-step energy minimization, 60 ps heating, 1 ns energy equilibration, and 1 ns production. ⟨dU/dλ⟩ was obtained from the production T1 simulation. The free energy change ΔGM 0 →M2+ was then calculated using the trapezoidal integral rule by MATLAB.

q2 ξLS q2 1 − εs−1 · + · · L 8πε0εs L 8πε0 ⎡ 4π ⎛ R ⎞2 16π 2 ⎛ R ⎞5⎤ I⎟ ⎜ ⎢ ⎜ I⎟ − ⎥ 45 ⎝ L ⎠ ⎦ ⎣ 3 ⎝L⎠

PBC ΔGcorr =−

(9)

where q is the ionic charge (+2 e); εs is the solvent permittivity (78.0 for water); L is the box boundary length (31.9 Å for all the TI simulation systems); RI is the ionic radius, which can be estimated by using the radial ion-oxygen distribution function;40 ξLS is a constant and approximately equals −2.837 297.38−40 The second term ΔGϕcorr is the correction for the error in the potential at the ionic site due to its evaluation using an improper summation scheme, and due to the possible presence of an explicit liquid-vacuum interface.39 In this study, we used the formula provided by Kastenholz et al.39 to estimate this term: 3 ⎛ 4πRI 3 ⎞ ργ 1 4πRI ϕ ⎟ q ΔGcorr = −q⎜1 − − 4.18 3L3 3L3 ⎠ 6ε0 ⎝

(70.3 − 109/RI)

(10)

where ρ is the solvent number density (0.033 Å−3 for all the TI simulation systems); γ is the quadrupole-moment trace of the solvent molecule39 (0.77 e·Å2 for the TIP3P water molecule in Amber). Note that there is no need to add the correction for the truncated vdW interaction, because Amber has already added this correction into its calculation protocol (see the keyword “vdwmeth” in Amber’s manual, http://ambermd.org/doc12/ Amber14.pdf). Note that there is also no need to add the E


Article

Journal of Chemical Information and Modeling correction for standard state ΔG0→* here, because our revised experimental solvation free energy has been already corrected into the 1 M solution standard state. Therefore, the solvation free energy of Mg2+ in water is finally calculated by

were simulated in TIP3P water. Each simulation included 15000-step energy minimization, 60 ps heating, 0.5 ns energy equilibration, and 20 ns production. 2.7. Evaluating the Influence of Using Different Models in Different Force Fields. To investigate the reasons why different parametrization models have different simulation performances, a simple testing system containing only one Mg2+ ion (M) with one coordinated oxygen atom (O) was built (see Figure S5a). The metal ion was built by the Åqvist’s model, Duarte’s model, and our refined model, respectively; and the atom O was the oxygen atom of TIP3P water in AMBER force field. The comparison between different vdW and electrostatic potentials of these different models was then performed (see details in section 2.1 of the Supporting Information). To evaluate the influence of different force fields, we compared the total nonbonding potentials of atom O in OPLSnonbond AA force field (Unonbond OPLS ) and AMBER force field (UAMBER ) with the metal ion built in Duarte’s model (see details in section 2.2 of the Supporting Information). The relative error E of using Duarte’s model in AMBER force field can be defined as the following equation:

cal soft‐core PBC ϕ TI ΔGsol = ΔGappearing + ΔG M + ΔGcorr + ΔGcorr 0 → M2 +

(11)

2.6. Parameter Verification Procedure. 2.6.1. Comparing QM Optimized Geometry with MM Minimized Geometry. Two simplified metal center structures were used to perform the comparison between the QM optimized structures and the MM minimized structures. The first metal center structure was obtained from the yeast glutathione synthase (GS, PDBID: 1M0W41) including two Mg2+ ions and their ligand groups. The second metal center structure was obtained from the small GTPases OsRac1 (PDBID: 4U5X42) including one Mg2+ ions and its ligand groups. All the detials of the calculation method can be found in the previous work of Hu and Ryde.43 The comparison included the correlation coefficients (r2) between the Hessian matrix elements of QM and MM structures and the RMSd values between bond lengths, angles, dihedrals, and coordinates of QM and MM structures. The QM optimized structures were treated as the standards. 2.6.2. Comparing MD Simulated Structure with Crystal Structure. Three different metalloenzymes that contain Mg2+ ions were used to test our refined dummy model parameters. The first protein was the kinase GSK3β (PDBID: 1PYX44), whose active site contains two Mg2+ ions binding with AMPPNP (shown in Figure 3a). According to the study of Lu et al.,6 the first chain of GSK3β was retained and the missing residues were added. AMP-PNP was manually changed to ATP. To check if our refined model can solve the problem found by Lu et al.,6 the Mg2+ ions were parametrized by our refined dummy model and Duarte’s dummy model, respectively. The second protein was the yeast glutathione synthase (GS, PDBID: 1M0W41), whose active site also contains two Mg2+ ions binding with AMP-PNP (shown in Figure 3b). The missing residues were added and AMP-PNP was manually changed to ATP. To check if our refined model has advantages in simulating a double-metal center, the Mg2+ ions were parametrized by our refined dummy model, Duarte’s dummy model, Åqvist’s point charge model and Allnér’s point charge model, respectively. In addition, a restrained metal center with Åqvist’s point charge model was built and simulated to be compared with the performances of our dummy atom model. The force constants of the coordinate bonds were obtained from the frequency calculation using Gaussian 03.45 The harmonic restraints, not real bonds, were used to restrict the bond lengths to vibrate around the lengths obtained from the crystal structure. All the partial charges of the Mg2+ ions and the ligand oxygen atoms were replaced by the RESP charges derived from the QM optimization with B3LYP at 6-31G* level. The last protein was the plant Rho-type small GTPases OsRac1 (PDBID: 4U5X42), whose active site only contains one Mg2+ ion binding with GMP-PNP (shown in Figure 3c). The missing residues were added, and GMP-PNP was manually changed to GTP. To check if our model is applicable in the single-metal center, the Mg2+ ions were parametrized by our refined dummy model, Duarte’s dummy model, Åqvist’s point charge model and Allnér’s point charge model, respectively. The force field parameters of ATP and GTP were obtained from Meagher’s work.46 All these protein−cofactor complexes

E=

nonbond nonbond UOPLS − UAMBER nonbond UOPLS

= 1−

nonbond UAMBER nonbond UOPLS

(12)

2+

Here, four ligands of Mg (TIP3P-O, Glu-OE1, Thr-OG1 and Asp-OD1; parameters are shown in Table S3) that are common in metalloenzymes were selected to calculate the relative error E. All the above calculations were performed in MATLAB.

3. RESULTS 3.1. Refined Dummy Model Parameters of Mg2+. The refined dummy model parameters of Mg2+ (shown in Table 1) were generated by our parameter refining strategy described in section 2.3. By solving the optimization problem (eq 8), the parameters of the dummy atoms were refined (Ri = 1.3882 Å, εi = 1 × 10−8 kcal/mol). For the parameters of the metal core, total two iterations of the parameter refining procedure were finally used. The eligible models generated in these iterations were shown in Table 2. The threshold values of all the iteration steps were shown in Table S2. During our parameter refining procedure, the threshold values were decreased stepwise to improve the accuracy. The calculated free energies and their corrections described in section 2.5 were shown in Table 2. The final refined model, with RMC = 1.036 and εMC = 0.066, had the exp closest ΔGcal sol value (−466.37 kcal/mol) to ΔGsol . The dV/dλ plots of the final model were shown in Figure S1. The average of d(Mg2+−O) of the six ligands was 2.090 Å, which was also the closest calculated value to the experimental value. For Duarte’s model tested in the AMBER ff03 force field, the average value of d(Mg2+−O) exceeded 2.12 Å, which was 0.03 Å above the experimental value. In addition, the ΔGcal sol of Duarte’s model was about 13 kcal/mol higher than the experimental value (see Table 2). 3.2. Parameter Testing. The comparison of QM optimized structure and the MM minimized structure of the simplified metal center for the protein 1M0W was shown in Table S4, and that for the protein 4U5Xwas shown in Table S5. All of the dummy models showed a slightly lower correlation of Hessian matrixes to the QM structures. But for the doublemetal center (Table S4), the geometry produced by our refined F


Article


Figure 4. Molecular structure of ATP and GTP with atom name used in this work (a) and root-mean-square fluctuation (RMSF) of non-hydrogen atoms of ATP or GTP in the simulated protein systems for 1PYX (b), 1M0W (c), and 4U5X (d). RMSF calculations used the last 10 ns trajectories. The regions related to the broken hydrogen bonds discussed in the text (marked in Figure 5 as well) are covered by yellow rectangles.

better Mg−ligand distances that are close to the crystal structure than the point charge model. In addition, the difference between these distances produced by Duarte’s model and our refined model is small (about 0.04 Å). But the produced coordination mode is incorrect for all these models compared to the crystal structure, with Asp200 changing from bidentate to monodentate and instead a water molecule coming into the coordination sphere. For the hydrogen bond network formed between the protein 1PYX and the cofactor ATP, most of them had occupancies in excess of 80%, except the two hydrogen bonds: Lys183-NZ···ATPO2G and Gln185-O···ATP-O3* (see Figure 5a). The hydrogen bond occupancy of Lys183-NZ···ATP-O2G was 45% for our model but was zero for Duarte’s model. Note that, this hydrogen bond can cooperate with Mg2+ to properly orientate the ATP-O2G atom and further fix the correct orientation of the γ-phosphate group (see Figure S3a), which is considered to be crucial during the reaction procedure. With the disrupted hydrogen bond Lys183-NZ···ATP-O2G, the conformational changes and fluctuations of ATP in the system simulated by Duarte’s model were larger than that simulated by our model (see Table 3 and Figure 4b) and finally caused the large structural deviation of the ATP molecule (especially the γphosphate group) in Duarte’s model (shown in Figure S4). However, the hydrogen bond Gln185-O···ATP-O3* had no significant effect on the geometry of ATP, although its occupancy in Duarte’s model was much higher than that in our model. To verify the enhanced performance of our refined model, we compared the simulation results of another double-metal center system, the yeast glutathione synthase (1M0W), with the Mg2+ parametrized by Åqvist’s model, Allnér’s model, Duarte’s model, our refined model, and the restrained model,

model is the best one. For the single-metal center (Table S5), the geometries produced by both of the dummy models are not as good as those produced by point charge models. Note that the performances of all the present models are worse than those of the bonded models demonstrated by Hu and Ryde.43 Next, we performed MD simulations on three ATP/GTP-Mg2+protein complex systems to further test the parameter performance. 3.2.1. Parameter Testing in Double-Metal Center System. Lu et al.6 reported that the dummy model optimized by Oelschlaeger et al. was inadaptable in simulating the multimetal center of kinase GSK3β (PDBID: 1PYX). In that study, the crucial hydrogen bond formed between Lys85-NZ and ATPO2B (the atoms of ATP or GTP and their names are shown in Figure 4a), which is known to play an essential role in the enzymatic catalysis reaction, was broken during the 20 ns simulation using Oelschlaeger’s dummy model. In our test, the occupancy of this hydrogen bond in each system remained to above 80% (see Figure 5a). This indicates that both Duarte’s model and our refined model can overcome the problem caused by using Oelschlaeger’s model. In addition, the whole protein in each system simulated using dummy models underwent less conformational changes than that in Lu’s study simulated using Åqvist’s model (RMSDs of backbone of dummy model systems were both less than 2.00 Å, Table 3). The CN of each Mg2+ was kept to be 6 (especially for our refined model with a standard deviation of zero, Table 3). The distance between the two Mg2+ ions was also equilibrated but larger than that obtained from Åqvist’s model (see Table 3 and Figure S2a). All the three models produced larger d(Mg−Mg) values than the crystal structure (3.8 Å). The Mg−ligand distances averaged during 20 ns simulations were shown in Table S8. It is clear that dummy atom models can produce G


Article


respectively. The protein 1M0W had the lowest RMSD of backbone when simulated using our refined model (see Table 4). Similar to the protein 1PYX, CNs of the two Mg2+ were constant to be six (see Table 4) for our refined model, but inconstant for the other four models. The distance between the two Mg2+ ions was equilibrated in Duarte’s model, our refined model and the restrained model, but disequilibrated in Åqvist’s model and Allnér’s model (shown in Figure S2b). All the five models produced a smaller d(Mg−Mg) than the crystal structure (4.1 Å). The Mg−ligand distances averaged during 20 ns simulations were shown in Table S7. Similar to the results of 1PYX, dummy atom models can still produce better Mg− ligand distances than the point charge models. The restrained model shows the best Mg−liganddistances because of the restrained bond lengths. The difference between these distances produced by Duarte’s model and our refined model is still small (about 0.04 Å). This table also demonstrates that the point charge models, as well as the restrained model, produced a changing coordination mode of Glu386, causing an inconstant coordination number. But for the dummy atom models and the point charge models, Glu146 changed to monodentate and again an extra water molecule came in, like those observed for protein 1PYX. Only the restrained model gives the correct coordination mode of Glu146. For the fluctuation of each nonhydrogen atom in ATP (shown in Figure 4c), the whole ATP molecule showed a high fluctuation both in Åqvist’s model and in Allnér’s model. But low fluctuations with a little difference were observed in Duarte’s model, our refined model and the restrained model. The fluctuation of the ribose moiety of the ATP simulated by our refined model, as well as that simulated by the restrained model, was the lowest among all the five models. The fluctuation of the γ-phosphate group of the ATP simulated in our refined model was a little higher than that in Duarte’s model. But this tiny increase of the fluctuation is considered not to be significant to affect the whole geometry of the ATP-Mg2+-protein system. On the other hand, for the hydrogen bond network formed between the protein 1M0W and the cofactor ATP, only half of them could retain 80% occupancy (see Figure 5b) in each of the four systems. In spite of this, the hydrogen bond Lys324-NZ···ATP-O1B reached 75% occupancy when simulated by our model, but showed all less than 30% occupancy for Åqvist’s model, Allnér’s model and the restrained model, and even zero occupancy by Duarte’s model. It is worth noting that this hydrogen bond has a spatial orientation similar to the crucial hydrogen bond Lys85-NZ··· ATP-O2B in the protein 1PYX (see Figure S3b), indicating this hydrogen bond is also crucial for the enzymatic activity. That is, for the protein 1M0W, the ATP-Mg2+-protein systems obtained from the MD simulations using our refined model may be more reliable than those using other nonbonded models but the produced metal center geometry is worse than that produced by the restrained or bonded model. 3.2.2. Parameter Testing in Single-Metal Center System. To verify the performance of our refined model in a small and simple system with a single Mg2+ ion, the protein GTPases OsRac1 (PDBID: 4U5X) was selected and parametrized by using the Åqvist’s model, Allnér’s model, Duarte’s model, and our refined model. Similar to the above two systems, the point charge models, Åqvist’s model, and Allnér’s model, showed a more fluctuant CN of Mg2+ than other two dummy atom models (see Table 5). The point charge models still produced worse Mg−ligand distances than the dummy atom models (shown in Table S8) but all the models produced a constant

Figure 5. Hydrogen bond (HB) occupancy of the ATP/GTP−protein systems for (a) 1PYX, (b) 1M0W, and (c) 4U5X. The HB network is obtained from the crystal structure. All the occupancies are calculated using the last 10 ns trajectories. The HBs whose occupancies have large differences among different models are covered by yellow rectangles.

Table 3. Root Mean Square Deviation (RMSD) of Protein Backbone and ATP, the Coordination Number (CN) of Mg2+, and Mg2+−Mg2+ Distance in the Double-Metal Center for Protein 1PYXa dummy atom model 1PYX RMSD of protein (Å) RMSD of ATP (Å) CN of Mg2+ Mg I Mg II d(Mg−Mg) (Å)

point charge model Åqvist

b

Duarte

this work

2.00 ± 0.16

1.83 ± 0.13

1.97 ± 0.20

6.00 6.00 4.55 ± 0.22

0.72 5.99 6.00 4.70

± ± ± ±

0.10 0.05 0.00 0.12

0.48 5.97 6.00 4.79

± ± ± ±

0.10 0.18 0.00 0.11

a

All the values were calculated using the last 10 ns trajectory (5000 frames) and given as “average ± standard deviation”. bThe results with Åqvist’s model were taken from Lu’s work.6

H


Article


Table 4. RMSD of Protein Backbone and ATP, CN of Mg2+, and Mg2+−Mg2+ Distance in the Double-Metal Center for Protein 1M0Wa point charge model 1M0W RMSD of protein (Å) RMSD of ATP (Å) CN of Mg2+

Åqvist

Mg I Mg II

d(Mg−Mg) (Å) a

1.85 0.55 5.99 6.00 3.70

± ± ± ± ±

dummy atom model

Allnér

0.12 0.10 0.08 0.13 0.40

1.97 0.57 5.97 5.99 3.78

± ± ± ± ±

0.11 0.19 0.16 0.04 0.36

Duarte 1.92 0.62 6.00 6.00 3.69

± ± ± ± ±

0.11 0.12 0.00 0.04 0.08

restrained

this work 1.69 0.44 6.00 6.00 3.64

± ± ± ± ±

0.09 0.07 0.00 0.00 0.08

Åqvist 1.84 0.43 6.00 6.61 3.86

± ± ± ± ±

0.14 0.09 0.00 0.54 0.13

All the values were calculated using the last 10 ns trajectory (5000 frames) and given as average ± standard deviation.

Table 5. RMSD of Protein Backbone and ATP and CN of Mg2+ for Protein 4U5Xa point charge model

dummy atom model

4U5X

Åqvist

Allnér

Duarte

this work

RMSD of backbone (Å) RMSD of GTP (Å) CN of Mg2+

1.34 ± 0.13

1.46 ± 0.15

1.39 ± 0.13

1.38 ± 0.14

0.43 ± 0.14

0.38 ± 0.13

0.47 ± 0.14

0.43 ± 0.13

5.99 ± 0.03

5.99 ± 0.04

6.00 ± 0.00

6.00 ± 0.00

a

All the values were calculated using the last 10 ns trajectory (5000 frames) and given as average ± standard deviation.

Figure 6. Relative error (E) of using Duarte’s model with AMBER03 force field. Calculation details for E are given in the Supporting Information.

coordination mode of Mg2+, indicating that point charge models are also suitable. For the simulated GTP structures the main difference is in the γ-phosphate group of GTP, where the fluctuation was a little lower for the point charge models than that for the dummy atom models (see Figure 4d). Correspondingly, a hydrogen bond Gly67-N···GTP-O1G (shown in Figure S3c) was found to have only about 20% occupancy for both Duarte’s model and our refined model, but more than 80% for both Åqvist’s model and Allnér’s model (see Figure 5c). For Duarte’s model and our refined model, the subtle difference between the fluctuations of the two GTP structures was caused by the more stable hydrogen bond Lys165-N···GTP-O6 (shown in Figure S3c) observed in our model (see Figure 5c). All the other hydrogen bonds reached more than 80% occupancy in each system during our simulations. Therefore, for the single-metal center system, point charge models were verified to be better than dummy atom models. It is worth noting that our model can also produce a better result than Duarte’s model. In spite of that, all the four models produced similar GTP-Mg2+-protein conformations. 3.3. Influence of Using Different Models in Different Force Fields. Compared with the point charge model (Åqvist’s model), dummy atom models have both an enhanced vdW interaction and an enhanced electrostatic interaction (shown in Figure S5) between the metal ion and its ligands. For the influence of different force fields, the relative error E of using Duarte’s model in the AMBER ff03 force field was calculated for each four common ligands of Mg2+ (shown in Figure 6). Each of the relative error tends to zero when the distance r is less than 1.5 Å but tends to a specific value that less than 20% when r is larger than 2.4 Å. Moreover, when r is around 1.75 Å, E tends to infinity, because in this range the potential Unonbond OPLS tends to zero.

4. DISCUSSION According to the above results, our refined dummy atom model of Mg2+, which was optimized in the AMBER ff03 force field, had a better performance of reproducing the experimentally observed geometric and thermodynamic properties of Mg2+ in water than Duarte’s dummy atom models simulated in the AMBER ff03 force field. For a double-metal center, our refined model can produce the geometry closer to that obtained from a high-level QM optimization (Table S4) than previous dummy atom models. But the performance is worse than that of the standard bonded models. For simulating the metalloenzymes that contains a double-metal center, the advantage of our refined model compared with Duarte’s model was not significant. For the protein 1M0W, some of the Mg−ligand distances produced by our refined model were even worse than that produced by Duarte’s model and far away from that produced by the restrained model (Table S7). Dummy atom models favor a monodentate binding mode for carboxylate groups (Glu146 and Glu386). For the restrained model, the restricted bond between the carboxylate group of Glu386 and Mg 2+ (monodentate binding mode in the crystal structure) leads a bidentate binding mode. Such incorrect binding mode may be fixed by increasing the vdW radius of metal ions until one does not get a bidentate coordination. In addition, all the models, even including the restrained model, failed to reproduce the Mg−Mg distances in the crystal structures. But for the QM optimized structure of the metal center of protein 1M0W, the Mg−Mg distance is 3.63 Å, which is close to the simulated distance 3.64 Å produced by our refined model. This evidence, which has good and bad points, is not enough to judge which I


Article

Journal of Chemical Information and Modeling ⎧ ⎪ lim E = 1 − ⎪ ⎪ r→0 ⎨ ⎪ ⎪ lim E = 1 − ⎪ r →∞ ⎩

model is better. A good simulated structure of an activated enzyme must retain the geometry features related to the enzymatic activity. Thus, the hydrogen bonds that connect the cofactor ATP or GTP to the protein and the global conformational changes of metal center were compared between different models to finally evaluate their simulation performances. It is therefore credible that the refined dummy model can produce a more reliable ATP-Mg2+-protein conformation for the double-metal system than the point charge models and previous dummy atom models. For the metal center geometry, however, the refined dummy atom model is much worse than the bonded or restrained models. For a single-metal center, our refined model produced the geometry not as good as that produced by point charge models (Table S5). For simulating the metalloenzyme that contains a single-metal center, the results showed a very little difference between these models. The well-optimized point charge models, Åqvist’s model and Allnér’s model, showed a better performance for simulating the single-metal system. In spite of this, our refined model still had a better performance than Duarte’s model under the AMBER ff03 force field. Therefore, our parameter screening strategy can produce a high-performance dummy atom model for a specific force field. Compared with the point charge models, the enhanced performance of the dummy atom model is mainly because of the enhanced vdW and electrostatic interactions (shown in Figure S5). The increased repulsion of vdW interaction and the increased attraction of electrostatic interaction can make the CN of Mg2+ and the geometry of ATP-Mg2+ system steadier in the double-metal center, where the interactions are complex. However, the interaction in the single-metal center that is parametrized by the point charge model may be strong enough to keep its stability. Therefore, the enhanced interactions caused by the dummy atom model may not have significant advantages compared with the point charge model. The differences between the calculated geometric and thermodynamic properties of Duarte’s Mg2+ model in aqueous solution from the present work and those from Duarte’s work are caused by the following two factors. The main reason for this is that we used the revised experimental solvation free energy to optimize the parameters. Duarte et al.14 reported the calculated solvation free energy of their model that equals −454.4 kcal/mol, which is very close to their referred free energy (−454.2 kcal/mol), by using their free energy calculation protocol. In the present work, we also obtained a very similar calculated solvation free energy that equals −455.28 kcal/mol after using our free energy correction protocol (see section 2.2). That is, Duarte’s model is indeed a well-optimized model for reproducing Noyes’s solvation free energy. However, the simulative performance of Duarte’s model is worse than that of our model, which is optimized to reproduce the revised experimental solvation free energy, proving that the revised experimental solvation free energy is indeed more accurate and suitable. The other reason is that we used the same force field (the AMBER ff03 force field) in the parameter refining and testing, and this force field is different from the one used in Duarte’s model optimized for (the OPLSAA force field). We used the relative error E of using Duarte’s model in the AMBER ff03 force field to demonstrate this effect (see section 2.6). The limitation of E can be derived as follows:

A OAMBER A OOPLS Q OAMBER Q OOPLS

(13)

The limitation value at zero only depends on the vdW parameters A of atom O in the two different force fields. Similarly, the limitation value at infinity only depends on the atomic charges of atom O in the two different force fields. Therefore, the larger the difference between these two force fields is, the higher the relative error will be. Furthermore, when r is between 1.7 and 1.9 Å, the relative error E has a very steep slope (see Figure 6) indicating a large deviation between the two nonbonding potentials of each ligand. This deviation may influence the calculated geometry of Mg2+−O, when using Duarte’s model in the AMBER force field, but the influence is negligible for TIP3P water. Therefore, the force field where the model is optimized is the best force field to perform the related MD simulations and a change in the force field may cause inaccuracy to some extent. If the metal center needs to be simulated with a high-accurate geometry consistent with the crystal structure, the restrained or bonded model is recommended. Otherwise, for the monodentate binding mode of carboxylate groups in double metal center systems, the dummy atom model, which can yield similar results to the restrained model, is recommended. For the binodentate binding mode of carboxylate groups in a doublemetal center system, the dummy atom model with the restriction on the distance between the metal and the carboxylate groups will be fine. We also suggest using either the dummy model or the point charge model to simulate the simple single-metal system. In addition, force field, accurate experimental solvation free energy, and suitable free energy correction protocol are crucial factors influencing model development and simulation performance. Combining these factors, the parameter screening strategy proposed in this work can be applied to construct dummy atom models for use in diverse force fields containing Coulomb and Lennard-Jones potentials without a loss in accuracy.

5. CONCLUSIONS In this work we refined the dummy atom model of Mg2+ under the AMBER ff03 force field with TIP3P water model via a simple parameter screening strategy. The refined model accurately reproduces the experimental values of metal−ligand distance, coordination number and solvation free energy of Mg2+ in water. Compared with two point charge models (Åqvist’s and Allnér’s models), the refined model showed an improved performance for simulating the complex doublemetal systems. Compared with the previous dummy atom models (Oelschlaeger’s and Duarte’s models), our refined model produced a more reliable ATP/GTP-Mg2+-protein conformation in both single- and double-metal systems, in which the crucial hydrogen bonds were highly retained. Similar to other unbounded models, the refined model in this work failed to reproduce the Mg−Mg distance and favored a monodentate binding of carboxylate groups. If one really wants to reproduce a metal site, a bonded model is needed. A nonbonded model should be used only if you want a fast model of a site that is not of great interest in the simulation. The J


Article


Kinase Gsk3β: A Molecular Simulation Study. Proteins: Struct., Funct., Genet. 2013, 81, 740−753. (7) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (8) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (9) MacKerell, A. D.; Banavali, N.; Foloppe, N. Development and Current Status of the Charmm Force Field for Nucleic Acids. Biopolymers 2000, 56, 257−265. (10) Panteva, M. T.; Giambaşu, G. M.; York, D. M. Comparison of Structural, Thermodynamic, Kinetic and Mass Transport Properties of Mg2+ Ion Models Commonly Used in Biomolecular Simulations. J. Comput. Chem. 2015, 36, 970−982. (11) Pang, Y. P. Novel Zinc Protein Molecular Dynamics Simulations: Steps toward Antiangiogenesis for Cancer Treatment. J. Mol. Model. 1999, 5, 196−202. (12) Oelschlaeger, P.; Klahn, M.; Beard, W. A.; Wilson, S. H.; Warshel, A. Magnesium-Cationic Dummy Atom Molecules Enhance Representation of DNA Polymerase B in Molecular Dynamics Simulations: Improved Accuracy in Studies of Structural Features and Mutational Effects. J. Mol. Biol. 2007, 366, 687−701. (13) Saxena, A.; Sept, D. Multisite Ion Models That Improve Coordination and Free Energy Calculations in Molecular Dynamics Simulations. J. Chem. Theory Comput. 2013, 9, 3538−3542. (14) Duarte, F.; Bauer, P.; Barrozo, A.; Amrein, B. A.; Purg, M.; Åqvist, J.; Kamerlin, S. C. L. Force Field Independent Metal Parameters Using a Nonbonded Dummy Model. J. Phys. Chem. B 2014, 118, 4351−4362. (15) Åqvist, J.; Warshel, A. Free Energy Relationships in Metalloenzyme-Catalyzed Reactions. Calculations of the Effects of Metal Ion Substitutions in Staphylococcal Nuclease. J. Am. Chem. Soc. 1990, 112, 2860−2868. (16) Pang, Y. P. Successful Molecular Dynamics Simulation of Two Zinc Complexes Bridged by a Hydroxide in Phosphotriesterase Using the Cationic Dummy Atom Method. Proteins: Struct., Funct., Genet. 2001, 45, 183−189. (17) 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. (18) Åqvist, J. Ion-Water Interaction Potentials Derived from Free Energy Perturbation Simulations. J. Phys. Chem. 1990, 94, 8021−8024. (19) Allnér, O.; Nilsson, L.; Villa, A. Magnesium Ion−Water Coordination and Exchange in Biomolecular Simulations. J. Chem. Theory Comput. 2012, 8, 1493−1502. (20) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the Opls All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (21) Marcus, Y. Ionic Radii in Aqueous Solutions. Chem. Rev. 1988, 88, 1475−1498. (22) Marcus, Y. A Simple Empirical Model Describing the Thermodynamics of Hydration of Ions of Widely Varying Charges, Sizes, and Shapes. Biophys. Chem. 1994, 51, 111−127. (23) Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. H.; Earhart, A. D.; Coe, J. V.; Tuttle, T. R. The Proton’s Absolute Aqueous Enthalpy and Gibbs Free Energy of Solvation from Cluster-Ion Solvation Data. J. Phys. Chem. A 1998, 102, 7787−7794. (24) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Aqueous Solvation Free Energies of Ions and Ion-Water Clusters Based on an Accurate Value for the Absolute Aqueous Solvation Free Energy of the Proton. J. Phys. Chem. B 2006, 110, 16066−16081. (25) Camaioni, D. M.; Schwerdtfeger, C. A. Comment on “Accurate Experimental Values for the Free Energies of Hydration of H+, OH-, and H3O+. J. Phys. Chem. A 2005, 109, 10795−10797.

advantage of our model can be attributed to the use of a revised (more accurate) experimental solvation free energy and of a suitable free energy correction protocol. The refined model provides a good choice to study Mg2+ ions with ligand exchange in silico using the AMBER ff03 force field with the TIP3P water model. This work also proposes a simple parameter screening strategy for refinement of dummy atom models to achieve a better performance.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.5b00286. Test of Saxena’s dummy model, computational details, and supplementary tables and figures (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address †

H.Z.: Department of Biological Science and Engineering, School of Chemistry and Biological Engineering, University of Science and Technology Beijing, 100083 Beijing, China. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Basic Research Program of China (973 program) (2013CB733600), the National Nature Science Foundation of China (21436002, 21390202). All of the simulations were supported by CHEMCLOUDCOMPUTING.

■

ABBREVIATIONS MD, molecular dynamics; QM, quantum mechanics; vdW, van der Waals; MC, metal core; D, dummy atom; d(MC−O), distance between the metal core and coordinated oxygen atoms; d(Mg−O), distance between the Mg2+ and coordinated oxygen atoms; RMSD, root mean square deviation; RMSF, root mean square fluctuation; d(Mg−Mg), distance between the two Mg2+ ions; CN, coordination number; GS, glutathione synthase; TI, thermodynamic integration

■

REFERENCES

(1) Bertini, I.; Gray, H. B.; Stiefel, E. I.; Valentine, J. S. Biological Inorganic Chemistry: Structure and Reactivity; University Science Books: Mill Vallcy, CA, 2007. (2) Williams, R. J. Metal Ions in Biological Systems. Biol. Rev. 1953, 28, 381−412. (3) Suárez, D.; Merz, K. M. Molecular Dynamics Simulations of the Mononuclear Zinc-B-Lactamase from Bacillus Cereus. J. Am. Chem. Soc. 2001, 123, 3759−3770. (4) Pabis, A.; Geronimo, I.; York, D. M.; Paneth, P. Molecular Dynamics Simulation of Nitrobenzene Dioxygenase Using Amber Force Field. J. Chem. Theory Comput. 2014, 10, 2246−2254. (5) Widderich, N.; Pittelkow, M.; Höppner, A.; Mulnaes, D.; Buckel, W.; Gohlke, H.; Smits, S. H.; Bremer, E. Molecular Dynamics Simulations and Structure-Guided Mutagenesis Provide Insight into the Architecture of the Catalytic Core of the Ectoine Hydroxylase. J. Mol. Biol. 2014, 426, 586−600. (6) Lu, S. Y.; Huang, Z. M.; Huang, W. K.; Liu, X. Y.; Chen, Y. Y.; Shi, T.; Zhang, J. How Calcium Inhibits the Magnesium-Dependent K


Article

Journal of Chemical Information and Modeling (26) Marcus, Y. The Thermodynamics of Solvation of Ions. Part 4. Application of the Tetraphenylarsonium Tetraphenylborate (Tatb) Extrathermodynamic Assumption to the Hydration of Ions and to Properties of Hydrated Ions. J. Chem. Soc., Faraday Trans. 1 1987, 83, 2985−2992. (27) Rosseinsky, D. Electrode Potentials and Hydration Energies. Theories and Correlations. Chem. Rev. 1965, 65, 467−490. (28) Noyes, R. M. Thermodynamics of Ion Hydration as a Measure of Effective Dielectric Properties of Water. J. Am. Chem. Soc. 1962, 84, 513−522. (29) Matlab, version 2011b; The MathWorks, Inc.: Natick, MA, 2011. (30) Case, D.; Darden, T.; Cheatham III, T.; Simmerling, C.; Wang, J.; Duke, R.; Luo, R.; Walker, R.; Zhang, W.; Merz, K. Amber, version 11; University of California, San Francisco, 2010. (31) 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. The Amber Biomolecular Simulation Programs. J. Comput. Chem. 2005, 26, 1668−1688. (32) Pearlman, D. A.; Case, D. A.; Caldwell, J. W.; Ross, W. S.; Cheatham, T. E.; DeBolt, S.; Ferguson, D.; Seibel, G.; Kollman, P. Amber, a Package of Computer Programs for Applying Molecular Mechanics, Normal Mode Analysis, Molecular Dynamics and Free Energy Calculations to Simulate the Structural and Energetic Properties of Molecules. Comput. Phys. Commun. 1995, 91, 1−41. (33) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; et al. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24, 1999−2012. (34) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (35) Humphrey, W.; Dalke, A.; Schulten, K. Vmd: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (36) Steinbrecher, T.; Mobley, D. L.; Case, D. A. Nonlinear Scaling Schemes for Lennard-Jones Interactions in Free Energy Calculations. J. Chem. Phys. 2007, 127, 214108. (37) Hummer, G.; Pratt, L. R.; García, A. E. Ion Sizes and Finite-Size Corrections for Ionic-Solvation Free Energies. J. Chem. Phys. 1997, 107, 9275−9277. (38) Darden, T.; Pearlman, D.; Pedersen, L. G. Ionic Charging Free Energies: Spherical Versus Periodic Boundary Conditions. J. Chem. Phys. 1998, 109, 10921−10935. (39) Kastenholz, M. A.; Hünenberger, P. H. Computation of Methodology-Independent Ionic Solvation Free Energies from Molecular Simulations. Ii. The Hydration Free Energy of the Sodium Cation. J. Chem. Phys. 2006, 124, 224501. (40) Warren, G. L.; Patel, S. Hydration Free Energies of Monovalent Ions in Transferable Intermolecular Potential Four Point Fluctuating Charge Water: An Assessment of Simulation Methodology and Force Field Performance and Transferability. J. Chem. Phys. 2007, 127, 064509. (41) Gogos, A.; Shapiro, L. Large Conformational Changes in the Catalytic Cycle of Glutathione Synthase. Structure 2002, 10, 1669− 1676. (42) Kosami, K. i.; Ohki, I.; Nagano, M.; Furuita, K.; Sugiki, T.; Kawano, Y.; Kawasaki, T.; Fujiwara, T.; Nakagawa, A.; Shimamoto, K.; Kojima, C. The Crystal Structure of the Plant Small Gtpase Osrac1 Reveals Its Mode of Binding to Nadph Oxidase. J. Biol. Chem. 2014, 289, 28569−28578. (43) Hu, L.; Ryde, U. Comparison of Methods to Obtain Force-Field Parameters for Metal Sites. J. Chem. Theory Comput. 2011, 7, 2452− 2463. (44) Bertrand, J.; Thieffine, S.; Vulpetti, A.; Cristiani, C.; Valsasina, B.; Knapp, S.; Kalisz, H.; Flocco, M. Structural Characterization of the Gsk-3β Active Site Using Selective and Non-Selective Atp-Mimetic Inhibitors. J. Mol. Biol. 2003, 333, 393−407.

(45) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03; Gaussian, Inc., Wallingford CT, 2004. (46) Meagher, K. L.; Redman, L. T.; Carlson, H. A. Development of Polyphosphate Parameters for Use with the Amber Force Field. J. Comput. Chem. 2003, 24, 1016−1025.

L


Refined Dummy Atom Model of Mg2+

Recommend Documents