New Zinc Ion Parameters Suitable for Classical MD Simulations of

Article pubs.acs.org/jcim

Cite This: J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

New Zinc Ion Parameters Suitable for Classical MD Simulations of Zinc Metallopeptidases Antonija Tomic,́ †,§ Gordan Horvat,‡ Michael Ramek,§ Dejan Agic,́ ∥ Hrvoje Brkic,́ ⊥,# and Sanja Tomic*́ ,† †

Division of Organic Chemistry and Biochemistry, Ruđer Bošković Institute, Bijenička 54, 10 000 Zagreb, Croatia Department of Chemistry, Faculty of Science, University of Zagreb, Horvatovac 102A, 10 000 Zagreb, Croatia § Institute of Physical and Theoretical Chemistry, Graz University of Technology, Stremayrgasse 9, 8010 Graz, Austria ∥ Faculty of Agrobiotechnical Sciences Osijek, Josip Juraj Strossmayer University of Osijek, Petra Svačića 1d, 31 000 Osijek, Croatia ⊥ Faculty of Medicine, Josip Juraj Strossmayer University of Osijek, J. Huttlera 4, 31 000 Osijek, Croatia # Faculty of Dental Medicine and Health, Josip Juraj Strossmayer University of Osijek, Crkvena 21, 31 000 Osijek, Croatia

Downloaded via UNIV OF SOUTHERN INDIANA on July 17, 2019 at 14:49:15 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

S Supporting Information *

ABSTRACT: The main aim of this work was to find parameters for the zinc ion in human dipeptidyl peptidase III (DPP III) active site that would enable its reliable modeling. Since the parameters publicly available failed to reproduce the zinc ion coordination in the enzyme, we developed a new set of the hybrid bonded/nonbonded parameters for the zinc ion suitable for molecular modeling of the human DPP III, dynamics, and ligand binding. The parameters allowed exchange of the water molecules coordinating the zinc ion and proved to be robust enough to enable reliable modeling not only of human DPP III and its orthologues but also of the other zinc-dependent peptidases with the zinc ion coordination similar to that in dipeptidyl peptidases III, i.e., peptidases with the zinc ion coordinated with two histidines and one glutamate. The new parameters were tested on a set of 21 different systems comprising 8 different peptidases, 5 DPP III orthologues, thermolysin, neprilysin, and aminopeptidase N, and the results are summarized in the second part of the article.

■

INTRODUCTION Metal ions play a vital role in biological systems with estimates that over 50% of proteins contain metal ions.1 In proteins, metal ions play either a catalytic role, by participating directly in the chemical catalysis, or a structural role solely to maintain protein structure and stability.2 With computational tools increasingly becoming important in chemical research, methods have emerged to effectively face the challenge of modeling metal ions in the gas, aqueous, and solid phases. In their recent article, Li and Merz have reviewed both quantum and classical modeling strategies for metal ion containing systems that have been developed over the past few decades.3 Besides experiments,4,5 quantum mechanical (QM) studies of metal ion containing systems at the semiempirical, ab initio, and density functional levels of theory6,7 are a wealthy source of data used in classical modeling efforts. As a result there are numerous models for modeling metal ions in classical simulations nowadays, ranging from the simplest, 12-6 and 12-6-4 LJ-type nonbonded models,8−11 to the bonded models,12,13 which use predefined covalent bonds between the metal and its ligands while the partial atomic charges are obtained from RESP fitting14 or CMX models15 typically resulting in a noninteger charge on the metal ion unlike the formal charge typically used in the nonbonded model. Due to their rigidity the latter models often cause artifacts on the ligand (water, substrate, etc.) conformational sampling © XXXX American Chemical Society

and dynamics. A semibonded model proposed by Pang et al.16−18 for the zinc ion with four covalently bonded cationic dummy atoms (CaDA) around the zinc nucleus lessened rigidity of the metal ion ligands but overemphasized the tetrahedral coordination of the metal ion. Just recently Yang et al.19 have used the Metal Centre Parameter Builder (MCPB)20 to generate new hybrid bonded/nonbonded FF parameters to describe the zinc active site architecture in the M1 and M17 aminopeptidase from Plasmodium falciparum. In this approach, metal−protein interactions were modeled by bonding parameters, while leaving all other zinc ligands unrestrained (through nonbonded interactions). Finally, besides nonpolarized, there are significant efforts in the development of polarizable models (e.g., the fluctuating charge, Drude oscillator, and the induced dipole models)3 which explicitly treat polarization in MD simulations. In aqueous solution, it is well established that Zn2+ has the coordination number six and forms an octahedral structure with water. However, analysis of the protein X-ray structures revealed that among 126 structural proteins, a majority of the zinc sites (82%) are 4-coordinated, only 14% are 5coordinated, and 4% are 6-coordinated; of 147 catalytic zinc ion binding proteins, 58% are 4-coordinated, 31% are 5Received: March 18, 2019 Published: July 5, 2019 A

DOI: 10.1021/acs.jcim.9b00235 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling coordinated, and 11% are 6-coordinated.21 Therefore, the tetrahedral coordination is dominant in proteins although in a significant number of proteins zinc ions have higher coordination numbers. When treated in a nonbonded fashion, the four-coordinate (tetrahedral) zinc complex identified in the enzymes X-ray structures has often changed to a six-coordinate (octahedral) zinc complex during MD simulations.18,22−24 This strong preference for the octahedral coordination is a consequence of the force field parameters developed for a zinc ion coordinating to six water molecules and the assumption that the coordination geometry of the metal is determined by the repulsion among its ligands.18 Human dipeptidyl peptidase III (h.DPP III) is a two-domain metallopeptidase, for which the zinc cation participates in the peptide bond hydrolysis. It is involved in the intracellular protein catabolism, pain modulation, and defense against oxidative stress.25−27 Its broad specificity toward peptides of varying lengths and compositions is still not properly understood although it could be partially explained by its ability to fluctuate between an open form with a wide interdomain cleft and a more compact one referred to as closed (PDB codes 3FVY and 5EGY, respectively). In the experimentally determined h.DPP III crystal structures the zinc ion is mostly tetracoordinated, by two His (H450, H455) and Glu (E508) residues, from the conserved hexapeptide motifs H450EXXGH455 and E507ECXXE512, respectively, and water. One additional water molecule is found at 2.6 Å from the zinc ion in the structure of the open h.DPP III and therefore could be considered as a fifth metal ligand. This is in agreement with the h.DPP III catalyzed Leu-enkephalin hydrolysis, determined by quantum mechanics−molecular mechanics (QMMM) calculations, with the zinc ion coordination changing from four to five during the reaction.28 Distances between the zinc ion and its nearest residues in these complexes are given in Table 1. Coordination of the active site zinc ion by two histidines from the first and glutamate from the second conserved motif is common in all known DPP III orthologues. In the crystal structures of the yeast DPP III, y.DPP III (PDB code 3CSK), the zinc ion is tetracoordinated while in bacterial and mushroom DPP III orthologues (in

DPPs III from Caldithrix abyssi, PDB code 6EOM, Bacteroides thetaiotaomicron, PDB code 5NA7, and Armillaria tabescens, PDB code 5YFB, in the text below referred as Ca.DPP III, Bt.DPP III, and At.DPP III, respectively) it is pentacoordinated. Besides, in the structure of the Bt.DPP III complex with the TRIS buffer (PDB code 5NA6) there are two orientations of the ligand; in one TRIS coordinates the zinc ion via its nitrogen and in the other with two oxygen and nitrogen atoms resulting with tetra- and hexacoordinated zinc ion in this structure, respectively. In the above-mentioned structures, glutamate corresponding to Glu508 in human DPP III coordinates the metal ion monodentately, except in the structure of Ca.DPP III, where both carboxyl oxygens take part in the zinc ion coordination. A similar type of coordination (2His + Glu) was found in the structures of many other proteins,21,29 especially in peptidases, like for example, thermolysin, neprilysin, anthrax lethal factor, aminopeptidase N, .... In the experimentally determined structures of thermolysin, the zinc ion is mostly four- or five-coordinated, by two His residues (H142 and H146), Glu166, and either water or ligand, wherein the coordination by Glu and ligand is predominantly (>85%) monodentate (Table S1). In 6 out of 169 structures there are 2 ligand molecules coordinating the zinc ion, and in an additional 6 structures there are 2 water molecules within 2.7 Å of zinc. Also, Blumberger at al. found that during the thermolysin-catalyzed peptide hydrolysis the tetracoordinated zinc ion in enzyme−substrate complex transforms to pentacoordinated one in the first transition state.30 Exchange of the zinc ion coordination number from tetra to penta was found in some other enzyme-catalyzed reactions as well.31 Of 15 analyzed X-ray structures of neprilysin, the zinc ion is tetracoordinated in 10 (by two histidines, glutamate and ligand), and in the others it is pentacoordinated with the exception of one structure, where the substrate ([2(r,s)-2sulfanylheptanoyl]-Phe-Ala) is a bidentate ligand and the other carboxyl oxygen of glutamate is at 2.51 Å of zinc (Table S1). However, there is no experimentally determined structure of the ligand-free neprilysin. Among 73 analyzed structures of aminopeptidase N, 48 have a tetra- and 23 structures have a pentacoordinated zinc ion, while in 2 structures it is coordinated by 3 ligands (Table S1). During MD simulations of DPPs III utilizing nonbonding parameters for the zinc ion (charge 2.0 e, while the VdW radius and energy well are 1.22 Å and 0.250 kcal/mol, respectively)29,32−36 we had noticed a sudden change from tetrahedral to octahedral zinc coordination already at the optimization phase. Although in general the MD simulation results agreed well with the experimental results, i.e., provided accurate information on intrinsic protein dynamic and the relative ligand binding affinities, the strong preference for the octahedral zinc coordination established by one additional water molecule and carboxyl oxygen of the glutamate from the conserved HEXXGH motif22−24,37 makes the use of such parameters questionable. Therefore, our main aim was to determine (find) parameters for the zinc ion able to reproduce the reactive structure of the h.DPP active site during MD simulations with the Amber force field. In this structure, similarly to the structures of h.DPP III determined by X-ray diffraction, the zinc ion is monodentately coordinated by its ligands: two histidines, one glutamate, and one or two water molecules; i.e., the parameters should enable exchange of the coordination number (CN) between 4 and 5

Table 1. Distances between the Zinc Ion and Its Nearest Residues in the Michaelis Complex Structure Obtained by QMMM Optimization of Closed h.DPP III−Leu-enkephalin Complex, and the Average Distances (and Standard Deviations) in Crystal Structures of h.DPP III (with Zinc in the Active Site) Available in the Protein Data Bank (3FVY, 5E33, 5E3A, 5EGY, 5EHH, and 5E3C)

H450NE2−ZN H455NE2−ZN E508OE1−ZN E508OE2−ZN E451OE1−ZN E451OE2−ZN Os−ZNa Owater−ZN

QMMM-optimized DPP III−Leuenkephalin complex28

h.DPP III crystal structures

d/Å

(⟨d⟩ ± SD)/Å

2.12 2.12 2.04 3.07 3.70 3.96 2.89 2.02

2.02 2.07 2.07 3.10 3.90 4.40 3.10 2.50

± ± ± ± ± ± ± ±

0.04 0.05 0.05 0.32 0.00 0.00 0.95 0.08

a

Os is the carbonyl oxygen atom from the second peptide bond from the substrate N-terminus. B


Article

Journal of Chemical Information and Modeling

Figure 1. Definition of hybrid bonded/nonbonded models of the h.DPP III active site zinc ion. Schematic representation of the CH3-capped metal center (small model) used for the QM calculations: (a) 3-ligand model, (b) 4-ligand model, and (c) 4-ligand-extended model. Newly parametrized bonds are shown in thick black lines, while water and amino acid residues whose charges were fixed to values from the ff14sb FF during RESP fitting are colored blue.

EECXXE hexapeptide motifs, while in the Ca.DPP III (PDB code 6EOM) the first of them is replaced by the pentapeptide HEISH.

coupled exclusively to exchange of water molecules in zinc coordination. In present work we (a) summarize the results of tests performed for human DPP III using already published parameters for the different zinc ion models: available nonbonded, semibonded CaDA, and hybrid bonded/nonbonded, (b) derived three sets of new hybrid bonded/ nonbonded parameters to be used for the zinc ion in human DPP III and tested them on its different conformations and h.DPP III−ligand complexes, and (c) discuss the set of parameters with the best performances on h.DPP III, also tested on a set of 15 different systems with some other zincdependent peptidases. The following sets of parameters were derived: (1) 3-ligand, derived using exactly the same procedure as Yang et al.,19 (2) 4-ligand parameters which were obtained including the water coordinating the zinc ion in the procedure, and (3) the 4extended model parameters which were developed in a way that besides the zinc ligands (i.e., three amino acids from the first coordination sphere and water molecule), the amino acids residues from the second zinc ion coordination sphere were included in the parameters developing procedure as well, wherein the latest one proved to be the best. With this set of parameters the proper positioning of all active participants in the h.DPP III catalyzed reaction, the substrate, water molecule, and the active site amino acid residues, was obtained clearly showing the importance of including the second metal ion coordination sphere in the QM calculations based procedure for the parameters derivation. Testing was performed on 21 different systems comprising five DPP III orthologues, thermolysin, neprilysin, and aminopeptidase N, either ligand-free or in the complex with ligands. It should be noted that the sequence identity between human DPP III and the other four DPP III orthologues is rather low, ranging from only 14.29% for Ca.DPP III to 40.66% for At.DPP III.34,36,38,39 Further, human, yeast, Bt.DPP III, and mushroom At.DPPs III possess the conserved HEXXGH and

■

COMPUTATIONAL METHODS

Nonpolarizable Force Field MD Simulations. The zinc active site models were built using different approaches, from the one based on the simplest nonbonded model (i.e., the 126-4 LJ-type model) to the more complicated semibonded models, CaDA and hybrid bonded/nonbonded ones. Nonbonded Zinc Models. In order to consider the ioninduced dipole interactions during MD simulations, we added the r−4 term to the standard 12-6 LJ nonbonded zinc ion model using default parameters and combining the rule from Li and Merz.8 Specifically, by using the add12_6_4 command in the parmed program available within Amber16, we generated a topology file that contained the additional C4 term between the divalent zinc ion and every other atom type.8 For the C4 term between the metal ion and the oxygen atom in the water molecule, we used the C4 values stored in parmed for the TIP3P model40 of water. CaDA Model. The CaDA model of the zinc ion was built using the Amber FF parameters for the tetrahedron-shaped zinc divalent cation, which are available from the developer’s web page http://www.mayo.edu/research/labs/computeraided-molecular-design/projects/zinc-protein-simulationsusing-cationic-dummy-atom-cada-approach. In order to simulate the propensity of zinc to tetrahedral coordination geometry, the metal ion was replaced by a five-atom molecule termed tetrahedral zinc divalent cation, in which four identical dummy (pseudo) atoms representing the four vacant 4s4p3 orbitals of zinc18 were placed at the four apexes of the tetrahedron, with the zinc nucleus located in its center. The dummy atoms, covalently bonded to the zinc nucleus, interacted with the surrounding electrostatically (q = +0.5 e), while to the zinc nucleus experienced only van der Waals interactions (r = 3.1 Å, ε = 1 × 10−6 kcal/mol, and q = 0). C


Article

Journal of Chemical Information and Modeling Pang et al.17,18 have shown that in a case when the second zinc coordination shell Glu residue forms a hydrogen bond with a His residue from the first zinc coordination shell, the histidine imidazole ring will be deprotonated and the glutamates will be protonated. Therefore, we have built two types of h.DPP III structures, with histidines coordinated to the tetrahedronshaped zinc divalent cation, H450 and H455: (a) in their anionic form and the second zinc coordination shell Glu residues neutral, and (b) with neutral His residues and negatively charged Glu. Hybrid Bonded/Nonbonded Models. Four different sets of bonded plus nonbonded parameters were used to describe the architecture of the h.DPP III active site during our MD simulations that utilize the hybrid bonded/nonbonded approach. The first set of parameters, referred as M1 model, is the one developed by Yang et al.19 for the zinc ion in M1 aminopeptidase. In their parametrization procedure they considered three metal ion ligands, two His residues and one Glu, but not the fourth one, a water molecule. Further, using the MCPB.py20 program, we generated three sets of new parameters for the zinc ion in the h.DPP III active site: 3ligand model, 4-ligand model, and 4-ligand-extended model. MCPB.py uses two models to perform the parametrization: a smaller one to obtain the metal-associated bond and angle parameters and a larger one to parametrize the partial charges. The same amino acid and non amino acid ligands are included in both models. The amino acid residues that coordinate the metal ion are represented by the CH3R (R is an amino acid side chain) in the small model, while in the large one the amino acid residues are capped by acetyl (ACE) and N-methyl (NME) groups. In the 3-ligand model parameters (bonded and nonbonded) were calculated for three amino acid residues that monodentately coordinate the zinc ion (H450, H455, and E508) in the closed h.DPP III active site (Figure 1a). The B3LYP/6-31G(d)-optimized structure (using Gaussian 09)41 of a small model (see Figure 1a) was used to determine the bond (ZN−H450NE2, ZN−H455NE2, and ZN−E508OE1) and angle parameters applying the Seminario method built in the MCPB.py program. For dihedral angles, MCPB.py assigns zero torsion barriers.20 Partial atomic charges were determined on a larger model after performing hydrogen-only geometry optimization using the B3LYP/6-31G(d) level of theory. RESP (restrained electrostatic potential)20 fitted charges of the metal ion and the residues coordinating it (including their side chains) were obtained using charge model B (ChgModB) available in MCPB.py. Consequently, these residues were assigned new Amber internal atom types and new residue names to differentiate them from those in the ff14sb FF (Table S2).42 The VDW radius of the zinc ion was set at 1.395 Å for the RESP fits. The New Hybrid Bonded/Nonbonded Models Derived Using an Extended Metal Ion Environment. In the next two models, referred as 4-ligand model (Figure 1b) and 4-ligandextended model (Figure 1c), the water molecule coordinating the zinc ion was included in the parametrization procedure, and, in the case of 4-ligand-extended model the second coordination shell amino acid residues participating in stabilization of the metal-coordinated residues Y318, E451, E507, and E512 were also considered in the parametrization procedure. The active site geometry with the zinc ioncoordinated water molecule was obtained from the QMMMoptimized Michaelis complex structure.28 Hu and Ryde43 showed that the metal coordinating water molecules treatment

might be a problem when using the bonded models. In our model water molecule was represented by nonbonding parameters. The water and the second shell residues were allowed to move during the optimization of the small model, and the positions of their hydrogen atoms were optimized during the hydrogen-only minimization of the larger model. During the RESP charge fitting procedure the partial charges on Y318, E451, E507, E512, and the water molecule were fixed to values assigned from the ff14sb FF. In addition, the water molecule−zinc ion contact was treated as nonbonded, enabling a water exchange during the MD simulations. MD System Setup and Simulation Protocol. Initial structures of the open (oDPP III) and closed (cDPP III) h.DPP III in the ligand-free MD simulations were derived from the crystal structures of h.DPP III, PDB codes 3FVY and 5EGY, respectively. In each structure, the water molecule closest to the zinc ion was preserved. Similarly, the initial structures in the simulation of the other systems were also the experimental ones (single crystal diffraction). For the simulations of yeast, bacterial, and mushroom (Armillaria tabescens) DPP III orthologues the structures determined with the highest resolution were used as initial (PDB codes 5NA6, 6EOM, 3CSK, and 5YFB, respectively). Initial structures for simulation of thermolysin, neprilysin, and aminopeptidase N were the structures available in the PDB database under codes 5A3Y, 4ZR5, and 4XN4, respectively. The crystal structure of the h.DPP III−tynorphin complex (PDB code 3T6B) was used as a template for building complexes with synthetic substrate Arg-Arg-2-naphthylamide (RRNA) as well as with Leu-enkephalin. The latter one was used in QMMM calculations of the chemical reaction (for details, see Tomić et al.)28 since at that time the crystal structure of DPP III−Leu-enkephalin complex (PDB code 5E3A) was not available. In the QMMM-optimized h.DPP III−Leu-enkephalin structure, the water molecule coordinated to the zinc ion (ZN−WAT distance 2.0 Å) is hydrogen bonded to the catalytically important E451. There is no water molecule within 2.5 Å of the assumed zinc ion position in the crystal h.DPP III−tynorphin complex structure (PDB code 3T6B), while in the h.DPP III−Leu-enkephalin complex structure there is one at 3.7 Å from the zinc ion, and it was preserved. The aforementioned experimentally determined structures of the closed h.DPP III (PDB codes 5EGY, 3T6B, and 5E3A) represent its E451A mutant, so we mutated Ala451 back to Glu before the simulations using the program Leap which adds the missing side chain atoms according to the library template. The same procedure was used for preparation of the other systems for MD simulations as well; i.e., only the systems with the wild type proteins were simulated. In addition, in the 3T6B structure the missing zinc ion was added using the crystal structure of the open h.DPP III (PDB code 3FVY) as a template. Missing amino acid atoms and residues were added using the program Modeller9v2.44 All Arg and Lys residues in the structure were positively charged (+1 e), while Glu and Asp residues were negatively charged (−1 e), as expected at physiological (experimental) conditions. The protonation of histidines was checked according to their ability to form hydrogen bonds with neighboring amino acid residues or to coordinate the metal ion. Exception were the CaDA simulations, where the histidines coordinating the zinc divalent cation in h.DPP III, H450, and H455 were in their anionic form and the second-shell Glu residues (E512 and E507) were D


Article


Pang et al.17 that the imidazole ring of histidine coordinated to the zinc ion will exist as imidazolate anion, when it is hydrogen bonded to the carboxylate group of Glu (or Asp) from the second zinc coordination shell. The Glu (or Asp) carboxylate group serves as an effective proton acceptor for the zinccoordinated imidazole proton donor. For this purpose we started from the QMMM-optimized structure of the h.DPP III−Leu-enkephalin complex28 and included E507 and E512 in the QM region. We have constructed two systems with different protonation states of H450−E512 and H455−E507 hydrogen-bonded pairs; the first system, named H0−E− (Figure 2a, left), had neutral histidines in the first coordination

protonated. Nonpeptide ligands were parametrized using the generalized Amber force field (gaff2). The protein and protein−substrate complexes were placed into the truncated octahedral box filled with TIP3P water molecules (the minimum distance between any atom originally present in solute and the edge of the periodic box was 14 Å), and Na+ ions were added in order to neutralize the system. All simulations were carried out using the Amber16 MD package45,46 and the ff14sb FF.42 Initial structures were minimized in three cycles. In the first cycle of optimization (1500 steps), water molecules were relaxed, while the rest of the system was harmonically restrained with a force constant of 32 kcal mol−1 Å−2. In the second cycle (2500 steps), a force constant (1 kcal mol−1 Å−2) was applied to the zinc ion and the protein backbone, while the final optimization (1500 step) was conducted without any constrains. Minimizations were performed using 470 steps of steepest descent followed by the conjugated gradient for the remaining steps. Afterward the systems were heated from 0 to 300 K and equilibrated: DPP III systems during 100 ps and distant peptidases during 7.3 ns (using NVT ensemble with an integration step of 1 fs). This was followed by productive MD simulations at 300 K using NpT ensemble and an integration step of 2 fs in duration of at least 100 ns (the productive simulations ranged from 100 to 220 ns). The SHAKE algorithm was used to constrain covalent bonds involving hydrogens. The pressure was maintained by the Berendsen barostat at 1 atm with the pressure relaxation time of 1 ps, while the system temperature was held constant at 300 K using the Langevin thermostat with a collision frequency of 1 ps−1. A cutoff of 11 Å was used for nonbonded interactions while long-range electrostatic interactions were treated using the particle mesh Ewald (PME) method.47,48 Polarizable Force Field MD Simulations. In order to conduct Drude FF MD simulations, a modified version of Gromacs was used.49 The zinc ion was simulated using the parameters from the Drude force field described by Haibo Yu et al.50 Initial h.DPP III structures were the same as the ones used in nonpolarizable FF MD simulations with the water closest to the zinc ion being included. These structures were solvated in NPT preequilibrated boxes of 35 000−45 000 SWM4-NDP51 water molecules. During the minimization, the position of both protein and ions in the system was kept fixed by the harmonic restraint potential with a force constant of 2.39 kcal mol−1 Å−2. After the minimization, the system was heated from 1 to 300 K during 0.5 ns using the simulated annealing protocol and the constant volume conditions, followed by preproduction simulation in the NVT ensemble at 300 K during which the system equilibrated. Finally, a production simulation of 5 ns under the NVT conditions at 300 K was conducted. During the NVT simulations, the Verlet algorithm52 was used for the atom propagation, while Drude particles propagated by the extended Lagrangian dynamics method49,53 using a time step of 1 fs. The temperature of the atoms was controlled by means of the Nose−Hoover algorithm54,55 with a time constant of 0.1 ps; the temperature of Drude particles was kept at 1 K by using the same thermostat and a time constant of 0.005 ps. The cutoff radius for nonbonded van der Waals and short-range Coulomb interactions was 12 Å. Long-range Coulomb interactions were treated by the Ewald method as implemented in the PME procedure.47 H450 and H455 Protonation States. QMMM calculations were performed to test the hypothesis suggested by

Figure 2. First and second zinc coordination shell residues in the active site of the h.DPP III: (a) H0−E− system and (b) H−−E0 system, before (left) and after (right) B3LYP/[6-31G(d)+LanL2DZECP] optimization.

shell and negative glutamates in the second coordination shell, and the second one, named H−−E0 (Figure 2b, left), had histidines in their anionic form and neutral glutamates. The systems were optimized using the two-layer ONIOM method as implemented in Gaussian 09.41 For this optimization the system was divided into two layers which were handled at different levels of theory. Link atoms were placed where borders “cut” covalent bonds (Cα−Cβ). The QM layer consisting of Y318, H450, E451, H455, E507, E508, E512, and H568 side chains, Leu-enkephalin, the zinc ion, and a water molecule were treated by the B3LYP DFT method using the 631G(d) basis set for the H, N, C, and O atoms and the LANL2DZ-ECP for the zinc ion. The MM part, i.e., the rest of the protein and the water molecules from the first and second solvation sphere (3501 waters), was treated using the Amber FF (parm96).56 The initial cycle of geometry optimizations was performed using mechanical embedding, and the final geometry optimization was performed using electronic embedding.57 The net charge of the whole system was −23 e, while the QM region had a +1 charge. During the QMMM geometry optimization, the protein residues and water molecules within 8 Å from the substrate molecule were allowed to move. E


Article


Table 2. Results of 40 ns Long (Unless Otherwise Specified) MD Simulations of Ligand-Free (Open, oDPP III, and Closed, cDPP III) and Peptide-Bound (Tynorphin and Leu-enkephalin, L-en) Human DPP III Using the ff14sb Force Field (Unless Otherwise Specified) and Different Models To Describe Interactions between Zinc Ion and Its Surroundingsa model 12-6LJ-type nonbondedb

12−6−4 LJ-type nonbonded CaDAc H0−E− H−−E0

M1

3-ligand 4-ligand

MD simulation

⟨CN⟩d ± SD

⟨N(H20⟩ ± SD

H450

E451

H455

E508

oDPP III (3FVY, ff12sb, 100 ns) cDPP III (3T6B, ff12sb, 100 ns) oDPP III (3FVY, ff03, 100 ns) cDPP III−L-en (3T6B, ff03, 30 ns) cDPP III − tynorphin (3T6B, ff14sb, 100 ns) cDPP III (5EGY, 20 ns)

6.46 6.59 6.12 6.89 5.66 5.91

± ± ± ± ± ±

0.55 0.67 0.33 0.32 0.51 0.28

2.15 1.91 2.00 0.00 0.00 2.00

± ± ± ± ± ±

0.57 0.33 0.07 0.07 0.00 0.03

+ ± + + ± +

M/B M/B M B M −

± ± + + + +

M/B M/B M/B B B B

cDPP III (5EGY) cDPP III − L-en (QMMM) cDPP III (5EGY) cDPP III (5EGY) cDPP III−L-en (QMMM) cDPP III (3T6B) cDPP III (5EGY) cDPP III−L-en (3T6B) cDPP III−L-en (5E3A) cDPP III (5EGY) cDPP III−L-en (5E3A) cDPP III (5EGY) oDPP III (3FVY, 200 ns) cDPP III−L-en (5E3A) cDPP III−L-en (QMMM)

4.67 5.13 4.42 4.36 4.54 5.32 5.01 5.00 5.00 5.04 5.02 5.22 5.47 4.96 5.25

± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.70 0.66 0.50 0.49 0.50 0.47 0.07 0.00 0.00 0.28 0.27 0.58 0.57 0.49 0.57

1.88 1.00 0.49 0.80 0.98 1.34 1.01 0.00 1.00 1.05 1.03 1.64 1.71 1.15 1.44

± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.33 0.01 0.52 0.46 0.13 0.47 0.12 0.04 0.01 0.31 0.28 0.48 0.46 0.36 0.50

− ± + + + + + + + + + + + + +

M/B M/B M M M/B − − − − − − − − − −

± + + + + + + + + + + + + + +

M M/B M/B −/M −/M B B B B B B M/B M/B M/B M/B

LIG

+ +

−

−

+ − −

− −

a

Monodentat (M) or bidentate (B) coordination of the zinc ion by E451 and E508 during simulation is indicated, while for histidines and ligand (substrate/inhibitor) it is indicated whether they do (+) or do not coordinate (−) with the zinc ion during the simulation. We considered 2.5 Å to be a boundary distance that determines ion coordination. ⟨N(H2O)⟩ represents the average number of water molecules found within 2.5 Å from the zinc ion. For each simulation, the initial structure used for the system setup is indicated in parentheses. bSimulations with H568 in its neutral form. Simulations performed in our previous works.37,63,64 cCaDA parameters for the zinc ion combined with parameters for neutral H450 and H455 hydrogen bonded to negatively charged glutamates E512 and E507, respectively (marked H0−E−), and negatively charged H450 and H455 hydrogen bonded to neutral glutamates E512 and E507, respectively (marked H−−E0). d⟨CN⟩ is the average zinc ion coordination number.

the average metal ion coordination number ⟨CN⟩ and the average root-mean-square deviation, ⟨RMSD-1⟩ and ⟨RMSD2⟩, values. ⟨ RMSD-1⟩ and ⟨RMSD-2⟩ are defined as the average root-mean-square deviation of all of the heavy amino acid residues’ atoms from the first and the second metal coordination sphere, respectively. RMSD values were calculated against relevant X-ray structures or, in case of the closed h.DPP III, QMMM-optimized structure.

MM/PBSA. The substrate binding free energies were approximated by MM/PBSA energies calculated within the Amber16 suit. For each complex the MM/PBSA energies have been calculated on a set of 5 ns long intervals (50 evenly spaced snapshots per interval) sampled throughout simulations using a single trajectory approximation. The calculations were accomplished for the enzyme with a dielectric constant of 1.0 immersed into the solvent with a dielectric constant 80. The ion concentration was 0.15 M. The polar component of the solvation enthalpy was calculated using the AMBER built-in PBSA solver and atom-type/charge-based radii by Tan and Luo.58 The nonpolar component was determined by ΔHnonpol = γSASA + β, where the solvent accessible surface area (SASA) was calculated with the MolSurf program.59 The surface tension γ and the offset β were set to the standard values of 0.0378 kcal (mol Å2)−1 and −0.5692 kcal mol−1, respectively. As suggested by Maffucci and Contini,60 20 water molecules closest to the substrate molecule defining the ligand hydration shell were considered as part of the protein structure for the energy calculation. The number of 20 water molecules was based on the number of hydrogen bonds that the substrate molecule forms with water molecules in the enzyme binding site throughout the trajectory. Unless differently specified, default parameters were used.

■

RESULTS AND DISCUSSION A. Protonation State of Imidazoles Coordinated to the Zinc Ion in h.DPP III. The results of the QMMM optimization of the h.DPP III−Leu-enkephalin complexes with different protonation states of H450-E512 and H455-E507 hydrogen-bonded pairs, (H0−E− and H−−E0) clearly showed that both histidines coordinated to the central zinc ion are in their neutral form (see Figure 2a−b, right). While the QMMM-optimized structure (utilizing charge embedding) of the H0−E− system does not differ significantly from the starting structure, this is not so for the H−−E0 system (Figure 2b): the protons moved from glutamates to histidines during the electronic embedding optimization; i.e., the H−−E0 system transformed into the H0−E− system. B. Testing the Available Zinc Ion Parameters Using h.DPP III as a Test System. 12-6 LJ-Type Nonbonded Approach. Previous MD simulations utilizing only the electrostatic plus 12-6 LJ-type nonbonded parameters to describe interactions between the metal ion and its surroundings (some of them presented in Table 2) showed

■

STRUCTURE EVALUATION In order to quantify how good the parameters are for simulations of different zinc-dependent aminopeptidases, besides the average metal ion −ligands distances we utilized F


Article


Table 3. Results of 5 ns Long MD Simulations of Ligand-Free (Open, oDPP III, and Closed, cDPP III) and Leu-enkephalinBound Human DPP III Using the Drude-2013 Polarizable Force Fielda (⟨d⟩ ± SD)/Å

oDPP III (3FVY) cDPP III (5EGY) cDPP III−L-en (QMMM)

⟨CN⟩b ± SD

± SD

H450NE2− ZN

H455NE2− ZN

E508OE1−ZN

E508OE2−ZN

E451OE1−ZN

E451OE2−ZN

Osc−ZN

4.04 ± 0.19 4.06 ± 0.25 4.15 ± 0.36

1.00 ± 0.00 1.00 ± 0.00 1.00 ± 0.00

1.88 ± 0.03 1.89 ± 0.04 1.89 ± 0.04

1.88 ± 0.03 1.88 ± 0.03 1.89 ± 0.04

1.93 ± 0.06 2.39 ± 0.53 2.18 ± 0.44

2.98 ± 0.26 2.46 ± 0.53 2.63 ± 0.53

5.95 ± 0.89 6.01 ± 0.72 5.60 ± 0.73

5.91 ± 0.82 5.76 ± 1.03 5.61 ± 1.52

6.18 ± 0.90

a 2.5 Å was used as the boundary distance for zinc ion coordination. ⟨N(H20)⟩ is the average number of water molecules found within 2.5 Å from the zinc ion. For each simulation, the initial structure used for the system setup is indicated in parentheses. b⟨CN⟩ is the average zinc ion coordination number. cThe average distance (and standard deviation) between the ligand Os atom (carbonyl oxygen atom from the second peptide bond from the substrate N-terminus) and the zinc ion.

that the choice of the Amber force field does not influence the Zn2+ preference for higher (6 and above) coordination numbers during MD simulations. As a matter of fact, we performed a series of MD simulations for DPP III orthologues using either ff03, ff12SB, or ff14SB force field,22−24,28,33−37,61−63 and the results were consistent indicating that the zinc ion is beside by H450, H455, and E508, coordinated by E451 and two waters or substrate as well.37,63,64 Such a strong preference for the mostly octahedrally coordinated zinc is probably boosted by the large positive charge of the zinc ion (+2 e). In addition, MD simulations of the h.DPP III complexes with natural peptides revealed a strong coordination of the zinc ion by peptides throughout the simulation, thus preventing a water molecule from entering the zinc ion coordination sphere. This is in contrast with the theoretical QMMM calculations of the h.DPP III−Leu-enkephalin complex, which showed that the catalytic mechanism of this enzyme relies on the presence of a nucleophilic (zinc-coordinated) water to attack the carbonyl carbon of the substrate that is not coordinated to the zinc ion.28 12-6-4 LJ-Type Nonbonded Approach. MD simulations (20 ns long) of the closed h.DPP III structure performed with the 12-6-4 LJ-type nonbonded model also indicated a hexacoordinated zinc ion during the simulation, where both histidines, E508 (bidentately), and two water molecules coordinated the zinc ion throughout the simulations (Table 2). The average distances between the metal ion and its first coordination shell ligands determined during the MD simulation are, mostly, slightly smaller than those determined in either the X-ray or the QMMM-optimized h.DPP III structures (Table 1), with a slightly larger H450NE2−ZN average distance being an exception (Table S3). Importantly, the catalytically important glutamate residue E451 was not coordinated to the zinc ion during the simulations; rather, it was hydrogen bonded with one of the water molecules coordinated to the zinc ion. Apparently, the 12-6-4 LJ-type of parameters is better than the 12-6 one. CaDA Approach. In order to enable a free exchange of ligands in the metal coordination sphere, we utilized the zinc parameters available within the CaDA approach for MD simulations. Treating the zinc ion as a tetrahedron-shaped zinc divalent cation and histidines H450 and H455 as neutral residues resulted in a significant reduction of the average CN values (see Table 2, H0−E− section). Precisely, ⟨CN⟩ was 4.67 during MD simulations of the closed, ligand-free h.DPP III and 5.13 during the simulations of its complex with Leu-enkephalin (Table 2, H0−E− section). However, the low CN values obtained during the MD simulations with CaDA parameters

are mostly a consequence of the rather weak interactions between the zinc ion and histidines: H450 and H455 escaped from the zinc ion coordination sphere at the very beginning and after 2 ns, respectively, of the MD simulations of h.DPP III. The average ZN−H450NE2 distance of 2.45 ± 0.11 Å obtained during the h.DPP III−Leu-enkephalin MD simulations indicates the partial coordination of the metal ion by this histidine (see Tables 2 and S3, H0−E− section). Moreover, the catalytically important glutamate residue (E451) was coordinated to the zinc ion throughout the simulations, either monodentately or bidentately. The average ZN−Os distance (Os is a carbonyl oxygen from the second Gly residue from the substrate N terminus) of 4.22 ± 0.27 Å indicates that Leuenkephalin was not coordinated to the zinc ion during the MD simulations of the complex. E508 was monodentately coordinated to the zinc ion during the MD simulations of the ligand-free h.DPP III, while an interconversion between monodentate and bidentate coordination was observed in the case of the h.DPP III−Leu-enkephalin complex simulation. Although QMMM calculations have shown that both histidines coordinating to the central zinc ion are in their neutral form in h.DPP III, following the suggestion of Pang et al.17 we have also performed simulations with H450 and H455 negatively charged (results are given in Tables 2 and S3, H−− E0 section). During these simulations, both histidines were tightly bound to the zinc ion, and the average H450NE2−ZN and H455NE2−ZN distances (Table S3) are in agreement with those presented in Table 1. A further reduction of the zinc ion ⟨CN⟩ (4.42 and 4.36 in the case of two distinct ligand-free h.DPP III MD simulations and 4.54 in the case of the h.DPP III−Leu-enkephalin complex MD simulations, Table 2) is mostly the result of less water molecules that participate in the metal ion coordination (in average 0.49−0.98 water molecules) and a weaker ZN-E508 coordination bond during the simulations utilizing the H−−E0 parameters (Table 2). After 2 and 20 ns E508 escaped from the zinc ion coordination sphere during one of the MD simulations of the ligand-free h.DPP III as well as during the MD simulation of the h.DPP III−Leuenkephalin complex, respectively. E451 was monodentately coordinated to the zinc ion, with the carboxyl atom exchanging in coordination, during both MD simulations of the ligand-free h.DPP III with H−−E0 parameters (Tables 2 and S3). Similarly to the H0−E− simulations, the substrate molecule was not coordinated to the zinc ion during the H−−E0 simulation. However, the slightly larger average ZN−Os distance (4.84 ± 0.34 Å, see Table S3) was determined during the latter simulation. G


Article


Despite constant improvements in the applicability and tractability of simulations with polarizable potentials, there is a still a significant difference in the achievable simulation times, especially in the case of biomacromolecular systems. For example, for the closed ligand-free h.DPP III MD simulations utilizing the polarizable potential (altogether 77 645 atoms) on an Intel Xeon CPU E5-2630 v3 @ 2.40 GHz with a GeForce GTX 1070 GPU card, ∼0.5 ns/day of simulation time was achieved, while for the same system MD simulations utilizing a nonpolarizable potential (4-ligand-extended model parameters and altogether 72 610 atoms) an Intel i7-6700K CPU @ 4.00 GHz with a GeForce GTX 1070 GPU card generated ∼17 ns/ day of simulation time (using the same time step). Apparently, for studying the long-range properties of biomacromolecules as well as for comparative studies of a large number of protein complexes with different ligands, the nonpolarizable MD simulations are still a more convenient choice. C. Development of the New Zinc Ion Parameters for Peptidases with a Zinc Ion Coordination Similar to That in h.DPP III. 3-Ligand Model Parameters. We have applied the parametrization procedure reported by Yang et al.19 to h.DPP III and developed our own parameters, which we named 3-ligand model parameters. However, the zinc ion coordination geometry sampled during the simulation with these parameters only insignificantly differs from the geometries sampled with the M1 parameter MD simulations (see Tables 2 and S3). Again, the ligands of the pentacoordinated zinc ion were histidines, water, and E508 (bidentate) with the slightly larger average E508OE1−ZN distance than in the M1 model parameter simulations (see Table S3). Like in the M1 model parameter simulations, the zinc-coordinated water molecule was hydrogen bonded to E451 throughout the simulations (data not shown). 4-Ligand Model Parameters. In order to account for the water influence on the charge distribution around the metal ion, we included the zinc coordinating water molecule into the model and developed a new set of hybrid bonded/nonbonded parameters, named 4-ligand model parameters. As shown in Table 2, during the MD simulations with these 4-ligand model parameters the average value of CN (ranging from 4.96 to 5.47) was slightly larger than in the simulations with the M1 and the 3-ligand model parameters. The reason for this is the increased number of water molecules taking part in the zinc ion coordination (with the average number ranging from 1.15 to 1.71, Table 2). The average ZN−E508OE2 distance is larger than those determined in previous simulations, which indicates that OE2 is only partially coordinated to the zinc ion (see Table S3 and Figure S2). Like in the experimentally determined structures of the closed h.DPP III and in the QMMM-optimized Michaelis complex structure, 28 this E508OE2 oxygen atom (i.e., the oxygen not directly coordinated to the zinc ion) is hydrogen bonded to the Y318 hydroxide group during the MD simulations of the h.DPP III−Leuenkephalin complex with the 4-ligand model parameters. During the simulations of the ligand-free, open and closed, h.DPP III, E508 preferably interacted either with the H568 imidazole NH group or with the other water molecules found in the active site (see Figure S2). It should be noted that both zinc ion-coordinated histidines have the same partial charges (averaged over those calculated for each of them, H450 and H455) during the MD simulation. 4-Ligand-Extended Model Parameters. In view of the importance of the second coordination shell ligands for both

To summarize, the simulations utilizing the CaDA type of parameters, either H0−E− or H−−E0, were not able to reproduce the accurate h.DPP III active site structure. Hybrid Bonded/Nonbonded Approach. Using the so-called M1 model parameters in simulations of the ligand-free h.DPP III and its complexes, the average value of CN in a range from 5.0 to 5.32 (Table 2) has been obtained. During these simulations, the zinc ion was coordinated by both histidines, by E508, and either by water or the substrate molecule. The catalytically important glutamate (E451) has not participated in the zinc ion coordination during these simulations (see average E451−ZN distances in Table S3). Moreover, it has been hydrogen bonded to the water molecule coordinating to the metal ion in accord with the proposed reaction mechanism. However, a stable bidentate coordination of the zinc ion by E508 was not observed either in the QMMM-optimized Michaelis complex structure28 or in the crystal structures (see Table 1). MD Simulations with Polarizable Potential. The effect of charge polarization on the active site architecture was examined by performing 5 ns long MD simulations with the polarizable Drude-2013 potential. The simulations revealed a stable tetrahedral coordination geometry (⟨CN⟩ ranging from 4.04 to 4.15, Table 3) of the zinc ion made up by H450, H455, E508, and a water molecule. Coordination bonds to the histidine residues (H450 and H455) were slightly shorter than those in the reference structures (Tables 1 and 3), and E508 was monodentately coordinated to the zinc ion, only occasionally exchanging its carboxy oxygens in the metal coordination (indicated by large SD values in Table 3) during the simulation. However, the hydrogen bond between E451 and the zinc-coordinated water molecule, necessary for the catalysis, existed only sporadically. Rather, E451 formed hydrogen bonds with amino acids from the lower protein domain: with N406 and occasionally with S408 or N391 during the ligand-free h.DPP III MD simulations and with N391 and N406 during the cDPP III−Leu-enkephalin complex MD simulations (see Figure S1). As a result, slightly larger ⟨ ZN−E451⟩ distances were obtained (Table 3) than in the MD simulations with the hybrid parameters. It should be emphasized that the hydrogen bonds mentioned do not exist in the experimentally determined structures of h.DPP III and have not been noticed during the QMMM calculations of the reaction pathway. In addition, the average ZN−substrate distance measured during the simulations with the polarizable potential (6.18 ± 0.90 Å, Table 3) indicated that the substrate molecule was located a little further from the metal ion than in the crystal and QMMM-optimized structures (Table 1). The reason is a breakage of the hydrogen bonds between Leuenkephaline and the residues from the lower protein domain (i.e., N391, H568, and Y318) during the MD simulations with the polarizable potential. As a consequence, the substrate moved slightly deeper into the enzyme binding site and established a new hydrogen bond between the second amide group from the N-terminus and E508. However, one should be aware that the integration scheme, the nonbonded cutoff, as well as the temperature and the pressure coupling schemes used during the simulations with the polarizable potential differ from those used during the simulations with the Amber force fields, so they cannot be directly compared. Also the MD simulations with the polarizable potential were relatively short, so it might be that the equilibrium has not been reached during the simulation time. H


1.87 ± 0.35 1.92 ± 0.28 1.28 ± 0.45 1.00 ± 0.03 0.00 ± 0.01 1.52 ± 0.52

⟨CN⟩b ± SD

4.86 ± 0.35 4.91 ± 0.29 4.28 ± 0.45

4.01 ± 0.11 4.59 ± 0.50

4.55 ± 0.51

2.04 ± 0.05

2.04 ± 0.05 2.05 ± 0.05

2.06 ± 0.05 2.05 ± 0.05 2.04 ± 0.05

H450NE2− ZN

2.04 ± 0.05

2.03 ± 0.05 2.05 ± 0.05

2.06 ± 0.05 2.06 ± 0.05 2.04 ± 0.05

H455NE2− ZN

1.94 ± 0.04

1.93 ± 0.04 1.95 ± 0.04

1.94 ± 0.04 1.94 ± 0.04 1.93 ± 0.04

E508OE1−ZN

3.01 ± 0.13

2.96 ± 0.09 2.99 ± 0.09

3.16 ± 0.10 3.17 ± 0.12 3.02 ± 0.15

E508OE2−ZN

(⟨d⟩ ± SD)/Å

4.45 ± 0.32

4.35 ± 0.24 2.77 ± 0.74

5.07 ± 0.79 4.39 ± 0.23 4.30 ± 0.25

E451OE1−ZN

4.37 ± 0.22

4.23 ± 0.20 4.01 ± 0.70

4.67 ± 0.81 4.59 ± 0.26 4.35 ± 0.29

E451OE2−ZN

3.82 ± 0.61

3.93 ± 0.67 2.25 ± 0.15

4.27 ± 0.64

Osd−ZN

0.44 ± 0.09

0.41 ± 0.07 0.43 ± 0.05

0.93 ± 0.10 0.62 ± 0.06 0.53 ± 0.15

⟨RMSD-1⟩ ± SD

0.65 ± 0.06 (1.16 ± 0.14)

0.78 ± 0.04 (0.88 ± 0.09) 0.63 ± 0.05 (0.88 ± 0.14)

1.00 ± 0.11 (2.33 ± 0.49) 0.75 ± 0.04 (1.46 ± 0.21) 0.63 ± 0.08 (0.66 ± 0.09)

⟨RMSD-2⟩ ± SDe

a

2.5 Å was used as the boundary distance for zinc ion coordination. ⟨N(H20)⟩ is the average number of water molecules found within 2.5 Å from the zinc ion. For each simulation, the initial structure used for the system setup is indicated in parentheses. The average root-mean-square deviation of all of the heavy amino acid residues’ atoms in the first (H450, H455, and E508 residues) and the second (H450, E451, H455, E507, E508, E512 residues) coordination sphere is represented by ⟨RMSD-1⟩ and ⟨RMSD-2⟩ values, respectively. RMSD values were calculated against relevant, X-ray (pdb code: 3FVY), or QMMM-optimized structure28 in the case of the open DPP III and all other structures, respectively. b⟨CN⟩ is the average zinc ion coordination number. cThe initial cDPP III−RRNA1 complex does not contain the water molecule coordinated to the zinc ion, and the cDPP III−RRNA2 does. dThe average distance (and standard deviation) between the ligand Os atom (carbonyl oxygen atom from the second peptide bond from the substrate N-terminus) and the zinc ion. eIn parentheses average RMSD values obtained when Y318 was considered as a part of the zinc ion second coordination sphere.

oDPP III (3FVY, 220 ns) cDPP III (5EGY) cDPP III−L-en (QMMM) cDPP III−L-en (5E3A) cDPP III−RRNA1c (3T6B) cDPP III−RRNA2c (3T6B)

⟨N(H2O)⟩ ± SD

Table 4. Results of 100 ns Long (Unless Otherwise Specified) MD Simulations of Ligand-Free (Open, oDPP III, and Closed, cDPP III) and Ligand-Bound (Leu-enkephalin, L-en, and Arg-Arg-2-naphthylamide, RRNA) Human DPP III Using the ff14sb Force Field and Parameters Obtained for 4-Ligand-Extended Model of Zinc Coordination Spherea

Journal of Chemical Information and Modeling Article

I


Article


RMSD values for the metal first and second coordination spheres against relevant (QMMM-optimized28 or X-ray, PDB code: 3FVY) structures calculated for each MD trajectory (Figure S4) showed that the geometries of all metal ion coordination spheres were stable throughout the MD simulation with average RMSD-1 and RMSD-2 values for the closed enzyme form (0.412−0.618 and 0.626−0.780 Å in Table 4, respectively) slightly smaller than for the open form (0.932 and 0.998 Å in Table 4, respectively).The obtained results indicated the 4-ligand-extended model parameters as the most promising fixed-charge metal parameters for simulations of h.DPP III. D. Further Testing of the New 4-Ligand-Extended Model Zinc Ion Parameters. In order to examine performances of the 4-ligand-extended model parameters in detail, MD simulations of the systems with h.DPP III were additionally analyzed. Further, the parameters were tested on a series of other enzymes with His-His-Glu zinc coordination. Effect on Ligand Binding Affinities and Protein Dynamics. MM/PBSA calculations were used to compare the previous results, obtained utilizing the nonbonding, 12-6 LJ, parameters with the results obtained utilizing the new 4-ligand-extended model zinc ion parameters. The enthalpy contributions to the binding free energies calculated for Leu-enkephalin and ArgArg-2-naphthylamide bound to the closed h.DPP III structure during MD simulations are given in Figure 3. As indicated,

the charge distribution and positioning of the ligands in the first coordination shell, we developed a new 4-ligand-extended model parameters, where the second coordination shell ligands Y318, E451, E507, and E512 were also included in the parameter development. Their positions were optimized in the small model (see Computational Methods section Hybrid Bonded/Nonbonded Models) and their Amber FF partial charges were included into the RESP charge fitting procedure (Table S2, a detailed overview of the 4-ligand-extended model parameters is given in the the Supporting Information). In the 3-ligand and 4-ligand models the E508OE1 partial charge is 0.17 e and 0.2 e, respectively, more negative than the E508OE2 one. However, in the 4-ligand-extended model the oxygen bonded to the zinc ion (E508OE1) is only slightly more negative (0.03e) than the other one (E508OE2). The reason is probably the presence of Y318, hydrogen-bonded to E508OE2, in the RESP fitting model. This hydrogen bond remained preserved during the MD simulations of the h.DPP III−substrate complexes in agreement with the experimentally determined structures of the closed h.DPP III and QMMM-optimized Michaelis complex structure.28 According to the charge values obtained for different initial models (Table S2), it is clear that the RESPA charges are sensitive to the initial coordinates. It was shown that mutual relationship of the charges in a parameter model seems to be more important than differences of their values among the models. The slightly higher charge value at Glu508 OE2 than at Glu508 OE1 in the 4-ligandextended set of parameters resulted with proper positioning of Glu508 and enabled establishment of the catalytically active enzyme conformation. Like in the simulations with the 4ligand model parameters, the partial atomic charges at the zinccoordinated histidines were the same. The zinc ion was coordinated by three protein residues, i.e., H450, H455, and E508, while a water molecule mostly completed the tetrahedral coordination of the metal ion during the simulations (⟨CN⟩, ⟨NH20⟩, and the relevant ZN−ligand distances are given in Table 4). Entrance of the second water molecule into the zinc ion coordination sphere resulted in an exchange between the tetra- and pentacoordinated metal ion. Interestingly, the average coordination numbers determined during the simulations of the ligand-free h.DPP III structures (both open and closed, Table 4) are larger than those determined during the MD simulations of the h.DPP III complexes (with either Leuenkephalin or RRNA). Apparently, the substrate presence in the enzyme active site hinders water molecules to approach the zinc ion (Figure S3). However, when the zinc-coordinated water molecule was present in the initial systems, like in the cDPP III−RRNA2 and cDPP III−Leu-enkephalin complexes (in the accord with the QMMM-optimized Michaelis complex structure28), it remained coordinated to Zn2+, and the RRNA molecule did not enter the zinc coordination sphere during the MD simulations making the hydrolytic reaction feasible (Figure S3). The mean coordination numbers range from 4.01 to 4.90 (Table 4) and correspond to those in experimentally determined and QMMM-generated h.DPP III structures. Although the substrates were stabilized by a number of hydrogen bonds during all MD simulations (Tables S4 and S5), their binding was not productive in all complexes since the binding mode, in which the zinc-bound water molecule is not present, is considered to be nonproductive.28,65 Such a zinc coordination was observed in the h.DPP III−endomorphin-2 crystal structure (PDB code 5EHH).

Figure 3. MM/PBSA binding energies calculated for the closed DPP III complexes with Arg-Arg-2-naphthylamide (RRNA) and Leuenkephalin (L-en) on 5 ns long intervals during MD simulations utilizing the 4-ligand-extended model parameters for the zinc coordination shell. Average values calculated on 50 frames per interval are given. The 20 water molecules closest to the substrate molecule were considered as a part of the protein structure.

binding of the peptide substrate is less exothermic than binding of RRNA. This result agrees with the trend in binding free energies obtained during MD simulations with 12-6 LJ-type nonbonding parameters (see Table 1 in Tomić and Tomić;64 the cWT-RRNA and cWT-L-EN structures are equivalent to the structures used in this paper). During 220 ns long MD simulation of the open ligand-free h.DPP III, the protein radius of gyration (Rg) significantly decreased (min and max values are 25.66 and 28.12 Å, respectively; see Figure S5 upper left) resulting with the protein closure. The MD simulations started from the closed h.DPP III forms (ligand-free or ligand-bound) did not reveal J


Article


Table 5. Average Distances between the Zinc Ion and Its Nearest Residues in Crystal Structures of the Enzymes Used To Test the Newly Derived 4-Ligand-Extended Zinc Ion Parameters: DPP III Orthologues (Other than Human) Available in the Protein Data Bank, Thermolysin, Neprilysin, and Aminopeptidase N (⟨d⟩ ± SD)/Å zinc peptidase

H1−ZN

DPP III orthologues thermolysin neprilysin aminopeptidase N

a

2.05 2.04 1.98 2.06

± ± ± ±

0.02 0.06 0.09 0.06

H2−ZN 2.03 2.05 2.04 2.09

± ± ± ±

0.06 0.06 0.06 0.07

EOE1−ZN 2.03 2.04 1.98 1.99

± ± ± ±

0.09 0.15 0.07 0.07

EOE2−ZN 2.90 2.86 2.70 2.82

± ± ± ±

0.56 0.26 0.30 0.17

Lig(O)−ZNb

Lig(N)−ZNb

Lig(S)-ZNb

O1water−ZNc

O2water−ZNd

± ± ± ±

2.15 ± 0.00 − − −

− 2.34 ± 0.05 2.11 ± 0.08 −

2.15 ± 0.12 2.10 ± 0.26 − 2.14 ± 0.06

2.31 ± 0.25 2.45 ± 0.21 − 2.61 ± 0.01

2.13 2.01 2.00 2.03

0.00 0.12 0.12 0.16

a Other than human. Structures with low resolution (above 3 Å) are excluded from the analysis, see Table S1. bAveraged over the distances with ligands closest to the zinc ion (if they are in the zinc coordination sphere) for the specify element, see Table S1. cAveraged over the distances to the closest water molecules coordinating the zinc ion, see Table S1. dAveraged over the distances to the second closest water molecule coordinating the zinc ion, see Table S1.

Table 6. Results of 150 ns Long MD Simulations of Several Ligand-Free DPP III Orthologues (Yeast DPP III, yDPP III; Bacteroides thetaiotaomicron DPP III, BtDPP III; Caldithrix abyssi DPP III, CaDPP III, and Mushroom Armillaria tabescens, AtDPP III) and Ligand-Bound (Arg-Arg-2-naphthylamide, RRNA) Bacteroides thetaiotaomicron DPP III, BtDPP III-RRNA Using the ff14sb Force Field and 4-Ligand-Extended Model Parameters of Zinc Coordination Spherea (⟨d⟩ ± SD)/Å DPP III orthologues yDPP III (3CSK) BtDPP III (5NA6) CaDPP III (6EOM) AtDPP III (5YFB) BtDPP IIIRRNA (5NA6)

run 1 2 1 2 1 2 1 2 1 2

⟨CN⟩b ± SD 4.59 4.64 4.55 4.53 4.80 4.75 4.87 4.80 4.66 4.56

± ± ± ± ± ± ± ± ± ±

0.52 0.51 0.57 0.56 0.49 0.52 0.39 0.39 0.49 0.52

⟨N(H2O)⟩ ± SD 1.54 1.59 1.55 1.53 1.80 1.75 1.87 1.80 1.66 1.56

± ± ± ± ± ± ± ± ± ±

0.52 0.51 0.57 0.56 0.49 0.52 0.39 0.39 0.49 0.52

H1−ZN 2.05 2.06 2.06 2.06 2.06 2.06 2.09 2.04 2.05 2.06

± ± ± ± ± ± ± ± ± ±

0.05 0.05 0.04 0.05 0.04 0.05 0.05 0.05 0.05 0.05

H2−ZN 2.05 2.05 2.04 2.04 2.04 2.05 2.07 2.03 2.05 1.95

± ± ± ± ± ± ± ± ± ±

0.05 0.05 0.05 0.05 0.04 0.05 0.05 0.05 0.05 0.05

EOE1−ZN 1.94 1.95 1.95 1.95 1.96 1.96 1.98 1.93 1.95 1.95

± ± ± ± ± ± ± ± ± ±

0.04 0.05 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04

EOE2−ZN 3.11 3.19 3.03 3.02 3.24 3.21 3.60 2.99 3.05 3.00

± ± ± ± ± ± ± ± ± ±

0.14 0.24 0.09 0.09 0.13 0.15 0.15 0.10 0.10 0.09

Osc−ZN

6.02 ± 0.92 7.87 ± 2.61

(⟨RMSD-1⟩ ± SD)/Å 0.73 0.76 0.42 0.42 0.64 0.60 0.76 0.66 0.76 0.66

± ± ± ± ± ± ± ± ± ±

0.18 0.14 0.06 0.06 0.08 0.10 0.08 0.04 0.08 0.04

(⟨RMSD-2⟩ ± SD)/Å 0.90 1.22 0.73 0.71 0.82 0.83 0.96 0.76 0.96 0.76

± ± ± ± ± ± ± ± ± ±

0.16 0.19 0.15 0.15 0.07 0.08 0.09 0.08 0.09 0.08

2.5 Å was used as the boundary distance for zinc ion coordination ⟨CN⟩ while ⟨N(H20)⟩ represents the average number of water molecules found within 2.5 Å from the zinc ion. For each simulation, the initial structure used for the system setup is indicated in parentheses. ⟨RMSD-1⟩ and ⟨RMSD-2⟩ represent the average root-mean-square deviation of all of the heavy amino acid residues’ atoms in the first and the second coordination sphere, respectively. The first zinc ion coordination sphere consists of the following residues: H460, H465, and E517 in yDPP III; H448, H453, and E476 in BtDPP III and BtDPP III-RRNA; H379, H383, and E412 in CaDPP III; and H443, H448, and E501 in the AtDPP III. The second zinc ion coordination sphere consists of the residues: H460, E461, H465, E516, E517, and E520 in yDPP III; H448, E449, H453, E475, E476, and D480 in BtDPP III and BtDPP III-RRNA; H379, E380, H383, E411, E412, and D416 in CaDPP III; and H443, E444, H448, E500, E501, and E505 in AtDPP III structure. RMSD values were calculated against the relevant crystal structures. b⟨CN⟩ is the average zinc ion coordination number. c Average distance (and standard deviation) between the ligand Os atom (carbonyl oxygen atom from the second peptide bond from the substrate N-terminus) and the zinc ion. a

major changes of the initial Rg values, indicating that the protein globularity has not changed during these simulations (Figure S5 upper left). A similar behavior of the protein was observed previously during MD simulations with the 12-6 LJtype nonbonded parameters for the zinc ion (see Figures S6 and S7).24,37 To summarize, the 4-ligand-extended model parameters enabled sampling of the structures of both the protein and its complexes, with the active site architecture close to real ones preserving the long-range molecular motions. MD Simulations of Other Peptidases. In order to find out if the newly developed hybrid bonded/nonbonded, 4-ligandextended model parameter is also appropriate for MD simulations of the other peptidases (containing the zinc ion coordinated with two His and one Glu residues), we studied the four other DPP III orthologues as well as thermolysin, neprilysin, and aminopeptidase N. The average distances between the zinc ion and its ligands in the structures of these enzymes available in the Protein Data Bank are given in Table 5 (see also Table S1 for details). These values and the average

values of RMSD-1 and RMSD-2 determined during MD simulations are considered as referent when evaluating the simulation results obtained for those peptidases. The simulation results are summarized in Tables 6 and 7. The distances determined during these MD simulations are close to those in the experimental structures, with discrepancies within the range of standard deviations. The glutamate from the first zinc coordination sphere coordinates the zinc ion mostly monodentately in both crystal and simulated structures. During the simulations the zinc ion was predominantly tetraand pentacoordinated (⟨CN⟩ ranging from 4.55 to 4.87 for DPPs III and from 4.47 to 5.25 for the other three peptidases) with one or two water molecules present in its first coordination sphere (⟨N(H2O)⟩ from 1.53 to 1.87 for DPPs III and from 1.47 to 2.25 for the other three peptidases). The exception is the aminopeptidase N in the complex with Lmethionine, where the zinc-coordinated water leaves the zinc coordination sphere and the bidentate substrate coordination of the metal ion is established during the MD simulation. The highest metal ion coordination was determined for thermolysin K


0.53 1.68 1.14 0.55 0.56 0.64 0.52 0.91 0.57 0.4

± ± ± ± ± ± ± ± ± ±

2.25 2.08 1.67 1.63 1.49 1.59 1.58 1.47 1.79 0.08

5.25 ± 1.24 5.17 ± 1.72 4.71 ± 1.13 4.63 ± 0.55 4.49 ± 0.56 4.59 ± 0.64 4.58 ± 0.52 4.47 ± 0.91 4.8 ± 0.54 4.93 ± 0.44 2.06 2.05 2.06 2.05 2.05 2.05 2.05 2.06 2.06 2.05

± ± ± ± ± ± ± ± ± ± 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05

H1−ZN 2.05 2.04 2.05 2.03 2.03 2.03 2.03 2.03 2.05 2.05

± ± ± ± ± ± ± ± ± ± 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04 0.05

H2−ZN 3.05 3.02 3,14 3.05 3.05 3.05 3.04 2.99 3.12 3.05

± ± ± ± ± ± ± ± ± ± 0.10 0.05 0,17 0.10 0.10 0.10 0.09 0.10 0.12 0.09

EOE1−ZN 1.94 1.94 1.94 1.94 1.94 1.94 1.94 1.94 1.94 1.95

± ± ± ± ± ± ± ± ± ± 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04

EOE2−ZN

2.19 ± 0.12 2.41 ± 0.59

Osc−ZN 0.69 0.56 0.53 0.45 0.49 0.45 0.41 0.44 0.83 0.66

± ± ± ± ± ± ± ± ± ± 0.11 0.11 0.07 0.06 0.07 0.08 0.06 0.05 0.16 0.06

( ± SD)/Å 0.89 0.75 0.70 0.63 0.63 0.53 0.56 1.06 1.26 0.96

± ± ± ± ± ± ± ± ± ±

0.13 0.16 0.11 0.07 0.13 0.09 0.08 0.12 0.17 0.30

( ± SD)/Å

a

2.5 Å was used as the boundary distance for zinc ion coordination. ⟨N(H20)⟩ is the average number of water molecules found within 2.5 Å from the zinc ion. For each simulation, the initial structure used for the system setup is indicated in parentheses (and the simulation time). ⟨RMSD-1⟩ and ⟨RMSD-2⟩ represent the average root-mean-square deviation of all of the heavy amino acid residues’ atoms in the first and the second coordination sphere, respectively. The first zinc ion coordination sphere consists of heavy atoms of the following residues: H142, H146, and E166 in the thermolysin; H297, H301, and E320 in the aminopeptidase N; and H584, H588, and E647 in the neprilysin. The second zinc ion coordination sphere consists of heavy atoms of the following: H142, H146, E166, D170, D150, G162, N165, and D138 in the thermolysin; H297, H301, E320, D327, V324, T304, L316, N306, and K319 in the aminopeptidase N; and H584, H588, E647, D591, D651, M580, and T644 in the neprilysin. RMSD values were calculated against relevant crystal structures. b⟨CN⟩ is the average zinc ion coordination number. cAverage distance (and standard deviation) between the ligand Os atom (carbonyl oxygen atom from the second peptide bond from the substrate N-terminus) and the zinc ion. dSimulations of the protein without methionine in the enzyme active site. eSimulations of the protein with methionine bound.

thermolysin (5a3y) 150 ns thermolysin (5a3y) 150 ns thermolysin (5a3y) 150 ns neprilysin (4zr5) chain A 100 ns chain B neprilysin (4zr5) 150 ns chain A chain B aminopeptidase Nd (4xn4) 120 ns aminopeptidase Nd (4xn4) 150 ns aminopeptidase Ne (4xn4) 110 ns

⟨N(H2O)⟩ ± SD

⟨CN⟩ ± SD b

( ± SD)/Å

Table 7. Results of MD Simulations of Thermolisyn, Neprilysin, and Aminopeptidase N (PDB codes 5A3Y, 4ZR5, and 4XN4, Respectively) Using the ff14sb Force Field and 4-Ligand-Extended Model Parameters of Zinca

Journal of Chemical Information and Modeling Article

L


Article


Figure 4. Overlay of the first and the second zinc coordination shell in crystal (carbon and zinc atoms are colored green) and optimized structures (carbon and zinc atoms are colored cyan) obtained after 150 ns of MD simulation of (a) y.DPP III, (b) At.DPP III, (c) Ca.DPP III, and (d) Bt.DPP III. Water molecules in crystal structures shown as red spheres only. Backbone and hydrogen atoms are not shown for clarity.

(⟨5.04⟩ ± 1.39, the value averaged during three independent MD simulations) where mostly two, and occasionally even three, water molecules can be considered as zinc ion ligands. However, the pentacoordinated zinc ion is also quite frequent in the crystal structures of thermolysin (in 43 of 167 PDB structures). It should be noted that the ligand molecules are, in general, closer to the zinc ion in the structures from the test set than in h.DPP III in accord with the experimental crystal structures (Tables 5−7). The conformation of both the first and the second coordination spheres remained conserved during the MD simulations (Figure S8 and S9), with the average RMSD values ranging from 0.41 to 0.83 Å (RMSD-1) and from 0.53 to 1.22 Å (RMSD-2), respectively (see Tables 6 and 7 and Figures 4 and 5). These RMSD values are similar to those determined for the systems with h.DPPs III and depend on the initial structure compactness. The range and type of the long-range motions characteristic for DPP III orthologues significantly

depended on the initial structure; if it was open, like the yeast DPP III (PDB code 3CSK) and human DPP III (PDB code 3FVY), domains have approached each other during the MD simulation, resulting in a decrease of the proteins’ radius of gyration (see Figures S5 and S6). However, MD simulations starting from the more compact structures revealed either very limited extension of the domain movement (Ca.DPP III, At.DPP III, and Bt.DPP III in the complex with RRNA) or subsequent opening and closing (ligand-free Bt.DPP III). These results are in accord with the previously performed simulations with the nonbonding zinc parameters (Figure S5) and in agreement with the experimental data.33−35,62 Nevertheless, the relative orientation of the amino acid residues important for the substrate stabilization and the enzyme activity was preserved during the simulations, as can be seen from the overlap of the zinc ion coordination spheres (the first and the second) in the initial (crystal) and final (MD simulated) structures (Figures 4 and 5). M


Article


Figure 5. Overlay of the first and the second zinc coordination shell in crystal (carbon and zinc atoms are colored green) and optimized structures (carbon and zinc atoms are colored cyan) obtained after 150 ns of MD simulation of (a) thermolysin, (b) neprilysin, and (c) aminopeptidase N. Water molecules in crystal structures shown as red spheres only. Backbone and hydrogen atoms are not shown for clarity.

approach, resulted in the relatively large average zinc CN values (up to 5.47, Table 2). This was the consequence of a rather stable bidentate coordination by E508, which is not present in the crystal structures of h.DPP III nor in the QMMM-optimized Michaelis complex structure.28 A stable tetrahedral coordination geometry of the zinc ion made up by H450, H455, E508, and a water molecule was observed during MD simulations utilizing a polarizable potential. However, during these simulations, the catalytically important hydrogen bond between E451 and the zinc-coordinated water molecule existed only sporadically, and also the substrate molecule moved slightly further from the metal ion. Finally, we developed the new hybrid 4-ligand-extended model parameters in a way that in the parametrization procedure the zinc ligands (H450, H455, E508, and water) and the amino acid residues from the second coordination shell (Y318, E451, E507, and E512) in h.DPP III have been taken into account. During MD simulations with these parameters, the relative orientations of the catalytically important amino acids were preserved, while exchange between 4- and 5-coordinated zinc ion occurred in accord with the proposed reaction mechanism for h.DPP III (Table 4). During MD simulations with these parameters, the relative

Overall, 4-ligand-extended model parameters were able to ensure proper (catalytically potent) active site conformation in all simulated systems.

■

CONCLUSION In this work we tested different strategies for modeling the zinc ion in h.DPP III in classical MD simulations: some developed to be used with nonpolarizable potential (nonbonded, hybrid bonded/nonbonded, and cationic dummy atom models) and the others suitable for the polarizable (Drude-2013) potential. Also, we developed new hybrid bonded/nonbonded zinc ion parameters. According to the calculations, the hybrid models that combine bonded and nonbonded approaches in a way that new bonds between the metal ion and protein residues are defined while leaving all other ligands unrestrained have shown the most promising results. On the other hand, the existing tetrahedron-shaped cationic dummy atom model has not proved suitable for simulating the h.DPP III active site geometry. The MD simulations utilizing hybrid bonded/nonbonded metal ion parameters of Yang et al.,19 and the ones we have developed specifically for the h.DPP III active site using their N


Article

Journal of Chemical Information and Modeling orientations of the catalytically important amino acids were preserved. For example, H568−substrate, Y318−substrate, and E508−Y318 hydrogen bonds, shown to be important for the enzyme activity, were present in the h.DPP III−substrate (Leuenkephalin or RRNA) complexes throughout MD simulations with the 4-ligand-extended model parameters. In addition, E508 was monodentately coordinated to the zinc ion, and the exchange between 4- and 5-coordinated zinc ion, mostly due to one or two water molecules coordinated to the metal ion during MD simulations, occurred in accord with the proposed reaction mechanism for h.DPP III (Table 4). We also noticed that it might be important to include a water molecule as the zinc ion ligand in the structure when simulating the enzyme− substrate complexes. Depending on the presence of a water molecule in the initial complexes, either nonproductive or productive substrate binding modes were sampled during the simulations of the h.DPP III−Arg-Arg-2-naphthylamide complex. The nonproductive binding mode was observed when the substrate (and occasionally E451), but not water, coordinated the zinc ion (Table 4) in the starting structure. On the other hand, the structural requirements around the zinc ion needed for the reaction to occur were fulfilled when simulations started from the complex with a water molecule coordinated to the zinc ion. During the simulations, water remained coordinated to the zinc ion, while the substrate and E451 never entered the zinc coordination shell; instead E451 was hydrogen bonded to the zinc-coordinated water molecule. Performances of our newly developed hybrid bonded/ nonbonded, 4-ligand-extended model parameters were tested in MD simulations of other systems as well: four other DPP III orthologues, thermolysin, neprilysin, and aminopeptidase N. The parameters enabled realistic modeling of the architecture around the catalytically important zinc ion in all cases (Tables 5−7). Furthermore, the MD simulations of the ligand-free DPP III orthologues revealed the interdomain movement in agreement with the available experimental results as well as with the results of the MD simulations obtained previously with the nonbonded, 12-6 LJ-type parameters for the zinc ion indicating rather weak influence of metal ion parameters on the long-range protein dynamic. In conclusion, the new 4-ligand-extended model parameters developed for human DPP III have enabled reliable modeling of the active site zinc coordination geometry in the other DPP III orthologues, as well as in distant hydrolases, and represent a value in molecular modeling of hydrolases. The credible active site architecture is indispensable for understand the substrate specificity of these enzymes and, so, important for development of novel ligands.

■

(Figures S4, S8, and S9); radius of gyration profiles obtained during the MD simulations utilizing nonbonded parameters (Figures S5 and S6); cartoon representation of closed h.DPP III structures (Figure S7); distances between the active site zinc ion and its ligands in structures of DPP III orthologues, aminopeptidase N, thermolysin, and neprilysin available in PDB (Table S1); RESP-calculated charge distribution, atom name, and atom type for atoms within selected amino acids for different hybrid bonded/nonbonded models (Table S2); selected average distances obtained during MD simulations utilizing different models to describe interactions between zinc ion and its surroundings (Table S3); hydrogen bond occupancies for different h.DPP III complexes (Tables S4 and S5); complete set of 4-ligand-extended model parameters (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +385 1 457 1251. Fax: +385 1 468 0245. ORCID

Antonija Tomić: 0000-0003-0109-3547 Author Contributions

A.T. and S.T. conceived and designed the study; A.T., G.H., D.A., and H.B. performed MD simulations and analyzed the data; M.R. and S.T. provided funding. All authors wrote the paper. Funding

Authors A.T., D.A., H.B., and S.T. received funding from Croatian Science Foundation under the project IP-2018-012936. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Authors A.T., D.A., H.B., and S.T. acknowledge Croatian Science Foundation (project number IP-2018-01-2936) for financial support and the Isabella cluster (http://www.srce. unizg.hr/en/isabella/) for computer time. We acknowledge Marko Tomin for his help in preparation of the initial structures of bacterial and mushroom DPPs III.

■

ABBREVIATIONS DPP III, dipeptidyle-peptidase III; h.DPP III, human dipeptidyle-peptidase III; cDPP III, closed human dipeptidyle-peptidase III; oDPP III, open human dipeptidylepeptidase III; y.DPP III, yeast dipeptidyle-peptidase III; At.DPP III, Armillaria tabescens dipeptidyle-peptidase III; Bt.DPP III, Bacteroides thetaiotaomicron dipeptidyle-peptidase III; Ca.DPP III, Caldithrix abyssi dipeptidyle-peptidase III; MD, molecular dynamic; QM, quantum mechanics; QMMM, quantum mechanics−molecular mechanics; FF, force field; RRNA, Arg-Arg-2-naphthylamide; L-en, Leu-enkephalin; LJ, Lennard-Jones; PDB, Protein Data Bank; MM/PBSA, Molecular Mechanics Poisson−Boltzmann Surface Area; CN, coordination number; MCPB, Metal Centre Parameter Builder

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.9b00235. Selected distances during MD simulations utilizing the polarizable Drude-2013 force field (Figure S1); selected distances along MD simulations utilizing 4-ligand model parameters (Figure S2); change in the number of water molecules during the MD simulations utilizing the 4ligand-extended model parameters (Figure S3); RMSD profiles of the zinc first and second coordination sphere monitored along the MD trajectory for human DPP III, DPP III ortholouges, and other zinc metallopeptidases O


Article


■

(23) Š poljarić, J.; Tomić, A.; Vukelić, B.; Salopek-Sondi, B.; Agić, D.; Tomić, S.; Abramić, M. Human Dipeptidyl Peptidase III: The Role of Asn406 in Ligand Binding and Hydrolysis. Croat. Chem. Acta 2011, 84 (2), 259−268. (24) Tomić, A.; González, M.; Tomić, S. The Large Scale Conformational Change of the Human DPP III-Substrate Prefers the “Closed” Form. J. Chem. Inf. Model. 2012, 52 (6), 1583−1594. (25) Chiba, T.; Li, Y.-H.; Yamane, T.; Ogikubo, O.; Fukuoka, M.; Arai, R.; Takahashi, S.; Ohtsuka, T.; Ohkubo, I.; Matsui, N. Inhibition of Recombinant Dipeptidyl Peptidase III by Synthetic Hemorphinlike Peptides. Peptides 2003, 24 (5), 773−778. (26) Barsun, M.; Jajcanin, N.; Vukelić, B.; Spoljarić, J.; Abramić, M. Human Dipeptidyl Peptidase III Acts as a Post-Proline-Cleaving Enzyme on Endomorphins. Biol. Chem. 2007, 388 (3), 343−348. (27) Liu, Y.; Kern, J. T.; Walker, J. R.; Johnson, J. A.; Schultz, P. G.; Luesch, H. A Genomic Screen for Activators of the Antioxidant Response Element. Proc. Natl. Acad. Sci. U. S. A. 2007, 104 (12), 5205−5210. (28) Tomić, A.; Kovačević, B.; Tomić, S. Concerted Nitrogen Inversion and Hydrogen Bonding to Glu451 Are Responsible for Protein-Controlled Suppression of the Reverse Reaction in Human DPP III. Phys. Chem. Chem. Phys. 2016, 18 (39), 27245−27256. (29) Dokmanić, I.; Sikić, M.; Tomić, S. Metals in Proteins: Correlation between the Metal-Ion Type, Coordination Number and the Amino-Acid Residues Involved in the Coordination. Acta Crystallogr. D. Biol. Crystallogr. 2008, 64 (3), 257−263. (30) Blumberger, J.; Lamoureux, G.; Klein, M. L. Peptide Hydrolysis in Thermolysin: Ab Initio QM/MM Investigation of the Glu143Assisted Water Addition Mechanism. J. Chem. Theory Comput. 2007, 3 (5), 1837−1850. (31) Zhang, C.; Liu, X.; Xue, Y. A General Acid−general Base Reaction Mechanism for Human Brain Aspartoacylase: A QM/MM Study. Comput. Theor. Chem. 2012, 980, 85−91. (32) Bertosa, B.; Kojić-Prodić, B.; Wade, R. C.; Tomić, S. Mechanism of Auxin Interaction with Auxin Binding Protein (ABP1): A Molecular Dynamics Simulation Study. Biophys. J. 2008, 94 (1), 27−37. (33) Hromić-Jahjefendić, A.; Jajčanin Jozić, N.; Kazazić, S.; Grabar Branilović, M.; Karačić, Z.; Schrittwieser, J. H.; Das, K. M. P.; Tomin, M.; Oberer, M.; Gruber, K.; Abramić, M.; Tomić, S. A Novel Porphyromonas Gingivalis Enzyme: An Atypical Dipeptidyl Peptidase III with an ARM Repeat Domain. PLoS One 2017, 12 (11), No. e0188915. (34) Sabljić, I.; Tomin, M.; Matovina, M.; Sučec, I.; Tomašić Paić, A.; Tomić, A.; Abramić, M.; Tomić, S. The First Dipeptidyl Peptidase III from a Thermophile: Structural Basis for Thermal Stability and Reduced Activity. PLoS One 2018, 13 (2), No. e0192488. (35) Tomin, M.; Tomić, S. Dynamic Properties of Dipeptidyl Peptidase III from Bacteroides Thetaiotaomicron and the Structural Basis for Its Substrate Specificity − a Computational Study. Mol. BioSyst. 2017, 13 (11), 2407−2417. (36) Tomin, M.; Tomić, S. Oxidase or Peptidase? A Computational Insight into a Putative Aflatoxin Oxidase from Armillariella Tabescens. Proteins: Struct., Funct., Bioinf. 2019, 87, 390. (37) Tomić, A.; Berynskyy, M.; Wade, R. C.; Tomić, S. Molecular Simulations Reveal That the Long Range Fluctuations of Human DPP III Change upon Ligand Binding. Mol. BioSyst. 2015, 11 (11), 3068− 3080. (38) Baral, P. K.; Jajcanin-Jozić, N.; Deller, S.; Macheroux, P.; Abramić, M.; Gruber, K. The First Structure of Dipeptidyl-Peptidase III Provides Insight into the Catalytic Mechanism and Mode of Substrate Binding. J. Biol. Chem. 2008, 283 (32), 22316−22324. (39) Bezerra, G. A.; Dobrovetsky, E.; Viertlmayr, R.; Dong, A.; Binter, A.; Abramic, M.; Macheroux, P.; Dhe-Paganon, S.; Gruber, K. Entropy-Driven Binding of Opioid Peptides Induces a Large Domain Motion in Human Dipeptidyl Peptidase III. Proc. Natl. Acad. Sci. U. S. A. 2012, 109 (17), 6525−6530.

REFERENCES

(1) Thomson; Gray. Bio-Inorganic Chemistry. Curr. Opin. Chem. Biol. 1998, 2 (2), 155−158. (2) Lee, Y.; Lim, C. Physical Basis of Structural and Catalytic ZnBinding Sites in Proteins. J. Mol. Biol. 2008, 379 (3), 545−553. (3) Li, P.; Merz, K. M. Metal Ion Modeling Using Classical Mechanics. Chem. Rev. 2017, 117 (3), 1564−1686. (4) Krężel, A.; Maret, W. The Biological Inorganic Chemistry of Zinc Ions. Arch. Biochem. Biophys. 2016, 611, 3−19. (5) Binding, Transport and Storage of Metal Ions in Biological Cells: Metallobiology; Maret, W., Wedd, A., Eds.; Royal Society of Chemistry: Cambridge, U.K., 2014. (6) Aromaticity and Metal Clusters; Chattaraj, P. K., Ed.; CRC Press: Boca Raton, FL, 2010. (7) Concepts and Methods in Modern Theoretical Chemistry, 1st ed.; Ghosh, S. K., Chattaraj, P. K., Eds.; CRC Press: Boca Raton, FL, 2013. (8) Li, P.; Merz, K. M. Taking into Account the Ion-Induced Dipole Interaction in the Nonbonded Model of Ions. J. Chem. Theory Comput. 2014, 10 (1), 289−297. (9) Li, P.; Song, L. F.; Merz, K. M. Parameterization of Highly Charged Metal Ions Using the 12−6-4 LJ-Type Nonbonded Model in Explicit Water. J. Phys. Chem. B 2015, 119 (3), 883−895. (10) Li, P.; Roberts, B. P.; Chakravorty, D. K.; Merz, K. M. Rational Design of Particle Mesh Ewald Compatible Lennard-Jones Parameters for + 2 Metal Cations in Explicit Solvent. J. Chem. Theory Comput. 2013, 9 (6), 2733−2748. (11) Li, P.; Song, L. F.; Merz, K. M. Systematic Parameterization of Monovalent Ions Employing the Nonbonded Model. J. Chem. Theory Comput. 2015, 11 (4), 1645−1657. (12) Yu, Z.; Li, P.; Merz, K. M. Extended Zinc AMBER Force Field (EZAFF). J. Chem. Theory Comput. 2018, 14 (1), 242−254. (13) Peters, M. B.; Yang, Y.; Wang, B.; Füsti-Molnár, L.; Weaver, M. N.; Merz, K. M. Structural Survey of Zinc-Containing Proteins and Development of the Zinc AMBER Force Field (ZAFF). J. Chem. Theory Comput. 2010, 6 (9), 2935−2947. (14) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97 (40), 10269−10280. (15) Li, J.; Zhu, T.; Cramer, C. J.; Truhlar, D. G. New Class IV Charge Model for Extracting Accurate Partial Charges from Wave Functions. J. Phys. Chem. A 1998, 102 (10), 1820−1831. (16) https://www.mayo.edu/research/labs/computer-aidedmolecular-design/projects/zinc-protein-simulations-using-cationicdummy-atom-cada-approach; https://www.mayo.edu/research/labs/ computer-aided-molecular-design/projects/zinc-protein-simulationsusing-cationic-dummy-atom-cada-approach. (17) Pang, Y. P.; Xu, K.; Yazal, J. E.; Prendergas, F. G. Successful Molecular Dynamics Simulation of the Zinc-Bound Farnesyltransferase Using the Cationic Dummy Atom Approach. Protein Sci. 2000, 9 (10), 1857−1865. (18) Pang, Y.-P. Novel Zinc Protein Molecular Dynamics Simulations: Steps Toward Antiangiogenesis for Cancer Treatment. J. Mol. Model. 1999, 5 (10), 196−202. (19) Yang, W.; Riley, B. T.; Lei, X.; Porebski, B. T.; Kass, I.; Buckle, A. M.; McGowan, S. Generation of AMBER Force Field Parameters for Zinc Centres of M1 and M17 Family Aminopeptidases. J. Biomol. Struct. Dyn. 2018, 36, 2595−2604. (20) Li, P.; Merz, K. M. MCPB.Py: A Python Based Metal Center Parameter Builder. J. Chem. Inf. Model. 2016, 56 (4), 599−604. (21) Patel, K.; Kumar, A.; Durani, S. Analysis of the Structural Consensus of the Zinc Coordination Centers of Metalloprotein Structures. Biochim. Biophys. Acta, Proteins Proteomics 2007, 1774 (10), 1247−1253. (22) Tomić, A.; Abramić, M.; Spoljarić, J.; Agić, D.; Smith, D. M.; Tomić, S. Human Dipeptidyl Peptidase III: Insights into Ligand Binding from a Combined Experimental and Computational Approach. J. Mol. Recognit. 2011, 24 (5), 804−814. P


Article

Journal of Chemical Information and Modeling (40) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926−935. (41) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; Li, X.; Caricato, M.; Marenich, A.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Ortiz, J. V.; Izmaylov, A. F.; Sonnenberg, J. L.; Williams-Young, D.; Ding, F.; Lipparini, F.; Egidi, F.; Goings, J.; Peng, B.; Petrone, A.; Henderson, T.; Ranasinghe, D.; Zakrzewski, V. G.; Gao, J.; Rega, N.; Zheng, G.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.; Fox, D. J. Gaussian 09, revision D.01; Gaussian, Inc: Wallingford, CT, 2009. (42) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. Ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from Ff99SB. J. Chem. Theory Comput. 2015, 11 (8), 3696−3713. (43) Hu, L.; Ryde, U. Comparison of Methods to Obtain ForceField Parameters for Metal Sites. J. Chem. Theory Comput. 2011, 7 (8), 2452−2463. (44) Sali, A.; Blundell, T. L. Comparative Protein Modelling by Satisfaction of Spatial Restraints. J. Mol. Biol. 1993, 234 (3), 779−815. (45) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R. C.; Walker, W.; Zhang, K. M.; Merz, B.; Roberts, S.; Hayik, A.; Roitberg, G. S.; Swails, J.; Götz, A. W., Kolossváry, I., Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M.-J.; Cui, D. R.; Roe, D. H.; Mathews, M. G.; Seetin, R.; Salomon-Ferrer, C.; Sagui, V.; Babin, T.; Luchko, S.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 12; University of California: San Francisco, CA, 2012. (46) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber Biomolecular Simulation Programs. J. Comput. Chem. 2005, 26 (16), 1668−1688. (47) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98 (12), 10089. (48) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103 (19), 8577−8593. (49) Lemkul, J. A.; Roux, B.; van der Spoel, D.; MacKerell, A. D. Implementation of Extended Lagrangian Dynamics in GROMACS for Polarizable Simulations Using the Classical Drude Oscillator Model. J. Comput. Chem. 2015, 36 (19), 1473−1479. (50) Yu, H.; Whitfield, T. W.; Harder, E.; Lamoureux, G.; Vorobyov, I.; Anisimov, V. M.; MacKerell, A. D.; Roux, B. Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field. J. Chem. Theory Comput. 2010, 6 (3), 774− 786. (51) Lamoureux, G.; Harder, E.; Vorobyov, I. V.; Roux, B.; MacKerell, A. D. A Polarizable Model of Water for Molecular Dynamics Simulations of Biomolecules. Chem. Phys. Lett. 2006, 418 (1−3), 245−249. (52) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. A Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules: Application to Small Water Clusters. J. Chem. Phys. 1982, 76 (1), 637−649. (53) Lamoureux, G.; Roux, B. Modeling Induced Polarization with Classical Drude Oscillators: Theory and Molecular Dynamics Simulation Algorithm. J. Chem. Phys. 2003, 119 (6), 3025−3039.

(54) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52 (2), 255−268. (55) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31 (3), 1695− 1697. (56) 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 (19), 5179−5197. (57) Vreven, T.; Morokuma, K.; Farkas, O.; Schlegel, H. B.; Frisch, M. J. Geometry Optimization with QM/MM, ONIOM, and Other Combined Methods. I. Microiterations and Constraints. J. Comput. Chem. 2003, 24 (6), 760−769. (58) Tan, C.; Yang, L.; Luo, R. How Well Does Poisson−Boltzmann Implicit Solvent Agree with Explicit Solvent? A Quantitative Analysis. J. Phys. Chem. B 2006, 110 (37), 18680−18687. (59) Connolly, M. L. Analytical Molecular Surface Calculation. J. Appl. Crystallogr. 1983, 16 (5), 548−558. (60) Maffucci, I.; Contini, A. Explicit Ligand Hydration Shells Improve the Correlation between MM-PB/GBSA Binding Energies and Experimental Activities. J. Chem. Theory Comput. 2013, 9 (6), 2706−2717. (61) Matovina, M.; Agić, D.; Abramić, M.; Matić, S.; Karačić, Z.; Tomić, S. New Findings about Human Dipeptidyl Peptidase III Based on Mutations Found in Cancer. RSC Adv. 2017, 7 (58), 36326− 36334. (62) Jajčanin-Jozić, N.; Tomić, S.; Abramić, M. Importance of the Three Basic Residues in the Vicinity of the Zinc-Binding Motifs for the Activity of the Yeast Dipeptidyl Peptidase III. J. Biochem. 2014, 155 (1), 43−50. (63) Kazazić, S.; Karačić, Z.; Sabljić, I.; Agić, D.; Tomin, M.; Abramić, M.; Dadlez, M.; Tomić, A.; Tomić, S. Conservation of the Conformational Dynamics and Ligand Binding within M49 Enzyme Family. RSC Adv. 2018, 8 (24), 13310−13322. (64) Tomić, A.; Tomić, S. Hunting the Human DPP III Active Conformation: Combined Thermodynamic and QM/MM Calculations. Dalt. Trans. 2014, 43 (41), 15503−15514. (65) Kumar, P.; Reithofer, V.; Reisinger, M.; Wallner, S.; PavkovKeller, T.; Macheroux, P.; Gruber, K. Substrate Complexes of Human Dipeptidyl Peptidase III Reveal the Mechanism of Enzyme Inhibition. Sci. Rep. 2016, 6, 23787.

Q


New Zinc Ion Parameters Suitable for Classical MD Simulations of

Recommend Documents