MM Linear-Interaction Energy Model for Substrate

Subscriber access provided by UNIV OF NEW ENGLAND

Article

Improved QM/MM Linear Interaction Energy Model for Substrate Recognition in Zinc-Containing Metalloenzymes Williams Ernesto Miranda, Van A Ngo, Pedro Alberto Valiente, and Sergei Yu. Noskov J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b05628 • Publication Date (Web): 22 Jul 2016 Downloaded from http://pubs.acs.org on July 23, 2016

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

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

Page 1 of 40

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

The Journal of Physical Chemistry

Improved QM/MM Linear Interaction Energy Model for Substrate Recognition in Zinc-Containing Metalloenzymes

Williams E. Miranda†‡, Van A. Ngo‡, Pedro A. Valiente†*, Sergei Yu. Noskov‡*

†

Computational Biology and Biomolecular Dynamics Laboratory, Center for Protein Studies, Faculty of Biology, University of Havana, Havana, Cuba. ‡ Centre for Molecular Simulations and Department of Biological Sciences, University of Calgary, Calgary, Alberta, Canada

* Corresponding Authors Sergei Noskov, CMS, Department of Biological Sciences, University of Calgary, 2500 University Drive, BI-449, Calgary, Canada, T2N 1N4 E-mail:[email protected] Phone: 402-210-7971 Pedro A. Valiente, CPPDL, Center for Protein Studies, Faculty of Biology, University of Havana, Havana, Cuba E-mail: [email protected]

ACS Paragon Plus Environment

1


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

Page 2 of 40

ABSTRACT: One of the essential challenges in the description of receptor-drug interactions in presence of various polyvalent cations (such as zinc, magnesium or iron) is the accurate assessment of the electronic effects due to cofactor binding. The effects can range from partial electronic polarization of the proximal atoms in a receptor and bound substrate to long-range effects related to partial charge transfer and electronic delocalization effects between the cofactor and the drug. Here, we examine the role of the explicit account for electronic effects for a panel of small-molecule inhibitors binding to the zinc-aminopeptidase PfA-M1, an essential target for anti-malarial drug development. Our study on PfA-M1:Inhibitor interactions at the QM level reveals that the partial charge and proton transfer due to bound zinc ion are important mechanisms in the inhibitors recognition and catalysis. The combination of classical MD simulations with a-posteriori QM/MM corrections with novel DFTB parameters for zinc cation and the Linear-Interaction Energy (LIE) approach offers by far the most accurate estimates for the PfA-M1:Inhibitor binding affinities, opening the door for future inhibitor design. ABBREVIATIONS LIE, linear interaction energy; FEP, free energy perturbation; TI, thermodynamic integration; MD, molecular dynamics; LR, linear response; QM/MM, quantum mechanics-molecular mechanics; Ki, inhibition constant; Ka, affinity constant; Kd, dissociation constant; ∆G, binding free energy difference; , mean absolute error; PfA-M1, alanyl-aminopeptidase 1 from Plasmodium falciparum; MEP, Molecular electrostatic potential; SCC, self-consistent charges, DFTB, Density Functional Tight Binding.


2

Page 3 of 40

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


1. Introduction The calculation of absolute binding affinities for protein-ligand complexes remains as one of the major challenges in computational Structure-Based Drug Design (SBDD) approaches, thus limiting its application in virtual screening or pre-clinical drug discovery. An impressive progress in the development and application of modeling techniques to evaluate binding affinities has been achieved for the last two decades.1-3 Many of the traditional challenges related to conformational-space sampling for receptor-ligand complexes appear to be solvable with the emergence of powerful computers and molecular dynamics (MD) simulations that reach physiologically-relevant time-scales.4 Although binding simulations are coming to an age, the factors limiting further improvements in the computed binding affinities relate to the key assumption hard-wired in all of the classical force-fields. The local polarization is implicitly accounted via static partial charges, and hence the receptor site is unable to account for dynamical changes in the proximal and distal electronic densities upon ligand and/or ion binding.5 The problem is aggravated by the fact that many enzymes contain cofactors such as zinc, magnesium or iron ions as part of the binding pockets that are targeted by substrates and small molecule-inhibitors. The charge transfer and local polarization effects are expected to play an important role in the molecular recognition, catalysis and inhibition phenomena in these systems.6-7 Therefore, accurate modeling of cofactor dependent interactions warrants the use of Quantum-Mechanical (QM) calculations.8 The accuracy comes at a cost because sampling of large QM-treated systems is often a challenging problem.9-10 The two main strategies adopted to address this issue rely on the use of i) QM for active site models as reduced representations of the ion binding site (also known as the cluster approach)5,11-14 and ii) QM/MM for the whole system, where the ion binding site is represented at QM level and the rest of the protein at MM


3


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

Page 4 of 40

level.9-11,14-16 The cluster approach has been the most used to study electronic effects in metalloenzymes, especially in the first stages of mechanistic studies, because of its lower cost compared to the QM/MM scheme.11-14,17 However, the latter allows taking into account the protein environment in a more accurate way.9-11,14-16 Furthermore, treating part of the system classically (MM) allows an enhancement of the conformational-space sampling with realistic treatment of the environmental effects.9-11,14-16 An important gain in QM calculations’ performance can be achieved by using the SelfConsistent-Charge

Density-Functional

Tight-Binding

(SCC-DFTB)

approximation.18-20

Recently, SCC-DFTB scheme has been extended to include third-order contributions, known now as DFTB3 and has been thoroughly benchmarked for cation binding to various metalloproteins.11,21 Compared to the second-order formulation, DFTB3 considers the local electrical field dependence of the atomic Hubbard parameters. Therefore, DFTB3 provides an improved description of chemical properties for biomolecular systems.11,20-25 The present study strives to assess the role of the quantum effects in the binding of seven inhibitors to the zinc-modulated PfA-M1 malarial enzyme by using the recently developed DFTB3 parameters for Zn2+.11 In the first stage, we used the cluster approach to study the charge transfer and polarization processes between the zinc cofactor and its ligating atoms from the protein and the inhibitors. In the second stage, we assessed the contribution of the electronic effects to the binding free energy of the seven PfA-M1:Inhibitor complexes by using the QM/MM-LR approach.9-10 Our findings show that the explicit account for quantum mechanics effects is useful to understand the cofactor-mediated binding of small molecules to the enzyme, the involvement of catalytic residues and also to improve the predictive power of the LIE method for binding free energy calculations. Furthermore, free energy calculations become more


4

Page 5 of 40

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


accurate when classical MD simulations are used to generate ensembles of structures for proteinligand complexes and then quantum corrections are added a-posteriori.20,26-27 Therefore, the information obtained from QM approaches in this study could be useful to guide the discovery and design of highly potent and specific inhibitor of PfA-M1 as novel anti-malarial therapeutic candidates and may also be adopted as a general strategy for the design of inhibitors targeting zinc-bound enzymes.


5


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

Page 6 of 40

2. Methods 2.1 Selection of PfA-M1:Inhibitor complexes In this work we selected seven PfA-M1:Inhibitor complexes (Tables 1 and 2 and Fig. 1). They were chosen based on the following criteria: i) availability of high-resolution 3D structures;28 ii) protein-ligand complexes formed by a non-covalent interactions and iii) availability of experimental inhibition (Ki) constants. 2.2 Building of active site models for the calculation of electronic effects. To understand the impact of the cofactor on the distribution of partial charges in the coordinating atoms (from the protein and the ligand), we computed Mulliken charges using DFTB3/3OB11 for models of the active site built from each complex (Table 1 and Fig. 1). Each model contains the zinc ion, the inhibitor and the residues E463, H496, E497, H500, and E519 (Fig. 1). The histidines and glutamates were modeled as methyl-imidazol and acetate groups, respectively.11,13-14 Mulliken charges from DFTB3/3OB were compared to the corresponding values from full DFT computations reported previously for other zinc-metalloenzymes.9-10 All DFTB3/3OB calculations were performed with DFTB+ v1.2.2,29 using the Slater-Koster parameter files recently developed by Elstner and colleagues.21-23 Electronic smearing was performed to obtain the convergence of the self-consistent charges (SCC) inside DFTB3/3OB. The Fermi-filling maximum temperature factor for producing the smearing was 300K. Jmol software (http://www.jmol.org/) was used to analyze the molecular electrostatic potential (MEP) maps. 2.3. MD simulations All systems were prepared with CHARMM-GUI web server.30 Protonation states to residues were assigned based on pKa computations, (pKa’s were calculated using PROPKA31 web server


6

Page 7 of 40

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


at http://propka.ki.ku.dk/). The protein-ligand complexes were centered in a cubic simulation box (15 Å distance between the solute and the edges) and solvated in 150 mM of neutralizing NaCl solution. The TIP3P water model was used32 while protein and inhibitor molecules were modeled with the CHARMM3633 and the /General Force Field (CGenFF),34 respectively. Production MD simulations were carried out for 50 ns using the GROMACS (version 4.5.6)35 software package. The LINCS algorithm was used to constrain all the covalent bonds in non-water molecules,36 while the SETTLE algorithm was used to constrain bond lengths and angles in water.37 Before the production phase, all systems were minimized by running 1000 steps of Steepest Descent, followed by 1000 steps by Conjugate Gradient minimization. The heavy atoms in protein were constrained with 10 kcal/(mol Å2) harmonic constraint. Initial velocities were randomly generated for Maxwell distribution at T= 1 oK. Next, the simulation systems were gradually heated for 300 ps to reach the experimentally reported assay temperatures. The following protocol was used for system’s equilibration and thermalization. The Leapfrog scheme,38 with an integration time step of 2 fs was used to integrate the equations of motion. Temperature was controlled using a weak coupling to a bath with a time constant of 0.1 ps.39 For pressure control the Berendsen coupling algorithm with a time constant of 1.0 ps was used.39 Next, Langevin dynamics40 with an integration time step of 2 fs scheme was used for all production runs. Temperature was controlled using a weak-coupling to a bath with a time constant of 2.0 ps. For pressure control a Parrinello-Rahman41 coupling algorithm with a time constant of 5.0 ps was used. Long range electrostatic interactions were handled by Particle Mesh Ewald (PME) summation.42-43 The van der Waals interactions were modeled by Lennard-Jones (6-12) potentials.35 To avoid Zn2+ ion unbinding from the active site and/or incorrect pseudo-valence bond formation (due to the use of a classical force field), we applied distance restraints to


7


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

Page 8 of 40

maintain the Zn2+ ion in its correct ligation state.9-10,44 The distances between the Zn2+ ion and the atoms Nε (H496 and H500) and Oδ (E519) were restrained with a harmonic restraint potential using a force constant equal to 400 kcal/(mol*Å2). The (respective) crystal structure was used as a starting point for each complex. The tool g_rms was used to calculate the Root-Mean Square Deviation (RMSD) for the protein backbone, the inhibitors and the zinc ion. The tool g_hbond was used to calculate the number of hydrogen bonds between the ligands and the protein. 2.4 QM/MM calculations The DFTB3/MM calculations were performed using CHARMM version c38b1.45 To develop an appropriate QM/MM-LR model for PfA-M1:Inhibitors complexes, we assessed three different approaches by computing the single-point QM/MM energies and SASA terms in eqs. 1 and 2 (see next section) for each complex from: i) The X-Rays structure (model A) ii) The time-averaged structure from MD simulations (model B) and iii) Ensemble averaged energies, using 2000-5000 frames (model C). In all three cases, the QM region includes the side chains of residues E463, H496, E497, H500, and E519, the inhibitor, the zinc ion and those water molecules that remained at the binding site during the simulation. The rest of the protein (MM part) was treated with a classical force field.33 The QM truncated side chains from the residues mentioned above were capped with hydrogen atoms (link atom method), placed between Cα-Cβ atoms. The QM and MM regions interact via electrostatic interactions between the MM point charges and the QM wave function, and van der Waals interactions between QM and MM atoms.46 For the MM and QM parts of the QM/MM calculations, CHARMM36 force field and DFTB3/3OB method11,21-23,47 were used, respectively. The convergence criterion for the self-consistent field (SCF) cycle was set to 10-7 kcal/mol.


8

Page 9 of 40

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


2.5

Evaluation of binding free energies

Several approaches to calculate the binding affinities of protein-ligand complexes from molecular modeling have been developed, ranging from empirical and ‘‘knowledge-based’’ scoring functions to those based on free energies calculations, such as the rigorous free energy perturbation (FEP) and thermodynamic integration (TI) methods.48,49 Although FEP and TI approaches have shown to be accurate for binding affinity calculations, they are exceedingly demanding in computational power if being applied to large molecular datasets.50-51 Thus, further development of fast and accurate methods for structure-based drug design is critical. On the other hand, methods based on the Linear Response (LR) approximation rely on MD conformational sampling for receptor-bound and unbound ligand states (known as “end-point” methods). Nowadays, LR-based methods offer an attractive computational cost/accuracy ratio for drug design purposes, as they are faster than FEP or TI simulations and more accurate than scoring function-based methods. In LR-based methods, binding free energy is calculated as the 〉, van der Waals linear combination of the differences ∆ in the electrostatic energies 〈 52 〉, and the solvent-accessible surface areas (SASA)53-54 obtained from MD or energies 〈

MC simulations (eq. 1):

〉 + ∆〈 〉 + 〈〉 + ∆ = ∆〈

(1)

The subscripts l and s denotes the ligand and its surrounding environment, respectively. The coefficients α and β are scaling factors for these energy terms, while γ is an empirical constant whose value has been related to the binding site hydrophobicity and more recently to the balance between electrostatic and van der Waals energies in protein-ligand complexes.55-57 The LR approximation theory provides a physical basis for treating the electrostatics contribution to the binding free energy, using a value of β = 0.5. The van der Waals interactions between a ligand


9


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

Page 10 of 40

and its environment, represented by Lennard-Jones potentials, are commonly used to calculate the non-polar contribution to binding free energy.52 Another LR-based approach replaces the van der Waals and electrostatic energies in the LIE by a single QM/MM interaction energy value calculated for protein-ligand complexes scaled by (eq. 2).10 ∆ = ∆/ + ∆ + (2) This approach, coined as QM/MM-LR, also considers a SASA term scaled by γ coefficient and a κ constant term (eq. 2).10 The replacement of the MM-based van der Waals and electrostatic terms in eq. 1 by the QM/MM term10 allows simple extension of the LR approach to zinc-mediated inhibitor binding. 2.6 Development of QM/MM-LR models for PfA-M1:Inhibitor complexes We used eq. 2 to assess three different QM/MM-LR models (models A-C, see section 2.4) with ∆/ and terms calculated as: %&

! #$ ∆/ = / − / − /

∆ =

!

−

#$

− %&

(3) (4)

To assess the influence of the electronic effects in the performance of our developed QM/MM-LR models, we also used classical methods based on docking algorithm (Autodock4Zn software)58 and eq. 1 (coined as MM-LR model in this work). For MM-LR model, we also performed MD simulations for the ligands in water as outlined in section 2.3 for the complexes, 〉 in eq. 1 as but for 5 ns. Furthermore, we considered intra-ligand electrostatic interactions 〈

suggested in a previous work for highly charged ligands.59 2.7 Statistical Analysis


10

Page 11 of 40

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


To measure the parameterization models performance relative to experimental binding free energies (∆Gcalc vs. ∆Gexp), we calculated the correlation and multiple determination coefficients (r and R2 respectively). The former is the squared Pearson product moment (r) to measure linear dependence between two data sets as: '=

∑(∆ ! − 〈∆ ! 〉)(∆ & − 〈∆ & 〉) +∑(∆ ! − 〈∆ ! 〉), ∑(∆ & − 〈∆ & 〉),

(5)

R2 coefficient reports variance proportion from the experimental binding data that is ‘explained’ as model response to the simulation input data (eq. 6). ., = 1 −

(6) 0

where SSE = Σi(∆Gicalc-∆Giexp)2 and SST = Σi(∆Giexp-‹∆Giexp›)2 ; ∆Gicalc and ∆Giexp are the calculated and experimental binding free energies for complex i, respectively. The average absolute error (‹|error|›), the standard deviation of the error (SDerror), and the maximum error |Emax| were calculated as measures of deviations, in kcal/mol, between the optimized models and the experimental values.


11


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

Page 12 of 40

3. Results and Discussion 3.1 Inhibitors of PfA-M1 are involved in partial charge and proton transfer to zinc-cofactor and the residues E463 and E497, respectively. The routine use of quantum mechanics (QM) in all phases of in silico drug design (DD) to account for electronic effects is the logical next step in the evolution of this field.8 For instance, in a study that involved 49 metalloenzyme complexes, the authors found an average of 0.6 electron units of charge (e.u.) transferred between the protein and the ligand.60 In our work, we found that the effective charge on zinc was ≈ 0.67 e.u. for the seven active site models (Fig. 2). This means ≈ 1.33 e.u (2.00 - 0.67) transferred to the Zn2+ ion from its ligating atoms (from the inhibitor and also the protein residues H496, H500 and E519). In fact, the calculated molecular electrostatic potentials (MEPs) for the seven active site models show significant charge transfer and polarization between the bound cofactor (zinc) and the atoms directly ligating it (Fig. 3). Interestingly, the calculated effective charge on zinc seems to depend on the level of theory used for its calculation. For example, studies on other metalloenzymes with different QM-based methods display a range of effective charges for zinc cofactor from 0.55 to 1.00 e.u. for HartreeFock (HF) and DFT, respectively (Table 3).9-10,61 The charge derived for zinc using DFTB3/3OB is different from the calculations using DFT (Table 3). This result could be influenced by the fact that in DFTB3 the electron-electron interaction (Hartree and exchange-correlation term) is approximated by a simple Coulomb term that describes the interaction of atomic partial charges.21 Our study on the active site models at QM level can also be a step forward to understand the catalytic mechanism of the enzyme. In inhibitor-bound crystal structures, E497 has been reported to be within hydrogen bond distance of the ligand’s oxygen (expected to mimic that of


12

Page 13 of 40

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


nucleophilic water), and is postulated to act as the catalytic base of the hydrolysis reaction in other enzymes of the same aminopeptidase superfamily.62-66 Another residue, E463, has been found to interact with the ligand N-terminal nitrogen in crystal structures, although its implication in the catalytic mechanism of PfA-M1 remains unclear.62,66-67 Recently, classical MD simulations of PfA-M1 in complex with a substrate mimetic shed light on the hypothetical catalytic mechanism and the roles of active-site residues.68 The residues E463 and E497 displayed hydrogen-bond distances (≈ 2.5 Å) from the nucleophilic water during the simulations, although E463 also shown to interact with the N-terminal nitrogen of the ligand. Then, the authors proposed a catalytic mechanism for PfA-M1 where E497 and E463 interacted with the hydroxyl and N-terminal hydrogens respectively during the progression of the catalytic cycle.68 Nevertheless, as the proposed mechanism was based on classical MD simulations, relevant information on quantum effects such as proton transfer between the ligand and the residues E497 and E463 couldn’t be taken into account. As the inhibitors used in our work are transition-state analogs, the QM simulations could deepen into the role of E497 and E463 residues in the catalytic mechanism. We compared the distances between the ligands’ hydroxyl hydrogen, and the nearest carboxyl oxygen atom from the residues E497 and E463 before (crystal structure) and after QM minimization, to identify proton transfer in our systems. For the bestatin-based inhibitors, both E497 and E463 are within hydrogen bond distance (< 3 Å) with respect to the ligands’ hydroxyl hydrogen (Table 4 and Fig. 4A). As E497 is closer than E463 to the hydroxyl group (1.61 and 2.77 Å on average, respectively), the hydroxyl hydrogen was transferred to the former residue after the QM minimization (Table 4 and Fig 4B). However, in the case of PfA-M1:Bestatin system (3EBH) the distance is within the hydrogen bond range for E497, although 0.7 Å closer to the ligand’s


13


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

Page 14 of 40

hydroxyl hydrogen compared to the crystal structure (Table 4). On the other hand, E463 is closer on average to the N-terminal hydrogen than to the hydroxyl hydrogen (1.95 and 2.77Å, respectively) in the four bestatin-based inhibitors’ crystals (Tables 4 and 5 and Fig. 4A). This could explain why the N-terminal hydrogen is transferred to E463 instead of E497 after the QM minimization (Table 5 and Fig. 4B). For the hydroxamate-based ligands, the proton transfer from the hydroxyl-hydrogen to E497 or E463 residues will also depend on the distances in the crystal structures. In 4R5T and 4R5V systems, the proton was transferred to E463 after the QM minimization. This residue is closer to the hydroxyl-hydrogen in the crystal structures than residue E497 (1.71 and 2.71 Å on average, respectively, see Table 4 and Fig. 5A and B). On the contrary, E497 is closer to the hydroxyl hydrogen than E463 for 4R5X system (1.97 and 2.29 Å, respectively, see Table 5). Then, the proton is transferred to the former residue. To sum up, for bestatin analogs the QM minimization predicts two proton transfer processes: i) from the inhibitor’s hydroxyl hydrogen to E497 and ii) from the N-terminal hydrogen to E463. However, in the case of hydroxamate analogs, the hydroxyl hydrogen can be transferred to either E497 or E463 depending on the distance. Our QM studies reveal that the residues E497 and E463 could indeed play a very important role during the catalysis. Depending on the distance, both of them are able to activate the nucleophilic water molecule (represented here by the inhibitors’ hydroxyl group) by accepting a proton. It is known that the interaction of hydroxyl groups with a zinc cofactor lowers their pKa, which promotes proton transfer reactions.11,13,69 Furthermore, E463 is also able to accept a proton from the ligand N-terminal. Although in other works the residues E497 and E463 have been postulated to act as a catalytic base through experimental and theoretical approaches,62-66,68 this is the first time a QM study is done for Pf-


14

Page 15 of 40

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


AM1 enzyme to model the role of both residues as catalytic bases. The observed partial charge and proton transfer in PfA-M1:Inhibitor complexes indicates that such quantum effects play a very important role in ligand binding. Therefore, these effects must be incorporated into binding free energy calculation methodologies (i.e. the LIE method)52 to be used for ligand optimization in drug design schemes. 3.2 Including quantum effects and conformatinal sampling improves the accuracy of binding free energy calculations for PfA-M1:Inhibitor complexes. The zinc-mediated binding of inhibitors to the metalloenzymes is challenging to describe with the standard arsenal of modern biomolecular force fields. The apparent need to introduce a pseudo-bonding, adjust partial charge or connectivity at the zinc-binding site are a few of many strategies developed in the past.7-8,70-71 Our QM studies on the active site models suggest the occurrence of partial charge and proton transfer in the active site of PfA-M1:Inhibitor complexes (Figs. 2-5). However, to assess the influence of the electronic effects on the binding strength of the seven ligands, we need to go beyond the simplistic active site models and consider full protein-ligand interactions. So far, the most viable strategy for binding free energy calculation of ligands to metalloenzymes combines the QM/MM multiscale hybrid scheme with MD simulations and the linear response (LR) assumption (aka QM/MM-LR method).9-10 The inclusion of the electronic effects in the study of the zinc-dependent matrix metalloproteases MMP-9 and MMP-3, significantly improved the prediction ability of the LR approach.9-10 Even though both enzymes have quite similar binding sites,72-73 the QM/MM-LR approach was able to capture subtle differences (about one order of magnitude) in the binding affinity for the most potent inhibitors.9 Then, to illustrate the extent both the electronic and conformational sampling effects can influence binding free energy calculations in our PfA-M1:Inhibitor complexes, we


15


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

Page 16 of 40

computed binding affinities using the hybrid method QM/MM-LR (eq. 2) and also classical methods like molecular docking and the MM-LR scheme (eq. 1). Our results show that from among the QM/MM-LR models, the one that uses the ensemble average of the energy terms in eq. 2 (model C) provides the best correlation to experimentally measured binding free energies for PfA-M1:Inhibitor complexes (‹|error|› = 0.30 kcal/mol, SDerror= 0.22 kcal/mol, |Emax| = 0.58 kcal/mol, r2= 0.91, and R2= 0.82, see Fig. 6 and Table 6). On the other hand, a lower (but still good) performance is shown by those QM/MM-LR models based on a single-structure: Model A (X-rays structure): ‹|error|› = 0.44 kcal/mol, SDerror= 0.32 kcal/mol, |Emax| = 0.90 kcal/mol, r2= 0.78, R2= 0.61; and Model B (average structure): ‹|error|› = 0.44 kcal/mol, SDerror= 0.32 kcal/mol, |Emax| = 0.99 kcal/mol, r2= 0.79, R2= 0.62. The use of a single average structure from MD-generated sampling was proposed to provide similar results to ensemble averaging in a number of studies aiming at the binding affinity correlations74 and in protein pKa calculations.75 For example, the QM/MM-LR parameterization models developed for MMP-9 and MMP-3 enzymes based on MD timeaveraged structures showed high accuracy for the binding free energy calculations of metalloenzyme:inhibitor complexes.9-10 Our results, however, show that the ensemble averaging (model C) led to a better accuracy compared to single structure-based QM/MM-LR models A and B (Table 6 and Fig. 6). The role of the local dynamics in the assessment of the binding strength can be illustrated by comparing high- (D66) and low-(R5X) affinity inhibitors and the corresponding RMSD of the binding pockets (Fig. 7 and Table S1). The position restraints applied to the zinc ion stabilized the side chains of residues H496, H500 and E519 in their starting positions (RMSD < 1Å, see Table S1). The unrestrained backbone atoms fluctuate at around 1.5 Å. The high-affinity


16

Page 17 of 40

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


inhibitor D66 shows an average fluctuation of 0.65 Å, while the low-affinity inhibitor R5X displayed RMSD of 1.19 Å (Table S1). The inhibitors D66 and R5X establish 7-8 and 3-4 Hbonds, respectively, during the MD simulations (Table S2). The sampling of the conformational space visited by both bound ligands leads to a considerable improvement in the predicted binding affinity compared to QM/MM single-structure approaches: Both inhibitors show 2∆∆ ! 2 = 2.99 kcal/mol (Table 1) and the closest |∆∆ & | value was obtained by using QM/MM-LR model C (2.07 kcal/mol, see Table S3), while models A and B were much less accurate (|∆∆ & | = 1.35 and 1.32 kcal/mol respectively, see Table S3).The bifurcation in hydrogenbonding is an essential factor in accuracy of the computed binding affinities. On the other hand, the docking results from AutoDock4ZN58 (for each ligand’s best pose) “as is” positively correlate with available experimental data (‹|error|› = 1.22 kcal/mol, SDerror= 0.61 kcal/mol, |Emax| = 1.95 kcal/mol, r2= 0.64, R2= 0.41), but remains clearly inferior when compared to the also single-structure based QM/MM-LR models A and B (Fig. 6, Tables 6 and S4). The scoring function includes a specialized potential for describing the interactions of zinccoordinating ligands. This potential considers both energetic and geometric components and was calibrated on a data set of 292 crystal complexes containing zinc. Re-docking experiments in the original work showed a significant improvement for ligand pose prediction and binding free energy calculations.58 Actually, the re-docking experiments performed on PfA-M1:Inhibitor complexes in our work shows a low average RMSD of 1.13 Å relative to the pose observed in the crystal structure (Table S4). These results support the use of AutoDock4ZN for ligand database virtual screening and the use of QM/MM-LR for ligand optimization in drug design schemes aimed to design inhibitors for PfA-M1 and other zinc-metalloenzymes.


17


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

Page 18 of 40

The MM-LR model is intermediate in speed and accuracy between docking and QM/MM-LR methods. This model considers structures generated by MD simulations, leading to ~0.64 kcal/mol improvement with respect to AutoDock4Zn results (‹|error|› = 0.58 kcal/mol, SDerror= 0.61 kcal/mol, |Emax| = 1.57 kcal/mol see Fig. 6, Tables 6 and S5). The improvement again highlights the advantage of considering binding pocket flexibility. However, although the MMLR model and the QM/MM-LR model C are both based on the same MD generated ensemble for each complex, considering the QM effects in the latter improves the ‹|error|› and |Emax| values in 0.28 and 0.99 kcal/mol respectively (Table 6). 3.3 The contributions to the binding free energy in QM/MM-LR models. The energy term ∆/ that describes protein-ligand interactions (including the electronic effects) contributes favorably to the binding of the inhibitors to the enzyme (Tables 6, S3 and Fig 8). Furthermore, the ∆/ -term is also a reasonably accurate predictor of the binding affinities in the studied panel of inhibitors (R2 = 0.50 for ∆EQM/MM vs. ∆Gexp in QM/MM-LR model C, see Fig. 9 A). These results reflects that the interactions established between the inhibitors and the binding site of the enzyme, along with the electronic effects that take place by the presence of the zinc cofactor, favor the intermolecular association. On the other hand, the ∆-term (that represents the desolvation energy) opposes to the formation of the complex (Tables 6, S3 and Fig. 8). In fact, the correlation between the ∆term vs. ∆Gexp in QM/MM-LR model C is almost negligible (R2 = 0.16*10-3, see Fig. 9 B). It is known that the desolvation of polar groups is energetically more expensive than for non-polar ones.76-78 Then, the unfavorable contribution of the ∆-term in the QM/MM-LR models may be reflecting the desolvation of polar groups from both the ligands and the binding site of


18

Page 19 of 40

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


the protein, suggesting that the lost interactions between those polar groups and the solvent are not properly compensated when the complex is formed. The κ-term in our QM/MM-LR models makes the largest favorable contribution to the formation of the complexes by compensation of the unfavorable desolvation energy (Table 6 and Fig. 8). As the κ-term in LIE-based models is associated with the hydrophobicity of the binding site in the protein,55-57 our results suggest that the binding of the inhibitors to PfA-M1 is driven by the hydrophobic effect. However, the binding specificity towards PfA-M1 appears to depend on the interactions captured by the QM/MM energy term. For example, bestatin-based inhibitors show QM/MM energy values between -207.21 and -279.75 kcal/mol, while hydroxamate-based inhibitors display values ranging from -103.40 to -111.96 kcal/mol (Table S3, model C). Scaling the QM/MM energy term averages separately for bestatin-based and hydroxamate-based inhibitors (-251.22 and -113.24 kcal/mol respectively) by α coefficient (Table 6, model C) results in ∆/ term values equal to -3.25 and -1.54 kcal/mol respectively (1.71 kcal/mol difference). On the other hand, scaling SASA term averages for bestatin-based and hydroxamatebased inhibitors (-891.52 and -814.17 Å2 respectively, see Table S3 model C) by γ coefficient (Table 6, model C) results in ∆ term values equal to +5.86 and +5.35 kcal/mol (only 0.51 kcal/mol difference). So, the inhibitors (relative) binding trend is reflected by the QM/MM energy term (R2 = 0.49 for ∆∆EQM/MM vs. ∆∆Gexp in QM/MM-LR model C, see Fig. 9 C), that may be influenced by the different hydrogen bond forming capabilities of bestatin and hydroxamate-based inhibitors during the MD simulations (Table S3). This is in agreement with experimental data that shows bestatin-based inhibitors as stronger PfA-M1 binders with respect to hydroxamate-based inhibitors (Table 1). It is also worth mentioning that there is a weak correlation between QM/MM energy and SASA terms (R2 = 0.36 for model C, see Fig. S1). This


19


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

Page 20 of 40

fact also supports that both terms are reflecting different aspects of the binding process. In the case of the MM-LR model, the contributions from hydrogen bonds to the binding affinity are 52 〉. However, the prediction of the binding trend of the inhibitors considered in the term ∆〈 〉 vs. ∆∆Gexp, see Fig. S2), by this term is considerably less accurate (R2 = 0.30 for ∆∆〈

compared to ∆∆EQM/MM (Fig. 9 C). Likely, the main reason for this result is because the MMLR model doesn’t account for the electronic effects involved in the ligand binding. Although the training set used is small, the QM/MM-LR models (especially model C) show a strong sensitivity to the chemical features of the seven ligands by predicting their affinities accordingly. This is important considering that there are ligands with very similar chemical structure (Table S1). This fact suggests that our QM/MM-LR models (especially model C) could be used as tools for the optimization phase of bestatin and hydroxamate-like inhibitors of PfA-M1. Such specificity is expected to be useful for drug design schemes that target the inhibition of this particular enzyme. Previous reports on classical-LR models show that the training sets must include a diverse set of protein-ligand complexes in order to be transferable to other systems.5657,79-80

Therefore, our parameters should not be readily transferable, but the methodology should

be very useful for other zinc-metalloenzymes. 4. CONCLUDING REMARKS We have assessed the QM effects for seven PfA-M1:Inhibitor complexes with the novel zinc parameters for DFTB3 method. The cluster approach provides evidence of partial charge and proton transfer in zinc-mediated inhibitor binding. Furthermore, the calculations support the reports on the involvement of the residues E463 and E497 in the catalytic mechanism of the enzyme. On the other hand, considering the electronic effects together with conformational sampling in the LR method (QM/MM-LR) leads to an excellent performance in predicting the


20

Page 21 of 40

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


binding affinities for the seven complexes. Compared to the classical approaches that neglect electronic effects, the QM/MM-LR models are more accurate. Finally, we show conclusively that the novel DFTB3/3OB parameter set provides a reasonable level of accuracy at modest computational cost. We believe that our work is a step toward future research that aims to rationally design small molecule inhibitors for malarial zinc-containing enzymes. ACKNOWLEDGMENTS This work was supported by grants from the Natural Sciences and Engineering Research Council (Canada) to SYN (RGPIN-315019) and the Alberta Innovates Technical Futures Strategic Chair in (Bio)Molecular Simulations. WEM was also supported by the Emerging Leaders in the Americas Program (ELAP). All of the computations were performed on the WestGrid/ Compute Canada facilities, and the University of Calgary TNK cluster acquired with direct support by the Canada Foundation for Innovation. VN is supported by Alberta Innovates Health Solutions and Canadian Institute for Health Research postdoctoral fellowships. We would like to thank Drs. Qiang Cui and Xiya Lu for sharing latest parameters for Zn2+ in biological systems and Dr. Mauricio Da Silva for his valuable advice on using DFTB3 method. Supporting Information. The tables for QM/MM and SASA terms, RMSD values and H-Bonding from MD simulations for the data set complexes. Sample input files for QM/MM, QM and MD simulations. This material is available free of charge via the Internet at http://pubs.acs.org


21


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

Page 22 of 40

REFERENCES 1. Gilson, M. K.; Zhou, H. X. Calculation of Protein-Ligand Binding Affinities. Annu. Rev. Biophys. Biomol. Struct. 2007, 36, 21-42. 2. Jiao, D.; Golubkov, P. A.; Darden, T. A.; Ren, P. Calculation of Protein-Ligand Binding Free Energy by Using a Polarizable Potential. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 62906295. 3. Singh, N.; Warshel, A. Absolute Binding Free Energy Calculations: On the Accuracy of Computational Scoring of Protein-Ligand Interactions. Proteins: Struct., Funct., Bioinf. 2010, 78, 1705-1723. 4. Dror, R. O.; Dirks, R. M.; Grossman, J. P.; Xu, H. F.; Shaw, D. E. Biomolecular Simulation: A Computational Microscope for Molecular Biology. Annu. Rev. Biophys. Biomol. Struct. 2012, 41, 429-452. 5. Ngo, V. A.; Da Silva, M. C.; Kubillus, M.; Li, H.; Roux, B.; Elstner, M.; Cui, Q.; Salahub, D. R.; Noskov, S. Y. Quantum Effects in Cation Interactions with First and Second Coordination Shell Ligands in Metalloproteins. J. Chem. Theory Comput. 2015, 11, 4992−5001. 6. Grater, F.; Schwarzl, S. M.; Dejaegere, A.; Fischer, S.; Smith, J. C. Protein/Ligand Binding Free Energies Calculated with Quantum Mechanics/Molecular Mechanics. J. Phys. Chem. B 2005, 109, 10474-10483. 7. Raha, K.; Merz, K. M. A Quantum Mechanics-Based Scoring Function: Study of Zinc Ion-Mediated Ligand Binding. J. Am. Chem. Soc. 2004, 126, 1020-1021. 8. Raha, K.; Peters, M. B.; Wang, B.; Yu, N.; Wollacott, A. M.; Westerhoff, L. M.; Merz, K. M. J. The Role of Quantum Mechanics in Structure-Based Drug Design. Drug Discovery Today 2007, 12, 725-731. 9. Khandelwal, A.; Balaz, S. QM/MM Linear Response Method Distinguishes Ligand Affinities for Closely Related Metalloproteins. Proteins: Struct., Funct., Bioinf. 2007, 69, 326339. 10. Khandelwal, A.; Lukacova, V.; Comez, D.; Kroll, D. M.; Raha, S.; Balaz, S. A Combination of Docking, QM/MM Methods, and MD Simulation for Binding Affinity Estimation of Metalloprotein Ligands. . J. Med. Chem. 2005, 48, 5437–5447. 11. Lu, X.; Gaus, M.; Elstner, M.; Cui, Q. Parametrization of DFTB3/3OB for Magnesium and Zinc for Chemical and Biological Applications. J. Phys. Chem. B 2015, 119, 1062−1082. 12. Li, H.; Ngo, V. A.; Da Silva, M. C.; Callahan, K. M.; Salahub, D. R.; Roux, B.; Noskov, S. Y. Representation of Ion−Protein Interactions Using the Drude Polarizable Force-Field. J. Phys. Chem. B 2015, 119, 9401−9416. 13. Elstner, M.; Cui, Q.; Munih, P.; Kaxiras, E.; Frauenheim, T.; Karplus, M. Modeling Zinc in Biomolecules with the Self Consistent Charge-Density Functional Tight Binding (SCCDFTB) Method: Applications to Structural and Energetic Analysis. J. Comput. Chem. 2003, 24, 565-581. 14. Blomberg, M. R. A.; Borowski, T.; Himo, F.; Liao, R.-Z.; Siegbahn, P. E. M. Quantum Chemical Studies of Mechanisms for Metalloenzymes. Chem. Rev. 2014, 114, 3601−3658. 15. Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem. Int. Ed. 2009, 48, 1198-1229.


22

Page 23 of 40

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


16. Shaik, S.; Cohen, S.; Wang, Y.; Chen, H.; Kumar, D.; Thiel, W. P450 Enzymes: Their Structure, Reactivity, and Selectivity Modeled by QM/MM Calculations. Chem. Rev. 2010, 110, 949-1017. 17. Ngo, V.; da Silva, M. C.; Kubillus, M.; Li, H.; Roux, B.; Elstner, M.; Cui, Q.; Salahub, D. R.; Noskov, S. Y. Quantum Effects in Cation Interactions with First and Second Coordination Shell Ligands in Metalloproteins. J. Chem. Theory Comput. 2015, 11, 4992-5001. 18. Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge Density-Functional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260-7268. 19. Elstner, M. SCC-DFTB: What is the Proper Degree of Self-Consistency? J. Phys. Chem. A 2007, 111, 5614−5621. 20. Gaus, M.; Cui, Q.; Elstner, M. Density Functional Tight Binding (DFTB): Application to Organic and Biological Molecules. WIREs Comput. Mol. Sci. 2014, 4, 49-61. 21. Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the Self-Consistent-Charge DensityFunctional Tight-Binding Method (SCC-DFTB). J. Chem. Theory Comput. 2011, 7, 931-948. 22. Gaus, M.; Goez, A.; Elstner, M. Parametrization and Benchmark of DFTB3 for Organic Molecules. J. Chem. Theory Comput. 2013, 9, 338−354. 23. Gaus, M.; Lu, X.; Elstner, M.; Cui, Q. Parameterization of DFTB3/3OB for Sulfur and Phosphorus for Chemical and Biological Applications. J. Chem. Theory Comput. 2014, 10, 1518−1537. 24. Elstner, M. SCC-DFTB: What is the Proper Degree of Self-Consistency? . J. Phys. Chem. A 2007, 111, 5614−5621. 25. Yang, Y.; Yu, H.; York, D.; Cui, Q.; Elstner, M. Extension of the Self-Consistent-Charge Tight-Binding-Density-Functional (SCCDFTB) Method: Third Order Expansion of the DFT Total Energy and Introduction of a Modified Effective Coulomb Interaction. J. Phys. Chem. A 2007, 111, 10861−10873. 26. Seifert, G.; Joswig, J. O. Density-Functional Tight Binding-An Approximate DensityFunctional Theory Method. WIREs Comput. Mol. Sci. 2012, 2, 456-465. 27. Cui, Q.; Elstner, M. Density Functional Tight Binding: Values of Semi-Empirical Methods in an Ab Initio Era. Phys. Chem. Chem. Phys. 2014, 16, 14368−14377. 28. Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. The Protein Data Bank. Nucleic Acids Res. 2000, 28, 235-242. 29. Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a Sparse Matrix-Based Implementation of the DFTB Method. J. Phys. Chem. A 2007, 111, 5678−5684. 30. Jo, S.; Kim, T.; Iyer, V. G.; IM, W. CHARMM-GUI: A Web-Based Graphical User Interface for CHARMM. J. Comput. Chem. 2008, 29, 1859-1865. 31. Li, H.; Robertson, A. D.; Jensen, J. H. Very Fast Empirical Prediction and Rationalization of Protein pKa Values. Proteins: Struct., Funct., Bioinf. 2005, 61, 704 -721. 32. Jorgensen, W.; Chandrasekhar, J.; Madura, J.; Klein, M. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926-935. 33. Huang, J.; MacKerell, A. D. CHARMM36 All-Atom Additive Protein Force Field: Validation Based on Comparison to NMR Data. J. Comput. Chem. 2013, 0, 00-00. 34. Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I., et al. CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31, 671-690.


23


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

Page 24 of 40

35. van der Spoel, D.; Lindahl, E.; Hess, B.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A. L. T. M.; Feenstra, K. A.; van Drunenand, R., et al., Gromacs User Manual version 4.5.6. Royal Institute of Technology and Uppsala University: Sweden, 2010. 36. Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463-1476. 37. Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952-962. 38. Verlet, L. Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev. 1967, 159, 98-103. 39. Berendsen, H. J. C.; Postma, J. P. M.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684-3690. 40. van Gunsteren, W. F.; Berendsen, H. J. C. A Leap-Frog Algorithm for Stochastic Dynamics. Mol. Simul. 1988, 1, 173-185. 41. Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182-7190. 42. Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An W log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089-10093. 43. 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, 8577-8592. 44. Manzetti, S.; McCullocha, D. R.; Heringtona, A. C.; van der Spoel, D. Modeling of Enzyme-Substrate Complexes for the Metalloproteases MMP-3, ADAM-9 and ADAM-10. J. Comput.-Aided Mol. Des. 2003, 17, 551-565. 45. Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S., et al. CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545-1614. 46. Cui, Q.; Elsner, J.; Kaxiras, E.; Frauenheim, T.; Karplus, M. A QM/MM Implementation of the Self-Consistent Charge Density Functional Tight Binding (SCC-DFTB) Method. J. Phys. Chem. B 2001, 105, 569-585. 47. Kubillus, M.; Kubar, T.; Gaus, M.; Rezac, J.; Elstner, M. Parameterization of the DFTB3 Method for Br, Ca, Cl, F, I, K, and Na in Organic and Biological Systems. J. Chem. Theory Comput. 2015, 11, 332-342. 48. Kollman, P. A. Free Energy Calculations: Applications to Chemical and Biochemical Phenomena. Chem. Rev. 1993, 93, 2395-2417. 49. Brandsdal, B. O.; Osterberg, F.; Almlof, M.; Feierberg, I.; Luzhkov, V. B.; Aqvist, J. Free Energy Calculations and Ligand Binding. Adv. Protein Chem. 2003, 66, 123-158. 50. Deng, Y.; Roux, B. Computations of Standard Binding Free Energies with Molecular Dynamics Simulations. J. Phys. Chem. B 2009, 113, 2234-2246. 51. Zhao, C.; Caplan, D. A.; Noskov, S. Y. Evaluations of the Absolute and Relative Free Energies for Antidepressant Binding to the Aminoacid Membrane Transporter LeuT with Free Energy Simulations. J. Chem. Theory Comput. 2010, 6, 1900-1914. 52. Aqvist, J.; Medina, C.; Samuelsson, J. E. A New Method for Predicting Binding Affinity in Computer-Aided Drug Design. Protein Eng. 1994, 7, 385-391. 53. Carlson, H. A.; Jorgensen, W. L. An Extended Linear Response Method for Determining Free Energies of Hydration J. Phys. Chem. 1995, 99, 10667-10673 54. Lamb, M. L.; Tirado-Rives, J.; Jorgensen, W. L. Estimation of the Binding Affinities of FKBP12 Inhibitors Using a Linear Response Method. Bioorg. Med. Chem. 1999, 7, 851-860.


24

Page 25 of 40

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


55. Almlof, M.; Brandsdal, B. O.; Aqvist, J. Binding Affinity Prediction with Different Force Fields: Examination of the Linear Interaction Energy Method. J. Comput. Chem. 2004, 25, 12421254. 56. Valiente, P. A.; Gil, A.; Batista, P. R.; Caffarena, E. R.; Pons, T.; Pascutti, P. G. New Parameterization Approaches of the LIE Method to Improve Free Energy Calculations of PlmIIInhibitors Complexes. J. Comput. Chem. 2010, 31, 2723-2734. 57. Miranda, W. E.; Noskov, S. Y.; Valiente, P. A. Improving the LIE Method for Binding Free Energy Calculations of Protein−Ligand Complexes. J. Chem. Inf. Mod. 2015. 58. Santos-Martins, D.; Forli, S.; Ramos, M. J.; Olson, A. J. AutoDock4Zn: An Improved AutoDock Force Field for Small-Molecule Docking to Zinc Metalloproteins. J. Chem. Inf. Mod. 2014, 54, 2371−2379. 59. Almlof, M.; Carlsson, J.; Aqvist, J. Improving the Accuracy of the Linear Interaction Energy Method for Solvation Free Energies. J. Chem. Theory Comput. 2007, 3, 2162-2175. 60. Raha, K.; Merz, K. M. Large-Scale Validation of a Quantum Mechanics Based Scoring Function: Predicting the Binding Affinity and the Binding Mode of a Diverse Set of Protein– Ligand Complexes. J. Med. Chem. 2005, 48, 4558-4575. 61. Hou, T. J.; Zhang, Y.; Xu, X. J. Binding Affinities for a Series of Selective Inhibitors of Gelatinase-A Using Molecular Dynamics with a Linear Interaction Energy Approach. J. Phys. Chem. B 2001, 105, 5304-5315. 62. Ito K; Nakajima Y; Onohara Y; Takeo M; Nakashima K Crystal Structure of Aminopeptidase N (Proteobacteria Alanyl Aminopeptidase) from Escherichia coli and Conformational Change of Methionine 260 Involved in Substrate Recognition. J. Biol. Chem. 2006, 281, 33664-33676. 63. Kyrieleis O.J.; Goettig P.; Kiefersauer R.; Huber R.; H.., B. Crystal Structures of the Tricorn Interacting Factor F3 from Thermoplasma acidophilum, a Zinc Aminopeptidase in Three Different Conformations. J. Mol. Biol. 2005, 349, 787–800. 64. Luciani N.; Marie-Claire C.; Ruffet E.; Beaumont A.; Roques B.P. Characterization of Glu350 as a Critical Residue Involved in the N-Terminal Amine Binding Site of Aminopeptidase N (EC 3.4.11.2): Insights Into Its Mechanism of Action. Biochemistry 1998, 37, 686–692. 65. Hangauer D.G.; Monzingo A.F.; Matthews B.W. An Interactive Computer Graphics Study of Thermolysin-Catalyzed Peptide Cleavage and Inhibition by N-Carboxymethyl Dipeptides. Biochemistry 1984, 23, 5730–5741. 66. Rudberg P.C.; Tholander F.; Thunnissen M.M.; Haeggstrom J.Z. Leukotriene A4 Hydrolase/Aminopeptidase. Glutamate 271 is a Catalytic Residue with Specific Roles in Two Distinct Enzyme Mechanisms. . J. Biol. Chem. 2002, 277, 1398-1404. 67. McGowan, S.; Porter, C. J.; Lowther, J.; Stack, C. M.; Golding, S. J.; Skinner-Adams, T. S.; Trenholme, K. R.; Teuscher, F.; Donnelly, S. M.; Grembecka, J., et al. Structural Basis for the Inhibition of the Essential Plasmodium falciparum M1 Neutral Aminopeptidase. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 2537−2542. 68. Jones, P. M.; Robinson, M. W.; Dalton, J. P.; George, A. M. The Plasmodium falciparum Malaria M1 Alanyl Aminopeptidase (PfA-M1): Insights of Catalytic Mechanism and Function from MD Simulations. PloS One 2012, 6, e28589. 69. de Courcy, B.; Gresh, N.; Piquemal, J.-P. Importance of Lone Pair Interactions/Redistribution in Hard and Soft Ligands within the Active Site of Alcohol Dehydrogenase Zn-Metalloenzyme: Insights from Electron Localization Function. Interdiscip. Sci. Comput. Life Sci. 2009, 1, 55-60.


25


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

Page 26 of 40

70. Sakharov, D. V.; Lim, C. Zn Protein Simulations Including Charge Transfer and Local Polarization Effects. J. Am. Chem. Soc. 2005, 127, 4921–4929. 71. Riahi, S.; Roux, B.; Rowley, C. N. QM/MM Molecular Dynamics Simulations of the Hydration of Mg(II) and Zn(II) Ions. Can. J. Chem. 2013, 91, 552-558. 72. Lukacova, V.; Zhang, Y.; Mackov, M.; Baricic, P.; Raha, S.; Calvo, J. A.; Balaz, S. Similarity of Binding Sites of Human Matrix Metalloproteinases. J. Biol. Chem. 2004, 279, 14194–14200. 73. Lukacova, V.; Zhang, Y.; Kroll, D. M.; Raha, S.; Comez, D.; Balaz, S. A Comparison of the Binding Sites of Matrix Metalloproteinases and Tumor Necrosis Factor-a Converting Enzyme: Implications for Selectivity. J. Med. Chem. 2005, 48, 2361–2370. 74. Zoete, V.; Michielin, O.; Karplus, M. Protein-Ligand Binding Free Energy Estimation Using Molecular Mechanics and Continuum Electrostatics. Application to HIV-1 Protease Inhibitors. J. Comput.-Aided Mol. Des. 2003, 17, 861-880. 75. Van Vlijmen, H. W. T.; Schaefer, M.; Karplus, M. Improving the Accuracy of Protein pKa Calculations: Conformational Averaging versus the Average Structure. Proteins: Struct., Funct., Bioinf. 1998, 33, 145-158. 76. Mobley, D. L.; Dill, D. A. Binding of Small-Molecule Ligands to Proteins: "What You See" Is Not Always "What You Get". Structure 2009, 17, 489-498. 77. Kawasaki, Y.; Freire, E. Finding a Better Path to Drug Selectivity. Drug Discovery Today 2011, 16, 985-990. 78. Barrat, E.; Bronowska, A.; Vondrasek, J.; Cerny, J.; Bingham, R.; Phillips, S.; Homans, S. W. Thermodynamic Penalty Arising from Burial of a Ligand Polar Group within a Hydrophobic Pocket of a Protein Receptor. J. Mol. Biol. 2006, 362, 994-1003. 79. Linder, M.; Ranganathan, A.; Brinck, T. Adapted Linear Interaction Energy: A StructureBased LIE Parameterization for Fast Prediction of Protein-Ligand Affinities. J. Chem. Theory Comput. 2012, 9, 1230-1239. 80. Wang, W.; Wang, J.; Kollman, A. P. What Determines the van der Waals Coefficient β in the LIE Method to Estimate Binding Free Energies Using Molecular Dynamics Simulations? Proteins: Struct., Funct., Genet. 1999, 34, 395-402. 81. Velmourougane, G.; Harbut, M. B.; Dalal, S.; McGowan, S.; Oellig, C. A.; Meinhardt, N.; Whisstock, J. C.; Klemba, M.; Greenbaum, D. C. Synthesis of New (-)-Bestatin-Based Inhibitor Libraries Reveals a Novel Binding Mode in the S1 Pocket of the Essential Malaria M1 Metalloaminopeptidase. J. Med. Chem. 2011, 54, 1655–1666. 82. Harbut, M. B.; Velmourougane, G.; Dalal, S.; Reissa, G.; Whisstock, J. C.; Onderd, O.; Brissond, D.; McGowan, S.; Klemba, M.; Greenbaum, D. C. Bestatin-Based Chemical Biology Strategy Reveals Distinct Roles for Malaria M1- and M17-Family Aminopeptidases. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 526-534. 83. Mistry, S. N.; Drinkwater, N.; Ruggeri, C.; Sivaraman, K. K.; Loganathan, S.; Fletcher, S.; Drag, M.; Paiardini, A.; Avery, A. M.; Scammels, P. J., et al. Two-Pronged Attack: Dual Inhibition of Plasmodium falciparum M1 and M17 Metalloaminopeptidases by a Novel Series of Hydroxamic Acid-Based Inhibitors. J. Med. Chem. 2014, 57, 9168-9183.


26

Page 27 of 40

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


TABLES

Table 1. Protein-ligand complexes used in this work and their experimental binding free energies.

a

Ki(nM)

: ∆G789 [kcal/mol]

D66

4381

-10.46

1.8

D50

49081

-8.96

3EBH67

1.6

BES

19067

-9.55

3T8V82

1.8

BTJ

26082

-9.35

4R5T83

2.0

R5T

80083

-8.66

4R5V83

2.1

R5V

70083

-8.74

4R5X83

1.8

R5X

550083

-7.47

PDB ID28

Resolution (Å)

3Q4381

1.8

3Q4481

Ligand IDa

Experimental binding free energies were obtained considering ∆Gexp= -RT ln (1/Ki) where Ki is the inhibition constant.


27


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

Page 28 of 40

Table 2. Chemical structures of PfA-M1 inhibitors used in the training data set. Scaffold

R1

R2

Bestatin-based

-CH3

LigandID

IUPAC name

D66

N-[(2S,3R)-3-AMINO-2HYDROXY-4-(4METHOXYPHENYL) BUTANOYL]-L-LEUCINE

D50

N-{(2S,3R)-3-AMINO-4-[4(BENZYLOXY)PHENYL]-2HYDROXYBUTANOYL}-LLEUCINE

BES

2-(3-AMINO-2-HYDROXY4-PHENYLBUTYRYLAMINO)-4METHYL-PENTANOIC ACID

BTJ

N-[(2-{2-[(N-{(2S,3R)-3AMINO-4-[4(BENZYLOXY)PHENYL]2-HYDROXYBUTANOYL}LALANYL)AMINO]ETHOXY }ETHOXY) ACETYL]-4-BENZOYL-LPHENYLALANYL-N~6~HEX-5YNOYLLYSINAMIDE

-NHOH

R5T

TERT-BUTYL (1-(4-(1HPYRAZOL-1-YL)PHENYL)2-(HYDROXYAMINO)-2OXOETHYL)CARBAMATE

Hydroxamate-based

-NHOH

R5V

N-(1-(4-(1H-PYRAZOL-1YL)PHENYL)-2(HYDROXYAMINO)-2OXOETHYL)PIVALAMIDE

-NHOH

R5X

N-(1-(4-(1H-PYRAZOL-1YL)PHENYL)-2(HYDROXYAMINO)-2OXOETHYL)4-FLUOROBENZAMIDE


28

Page 29 of 40

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


Table 3. Calculated charge for zinc ion in metalloproteins using different QM methods Protein

Ligand type

QM method

qZn (e.u.)

Reference

MMP-2

Hydroxamate

HF/6-31G*

0.55

61

MMP-3

Hydroxamate

B3LYP/6-31G*

1.16

9

MMP-9

Hydroxamate

B3LYP/6-31G*

1.06

9-10

PfA-M1

Hydroxamate/Bestatin

DFTB3/3OB

≈0.67

(this work)

Table 4. Distances (Å) between the ligands’ hydroxyl hydrogen and the closest carboxyl oxygen from residues E497 and E463 before (crystal) and after quantum mechanics minimization (QMmin). E497

E463

PDBID

Inhibitor type

Crystal

QM min

Crystal

QM min

3Q43

Bestatin

1.68

1.06

2.87

*

3Q44

Bestatin

1.54

1.05

2.89

*

3T8V

Bestatin

1.65

1.01

2.69

*

3EBH

Bestatin

1.57

1.50

2.61

*

4R5T

Hydroxamate

2.80

*

1.80

0.99

4R5V

Hydroxamate

2.62

*

1.62

0.99

4R5X

Hydroxamate

1.97

1.03

2.29

*

*Distance > 3Å


29


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

Page 30 of 40

Table 5. Distances (Å) between the bestatin-based ligands’ N-terminal hydrogen and the closest carboxyl oxygen from the residue E463 before (crystal) and after quantum mechanics minimization (QMmin)

PDBID

Inhibitor type

Crystal

QM min

3Q43

Bestatin

1.78

1.01

3Q44

Bestatin

1.65

1.01

3T8V

Bestatin

2.63

1.01

3EBH

Bestatin

1.73

1.00

Table 6. Statistical Metrics for Various Binding Free Energy estimates Method

†

SDerror†

|Emax|†

r2

R2

β (*10-2)

α (*10-2)

γ(10-3)‡

κ†

AutoDock4Zn

1.22

0.61

1.95

0.64

0.41

-

-

-

-

MM-LR

0.58

0.61

1.57

0.32

0.10

-3.22

-11.92

0.21

-8.96

QM/MM-LRA

0.44

0.32

0.90

0.78

0.61

-

0.97

-6.53

-12.40

QM/MM-LRB

0.44

0.32

0.99

0.79

0.62

-

1.22

-3.79

-9.90

QM/MM-LRC

0.30

0.22

0.58

0.91

0.82

-

1.35

-6.58

-12.08

†Reported values are in kcal/mol. ‡Reported values are in kcal/mol*Å2


30

Page 31 of 40

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


Figure captions Figure 1. Inhibitors of PfA-M1 bound to the active site of the enzyme. A) D66 (PDBID: 3Q43), B) D50 (PDBID: 3Q44), C) BES (PDBID: 3EBH), D) BTJ (PDBID: 3T8V), E) R5T (PDBID: 4R5T), F) R5V (PDBID: 4R5V), and G) R5X (PDBID: 4R5X). The gray sphere represents the zinc cofactor. The ligand is shown in licorice representation. The zinc-coordinating residues H496, E497, and H500 and the acidic residues E463 and E497 that interact with the ligand are shown in ball and sticks representation. The cyan, blue and red spheres represent C, N and O atoms respectively. The rest of the protein is shown in new cartoon representation, where the secondary structure elements like α-helices, β-sheets, and loops are coloured in purple, yellow and white respectively.

Figure 2. DFTB3/3OB calculated charges (Mulliken) for the zinc ion and its coordination atoms. In PfA-M1:Inhibitor complexes, the zinc ion is coordinated by three oxygen atoms (two from the ligand, whether bestatin or hydroxamate, and Oδ from E519) and two Nε atoms from H496 and H500, respectively. An hydroxamate-based ligand is shown as example.

Figure 3. Molecular electrostatic potential (MEP) map in electronic units (e.u.) for the active site model of PfA-M1:R5X complex. A) The MEP map plane is defined by the zinc ion and the two coordinating oxygen atoms from the ligand (shown in balls and sticks model). The residues H496 and E519 are above the plane. B) The positively charged zinc ion (blue spot in the center) is surrounded by a high electronic density environment. The MEP maps for other six PfAM1:Inhibitor complexes are similar (not shown).


31


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

Page 32 of 40

Figure 4. Proton transfer between a bestatin-based ligand and the residues E497 and E463. A) Distances from the ligand’s hydroxyl and N-terminal hydrogens to the closest carboxyl oxygen from residues E497 and E463 respectively, in the crystal structure (from PDBID: 3Q43). The arrows represent the direction of proton capture. B) The residues E497 and E463 became protonated after the QM minimization. C, H, O, and N atoms are represented as black, white, red and blue balls, respectively. The side-chains of residues E497 and E463 are represented as acetate groups. The zinc ion and its coordinating residues H496, H500 and E519 were omitted for clarity.

Figure 5. Proton transfer between a hydroxamate-based ligand and the residue E463. A) Distances from the ligand’s hydroxyl hydrogen to the closest carboxyl oxygen from residues E497 and E463 respectively, in the crystal structure (from PDBID: 4R5V). The arrows represent the direction of proton capture. B) The residue E463 became protonated after the QM minimization. C, H, O, and N atoms are represented as black, white, red and blue balls, respectively. The side-chains of residues E497 and E463 are represented as acetate groups. The zinc ion and its coordinating residues H496, H500 and E519 were omitted for clarity.

Figure 6. Scatter diagram of calculated (∆Gcalc) versus experimental (∆Gexp) binding free energies for the seven training data set complexes using AutoDock4Zn (gray dots), MM-LR model (blue dots), and QM/MM-LR models A (green dots), B (orange dots) and C (brown dots).

Figure 7. Comparison of the (single) adopted pose for ligands D66 and R5X from X-Rays derived structures (panels A and B respectively) with respect to their conformational changes in


32

Page 33 of 40

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


PfA-M1 binding site during 50 ns MD simulations (panels C and D respectively). Protein backbone atoms at 4 Å distances from the ligands are shown in green. Carbon, Oxygen and Nitrogen atoms in the ligands are colored in cyan, red and blue, respectively. The zinc ion was omitted for clarity.

Figure 8. Binding free energy contributions from QM/MM-LR model C: QM/MM energies (yellow), SASA term (green) and κ coefficient (brown). The calculated and experimental binding free energies are represented in orange and blue respectively.

Figure 9. A) Scatter diagram of calculated α∆EQM/MM-term values(from model C) versus experimental binding free energies (∆Gexp) for the seven PfA-M1:Inhibitor complexes.B) Scatter diagram of calculated ∆-term values (from model C) versus experimental binding free energies (∆Gexp) for the seven PfA-M1:Inhibitor complexes. C) Scatter diagram of calculated α∆∆EQM/MM-term values (from model C) versus experimental binding free energies (∆∆Gexp) for the seven PfA-M1:Inhibitor complexes relative to the strongest inhibitor (D66).


33


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

Page 34 of 40

Figures Figure 1


34

Page 35 of 40

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


Figure 2

Figure 3


35


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

Page 36 of 40

Figure 4

Figure 5


36

Page 37 of 40

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


Figure 6

Figure 7


37


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

Page 38 of 40

Figure 8


38

Page 39 of 40

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


Figure 9


39


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

Page 40 of 40

For Table of Contents use only Graphic

Improved QM/MM Linear Interaction Energy Model for Substrate Recognition in Zinc-Containing Metalloenzymes. WILLIAMS E. MIRANDA,VAN A. NGO, PEDRO A. VALIENTE, and SERGEI YU. NOSKOV

The zinc-mediated binding was assessed for seven PfA-M1:Inhibitor complexes with a combination of recently developed DFTB3/3OB parameters and Molecular Dynamics Simulations. Charge and proton transfer processes were found to be important for accurate description of binding affinity and specificity in presence of zinc co-factor. The combination of QM/MM-Linear Response models with sampling from MD simulations leads to an excellent performance in predicting binding affinities of the seven inhibitors compared to classical methods. Our results could be useful for designing more potent and selective inhibitors for zinccontaining metalloenzymes.


40

MM Linear-Interaction Energy Model for Substrate

Recommend Documents