Refined models of New Delhi metallo-beta-lactamase-1 with inhibitors: an QM/MM modeling study

New Delhi metallo-beta-lactamase 1 (NDM-1) has been identified as a potential target for the treatment of multi-drug resistance bacterial infections. We used molecular docking, normal MD, SIE, QM/MM MD simulations, QM/MM GBSA binding free energy, and QM/MM GBSA alanine-scanning mutagenesis techniques to investigate interactions of the NDM-1 with 11 inhibitors (Tigecycline, BAL30072, D-captopril, Penicillin G, Ampicillin, Carbenicillin, Cephalexin, Cefaclor, Nitrocefin, Meropenem, and Imipenem). From our normal MD and QM/MM simulations, the correlation coefficients between the predicted binding free energies and experimental values are .88 and .93, respectively. Then simulations, which combined QM/MM/GBSA and alanine-scanning mutagenesis techniques, were performed and our results show that two residues (Lys211 and His250) have the strongest impact on the binding affinities of the 11 NDM-1/inhibitors. Therefore, our approach theoretically suggests that the two residues (Lys211 and His250) are responsible for the selectivity of NDM-1 associated inhibitors. Graphical abstract Difference in the predicted binding free energy of the virtual alanine mutations of the NDM-1 with the 11 inhibitors.


Introduction
Antimicrobial resistance has been recognized since the beginning of the antibiotic era in the mid-twentieth century. The history of β-lactamases, starting from the discovery of the TEM in 1965, has demonstrated how the microbial world keeps pace with technical advances. The metallo β-lactamases are more versatile enzymes that can convert the host bacteria into almost total β-lactam insusceptibily. Currently, multidrug-resistant gram-negative bacteria cannot be effectively treated with proven therapeutics (Giamarellou & Poulakou, 2009;Giamarellou, 2010). This problem is a great risk to human health. The function of beta-lactamases is to cleave the amide bond of the beta-lactam ring so as to inactivate the beta-lactam inhibitors (Livermore, 2009). Thus, beta-lactamases have been classified as the point of resistance for beta-lactam inhibitors. Beta-lactamases are divided into four classes (A, B, C, and D). Classes A, C, and D enzymes contain serine groups in their active sites. Class B (metallo-betalactamases) enzymes require one or two zinc ions for their activity. In 2009, a Swedish patient of Indian origin traveled to New Delhi and acquired a urinary tract infection caused by a carbapenem-resistant Klebsiella pneumoniae strain (Yong et al., 2009). Since then, thousands of infections in humans by the bacteria have providing new inhibitors for NDM-1 is an important issue. The NDM-1 structures have been determined by X-ray experiments, Kim et al., 2011) and we chose the structures (PBD ID: 4RBS, 4RAW, and 4HL2) as the target structures for our studies. From the incomplete crystal information of NDM-1/inhibitors complex structures (Hydrolysis of Penicillin G and Ampicillin), the active site residues of NDM-1 are the two zinc atoms, Leu65, Met67, Gln123, Asp124, Lys211, and Asn220 (Table 1). The inhibiting NDM-1 abilities of 11 inhibitors (Tigecycline, BAL30072, Dcaptopril, Penicillin G, Ampicillin, Carbenicillin, Cephalexin, Cefaclor, Nitrocefin, Meropenem, and Imipenem) have been reported Kumarasamy et al., 2010;Walsh, 2011). Thus, we used these reported results to verify our simulations.
Computational alanine-scanning mutagenesis (Chiappori, D'Ursi, Merelli, Milanesi, & Rovida, 2009) methods can predict important residues in the inhibitors/ antigen/antibody/peptides-associated protein interactions (Tokunaga, Arakawa, & Tokunaga, 2008). Therefore, our study carried out quantum mechanics/molecular mechanics (QM/MM) GBSA computational alanine-scanning mutagenesis simulations for key residues of the inhibitors/NDM-1 complex structures. Researchers have successfully used the computational alanine-scanning mutagenesis method to predict the D120 N mutation effect and the hydrolysis mechanism of NDM-1 (Chen et al., 2013;Zhu et al., 2013). To date, detailed information of the antibiotic/key residues of NDM-1 is still lacking. Accordingly, this method was used to study the mutation effects of the inhibitors/NDM-1 complex structures.
Due to the limitation of molecular mechanics methods alone, combined QM/MM methods can more accurately predict the enzyme-inhibitor interactions (Barnett & Naidoo, 2010;Mulholland, 2005Mulholland, , 2007Mulholland, Lyne, & Karplus, 2000;Solt et al., 2009;Warshel & Levitt, 1976) by allowing those enzyme-inhibitor interactions to be modeled: a small region (such as inhibitors, ions, and active site residues) is simulated by a quantum mechanical method; the other regions (such as protein and solvent molecules) are performed with an empirical force field (Claeyssens, Ranaghan, Manby, Harvey, & Mulholland, 2005;Friesner & Guallar, 2005;Lyne, Mulholland, & Richards, 1995;Mulholland, 2007;Mulholland et al., 2000;Ridder, Harvey, Rietjens, Vervoort, & Mulholland, 2003). QM theory generally limits QM/MM investigations. While semi-empirical methods provide fast computation speed for QM/MM molecular dynamics simulations, those methods are error-prone. DFT/HF methods can offer improved accuracy, however, and have opened new classes of enzymes to computational investigation (Bathelt, Zurek, Mulholland, & Harvey, 2005;Himo & Siegbahn, 2003). The original DFT/HF methods, however, are rather time-consuming and require significant computing power (Rathore et al., 2013). A modified DFT approach entitled SCC-DFTB has been developed by Prof. Cui, Elstner, Kaxiras, Frauenheim, and Karplus (2000), and Prof. Elstner et al. (2003) for reducing computing time and power (Elstner et al., 1998). It is based on a second-order expansion of the Kohn-Sham total energy with respect to charge density variations in the LCAO framework. This method has been applied to calculate the energies, geometries, and vibrational frequencies of small organic molecules, and the resulting mean average deviations from the experimental values were comparable to those of full density functional theory calculations with a double-ζ plus polarization basis set. This method has been tested for biologically relevant molecules, including structures of enzyme-inhibitor interactions (Elstner, Frauenheim, Kaxiras, Seifert, & Suhai, 2000).
Therefore the SCC-DFTB method is used for our QM/MM MD simulations. The approach consists of the following steps: Step 1. Perform molecular docking and normal 100-ns MD simulations.
Step 2. Performance SIE (solvated interaction energies) alanine-scanning mutagenesis free energies calculations to find out the binding hot-spots residues of the 11 complex binding modes. Step 3. Perform 1-ns QM/MM MD simulations and define the binding modes between the 11 NDM-1/ inhibitors.

Molecular docking and initial NDM-1/inhibitors complex structures
Eight inhibitors (Tigecycline, BAL30072, D-captopril, Carbenicillin, Cephalexin, Cefaclor, Nitrocefin, and Imipenem) were optimized using SYBYL (Tripos, 2009) software. The NDM-1 active site residues Kim et al., 2011) is a deep cavity that is formed by the loop regions between β5-α2 and β10-α4, in which the catalytic zinc ions are deeply buried. This cavity is clamped by the loop regions between the three regions (Thr119-Met126, Ser217-Asp225, and Met67-Gly71). These residues were chosen for molecular docking simulations. AutoDock Vina (Huey, Morris, Olson, & Goodsell, 2007) docking software was used to dock these eight inhibitors into the active site residues of the NDM-1 structure (PDB ID: 4RBS). For every inhibitor, the 2000 conformations of inhibitors were generated and then scored by the AutoDock score. The binding modes of the best AutoDock score were selected for subsequent 100-ns MD simulations. For further three inhibitors (Penicillin G, Ampicillin, and Meropenem), the crystal structures (PDB ID: 4HL2 for Ampicillin; 4RBS for Meropenem; 4RAW for Penicillin G,) were chosen as the initial complex structures for computational simulations. For the Penicillin G, the crystal structure (PDB ID: 4RAW) was chosen as the initial complex structure. Because this NDM-1 has been mutated (M67 V) and the Penicillin G is hydrolyzed, we used SYBYL (Tripos, 2009) software to modify the complex structure (V67 M) and the Penicillin G. For the Ampicillin, the crystal structure (PDB ID: 4HL2) was chosen as the initial complex structure. Due to the hydrolyzed Ampicillin structure, we also used SYBYL software to modify the complex structure. For the Meropenem, the crystal structure (PDB ID: 4RBS) was chosen as the initial complex structure. Due to the hydrolyzed Ampicillin structure, we also used SYBYL software to modify the complex structure. Then, these complex structures were selected for subsequent the 100-ns normal MD. From the binding mode analysis of the 11 NDM-1/inhibitors, our simulations results indicate that the residues (Zn atoms, Ala74, Asn220, Asp124, Cys208, Gln123, His122, His189, His250, Ile35, Leu65, Lys211, Lys73, Met154, Met67, Phe70, Trp93, and Val73) were chosen to SIE alaninescanning mutagenesis free energies calculations.

SIE alanine-scanning mutagenesis free energies calculations
The SIE method is a kind of the linear interaction energy methods (Perdih, Bren, & Solmajer, 2009). The binding free energies of the 11 NDM-1/inhibitors were calculated for snapshot structures. The SIE (Naïm et al., 2007) function is written as: The terms E C and E vdw are the intermolecular Coulomb and van der Waals interaction energies, respectively. The two terms were obtained using the AMBER force field (FF99SB). The term DG R bind is the reaction field energy change between the free and bound states, and the term is estimated by solving the Poisson equation. The terms of BRI and BEM (Purisima, 1998;Purisima & Nilar, 1995) are the molecular surfaces by generating with a variable-radius solvent probe (Bhat & Purisima, 2006). The ΔMSA term is the change in the molecular surface area upon binding. The constant values of these parameters are α = .1048, D in = 2.25, ρ = 1.1, γ = .0129 kcal/(mol Å 2 ), and C = −2.89 kcal/mol (Perdih et al., 2009). Our results also indicate that the 8 residues (His120, His122, His189, His250, Lys211, Trp93, Val73, and Zn atoms) are hot-spots residues of NDM-1 proteins. Therefore, the 8 residues and the 11 inhibitors were chosen to QM regions in the QM/MM MD simulations.

Computational QM/MM GBSA alanine-scanning mutagenesis
The binding free energy is calculated according to The free energies of the three parts (inhibitors, NDM-1, and NDM-1/inhibitors) are estimated from The term E QM/MM is the molecular mechanics gas phase energy, the term G pol is the polar contribution to the solvation free energy (Bashford & Case, 2000), the term G np is the non-polar solvation free energy, T is the absolute temperature of the system, and S QM/MM is the entropy. The E QM/MM is shown below: The terms E cov , E elec , E vdw , and E QM are the covalent, electrostatic, van der Waals, and quantum interaction energies, respectively. The E cov includes terms representing bond stretching (E bond ), angle vibrational (E angle ), and the dihedral angle torsion energies (E dihedral ), according to The entropy (S QM/MM ) consists of the entropy of the quantum calculations system (S QM ), and the entropy of the molecular mechanism (S MM ). This equation is according to The G np is estimated from Sitkoff's contributions, as (Sitkoff, Sharp, & Honig, 1994) with γ = .00542 kcal/mol −1 Å −2 and b = 0. The term A is the solvent accessible surface area which was estimated with a fast linear combination (Weiser, Shenkin, & Still, 1999). The term ΔG bind is obtained by averaging over an ensemble of the complex structural configurations taken from a QM/MM MD simulation. The binding free energies are defined as: (9) The computational alanine-scanning mutagenesis approach can accurately predict the hot-spots residues of proteins. The method presented can achieve an overall success rate of 80% and an 82% success rate in residues for which alanine mutation causes an increase in the binding free energy >1.0 kcal/mol (warm-and hot-spots) (Ramos & Moreira, 2013). The Amber14 and AmberTools14 software (sander.MPI and MMPBSA.py) were used to determine the binding free energies/computational alanine-scanning mutagenesis free energies of the complex structures.

Molecular docking and initial NDM-1/inhibitors complex structures
The 11 inhibitors were docked into the NDM-1 structure. After the 100-ns MD trajectories of the NDM-1 with solvent molecules and the 11 inhibitors, the results were used to perform the 100-ns normal MD simulations.
frequently. Figure 1 shows that the RMSD profiles of the 11 NDM-1/inhibitors are equilibrated at 20-ns. From each of the 80-ns MD trajectories of the NDM-1/inhibitors, 800 snapshots were taken at regular intervals for the SIE binding energy analyses.

SIE alanine-scanning mutagenesis free energies calculations
The simulations revealed satisfactory correlations between the calculated and experimental binding affinities of all 11 inhibitors, as indicated by the obtained R 2 value of .88 (Table 3 and Supplementary Figure S12). The residues (Ala215, Asn193, Asn220, Asp124, Asp212, Asp223, Cys208, Gly192, Gly219, Hie162, Hie223, His120, His122, His162, His189, His250, His95, Ile35, Ile8, Leu218, Lys184, Lys211, Lys216, Ser217, Ser251, and Val73) were mutated to alanine for computational alaninescanning mutagenesis calculations. The binding free energy of each inhibitor was obtained from the equilibrated 80-ns MD simulation and the SIE method, with both processes using the same parameters. All the results were listed in Figure 2. From Figure 2 and Supplementary  (Table 3 and Supplementary Figure S12). Our QM/MM GBSA computational alanine-scanning mutagenesis results are shown in Figure 3. For the computational alanine-scanning mutagenesis of NDM-1/11 inhibitors, the ΔG Diff of the two residue mutations (Lys211 and His250) can affect the binding affinities of the NDM-1/inhibitors except Nitrocefin. Our results show that the two residue mutations (Val73 and Trp93) cannot affect the binding affinities. The three residues also can affect the binding affinities depend on inhibitors. Accordingly, the two residues (Lys211 and His250) most affect the binding interactions among the 11 inhibitors.

Conclusion
In order to gain insight into the effect of the 11 inhibitors (Tigecycline, BAL30072, D-captopril, Penicillin G, Ampicillin, Carbenicillin, Cephalexin, Cefaclor, Nitrocefin, Meropenem, and Imipenem)/NDM1 complex structures on binding at the atomic level, we used the AutoDock Vina program, normal MD simulations, SIE free energies prediction method, QM/MM MD simulations techniques, QM/MM GBSA alanine-scanning mutagenesis method to predict the binding modes for 11 inhibitors (Tigecycline, BAL30072, D-captopril, Penicillin G, Ampicillin, Carbenicillin, Cephalexin, Cefaclor, Nitrocefin, Meropenem, and Imipenem) interacting with NDM-1. From our normal MD and QM/MM simulations, the correlation coefficients between the predicted binding free energies and experimental values are .88 Note: Here, the difference in the binding free energy is: ΔG Binding (alanine-scanning mutagenesis) -ΔG Binding (Table 3).