Rational Design of Methodology-Independent Metal Parameters Using

Subscriber access provided by Temple University Libraries

Article

Rational Design of Methodology-Independent Metal Parameters Using a Nonbonded Dummy Model Yang Jiang, Haiyang Zhang, and Tianwei Tan J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00223 • Publication Date (Web): 16 May 2016 Downloaded from http://pubs.acs.org on June 5, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Rational Design of Methodology-Independent Metal Parameters Using a Nonbonded Dummy Model Yang Jiang†, Haiyang Zhang‡, and Tianwei Tan* †

†

Beijing Key Lab of Bioprocess, College of Life Science and Technology, Beijing University of

Chemical Technology, Beijing 100029, China

‡

Department of Biological Science and Engineering, School of Chemistry and Biological

Engineering, University of Science and Technology Beijing, 100083 Beijing, China

Corresponding author: *[email protected]

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT

A nonbonded dummy model for metal ions is highly imperative for the computation of complex biological systems with for instance multiple metal centers. Here we present nonbonded dummy parameters of 11 divalent metallic cations, namely, Mg2+, V2+, Cr2+, Mn2+, Fe2+, Co2+, Ni2+, Zn2+, Cd2+, Sn2+, and Hg2+, that are optimized compatible with three widely used water models (TIP3P, SPC/E, and TIP4P-EW). The three sets of metal parameters reproduce simultaneously the solvation free energies (∆Gsol), the ion−oxygen distance in the first solvation shell (IOD), and coordination numbers (CN) in explicit water with a relative error less than 1%. The main sources of errors to ∆Gsol that arise from the boundary conditions and treatment of electrostatic interactions are corrected rationally, which ensures the independence of the proposed parameters on the methodology used in the calculation. This work will be of great value for the computational study of metal-containing biological systems.

2


Page 2 of 39

Page 3 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. INTRODUCTION Metal ions in the biomolecular system have multiple functions and play an important role in the catalytic process.1 For computational studies of metalloenzymes, one of the most challenging tasks is to develop suitable force field parameters for metal ions, especially for divalent metal ions. Three strategies for metal modeling have been proposed: the bonded2, nonbonded point charge,3, 4 and nonbonded dummy atom model.5-9 The bonded model has predefined covalent bonds or harmonic restrains between the metal and ligands and can reproduce a highly accurate crystal structure. The force field parameters of the metal-ligand complex are usually generated and optimized by quantum mechanics (QM) calculations2. However, once the ligand is altered, the bonded model parameters need be reoptimized to match the new metal site. Moreover, it does not allow ligand exchange during the simulation. The point charge model is simply described as a point charge surrounded by van der Waals (vdW) potentials. Only the parameters of vdW potentials need to be optimized.3, 4, 10 Although the point charge model allows ligand exchange, it will be inadequate when the simulation system has a complex metal center (two or more metal ions in one active site) or transition metals are present.6, 11

Recently, Li and Merz have presented the Particle Mesh Ewald (PME)12 compatible point charge

model parameters for 16 divalent metal ions, and these parameters reproduce the experimental solvation free energy and the coordination geometry simultaneously by introducing a 12-6-4 Lennard-Jones (LJ) potential4. However, the coordination numbers (CNs) of V2+, Mn2+, Cd2+, Sn2+ and Hg2+ were not modeled accurately. The nonbonded dummy atom model makes a good balance between the former two models. It was originally presented by Åqvist and Warshel11, and then 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

improved by many researchers.5-9, 13. This model was originally designed for the divalent metal ions with CN equal to six. The metal ion has a metal core (MC) surrounded by six dummy atoms (D) to form an octahedron (Figure 1). There are covalent bonds (b0) between MC and the dummy atoms, but no covalent bond between the dummy atom and the ligand, which hence allows ligand exchange. Recently, Duarte et al.8 have optimized dummy model parameters of seven divalent metal ions, including Mg2+、Mn2+、Fe2+、Co2+、Ni2+、Zn2+、Ca2+, using the OPLS-AA force field and Liao et al.9 have provided the dummy atom model of Cu2+ that includes the Jahn-Teller effect. Till now, all the proposed dummy atom models, even most of the nonbonded models, were optimized to reproduce the experimental solvation free energies (∆Gsol) provided by Noyes14, Rosseinsky15 or Marcus.16 These three sets of ∆Gsol were derived from the so-called “conventional” values calculated using the tabulated thermodynamics data from National Bureau of Standards (NBS), followed by adding the proton hydration free energy and then converted to the “absolute” values. The proton hydration free energy cannot be directly measured by experiments, so different methods were used to obtain this value. Marcus derived the proton hydration free energy based on the extra-thermodynamic assumption “TATB hypothesis”,17 which is different from the methods reported by Rosseinsky15 and Noyes14, leading to deviations between the three sets of ∆Gsol. Later on, another proton hydration free energy was provided by Tissandier et al.18 using the cluster pair approximation method, which was considered to be accurate and recommended by several reports.19-21 However, no revised experimental solvation free energies for the divalent metal ions based on Tissandier’s proton hydration free energy have been reported and used as the reference for

4


Page 4 of 39

Page 5 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


force field development, except for our previous work13 where we demonstrated a remarkable advantage of using the revised experimental solvation free energy for the Mg2+ dummy model in the parameter refining. In the present work, the experimental solvation free energies of the divalent metal ions taken from Marcus16 are revised using Tissandier’s proton hydration free energy (details were shown in Table S1). Considering that Cu2+ has the Jahn-Teller effect and has been optimized well before9, we take the rest 11 ions, namely, Mg2+, V2+, Cr2+, Mn2+, Fe2+, Co2+, Ni2+, Zn2+, Cd2+, Sn2+, and Hg2+, to develop the dummy atom models for TIP3P, SPC/E and TIP4P-EW water models, respectively. In addition, a free energy correction is introduced into the parameter generation procedure, including the correction for use of periodic boundary conditions with the PME method and the correction for an improper summation scheme.22 More details on the free energy corrections are presented in the Section 2.5. Considering these corrections makes the calculated free energy independent on the simulation methodology such as the finite system, the lattice-sum and the cutoff-based methods.22 That is, the optimized models in this work can reproduce the experimental solvation free energies using any of the above three simulation methodologies just by adding the corresponding corrections.22 The major strength of our models is that they can reproduce accurate coordination geometries without high-level QM calculations and without additional modifications of the vdW potential such as the 12-6-4 LJ potential proposed by Li and Merz4, especially for V2+, Cd2+, Sn2+ and Hg2+, whose CN values cannot be reproduced correctly by a nonbonded model before. Unlike the point charge metal parameters presented by Li and Merz that are only compatible

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 39

with the PME simulations, the models in this work are in principle compatible with almost all the simulation systems.

2. MATERIALS AND METHODS 2.1. Force Field Used in This Work To generate the dummy atom model parameters of the eleven divalent metal ions, we used Duarte’s parameters as the reference8. The values of force constants (Kb for bonds and Kθ for angles), equilibrium bond length (b0) and angle (θ0), atomic mass and charge of Duarte’s model were still used here. The vdW parameters for each metal type were set to be the vdW radius (R) and the 12-6 LJ potential well depth (ε). The nonbonding interaction potential function of two atoms i and j in the AMBER ff03 force field used in the present study has the following form23:

=

+

=

− 2 + ,

(Eq. 1) where rij is the distance between atom i and j in Å, qi is the charge of atom i in e, and k is a constant equaling 332. Rij and εij satisfy the Lorentz/Berthelot mixing rules: =

1 + = + 2 .

= " = "

Note that Ri here is the same with Rmin/2 in many other papers. 2.2. Experimental Values Used in This Work

6


(Eq. 2)

Page 7 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The experimental values of the ion-oxygen distance (IOD) and the coordination number (CN) for each metal were obtained from several works24-27 (shown in Table S1). The experimental value of the solvation free energy ∆%& was obtained from Marcus’s work16 and then revised by the '(

accurate value of the proton hydration free energy recommended by many reports18-21. The proton

hydration energy (-265.9 kcal/mol) has been corrected by converting the 1-atm gas phase standard ∗ state ∆%& to the 1-M solution standard state ∆%& (by subtracting ∆% → ∗ = 1.9 kcal/mol)19.

Note that, Marcus’s solvation free energies also have been corrected to the 1-M solution standard state, but the correction seems to have a wrong direction in his work16 (∆% → ∗ was added to

∗ ∆%& ). Thus, for the absolute solvation free energy provided by Marcus ( ∆%+,-./& ), the

. conventional ionic solvation free energy ∆%& can be derived by

. ∗ ∆%& = ∆%+,-./& − ∆% → ∗ − 0∆%1 ,

(Eq. 3) ∗ where ∆%+,-./& is under 1-M solution standard state, z is the net charge of the ion (2 for divalent

ions) and ∆%1 is the proton hydration free energy (-252.3 kcal/mol from Marcus). Then our revised metal ionic solvation free energy can be derived by using the accurate proton hydration free energy (∆%1 = -265.9 kcal/mol) and then corrected to the 1-M solution standard state: '( . ∆%& = ∆%& + 0∆%1 − ∆% → ∗ .

(Eq. 4)

The revised experimental solvation free energies for the divalent metal ions are shown in Table S1. 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.3. General Approach of Molecular Simulations All simulations were performed using the AMBER14 suite28 with the all-atom ff03 force field29. A 1000-step energy minimization was performed to eliminate improper contacts in each system. After energy minimization, each system was heated gradually from 0 to 298 K. During the heating steps, position restraints were imposed on the metal ions with a force constant of 10.0 kcal/mol/Å2. Each system was then equilibrated at the constant temperature (298 K) and pressure (1 bar) conditions via the Langevin dynamics (the collision frequency is 1.0 ps-1), with a coupling constant of 0.2 ps for both parameters. All the production simulations were performed at 298 K and 1 bar using a time step of 2 fs. Electrostatic interactions were calculated using the PME algorithm. The cutoff distance for van der Waals interactions was 10.0 Å. The SHAKE algorithm30 was applied to the bonds involving hydrogen. 2.4. Two-step Thermodynamic Integration The solvation free energy ∆%& of metal ion M2+ was calculated using the two-step thermodynamic integration (TI) method. The first step was perturbing nothing to the metal atom M (with no charge) in water using the soft-core potential (“appearing” part) in Amber11 suite31. For the parameter generation, 〈

345 6

〉 was calculated from five independent simulations at λ = 0.04691,

0.23076, 0.50000, 0.76923 and 0.95308, including a 400-ps energy equilibration and a 100-ps production. For the model testing, 〈

345 6

〉 was calculated from nine independent simulations at λ =

0.01592, 0.08198, 0.19331, 0.33787, 0.50000, 0.66213, 0.80669, 0.91802 and 0.98408. The parameters of soft-core potential were set to default. 〈

345 6

〉 was obtained from the production

8


Page 8 of 39

Page 9 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


89 simulation. The free energy change ∆% was then calculated using the Gaussian quadrature

formulas. The 〈

345 6

〉 plots for the metal ions calculated using the final parameters for TIP3P water

were shown in Figure S1. Others are similar with these and not shown.

The second step was charging the metal M from 0 e to +2 e by a parameter λ. For the parameter generation, 〈 6 〉 was calculated from five independent simulations at λ = 0.04691, 0.23076, :

0.50000, 0.76923 and 0.95308, including a 400-ps energy equilibration and a 100-ps production. For the model testing, 〈 6 〉 was calculated from nine independent simulations at λ = 0.01592, 0.08198, :

0.19331, 0.33787, 0.50000, 0.66213, 0.80669, 0.91802 and 0.98408. 〈 〉 was obtained from the 6 :

>? production simulation. The free energy change ∆%+ ; →+

Rational Design of Methodology-Independent Metal Parameters Using

Recommend Documents