3D-QSAR and molecular recognition of Klebsiella pneumoniae NDM-1 inhibitors

ABSTRACT New Delhi metallo-β-lactamase-1 (NDM-1) as a target for the development of anti-superbug agents, plays an important role in the resistance of β-lactam antibiotics and has received worldwide attention. Sulfhydryl propionic acid derivatives can effectively inhibit the catalytic activity of NDM-1, but the quantitative structure–activity relationship (QSAR) and inhibitor-target recognition mechanism both remain unclear. In this work, CoMFA and CoMSIA models of sulfhydryl propionic acid inhibitors with high predictive ability were obtained, from which the effect of different substituents on the inhibitory activity against NDM-1 were revealed at the molecular level. Then, two 120-nanosecond comparative molecular dynamics (MD) simulations for NDM-1 enzyme and NDM-1-inhibitor complex systems were performed to study the recognition and inhibition mechanism of sulfhydryl propionic acid derivatives. The binding of inhibitors alters the entire H-bond network of the NDM-1 system accompanied by the formation of strong interactions with I35, W93, H120, H122, D124, H189 and H250, further weakens the recognition of NDM-1 by its endogenic substrates. Finally, the results of free energy landscape and conformation cluster analyses show that NDM-1 underwent a significant conformational change after the association with sulfhydryl propionic acid inhibitors. Our findings can provide theoretical support and help for anti-superbug agents design based on the structures of NDM-1 and sulfhydryl propionic acid derivatives.


Introduction
β-Lactamase inhibitors are the most common type of antibiotics in clinical application, which are mainly used to treat pneumonia, urinary tract infection, abdominal cavity infection, blood infection, and neutropenia lack accompanied by fever diseases resulted from bacterial infection [1][2][3][4]. Due to antibiotics abuse and bacterial variation, drug resistance becomes more and more prominent [5,6]. So far, the infections caused by multidrug-resistant (MDR) bacteria have been a serious threat to patients' lives for lack of effective drug treatment [7,8].
In drug-resistant bacteria, β-lactamase can specifically and efficiently cleavage the β-lactam bond, which is one of the main drug-resistance mechanisms against β-lactam antibiotics [9][10][11][12][13]. According to Ambler classification, β-lactamase can be divided into four types [14] (i.e. A, B, C and D), among which the functional residue in the activity centre of A, C and D enzymes is serine, known as serine β-lactamase (SBL). As for B class, the functional one is bonding metal ion(s), thus also called metallo-β-lactamase (MBL) [15,16]. Based on the difference of amino acid sequence, B class can be further divided into three sub-categories: B1, B2 and B3. In the clinical treatment of the drug-resistance bacteria of β-lactam, the combined use of different β-lactamase inhibitors can partially restore the efficacy of β-lactam antibiotics. The research and development of β-lactamase inhibitors is quite fruitful, and series of drugs including avibactam, potassium clavulanic acid, sulbactam and tazobactam have been widely used in clinical practice. Nevertheless, most clinical agents only exhibit good inhibitory potency of SBLs, and still show no activity to MBLs [17,18], thus the development of novel β-lactamase inhibitors remains urgent.
New Delhi metallo-β-lactamase-1 (NDM-1) belongs to B1 class MBLs, has a αβ/βα interlayer structure and can hydrolyze a variety of β-lactam antibiotics including penems [15,16,19]. The active site of NDM-1 is located at the bottom of the small groove formed by two β-sheets. Compared with other MBLs, NDM-1 has a larger surface area and relies on two bonded zinc ions (i.e. Zn1 and Zn2) to exert catalytic functions [19]. Specifically, Zn1 is coordinated with H120, H122 and H189, then combines with a bridge OHto form a tetrahedral structure; Zn2 is coordinated with residues D124, C208 and H250, and connects to Zn1 with the help of the bridge OH -, maintaining a triangular cone configuration [20,21].
Based on structural features of NDM-1, some research groups have carried out the design and screening of NDM-1 inhibitors and reported many highly-active compounds with novel structure, such as sulfhydryl thiophene [22][23][24][25][26], benzopyran [27][28][29], covalently bonded [30,31], metal chelating agents [32][33][34] and other classes [35][36][37]. Among them, sulfhydryl propionic acid derivatives have gained more attention due to their structural diversity, good inhibitory activity, and high druggability [38][39][40]. Thomas et al. [41] used a flux screening model and obtained two NDM-1 inhibitors with sulphur alcohol groupp-chloromercuribenzoic acid (p-CMB) and nitroprussidefrom the LOPAC compound library. In addition, other 3-sulfonamides NDM-1 inhibitors were successfully picked out using high-throughput screening model by Wang's group [26]. In 2016, King and co-workers [23] conducted gene knockout experiment to single out the metalloβ-lactamase inhibitors -Aspergillomarasmine A (AMA). Li and Shen [24,25] designed and synthesised active 3-mercaptopropionic acid and 3-mercapto propionamide derivatives, which might be novel potential antimicrobial agents effectively inhibiting NDM-1. These previous theoretical and experimental works enhance our understanding of active data and screening models for sulfhydryl propionic acid derivatives (see fig.  S1). Nevertheless, the structure-activity relationship (SAR) analyses of this kind of compounds, serving as an important strategy in drug development, have not been carried out systemically using molecular simulation method. Other important scientific issues such as molecular recognition mechanism between inhibitors and NDM-1, as well as the details of conformational change remain to be further explored.
In this work, the structures and inhibitory activity data for a total of 23 sulfhydryl propionic acid inhibitors both were firstly collected to analyse the three-dimensional quantitative structure-activity relationship (3D-QSAR). Next, molecular docking, molecular dynamics (MD) simulation, MM/PBSA, free energy landscape and conformational cluster methods all were employed to reveal molecular recognition and inhibition mechanism of sulfhydryl propionic acid inhibitors (see Figure  1). These findings not only enhance our understanding of molecular recognition between NDM-1 and sulfhydryl propionic acids, but also provide a theoretical guidance for drug design and structural modification of NDM-1 inhibitors based on the structures of receptor and ligand.

Molecular docking
The 23 sulfhydryl propionic acid derivatives were constructed using ChemBio3D Ultra 12.0 [42], and roughly optimised under the MM2 molecular force field with the convergence of root mean standard (RMS) less than 0.0001 kcal·mol −1 ·Å −1 . Then, precisely geometric optimisation was carried out using the density functional theory (DFT) method at B3LYP/6-31G (d) [43] theoretical level with Gaussian 09 [44] software package to generate the initial conformations for the following molecular docking.
Molecular docking was conducted with AutoDock 4.2 software package [45]. In order to more accurately reflect the importance role of Zn (II) in molecular recognition, AutoDock4 Zn force field [46] was adopted. The coordinates of Zn (II) was set as the centre of docking box, where the lattice points were spaced at 3.75 × 10 −2 nm and the amounts of grid points were 40. The Lamarckian genetic algorithm (LGA) was used to generate and optimise the docking conformation. Based on LGA in combination with the energy of global and local search results, semi-empirical free energy evaluation function was adopted to track and assess molecular recognition between receptor and ligand. The semi-empirical free energy is composed of 4 terms, namely molecular internal energy, Van der Waals interactions, H-bonds and electrostatic interactions. In order to fully take into account the flexibility of small molecules, all the rotations of single bonds were fully sampled during molecular docking. Total 128 inhibitor conformations were collected for each docking and the snapshot with the lowest energy in the largest cluster was defined as the near-native docked result.

3D-QSAR
The 23 NDM-1 inhibitor compounds were randomly divided into a training set (17 molecules) and a test set (6 molecules), and the experimental IC 50 values [24,25] were converted into negative logarithmic form (pIC 50 ). Before constructing 3D-QSAR model, it is necessary to align molecules well. In general, the same type of inhibitor derivatives has similar binding modes with their receptor. In the superposition step of QSAR in this work, the compound 17 (Cmpd # 17) with the crystal structure was chosen as template molecule, and the common thiol substructure was defined as alignment rule.
All the 3D-QSAR operation was performed with SYBYL-X 2.0 package [47]. In the process of establishing competitive molecular field analysis (CoMFA) model, the aligned NDM-1 inhibitors were placed in the spatial grids, where sp 3 hybridised probe particles (such as: C + , H 2 O, CH 4 , H + , etc.) were used to calculate the interaction force between the probe and inhibitors. Each monitoring point was placed every 0.2 nm, then the data of the steric field (S) and the electrostatic field (E) were obtained by recording different spatial coordinates. Compared with CoMFA, comparative molecular similarity indices analysis (CoMSIA) model not only analyses S and E, but also detects the hydrophobic field (HD), H-bond acceptor field (A), and H-bond donor field (D), which helps to further explore the structure-activity relationship of ligand molecules.
Regression analysis of NDM-1 inhibitors in the training set was performed using partial least squares (PLS). First of all, the leave-one-out (LOO) method was adopted to cross validate, obtaining the optimal numbers of component (ONC) and the decision coefficient q 2 . Then, the non-cross validation was performed based on the ONC value, the 3D-QSAR model was eventually established, from which the correlation coefficient r 2 , estimated standard error Es, root mean square RMSE and F test value all were obtained.
The relevant parameters of PLS analysis can be used to assess model stability, as well as the potency to predict inhibitors' biological activity in the test set. Prediction correlation coefficients of the test set are calculated by the following equation: Among them, SD represents the sum of the standard deviation of experimental biological activity in test set and that in training set. PRESS is the sum of the squared error between the predicted biological activity value and the experimental ones in training set.

MD simulations
Molecular dynamics (MD) simulation is an important method to study the thermodynamic and kinetic characteristics of biomolecules at the molecular level. Two 120.5 ns comparative allatoms MD simulations were carried out for NDM-1 monomer and NDM-1-inhibitor complex, respectively. Firstly, MCPB.py [48] is used to build Zn (II) bond model, and the corresponding force field parameters were calculated with the Seminario method. The quantum mechanical calculation including geometric optimisation and Merz-Kollman charge acquisition was performed with Gaussian 09 software package [44], adopting B3LYP hybrid functional [43] and 6-31G* base group to define each atomic type.
MD simulations were implemented with AMBER 16 software package [49] using AMBER ff14SB full atomic force field. The simulation temperature was set to 300 K, and the solvent effect was characterised by TIP3P water model [50]. The minimum distance from the box boundary was set to 10.0 Å, then Na+ ions were added to neutralise system charge, finally the amounts of water molecules, total atoms, Na + ions in NDM-1 monomers and NDM-1-inhibitor complexes were 7849/ 27093/ 1 and 6131/ 21913/ 1, respectively.
Before the productive MD simulations, two-step energy minimisation was performed, i.e. the solute-constrained (the binding constant is set to 500 kcal·mol −1 ·Å −1 ) and unconstrained ones. Each-step geometry optimisation includes 5000-cycle steepest descent minimisations followed by 5000cycle conjugate gradient ones, with the convergence criterion of energy difference between two adjacent snapshots less than 0.01 kcal·mol −1 ·Å −1 . The productive MD simulations in this work were also composed of two steps, i.e. 0.5 ns solute-constrained and 120 ns unconstrained simulations. In solute-constrained stage, the constrained force constant was set to 10 kcal·mol −1 ·Å −1 ; the temperature gradually increased from 0 to 300 K. Then in the constant temperature (300 K) unconstrained MD stage, SHAKE algorithm [51] was used to constrain bond length; the truncation radius of non-bonded function was set to 10 Å; the integral step was 2 fs and the trajectories were sampled every 1ps; eventually, total 120,000 snapshots were collected and VMD was used to monitor conformational fluctuation.

Binding free energy calculation and energy decomposition
As an important physical and chemical parameter, binding free energy can be effectively used to judge molecular recognition between receptors and ligands, and is also a key criterion for evaluating the activity of drug molecules. Molecular mechanics/Poisson Boltzmann (MM/PBSA) method was used to calculate the binding free energy between NDM-1 and sulfhydryl propionic inhibitor Cmpd # 17. Based on the MD simulation of NDM-1-Cmpd # 17 complex, one snapshot was collected every 5 ns from trajectories, and a total of 24 conformations were collected, which were then optimised by 100,000 steps of energy minimisation with the convergence of energy gradient less than 0.0001 kcal·mol −1 ·Å −1 , and ultimately used for the calculation of average binding free energy.
Energy decomposition was carried out by Molecular mechanics/Generalised Born area (MM/GBSA) method [52]. According to MM/GBSA, the binding free energy values between receptors and ligands can first be divided into each residue, then be subdivided into internal energy in vacuum (including polar electrostatic interactions and non-polar Van der Waals interactions), polar solvation energy and non-polar solvation energy, and eventually be decomposed into the main chain and side chain of each residue. Here, the polar and non-polar solvation energy can be calculated with generalised Born model [53] and LCPO algorithm [54], respectively.

Conformational cluster
MMTSB package [55] was used to execute conformational cluster for the 120,000 snapshots obtained from MD simulation of NDM-1 hydrolase complexed with its inhibitor (i.e. Cmpd # 17). The conformational cluster analysis was based on the viewpoint proposed by Daura's group [56]: the RMSD value between each conformation is calculated, then the RMSD matrix (N × N, N is the conformation number) is established; A RMSD threshold is defined; if the RMSD values in the matrix are less than the threshold, the corresponding conformations are classified into one cluster, until all the conformations are grouped into a particular clusters, where the conformation with the lowest energy is defined as the representative conformation.

Free energy landscape
The free energy landscape (FEL) is used to analyse conformational distribution, free energy minimum and free energy barrier of the system. The exploration of conformational distribution probability helps to further reveal molecular recognition of receptor-ligand, protein folding and aggregation etc. The system free energy minimum is a series of steady state conformations that can be identified by endogenous ligand molecules under physiological conditions. The free energy barrier is used to represent the transient state between the various steady state conformations. In FEL analysis, the required conformational sampling is based on principal component analysis (PCA) of MD trajectories [39]. The expression for conformational free energy calculation is: Where X represents a particular principal component (PC), and P(X) is the conformational distribution probability of the system along this PC; kB and T are Boltzmann constant and absolute temperature, respectively.

Current status of NDM-1 inhibitors
In recent years, the study on novel structured antibiotics resistant to the raging superbug has been systematically carried out by scientists. Among them, the research and design strategies of NDM-1 inhibitors is a hotspot. Table 1 lists five representative NDM-1 inhibitors, respectively belonging to sulphur specie, benzopyrans, covalent binding, metal chelating agent, and other classes. Sulphur-containing inhibitors are novel NDM-1 inhibitors, and have some obvious advantages over other types of inhibitors, such as stronger inhibitory activity, small toxicity and etc. In 2007, captopril, as an angiotensin-converting enzyme inhibitor, was reported to show a strong inhibitory effect on NDM-1 [57], and the active site was formed by the chelation of thiol group with two zinc ions. Since then, sulphur-containing inhibitors have attracted more attention. Li and Shen [24,25] designed and synthesised a series of high-activity NDM-1 inhibitors by modifying proline and thiol side chains according to the inhibition mechanism of captopril.

CoMFA and CoMSIA models of sulfhydryl propionic acid
Currently, new drug development mainly includes the following four ideas: developing based on the existing antibiotics; designing from the molecular structure of known drugs; extracting from natural products; searching with computer-aided drug design (CADD) method. In fact, according to the existing inhibitors' structures and QSAR data, more active inhibitors are possibly obtained, which has become Table 1. Structures, inhibition mechanisms and characteristics of five classes of NDM-1 inhibitors.

Nos.
Class and structure pIC 50 (μM) Inhibition mechanisms Characteristics 1 Sulfhydryl propionic acid 1.5 The chelation of thiol with Zn1 hinders the hydrolysis reaction.
Diversity structure, good inhibitory activity, high probability of patent medicine modification 2 Benzopyrans 50 The negative hydrogen of benzopyrans is susceptible to a negative hydrogen transfer reaction, thus weakening its activity.
It is difficult to systematically study such compounds due to that there are few inhibitors, especially those with active data. It is not a common dominant inhibitor, and most of them are highly toxic.

others 4.0
The irreversible inhibition is achieved by depleting Zn 2+ in NDM-1 active site.
In spite of good inhibitory activity, too few inhibitor molecules lead to insufficient research.
the most important and mainstream drug development strategy.
The CADD operation relies on the structural information of drug targets and chemical information of small molecule drugs. Since NDM-1 was reported more than 10 years ago, the study on NDM-1 sulfhydryl propionic acid inhibitors has received wide attention. All available sulfhydryl propionic acid inhibitors were collected from literature reports and databases. Figure 2 shows molecular structures and pIC50 values (in μM) of 23 sulfhydryl propionic acid inhibitors used for subsequent 3D-QSAR analysis. To establish CoMFA and CoMSIA models, a training set was firstly created by randomly selecting 17 of 23 inhibitors, and then the rest 6 molecules would become elements of test set (see Figure 2). Figure 3 shows structural superimposition of sulfhydryl propionic acid inhibitors in training set. As shown from Figure 3, all the inhibitors were well superimposed, and the active sulfhydryl moiety was also located between two zinc ions, closer to the Zn1, which was basically consistent with the eutectic structure [58]. Good structural superimposition result is the basis of establishing a reasonable 3D-QSAR model. Table 2 lists the statistical parameters of CoMFA and CoM-SIA models. The cross-validation correlation coefficient (q 2 ) and non-cross-validation correlation coefficient (r 2 ) of CoMFA model are 0.589 and 0.998, respectively. RMSE is 0.120, and F-test value is 956.348. All of the above parameters show that the CoMFA model possesses good prediction ability.
The inhibitors' activity in the test set was predicted using CoMFA model, with the correlation coefficient of 0.959 between predicted and experimental values. In terms of contribution rate, stereo field (S) field and electrostatic filed (E) respectively account for 69.7% and 30.3%, showing the volume of sulfhydryl propionic acid inhibitor has an important influence on its inhibitory activity.
In the CoMSIA model, Five PLS statistics parameters (i.e. cross-validation correlation coefficient, non-cross-validation correlation coefficient, RMSE, F-test value and predicted correlation coefficient) are 0.602/ 0.993/ 0.112/ 280.205/ 0.964, also showing excellent prediction ability. In terms of each force field, the contribution rates of S, E, H, D and A are 11.6%, 24.6%, 30.6%, 20.4% and 12.8%, respectively, indicating that the hydrophobic field of sulfhydryl propionic acid inhibitors has great impact on its biological activity.
Based on the established CoMFA and CoMSIA models, Figure 4 shows the correlations between predicted pIC50 values and experimental ones for all the sulfhydryl propionic acid molecules. Excellent linear correlation (R2 is greater than 0.9) proves once again that both models have good prediction ability.
Table S1 lists all the predicted activity values for inhibitor molecules in training set and test set. The deviation (Res.) between the predicted value (Pred.) and the experimental value (pIC50) is less than 0.3 nM, and some data fit perfectly, such as Cmpd # 14. The CoMFA and CoMSIA models both show good stability and predictive ability again, which can be used to evaluate inhibitory potency of sulfhydryl propionic acid molecules. In comparison, the overall prediction ability of CoMFA model is slightly lower than that of CoMSIA model, still maintaining high prediction correlation coefficient. The combined use of two models can provide theoretical guidance for the design and modification of sulfhydryl propionic acid inhibitors.

Contour maps of CoMFA and CoMSIA models
In 3D-QSAR analyses, Cmpd # 17 with the strongest inhibitory activity was selected as the superimposition template, and the construction of contour map here was based on Cmpd # 17. Figure 5 shows the contour maps of CoMFA and CoMSIA models, where cutoff ratio of equipotential lines is 80%: 20%. As for molecular steric field distribution of CoMFA model    (see Figure 5(A)), the green and the yellow blocks around the representative inhibitor Cmpd # 17 respectively represent that the introduction of large volume groups in this region will be beneficial to and not conducive to inhibitory activity. The yellow section is mostly on the hydrophobic side and far away from sulfhydryl groups; the green section is mainly distributed on the opposite side. This explains why the activity of Cmpd # 18-23 is far lower than that of Cmpd # 17. In the molecular electrostatic field contour map of CoMFA model (see Figure 5(B)), Cmpd # 17 is surrounded by three red blocks and one blue block. Red means that the introduction of electronegative groups is conducive to increasing the inhibitory activity, and blue means the opposite. Compared with CoMFA model, CoMSIA model can provide extra hydrophobic field and H-bond donor/acceptor field information for molecular design. In Figure 5(C), almost all hydrophobic chains are surrounded by grey blocks, indicating that the introduction of hydrophobic groups at the site is not conducive to the improvement of inhibitory activity. On the contrary, a small amount of hydrophobic groups introduced on benzene ring side and methyl side (yellow block) can improve the activity of inhibitors. In the H-bond donor field (see Figure 5(D)), the blue block represents the H-bond donor region, and if H-bond donor atom-containing groups (such as hydroxyl, carboxyl, etc.) are introduced here, the inhibitor activity can be improved. Figure 5(E) shows the contour map of H-bond acceptor field. The red region indicates the favoured regions of the H-bond acceptor, which is mainly distributed on the opposite of the amide group, while the purple shows the adverse regions of H-bond acceptor. H-bond, as one of the important non-bonded interactions between drug molecules and receptors, plays a key role in molecular specific identification, binding and efficacy of drugs with therapeutic targets. According to the information of H-bond donor/acceptor fields, the introduction of more H-bond donors on the backside of sulfhydryl groups and the side of amide groups can effectively increase the catalytic activity of sulfhydryl propionic acid inhibitors. Figure 2 shows the structures of all 23 sulfhydryl propionic acids NDM-1 inhibitors, from which captopril (i.e. Cmpd # 17) is the first sulfhydryl propionic acid inhibitor to be discovered, showing relative higher NDM-1 inhibitory activity [24,25]. Figure 6 shows the binding model of Cmpd # 17 with NDM-1 via molecular docking method.

Molecular docking of Cmpd # 17
In the NDM-1-Cmpd # 17 complex model, the π-alkyl interaction between W93 benzene ring and methyl group of Cmpd # 17 was one of the important hydrophobic forces. It was consistent with previous experimental data that enzyme activity decreased by 40% ∼ 55% after W93A mutation, proving that W93 was crucial to the enzyme's catalytic function [20]. Two H-bonds of NDM-1 with K211 and N220 contribute to the specific identification between NDM-1 and Cmpd # 17. In the mechanism of NDM 1 catalysing hydrolysis cephalosporin proposed by Das [10], K211 is a key amino acid residue in proton transfer, and can stabilise the substrate in the active pocket  by salt bond [39,59]. This specific H-bond of K211 can affect the protonation step in the hydrolysis mechanism, and then inhibit NDM-1's enzymatic reaction. According to Chiou and Liang, the H-bonds of Cmpd # 17 with N220 amide can competitively inhibit the substrate identification by NDM-1 [57,59]. In addition, it was observed that the nitrogen atom and Zn2+ on the H250 imidazole ring also had some interaction with the Cmpd # 17 benzene ring respectively.

Total dynamic characteristics
In order to optimise the conformation of NDM-1 and NDM-1-Cmpd # 17, two 120.5 ns comparative MD simulations were conducted, which is also the premise for the following studies of conformational fluctuation and recognition mechanism. Figure 7 shows the total dynamics characteristics, including the change of potential energy, RMSD, gyration radius, as well as RMSF's distributions on Cα atoms.
As shown from Figure 7(A), the average potential energy and standard deviation of NDM-1-Cmpd # 17 and NDM-1 monomer were −8.035×105/ 0.02×105 and −6.406×105/ 0.01×105 kcal•mol −1 , respectively. Both systems were stable during the simul.6ation, especially NDM-1-Cmpd # 17, suggesting that the association of Cmpd # 17 helps to improve the stability of NDM-1. From RMSD analysis (see Figure 7(B)), NDM-1-Cmpd # 17 complex and NDM-1 monomer both tend to be stable after 500 ps. The average and standard deviation values of RMSD in NDM-1 complex and monomer are 0.134 ± 0.015 nm and 0.102 ± 0.011 nm, respectively. Although RMSD data occasionally fluctuated during MD simulations, they eventually became balanced and converged independently. Figure 7(C) shows the change of gyration radius over time, where the larger the value is, the more extensive the system becomes. In the whole simulation process, the gyration radius of NDM-1-Cmpd # 17 complex and NDM-1 monomer both remained stable, with the average and standard deviation values of 1.671 ± 0.011 and 1.653 ± 0.005 nm, respectively. In comparison, the gyration radius of NDM-1 complex slightly increased at 20 and 60 ns, nevertheless the system was still within the normal expansion range caused by full solvation. The stable potential energy, RMSD and gyration radius all guaranteed the reliability of the subsequent trajectory-based analyses of H-bond, energy and conformational change etc.
According to RMSF's distribution of Cα atoms (see Figure 7 (D)), seven residue fragments (i.e. S75∼L78, V89∼D90, V118∼T119, D124∼K125, T195∼V196, F205∼C208 and M248∼S249) in NDM-1-Cmpd # 17 complex have low RMSF values less than 0.7 nm. These fragments are mainly located inside of NDM-1 and are important to maintain overall system stability. For example, stable H-bonds were formed between K125 and the surrounding residues, which exhibits a certain steric hindrance effect on the residues bound to Zn1. G206∼C208 is the beginning of a flexible ring: α-helix (H228∼A238) and β-sheet (T195∼F205) junction, and low flexibility helps the two fragment maintain overall structural stability and participate in the formation of cysteine sites [60]. The association of Cmpd # 17 did not obviously affect the flexibility of the above fragment, suggesting that the inhibition mechanism of Cmpd # 17 may not be achieved by destabilising the key residues in the pocket, but by other ways.

H-bonds
H-bond is one of the most important non-bonded interactions to maintain structural stability of receptor-ligands and to ensure its proper functioning of biology. Table S2 shows the relatively stable H-bonds with occupancy over 60% in NDM-1 and NDM-1-Cmpd # 17 systems. Here, the formation of H-bond adopts the geometric criterion [45] the angle of donor-hydrogen-receptor is larger than 135°, and the distance of donor-receptor is less than 0.35 nm.
As shown from table S2, there are 10 and 13 stable H-bonds with occupancy over 60% in NDM-1-Cmpd # 17 complex and NDM-1 monomer. Although Cmpd # 17 did not directly form H-bonds with NDM-1, it greatly changed the original H-bonds network in NDM-1 monomer. For example, Five H-bonds (i.e. D254-S213, T90-H249, L160-Q146, T225-Y228 and S254-K213) were enhanced after the association with Cmpd # 17, while four different H-bonds (i.e. D267-R234, G76-S62, S212-A251 and Q76-Q37) were correspondingly weakened. Most of these enhanced H-bonds are located in the loop region connecting the internal α-helix and β-sheet, which reduces the flexibility of the loop region and stabilise the whole system conformations. In terms of the weakened H-bonds, most of them are distributed in the α-helix and β-sheet near the N-terminal and C-terminal outside NDM-Cmpd # 17 complex, far away from the binding pocket, which further indicates that the Hbond network of complex molecules is tightened in the metal catalytic centre.

Binding free energy
Binding free energy can effectively evaluate mutual recognition of receptor-ligands, and also is an important criterion to evaluate the activity of drug molecules. Table 3 lists the composition and contribution of each energy item to binding free energy between NDM-1 and Cmpd # 17. After the association of Cmpd # 17 with NDM-1, molecular internal energy significantly reduces, mainly due to a sharp decrease in molecular electrostatic energy (ELEIN, releasing 586.99 kcal•mol −1 ) and a slight increase in molecular Van der Waals internal energy (VDWIN, absorbing 3.55 kcal•mol −1 ). In addition, the solvated electrostatic energy (ELEPB) of NDM-1 rises and 558.14 kcal•mol −1 energy is absorbed, while the Van der Waals energy (VDWPB) partially reduces and 3.26 kcal•mol −1 energy is released. In conclusion, the polar electrostatic interactions are the main driving force for molecular recognition of NDM-1 by Cmpd # 17.
It can also be seen from Table 3, the changes of enthalpy, product of absolute temperature and entropy, and calculated binding free energy are −28.56/ −17.68/ −10.88 kcal•mol −1 , respectively. The results show that the recognition and association of Cmpd # 17 inhibitor with NDM-1 is spontaneous, and the NDM-1-Cmpd # 17 complex is sufficiently stable. What's more, the calculated binding free energy agree well with the experimental data (ΔGexp) of −7.30 kcal•mol −1 .

Energy decomposition
In the design of sulfhydryl propionic acid derivatives, it is necessary to determine the key residues of NDM-1 recognised by these kinds of inhibitors, so binding free energy decomposition was performed for the association of NDM-1 with Cmpd # 17 at residual level. Table 4 lists the key residues favouring the binding of NDM-1 to Cmpd # 17.
As shown from Table 4, the polar solvation (EGB) is the main factor to promote molecular association. The top 5 residues (i.e. D124, H250, H120, H189 and H122), which contribute the most to molecular recognition, are related to the composition of metal active site. In addition, the energy contribution of W93 outside of metal active site is also greater. Referring to the previous molecular docking (see Figure 6), W93's benzene ring can form a π-alkyl non-polar interaction with Cmpd # 17, and W93 is one of the key residues to maintain catalytic activity and molecular stability [20]. Similarly, I35 is also important for the stability of NDM-1 overall structure, molecular recognition of substrates, as well as enzymatic catalytic potency [30,61].
Based on the above energy decomposition results, fig. S2A shows the position of the top 7 key residues in metal active  centre. Centreing on Cmpd # 17, fig. S2B shows all the contact residues (13 in total) within a radius of 0.4 nm. Interestingly, the top 7 key residues obtained from energy decomposition were all within these 13 contact residues. Figure 8 shows FEL (A, C) and conformational clustering (B, D) of NDM-1 monomer and NDM-1-Cmpd # 17 complex at 300 K. According to Figure 8(A), there are four independent low-free energy regions (i.e. M1, M2, M3 and M4) in NDM-1-Cmpd # 17 complex, with two on each side and two in the middle. The four trajectories with lower free energy are mainly distributed in the time periods of 0∼25, 95∼120, 85∼110 and 25∼85 ns. As for NDM-1 monomer, there are only three lower free energy regions (i.e. M1, M2 and M3) located in the upper left and right corner (see Figure 8(c)). It is worth noting that the PC1 and PC2 fluctuation range of NDM-1-Cmpd # 17 complex is significantly larger than that of NDM-1 monomer, indicating that the binding of Cmpd # 17 to NDM-1 increases NDM-1's overall mobility, and makes the stable conformation more loose and disordered. The conformational flexibility results are consistent with that obtained from the previous RMSD analysis (see Figure 7). In order to further analyse the conformation with the lowest free energy at different simulation time, comparative conformational cluster analyses of NDM-1-Cmpd # 17 and NDM-1 were respectively performed, where the RMSD threshold was set to 0.03 nm (see Figure 8(B,D)). There are 3 and 4 clusters respectively in the two systems, and each of them represents a relatively independent representative conformation, which is consistent with FEL analysis. Specially, there were 4 clusters in NDM-1-Cmpd # 17, with their corresponding conformations respectively accounting for 48.1%, 7.4%, 10.2% and 34.3%. Nevertheless, NDM-1 has only 3 clusters and the conformational transition space is relatively small.

Possible recognition and inhibition mechanism of
Cmpd # 17 Figure 9 shows recognition and inhibition mechanisms of NDM-1 by Cmpd # 17. Based on the analyses of H-bond, free energy calculation, energy decomposition, the following molecular recognition process can be inferred: Cmpd # 17 firstly enters the active pocket through forming H-bonds between the phenolic hydroxyl group、amide oxygen of Cmpd # 17 and K211、 N220 of NDM-1; Then driven by electrostatic force, Cmpd # 17 tightly binds to I35, W93, H120, H122, D124, H189 and H250; Besides, the sulfhydryl is coordinated with Zn1, causing the NDM-1 H-bond network to redistribute to metal catalytic region and a reduction in binding free energy.
In addition, based on the FEL and conformational cluster analyses, it can be speculated that the association of Cmpd #17 may cause NDM-1's conformation to be disordered, leading to partial loss of catalytic activity, which may be one of the inhibition mechanisms of sulfhydryl propionic acid derivatives. In other words, the binding of Cmpd # 17 to the metal-active catalytic site of NDM-1weakened the recognition of NDM-1 by its endogenic substrates and competitively reduces penicillin interaction with NDM-1, which causes the penicillin to be free and effective.

Conclusions
After a theoretical exploration of 3D-QSAR for sulfhydryl propionic acid NDM-1 inhibitors, CoMFA and CoMSIA models with good predictive ability were established. According to thiol Figure 8. (Colour online) Free energy landscapes of NDM-1-Cmpd # 17 complex (A) and NDM-1 monomer (C) respectively, as well as the corresponding conformation cluster analyses (B and D). The colour depth in A/C represents conformational free energy level; the darker the colour is, the lower the free energy value of the corresponding conformation is.
inhibitors' contour maps, the introduction of a hydrophobic group at the α-methyl site, and a H-bond donor at the tail of the hydrophobic side chain both are helpful for increasing enzymatic activity, partly explaining the structure-activity relationship from biological experiences. To analyse molecular recognition of NDM-1 and sulfhydryl propionic acid NDM-1 inhibitors at atomic level, the binding modes between NDM-1 and Cmpd # 17 was obtained with molecular docking method, then refined with long-time MD simulations. The key residues and binding driving forces both were investigated by trajectory-based Hbond, binding free energy, energy decomposition analyses. The results show that I35, W93, H120, H122, D124, H189 and H250 are the key residues in molecular recognition of NDM-1 by Cmpd # 17, and the polar interactions are the main driving force. The association of thiol inhibitor Cmpd # 17 will lead to the redistribution of NDM-1 H-bond network, thus stabilising the inner metal catalytic centre. Free energy landscape and conformational cluster analyses both show that the magnitude of NDM-1 conformational motion become disorder after recognised by inhibitors, which may partially impair enzymatic catalytic function. The work can provide a theoretical basis for the design and modification of specific thiol inhibitors based on ligand structural characteristics and receptor-ligand recognition mechanism.

Disclosure statement
No potential conflict of interest was reported by the authors.