A multiscale approach to predict the binding mode of metallo beta-lactamase inhibitors

Antibiotic resistance is a major threat to global public health. β -lactamases, which catalyze breakdown of β -lactam antibiotics, are a principal cause. Metallo β -lactamases (MBLs) represent a particular challenge because they hydrolyze almost all β -lactams and to date no MBL inhibitor has been approved for clinical use. Molecular simulations can aid drug discovery, for example, predicting inhibitor complexes, but empirical molecular mechanics (MM) methods often perform poorly for metalloproteins. Here we present a multiscale approach to model thiol inhibitor binding to IMP-1, a clinically important MBL containing two catalytic zinc ions, and predict the binding mode of a 2-mercaptomethyl thiazolidine (MMTZ) inhibitor. Inhibitors were first docked into the IMP-1 active site, testing different docking programs and scoring functions on multiple crystal structures. Complexes were then subjected to molecular dynamics (MD) simulations and subsequently refined through QM/MM optimization with a density functional theory (DFT) method, B3LYP/6-31G(d), increasing the accuracy of the method with successive steps. This workflow was tested on two IMP-1: MMTZ complexes, for which it reproduced crystallographically observed binding, and applied to predict the binding mode of a third MMTZ inhibitor for which a complex structure was crystallographically intractable. We also tested a 12-6-4 nonbonded interaction model in MD simulations and optimization with a SCC-DFTB QM/MM approach. The results show the limitations of empirical models for treating these systems and indicate the need for higher level calculations, for example, DFT/MM, for reliable structural predictions. This study demonstrates a reliable computational pipe-line that can be applied to inhibitor design for MBLs and other zinc-metalloenzyme systems.


| INTRODUCTION
The spread of antimicrobial resistance (AMR) threatens our capability to treat infectious disease. 1 Microorganisms have developed several resistance mechanisms including efflux pump production, target alteration, decreased drug uptake and drug inactivation. 2 The latter mechanism is one of the most important and is most often seen in the form of the hydrolysis of β-lactam antibiotics, catalyzed by serine-β-lactamases (SBLs) and metallo-β-lactamases (MBLs). 3,4 One effective strategy to overcome this resistance mechanism is to inhibit β-lactamases. [5][6][7][8] However, while for the SBLs several compounds of different classes have already been approved for clinical use, no MBL inhibitor has yet been approved. 7 MBLs can be divided into three subclasses: B1, B2, and B3, with B1 as the most clinically significant because they catalyze the hydrolysis of a broad spectrum of β-lactams 9 and several family members, notably the NDM, IMP and VIM enzymes, are disseminating among Gram-negative bacterial pathogens (e.g., Escherichia coli, Pseudomonas aeruginosa) on mobile genetic elements. 4,10 MBLs contain two catalytic zinc ions within the active site which participate in the reaction mechanism by coordinating the substrate and activating a bound water/hydroxide molecule to undergo nucleophilic attack upon the carbonyl carbon of the scissile amide 4,[11][12][13] (Figure 1). Thus, efforts have been made to design inhibitors that coordinate one or both the zinc ions and compete with substrate for binding to the active site, 11 which could potentially be co-administered with β-lactams to restore their efficacy against MBL-producing strains.
In silico drug discovery to identify MBL inhibitors is more challenging than for SBLs because of the presence of zinc ions. Zinc ions possess d orbitals that participate in chemical bonding and, compared to s or p orbitals, have more complicated shapes. Furthermore, zinc ions can have multiple coordination numbers and different coordination geometries. 15 In molecular mechanics (MM) approaches, metal ions are often represented only as point charges with van der Waals interactions. 16,17 This simplistic representation cannot reproduce all the properties of zinc described above and consequently fails to successfully model the dynamic and versatile nature of zinc coordination in many biological systems. [18][19][20] Over the years there has been significant effort to overcome these limitations by introducing alternative MM models to treat metal ions, including bonded models, nonbonded models, ligand field molecular mechanics, cationic dummy atoms or polarizable force fields.
Bonded models add bonding potentials between the metal ions and ligands and can also include angles, torsions, electrostatic and van der Waals terms; drawbacks include the inability to model ligand exchange and dissociation, or changes in coordination. 20 The ligand field molecular mechanics model incorporates the ligand field stabilization energy directly into the potential energy expression of molecular mechanics. 21 The cationic dummy atom model places additional virtual atoms around the ion to mimic a covalent bond. 22 Polarizable force fields include terms to account for the polarization effect, but so far are computationally expensive. 23 Finally, in typical MM nonbonded approaches, interactions are evaluated by the sum of the Coulomb and Lennard-Jones terms (Equation 1). 24 However, this type of simple model tends to underestimate the strength of the coordination, and does not represent coordination geometry preferences and other electronic effects. 24 For this reason, other nonbonded interaction models have been proposed, such as the restrained nonbonded model or the LJ12-6-4 model (Equation 2), which has been tested for both Mg 2+25, 26 and Zn 2+ ions 27,28 and contains an r À4 term that is included to describe the ion-induced dipole interaction. 29 In this study, we investigate the binding of inhibitors to imipenemase (IMP)-1, a B1 MBL that is frequently found on mobile genetic elements that facilitate its spread through bacteria. 10,[30][31][32] The primary host, Pseudomonas aeruginosa, is a Gram-negative F I G U R E 1 Possible mechanism of β-lactam ring hydrolysis by subclass B1 metallo-β-lactamases. 12 of immunocompromised patients, and notably associated with infections of burn wounds and the lungs of individuals with cystic fibrosis. 33 In the continued absence of new antibacterial agents with anti-Gram-negative activity, infections by resistant P. aeruginosa strains are very difficult to treat. 34 Various compound classes have been investigated as potential MBL inhibitors (as reviewed 4,7,9,11,35 ); we here focus on thiol-containing compounds which exploit the high affinity of zinc for sulfur ligands and have been extensively investigated. 36,37 In consequence, relatively extensive structural data exist for complexes of IMP-1 with thiol compounds that may be exploited for development and assessment of approaches to model inhibitor binding. 36 Figure S1) to find the docking setting that best reproduced their interactions. In each case the inhibitor was modeled and docked as its thiolate anion. The region of the protein included in docking was selected by considering whole residues within a 10 Å radius around the two zinc ions of the IMP-1 active site, which comprises the whole catalytic pocket. Docking with PLANTS 41 (see below) applied distance restraints between the sulfur atom of the ligand and the zinc ions with a weight of À7.5 kcal/mol, and a maximal and minimal distance of  GAFF as force field, 44 for which benchmarking in agreement with quantum calculation has been reported, 45 while the protein was parameterized using the ff14SB force field. 46 Given the challenges associated with evaluating the interactions between the zinc ions and their coordinating atoms, we tested two different nonbonded models. Lennard-Jones equation 28 to represent the charge-induced dipole interaction, and for which no distance restraints were applied.

| Molecular dynamics simulations
We applied the SHAKE algorithm 47 and a time step of 2 fs. The simulation protocol included energy minimization of the hydrogen atoms followed by minimization of the water molecules, and finally of the whole complexes, using 100 steps of steepest descent and 2900 steps of the conjugate gradient algorithm. The system was then heated to 300 K using the Langevin thermostat for 20 ps. Equilibration was performed at a pressure of 1 atm with the Berendsen barostat for 240 ps. Finally, we computed the production run using the NVT ensemble over 100 ns.

| QM/MM optimization
For QM/MM calculations we tested an approximate DFT-based method, SCC-DFTB (self-consistent charge-density functional tight binding, abbreviated here as DFTB, in particular DFTB3 that is a third order expansion of DFTB 48 ) and a DFT method. 49  predicting accurate structures. 54 We considered the zinc ions as closed shells. Covalent C-C bonds between the MM and the QM regions were capped with hydrogen atoms.

| Method validation
We envisaged a workflow for developing models of IMP complexes with thiol-based inhibitors that involved successive steps (docking, MD, QM/MM), of higher levels of theory, to improve the accuracy of the final model. As described in Methods, for each step we tested multiple settings, in order to find the best combination to describe the structure and interactions with the metal ions.
We first tested docking options and scoring functions for thiolate anion-containing inhibitors (Table 1). Molecular docking predicts the binding mode of a ligand to a target molecule, usually a receptor or an enzyme. The resulting poses are generated by the combination of the conformational sampling algorithm, which explores the conformers of the ligand, and the scoring function that ranks those poses. When metal ions are included in the docking region, the scoring function needs to address a further level of complexity that can lead to the failure of the prediction of the correct binding mode. 55 We tested two docking tools that use two different search algorithms: the ant colony optimization (ACO) and the genetic algorithm (GA), namely PLANTS 45  Note: Equilibrium distances of the restrained nonbonded model between the zinc ions and the protein coordinating atoms are reported in Å.
Literature and database searches identified crystal structures of IMP complexes with seven different thiols ( Figure S1), to which we added our own recently determined structures of IMP complexes with two different MMTZ stereoisomers/diastereomers. 39   The results of the best docking runs (Opt3) for the two MMTZ compounds for which complex structures were available were then subjected to molecular dynamics (MD) simulations to improve the quality of the models obtained. These were run using either the restrained nonbonded (LJ12-6-R) model or the LJ12-6-4 model; averaged RMSD values for the production run frames are reported in Table 3. In the two MD simulations the average RMSD values (compared to crystal structures) of the binding site residues for the LJ12-6-4 model are consistently lower than those for the LJ12-6-R, although in the latter the inhibitor RMSD is slightly smaller, suggesting slightly greater stability of the complex. In general, for both treatments, the difference in coordination bond lengths compared to the crystal structure is small (Figure 3).
In the crystal structures of both MMTZ complexes the two zinc ions are each in a tetrahedral geometry, with the first (ZN1) coordinated by three histidine residues (His77, His79, His139) and by the thiol moiety of the inhibitor, while the second (ZN2) is coordinated by Cys158, Asp81, His197, and by the inhibitor thiol, which, as in other crystallographic complexes of IMP-1 with thiol-based inhibitors, 36 replaces the zinc-bridging hydroxide present in unliganded structures. 65 However, during MD simulations using either model, the coordination number of the zinc ions tends to increase, through addition of extra water molecules, up to a limit of six. This tendency is enhanced when the LJ12-6-4 model is used. It is worth noting that in general the distance between the two zinc ions did not change significantly (i.e., the largest difference was 0.77 Å, found in the L-anti-1b LJ12-6-4 simulation).
Finally, we refined the inhibitor complex models with QM/MM geometry optimization. In the refinement step we tried two different methods: DFTB3 and DFT, and in the latter tested two functionals (B3LYP and PBE0). The DFTB3/MM method was tested to establish whether it could improve structures from MM MD (e.g., whether it would decrease the RMSD and recover the tetrahedral geometry of the zinc ions that was observed in the crystal structures). Thus, the complexes underwent a DFTB3/MM geometry optimization. We found no or only slight improvement in RMSD values (a decrease of 0.43 Å in the RMSD of the binding site in D-syn-1b) (Table S1, Figure S6), but the geometry of the zinc ions was not recovered: they remained octahedrally coordinated ( Figure S7), even if DFTB3 was parameterized considering tetrahedral zinc-proteins as well. 66 Accordingly, DFTB3/MM was applied for other zinc tetrahedral complexes, but with further refinement with DFT. 67 DFT methods were then employed to provide a more accurate description of the metal ions, with application of the hybrid functionals B3LYP and PBE0 which are commonly used in QM/MM studies of zinc-containing proteins. 68 Table 3.

| Prediction of D-anti-1a MMTZ inhibitor binding mode to IMP-1
Previous study identified that a range of MMTZs, 39 including, but not limited to, D-syn-1b and L-anti-1b, act as micromolar inhibitors of IMP-1. However, in some cases crystal structures of inhibitor complexes with IMP-1 could not be obtained, even though biochemical data indicated that potency was comparable to that of compounds whose IMP-1 complexes could be structurally characterized. Accordingly, we sought to apply the approach developed here to model interactions of an additional MMTZ inhibitor, D-anti-1a, with IMP-1.

D-anti-1a is, respectively, a diastereoisomer and enantiomer of D-syn-
1b and L-anti-1b. As it was unclear which of these complexes was likely to better accommodate the binding of D-anti-1a, docking experiments were carried out using IMP-1 conformations derived from both the D-syn-1b and L-anti-1b crystal structures, and subsequent simulations were undertaken in parallel.
Based upon the results above, molecular docking of D-anti-1a into IMP-1 in both conformations (i.e., from the D-syn-1b and L-anti-1b complexes) used the Opt3 approach described above. Figure 5 shows the binding modes obtained using the two crystal structures.  Table 4. We obtained higher RMSD values (of both the ligand and the binding site) when docking was based on the L-anti-1b structure, suggesting that, in this structure, the protein was not in the optimal conformation to accommodate and establish interactions with D-anti-1a. In contrast, when based on the D-syn-1b structure, lower RMSD values of the binding site were obtained and the ligand remained more stable, with an RMSD along the trajectory that fluctuated around 0.5 Å ( Figure S8).
These values were little changed after QM/MM optimization using the B3LYP hybrid functional with the 6-31G(d) basis set ( Table 4). Inspection of coordination distances (Table S2)  During the MD simulation rearrangement of the loop L3, which contains Trp28, is observed, as well as of loop L10. This flexibility has been described in other studies and is consistent with the experimental evidence. However, in contrast to the behavior of these loops, the binding site of the enzyme remained stable for the entire trajectory, with average RMSD values below 2 Å. This is clearly evident from inspection of the RMSF values for these systems ( Figure S9) in which high RMSF values were obtained for loop L3 (residues 21-32) and loop L10 (residues 162-172), while low values were obtained for the binding site (residues: Phe51, His77, His79, Ser80, Asp81, His139, Thr140, Cys158, Ser196, His197).

| DISCUSSION
The drug discovery process often relies upon molecular docking and MD simulations, which can usefully be combined. 75  est. This indicates that the error in the final structure is low and that the approach is reasonable and suitable for application to simulate similar systems. However, we would also stress that, for zinc (and other metal) centers such as those explored here, RMSD values alone do not provide a sufficient criterion for accuracy (although they can be sufficient to demonstrate failure of a method) and, as detailed above, nonphysiological geometries such as octahedral zinc coordination and/or adoption of metal-bridging by Asp81 can occur even when overall RMSD values are low. These factors should be considered in general in assessing the quality of predicted zinc metalloprotein complexes.
Notwithstanding these considerations, our workflow ( Figure S10) produced a model for an inhibitor complex that is biochemically realistic both with respect to zinc coordination and interactions with active site residues (Lys161) known to be important to ligand binding by IMP-1. 61,86,87 Involvement of Lys161 in our model of D-anti-1a binding by IMP-1 is notable given that such interactions are not observed in the crystal structures of either the D-syn-1b or L-anti-1b complexes. We have proposed that S-π interactions involving the thiazolidine ring sulfur and conserved aromatic residues at position 51 contribute to the affinity of MMTZs for subclass B1 MBLs such as IMP-1. 39 However, the displacement of Trp28 and reorientation of the MMTZ thiazolidine observed in our model suggests that these may not be maintained in the IMP-1:D-anti-1a complex. In this case, the loss of S-π interactions is then compensated for by the addition of a hydrogen bond to Lys161, ensuring that comparable inhibition potency is retained. 39

| CONCLUSIONS
In this study, we tested several techniques and their ability to reproduce inhibitor binding modes observed in crystal structures. However, we also identified methods that failed to correctly handle the metal coordination sphere. In detail, docking without the inclusion of crystallographic water molecules was not able to produce ligand binding poses with RMSD values lower than 2 Å compared with the corresponding original crystal structure. The results also demonstrate the limitations of standard MM models for MD simulations of these systems. In particular, MD simulations using a LJ12-6-4 model increased the coordination number of the zinc ions to six; this model may be well suited for use with zinc sites with octahedral coordination geometry, 25 In summary, this study presents a validated, multiscale approach for the computational prediction of IMP-1 interactions with thiolate inhibitors. This approach is also suitable for application to studies of other MBL inhibitors, and other zinc-containing proteins.