3D-QSAR, Molecular Docking and Molecular Dynamics Simulation of Pseudomonas aeruginosa LpxC Inhibitors

As an important target for the development of novel antibiotics, UDP-3-O-(R-3-hydroxymyristoyl)-N-acetylglucosamine deacetylase (LpxC) has been widely studied. Pyridone methylsulfone hydroxamate (PMH) compounds can effectively inhibit the catalytic activity of LpxC. In this work, the three-dimensional quantitative structure-activity relationships of PMH inhibitors were explored and models with good predictive ability were established using comparative molecular field analysis and comparative molecular similarity index analysis methods. The effect of PMH inhibitors’ electrostatic potential on the inhibitory ability of Pseudomonas aeruginosa LpxC (PaLpxC) is revealed at the molecular level via molecular electrostatic potential analyses. Then, two molecular dynamics simulations for the PaLpxC and PaLpxC-inhibitor systems were also performed respectively to investigate the key residues of PaLpxC hydrolase binding to water molecules. The results indicate that orderly alternative water molecules can form stable hydrogen bonds with M62, E77, T190, and H264 in the catalytic center, and tetracoordinate to zinc ion along with H78, H237, and D241. It was found that the conformational transition space of PaLpxC changes after association with PMH inhibitors through free energy landscape and cluster analyses. Finally, a possible inhibitory mechanism of PMH inhibitors was proposed, based on our molecular simulation. This paper will provide a theoretical basis for the molecular design of LpxC-targeting antibiotics.


Introduction
Bacterial drug-resistance has become a public health problem worldwide and a serious threat to human health [1,2]. Recently, bacteria with the multiple drug-resistant (MDR) mechanism, also called "superbugs", have been reported all over the world [3][4][5], such as pan drug-resistant Acinectobacter baumannii, carbapenem-resistant Klebsiella pneumoniae, and Pseudomonas aeruginosa. Most of these "superbugs" are Gram-negative (G − ) bacteria possessing complex cell membranes, were performed comparatively. The difference of the motion patterns between PaLpxC and its complex with inhibitors were investigated using conformational cluster and free energy landscape (FEL) analyses (see Figure 1). These studies will provide a theoretical basis for the activity prediction, molecular design, and modification of PMH LpxC inhibitors. between PaLpxC and its complex with inhibitors were investigated using conformational cluster and free energy landscape (FEL) analyses (see Figure 1). These studies will provide a theoretical basis for the activity prediction, molecular design, and modification of PMH LpxC inhibitors.

Systems for Simulation
PMH LpxC inhibitors belong to a group of traditional hydroxamate molecules, which mainly suppress the activity of zinc ions at the bottom of LpxC's active pocket relying on the hydrophilic terminal hydroxamate moiety [10,[14][15][16][17][18][19][20][21][22]24]. Figure 2 shows the binding mode of Cmpd # 290 with PaLpxC and the molecular alignment of the PMH LpxC inhibitors. It is worth mentioning that the binding details will be analyzed below (see section on "molecular docking"). As shown in Figure 2, the public substructure of PMH molecules (i.e., pyridone methylsulfone hydroxamate) is aligned well, which maximizes the similarity with the spatial orientation of the molecules, and provides a good foundation for the subsequent generation of the comparative molecular field analysis (CoMFA) and comparative molecular similarity index analysis (CoMSIA) models.

CoMFA and CoMSIA Models
In this work, 31 PMH LpxC inhibitors (training set) were used for the establishment of the 3D-QSAR models, with the related parameters and results shown in Table S1. In the CoMFA model, the cross-validated correlation coefficient (q 2 ) and noncross-validated correlation coefficient (r 2 ) are 0.622 and 0.978, respectively. The root mean squared error (RMSE) is 0.102, and the combination with the high predicted correlation coefficient (r 2 p = 0.933) confirms the reasonability and reliability of this model. According to the CoMFA model, the contribution of the steric field (S) is 67.7%, and the electrostatic field (E) is 32.3%. The model indicates that the steric field surrounding the PMH LpxC inhibitors plays an important role in its inhibitory activity. The CoMSIA model also analyzes the hydrophobic field (H), hydrogen bond (H-bond) acceptor field (A), and H-bond donor field (D) of the training set molecules beyond the steric field and electrostatic field. In light of the CoMSIA model, the contribution of S is 35.3%, while that of E is 22.1%. Moreover, the hydrophobic filed portion occupies 30.0%, and the H-bond donor field and acceptor field hold 11.5% and 1.1%, Figure 1. Protocol of this work. 3D-QSAR: three-dimensional quantitative structure-activity relationships; CoMFA: comparative molecular field analysis; CoMSIA: comparative molecular similarity index analysis; MD: molecular dynamics.

Systems for Simulation
PMH LpxC inhibitors belong to a group of traditional hydroxamate molecules, which mainly suppress the activity of zinc ions at the bottom of LpxC's active pocket relying on the hydrophilic terminal hydroxamate moiety [10,[14][15][16][17][18][19][20][21][22]24]. Figure 2 shows the binding mode of Cmpd # 290 with PaLpxC and the molecular alignment of the PMH LpxC inhibitors. It is worth mentioning that the binding details will be analyzed below (see section on "molecular docking"). As shown in Figure 2, the public substructure of PMH molecules (i.e., pyridone methylsulfone hydroxamate) is aligned well, which maximizes the similarity with the spatial orientation of the molecules, and provides a good foundation for the subsequent generation of the comparative molecular field analysis (CoMFA) and comparative molecular similarity index analysis (CoMSIA) models.

CoMFA and CoMSIA Models
In this work, 31 PMH LpxC inhibitors (training set) were used for the establishment of the 3D-QSAR models, with the related parameters and results shown in Table S1. In the CoMFA model, the cross-validated correlation coefficient (q 2 ) and noncross-validated correlation coefficient (r 2 ) are 0.622 and 0.978, respectively. The root mean squared error (RMSE) is 0.102, and the combination with the high predicted correlation coefficient (r 2 p = 0.933) confirms the reasonability and reliability of this model. According to the CoMFA model, the contribution of the steric field (S) is 67.7%, and the electrostatic field (E) is 32.3%. The model indicates that the steric field surrounding the PMH LpxC inhibitors plays an important role in its inhibitory activity. The CoMSIA model also analyzes the hydrophobic field (H), hydrogen bond (H-bond) acceptor field (A), and H-bond donor field (D) of the training set molecules beyond the steric field and electrostatic field. In light of the CoMSIA model, the contribution of S is 35.3%, while that of E is 22.1%. Moreover, the hydrophobic filed portion occupies 30.0%, and the H-bond donor field and acceptor field hold 11.5% and 1.1%, respectively. The steric and hydrophobic fields of PMH LpxC inhibitors were shown to contribute greatly to their biological activities, followed by the electrostatic field and H-bond field. Based on the results of the CoMFA and CoMSIA models, it is speculated that changing the bulk and hydrophobicity of the molecules may be an important method to improve the biological activity of PMH LpxC inhibitors. respectively. The steric and hydrophobic fields of PMH LpxC inhibitors were shown to contribute greatly to their biological activities, followed by the electrostatic field and H-bond field. Based on the results of the CoMFA and CoMSIA models, it is speculated that changing the bulk and hydrophobicity of the molecules may be an important method to improve the biological activity of PMH LpxC inhibitors.   Figure 3, there is a significant linear correlation between the predicted pIC50 and the experimental values, which proves the reliability of the two models. Table S2 lists the predicted activity of molecules in the test set based on the established CoMFA model and CoMSIA model. The deviation of negative logarithm form of IC50 (pIC50) between the predicted and experimental values is less than 1, and some individual results are very consistent with the experimental values, such as Cmpd # 322. The results show that the CoMFA model established from this training set has good stability and predictive ability for the inhibitory activity of PMH LpxC molecules. By comparing the prediction results of the CoMSIA model with the CoMFA model, it is found that the overall predictive ability of the CoMSIA model is slightly inferior to that of the CoMFA model. Nevertheless, the CoMSIA model still has good predictive ability and high stability, and can complement the CoMFA model. The former provides more comprehensive   Figure 3, there is a significant linear correlation between the predicted pIC 50 and the experimental values, which proves the reliability of the two models. information for drug design based on 3D-QSAR, as well as for the modification of the lead compounds.  Figure S1 shows the contour maps of the CoMFA model for Cmpd # 326 with the highest activity; the truncation ratio of the contour line is 80%:20%. It can be seen from the steric field contour map (see Figure S1A) around the PMH LpxC inhibitor that there are a total of five yellow and one green contours. The yellow contours are mainly distributed in the hydrophobic side of the molecule and only one block is located in the methylsulfone terminus, suggesting that the introduction of large volume groups is not conducive to the increase of inhibitory activity. For example, the inhibitory activity of Cmpd # 292 with the introduction of methyl in the para-position of the pyridone ring is significantly lower than that of non-substituted Cmpd # 290. The green blocks at the end of the hydrophobic chain suggest that introducing large volume groups in that region will increase the inhibitory activity of the PMH LpxC inhibitors, which explains the higher inhibitory activity of Cmpds # 328-334 relative to Cmpd # 290 (see Table S2). However, in the electrostatic field contour map (see Figure S1B), there are a total of four red and one blue contours at the end of the phenyl ring, which indicate that an increase of the positive and negative charges in these regions are beneficial to improving inhibitory activity. Taking Cmpd # 326 and # 290 as an example, the introduction of two F atoms and one Cl atom with strong electron-withdrawing effects into the phenyl ring notably increases the electron cloud density. Therefore, the inhibitory activity of Cmpd # 326 is significantly higher than that of Cmpd # 290 (see Table S2). In addition, at the methylsulfone terminus, there is also a red contour, suggesting that increasing the number of negative charges of the methylsulfone moiety may enhance the inhibitory activity of the PMH LpxC inhibitors.

CoMFA and CoMSIA Contour Maps
The CoMSIA model can show five S, H, E, A, and D different molecular field results, which can provide more information for molecular design than the CoMFA model. Figure S2A-C, respectively, shows the hydrophobic field, H-bond acceptor field, and H-bond donor field contour maps of the CoMSIA model. The reference molecule in the drawing is Cmpd # 326. In Figure S2A, we can see that almost all of the hydrophobic chains are wrapped by grey contours, which indicates that the introduction of hydrophobic groups will not be conducive to the increase of molecular biological activity. In Figure S2B, cyan and purple contour maps represent the favored and disfavored regions of the H-bond donor, respectively. In Figure S2C, the red contour maps represent the disfavored regions of the H-bond acceptor. The H-bond is one of the important non-bonding interactions between drug molecules and receptors. According to these contour maps, the introduction of more H-bond donors at the oxime and the methylsulfone groups is unexpectedly not conducive to  Figure S1 shows the contour maps of the CoMFA model for Cmpd # 326 with the highest activity; the truncation ratio of the contour line is 80%:20%. It can be seen from the steric field contour map (see Figure S1A) around the PMH LpxC inhibitor that there are a total of five yellow and one green contours. The yellow contours are mainly distributed in the hydrophobic side of the molecule and only one block is located in the methylsulfone terminus, suggesting that the introduction of large volume groups is not conducive to the increase of inhibitory activity. For example, the inhibitory activity of Cmpd # 292 with the introduction of methyl in the para-position of the pyridone ring is significantly lower than that of non-substituted Cmpd # 290. The green blocks at the end of the hydrophobic chain suggest that introducing large volume groups in that region will increase the inhibitory activity of the PMH LpxC inhibitors, which explains the higher inhibitory activity of Cmpds # 328-334 relative to Cmpd # 290 (see Table S2). However, in the electrostatic field contour map (see Figure S1B), there are a total of four red and one blue contours at the end of the phenyl ring, which indicate that an increase of the positive and negative charges in these regions are beneficial to improving inhibitory activity. Taking Cmpd # 326 and # 290 as an example, the introduction of two F atoms and one Cl atom with strong electron-withdrawing effects into the phenyl ring notably increases the electron cloud density. Therefore, the inhibitory activity of Cmpd # 326 is significantly higher than that of Cmpd # 290 (see Table S2). In addition, at the methylsulfone terminus, there is also a red contour, suggesting that increasing the number of negative charges of the methylsulfone moiety may enhance the inhibitory activity of the PMH LpxC inhibitors.

CoMFA and CoMSIA Contour Maps
The CoMSIA model can show five S, H, E, A, and D different molecular field results, which can provide more information for molecular design than the CoMFA model. Figure S2A-C, respectively, shows the hydrophobic field, H-bond acceptor field, and H-bond donor field contour maps of the CoMSIA model. The reference molecule in the drawing is Cmpd # 326. In Figure S2A, we can see that almost all of the hydrophobic chains are wrapped by grey contours, which indicates that the introduction of hydrophobic groups will not be conducive to the increase of molecular biological activity. In Figure S2B, cyan and purple contour maps represent the favored and disfavored regions of the H-bond donor, respectively. In Figure S2C, the red contour maps represent the disfavored regions of the H-bond acceptor. The H-bond is one of the important non-bonding interactions between drug molecules and receptors. According to these contour maps, the introduction of more H-bond donors at the oxime and the methylsulfone groups is unexpectedly not conducive to enhancing the activity of PMH LpxC inhibitors. Similarly, the introduction of more H-bonds acceptors at the methylsulfone substituent is also not conducive to the improvement of the activity of these inhibitors theoretically.
In general, the CoMFA and CoMSIA models established in this work have good predictive ability. Additionally, the CoMFA and CoMSIA contour maps both also visually exhibit the relationship between the molecular field properties (including substituent groups, charge distribution, hydrophobicity, number of H-bond donor and acceptor atoms) and inhibitory activity of PMH LpxC inhibitors. On the basis of the contour maps, the introduction of bulky groups at the para-position of the carbonyl on the pyridone ring will enhance the negative charge of the phenyl terminus and also hinder the introduction of highly hydrophobic groups. Similarly, the avoidance of substituents containing H-bond donor or acceptor atoms near the end of the hydroxamate moiety may improve the inhibitory activity of PMH LpxC inhibitors.

MESP and Mulliken Charge Analyses
The strength of the interaction between the inhibitor molecule and the specific target is usually dependent on the unique structure and charge properties of the inhibitor molecule itself. Through calculation and analysis of the charge properties of inhibitors, it is helpful to investigate the relationship between the charge distribution and the inhibitory activity. In this work, Cmpd # 318, Cmpd # 290, and Cmpd # 326, which have similar molecular weights but different inhibitory activities, were chosen for molecular electrostatic potential (MESP) and Mulliken charge analyses. The results of the Mulliken charge calculation are listed in Table S3. Figure 4 shows the MESP distribution maps of Cmpd # 290 (A), Cmpd # 318 (B), and Cmpd # 326 (C). The red part indicates the region of the most negative charge, which shows the strongest electronic donation region. The blue part represents the region of the most positive charge, meaning the strongest electronic acceptance region. As can be seen from Figure 4, the negative charge region of Cmpd # 290 is mainly concentrated in the carbonyl oxygen of the pyridone ring and the oxygen of methylsulfone, suggesting these areas are H-bond acceptor regions; however, the positive charge region is mainly concentrated in the vicinity of the hydroximate moiety, which can coordinate with the lone electron pair of Zn 2+ , and inhibit its catalytic activity. Cmpd # 318 is obtained by introducing a methoxy group into the 2-position of the phenyl moiety in Cmpd # 290 (see Figure 4B). The resonance donating effect of methoxy increases the negative charge at the para-position of the conjugated system; while the negative charge at sulfone and the positive charge of the hydroxamate group decreases. Sulfone and hydroxamate both are the key groups of PMH LpxC inhibitors, and the change of the charge properties might be an important factor for the lower inhibitory activity of Cmpd # 318 compared to that of Cmpd # 290. Moreover, the hydrogen atoms at the 2,3-positon on the phenyl group of Cmpd # 290 are replaced by fluorine atoms. With the chlorine atom introduced into the 4-position of the phenyl group, Cmpd # 326 is obtained. The MESP distribution of Cmpd # 326 is shown in Figure 4C. Comparing it to Figure 4A, the halogen atoms have little effect on the charge properties of the above three sites, resulting in only a slight decrease. However, the polarization of the electrons in the molecular conjugate system changes obviously, which is beneficial to the enhancement of inhibitory activity. Thus, it is speculated that the increase of the inhibitory activity of Cmpd # 326 is related to the electrostatic potential change of the hydrophobic terminus.

Molecular Docking
Molecular docking was performed with the AutoDock4Zn force field, and the results are shown in Figure 5. As can be seen in Figure 5A, the 38 PMH LpxC inhibitors are well bound in the active pocket. The binding mode of the docking result is consistent with that of the co-crystal structure (3UHM.pdb): the hydroxamate moiety coordinates with Zn 2+ , and the phenyl pyridone chain is well combined in the hydrophobic channel of PaLpxC. Figure 5B shows the interactions between the inhibitor Cmpd # 290 and PaLpxC in the co-crystal structure. The purple dashed lines in the figure represent the hydrophobic interaction, which contains an amide-π packing interaction between the pyridone ring and G192, π-alkyl mixed hydrophobic interactions between F191 and the methyl of Cmpd # 290, and the π-σ mixed hydrophobic interactions formed between phenyl pyridone and L18, A206, A214, and V216. The green dashed lines indicate the H-bonds, which mainly exist between the sulfone oxygen atoms and the ω-amino group of K238, and between the β-hydroxyl of T190 and the carbonyl oxygen of the hydroxamate moiety. In addition, it is also observed that there is a certain electrostatic interaction between Zn 2+ and the carboxyl group of E77 (red dashed lines).  Red/blue color respectively indicates the highly electronegative/electropositive regions, while the regions with the other colors represent the potential between the two strongly electronegative and electropositive regions.

Molecular Docking
Molecular docking was performed with the AutoDock4 Zn force field, and the results are shown in Figure 5. As can be seen in Figure 5A, the 38 PMH LpxC inhibitors are well bound in the active pocket. The binding mode of the docking result is consistent with that of the co-crystal structure (3UHM.pdb): the hydroxamate moiety coordinates with Zn 2+ , and the phenyl pyridone chain is well combined in the hydrophobic channel of PaLpxC. Figure 5B shows the interactions between the inhibitor Cmpd # 290 and PaLpxC in the co-crystal structure. The purple dashed lines in the figure represent the hydrophobic interaction, which contains an amide-π packing interaction between the pyridone ring and G192, π-alkyl mixed hydrophobic interactions between F191 and the methyl of Cmpd # 290, and the π-σ mixed hydrophobic interactions formed between phenyl pyridone and L18, A206, A214, and V216. The green dashed lines indicate the H-bonds, which mainly exist between the sulfone oxygen atoms and the ω-amino group of K238, and between the β-hydroxyl of T190 and the carbonyl oxygen of the hydroxamate moiety. In addition, it is also observed that there is a certain electrostatic interaction between Zn 2+ and the carboxyl group of E77 (red dashed lines).

Molecular Docking
Molecular docking was performed with the AutoDock4Zn force field, and the results are shown in Figure 5. As can be seen in Figure 5A, the 38 PMH LpxC inhibitors are well bound in the active pocket. The binding mode of the docking result is consistent with that of the co-crystal structure (3UHM.pdb): the hydroxamate moiety coordinates with Zn 2+ , and the phenyl pyridone chain is well combined in the hydrophobic channel of PaLpxC. Figure 5B shows the interactions between the inhibitor Cmpd # 290 and PaLpxC in the co-crystal structure. The purple dashed lines in the figure represent the hydrophobic interaction, which contains an amide-π packing interaction between the pyridone ring and G192, π-alkyl mixed hydrophobic interactions between F191 and the methyl of Cmpd # 290, and the π-σ mixed hydrophobic interactions formed between phenyl pyridone and L18, A206, A214, and V216. The green dashed lines indicate the H-bonds, which mainly exist between the sulfone oxygen atoms and the ω-amino group of K238, and between the β-hydroxyl of T190 and the carbonyl oxygen of the hydroxamate moiety. In addition, it is also observed that there is a certain electrostatic interaction between Zn 2+ and the carboxyl group of E77 (red dashed lines).   Figure S3 shows the correlation between the calculated average binding free energies and the experimental pIC 50 values of the PMH LpxC inhibitors. As can be seen from Figure S3, the correlation coefficient (r 2 ) between them is 0.312. To some extent, the molecular docking method based on the AutoDock4 Zn force field can be used to predict the biological activity of PMH LpxC inhibitors roughly.

Stability and Flexibility of Simulation Systems
Throughout the MD simulations, the potential energy of the PaLpxC and PaLpxC-Cmpd # 290 systems are nearly kept stable at −1.16 × 10 5 kcal/mol, and the fluctuation rate is 0.172% (see Figure 6A). The root mean squared deviation (RMSD) analysis results of the systematic C α atoms (see Figure 6B) shows that the RMSD of the PaLpxC system in the process of simulation is basically maintained at 1.124 ± 0.081 Å; whereas that of the complex system formed by PaLpxC and Cmpd # 290 is basically maintained at 1.135 ± 0.090 Å. For systems with more than 4764 atoms, low potential energy and RMSD values indicate that the MD simulation process is stable and reliable.  Figure S3 shows the correlation between the calculated average binding free energies and the experimental pIC50 values of the PMH LpxC inhibitors. As can be seen from Figure S3, the correlation coefficient (r 2 ) between them is 0.312. To some extent, the molecular docking method based on the AutoDock4Zn force field can be used to predict the biological activity of PMH LpxC inhibitors roughly.

Stability and Flexibility of Simulation Systems
Throughout the MD simulations, the potential energy of the PaLpxC and PaLpxC-Cmpd # 290 systems are nearly kept stable at −1.16 × 10 5 kcal/mol, and the fluctuation rate is 0.172% (see Figure  6A). The root mean squared deviation (RMSD) analysis results of the systematic Cα atoms (see Figure 6B) shows that the RMSD of the PaLpxC system in the process of simulation is basically maintained at 1.124 ± 0.081 Å; whereas that of the complex system formed by PaLpxC and Cmpd # 290 is basically maintained at 1.135 ± 0.090 Å. For systems with more than 4764 atoms, low potential energy and RMSD values indicate that the MD simulation process is stable and reliable. The root mean squared fluctuation (RMSF) of systematic Cα atoms is exhibited in Figure 6C. The RMSF value of Cα atoms indicates the flexibility of the corresponding amino acid residues. According to Figure 6C, there are three amino acid fragments of the PaLpxC molecule with RMSF values less than 0.5 nm, that is, H78-L86, I239-L248, and G263-A265. These three fragments are located near the catalytic Zn 2+ , and the lower flexibility is beneficial to the specific recognition of the substrate. Generally, the formation of the complex did not affect the RMSF value of the PaLpxC molecules, which indicates that the binding of the PMH LpxC inhibitor Cmpd # 290 produced no The root mean squared fluctuation (RMSF) of systematic C α atoms is exhibited in Figure 6C. The RMSF value of C α atoms indicates the flexibility of the corresponding amino acid residues. According to Figure 6C, there are three amino acid fragments of the PaLpxC molecule with RMSF values less than 0.5 nm, that is, H78-L86, I239-L248, and G263-A265. These three fragments are located near the catalytic Zn 2+ , and the lower flexibility is beneficial to the specific recognition of the substrate. Generally, the formation of the complex did not affect the RMSF value of the PaLpxC molecules, which indicates that the binding of the PMH LpxC inhibitor Cmpd # 290 produced no significant change in the flexibility of the PaLpxC molecules, and it may be because of the inherent low flexibility of the PaLpxC molecule itself.
The B-factor can be converted to RMSF by Equation (1): According to the correlation analysis of the B-factor and RMSF in Figure 6D, with a significant correlation and a correlation coefficient of fitting curve r 2 = 0.30, the reliability of the MD simulation results is proven again.

Molecular Recognition Mechanism
The H-bond is one of the most important non-bond interactions that maintains the stability of complex molecules, playing a key role in the process of molecular recognition. In order to further explore the molecular recognition mechanism between Cmpd # 290 and PaLpxC, Table S4 lists the H-bonds formed between Cmpd # 290 and PaLpxC. The H-bond is defined by geometry criterion [26]: the distance between the donor atom (D) and acceptor atom (A) is less than 3.5 Å, and the angle D-H-A is more than 135 • .
As can be seen from Table S4, the stable H-bonds are formed between Cmpd # 290 and M62, T190, and H264 of PaLpxC. Miller's group [27] compared the co-crystal complex structure of Aquifex aeolicus LpxC (AaLpxC)-LPC-009 [28] with that of PaLpxC-LPC-009, proposing that the position of the carbonyl group in PaLpxC M62 is similar to that of AaLpxC H58, which is one of the key residues of the hydrolysis reaction. In the density functional theory (DFT) study [29] on the LpxC deacetylation mechanism, the formation of the H-bond between the hydroxyl at the side chain of T190 (residue number in PaLpxC) and the carbonyl oxygen atom of the substrate help to restrict it near the active site. In a series of studies of the AaLpxC hydrolysis mechanism [14,30,31], H264 (residue number in PaLpxC) is an important residue providing a proton for the amino group in hydrolysate. Thus, stable H-bonds between Cmpd # 290 and M62, T190, and H264 of PaLpxC play important roles in maintaining the stability of the complex, as well as the proton transfer process of the hydrolysis. It can also be seen from Table S4 that there is an H-bond formed between Cmpd # 290 and K238. Due to the rotation of the amino group, the H-bond between Cmpd # 290 and K238 is not as stable as the other H-bonds discussed above, with the total percentage of frames that showed the H-bond being only around 22.78%. Given that K238 is located in the uridine diphosphate (UDP) binding pocket [32,33], it is speculated that the H-bonds between Cmpd # 290 and K238 may be responsible for the high inhibitory capacity of Cmpd # 290.
Binding free energy is a dominant criterion in evaluating the activity of drug molecules. Molecular mechanics/Poisson-Boltzmann (MM/PBSA) algorithm was used to calculate the binding free energy. One snapshot was collected from the stable trajectory (20-100 ns) every 5 ns, and a total of 17 systemic conformations were collected from the MD trajectory for energy minimization. Each conformation was minimized to 100,000 steps with convergence of the energy gradient less than 0.0001 kcal·mol −1 ·Å −2 in entropy calculations. Table S5 lists the contributions to the binding free energy between PaLpxC and Cmpd # 290. In the form of each kind of energy item from Table S5, the free energy of PaLpxC decreased 13.79 kcal/mol after the formation of the complex with Cmpd # 290. The decrease of the systemic free energy is mainly due to two parts: the systemic van der Waals energy (VDW IN ) and the hydrophobic part of the systemic solvation energy (VDW PB ). From this point of view, the nonpolar interaction is the main driving force promoting the association of PaLpxC with the Cmpd # 290 inhibitor.
The energy decomposition results (see Figure 7 and Table S6) show that D241, F191, K238, M62, I197, T190, and H264 in the PaLpxC contribute greatly to the binding free energy. Combined with the analyses of the H-bond, it can be determined that M62, T190, K238, and H264 are the key residues in the molecular recognition between Cmpd # 290 and PaLpxC. According to the co-crystal complex structure of PaLpxC-Cmpd # 290 (see Figure 5B), D241 is also the key residue for fixing Zn 2+ and maintaining its pentacoordinate configuration. Moreover, the phenyl of F191 also formed the π-alkyl mixed hydrophobic interaction with the methyl of Cmpd # 290. The detailed binding mode of Cmpd # 290 in the PaLpxC active pocket with the lowest energy conformation is shown in Figure 7. of Cmpd # 290 in the PaLpxC active pocket with the lowest energy conformation is shown in Figure  7.

Conformational Change of the Systems
FEL is a kind of probability distribution function of free energy obtained by principal component analysis and statistics based on MD equilibrium trajectories. Figure 8A,C shows the free energy distribution of PaLpxC-Cmpd # 290 and the PaLpxC systems at 300 K, respectively. The depth of the color in the figure represents the free energy level of the corresponding regions. According to Figure 8A, there are two lower free energy independent regions in the complex PaLpxC-Cmpd # 290, with M1 at the upper left and M2 at the upper right. The minimum free energy region mainly corresponds to 20-35 ns, 50-65 ns and 75-120 ns. There are also two independent lower free energy regions in the PaLpxC protein corresponding to M1 at the lower right and M2 at the upper left during 65-120 ns. Based on the distribution of the two reaction coordinates (i.e., PC1 and PC2), it is found that, after binding with the inhibitor Cmpd # 290, the low free energy region of the PaLpxC does not significantly increase or decrease, but the conformational sampling space reduces. Specifically, PC1 reduces from the original 0.70 to −0.70 to 0 to −0.60, indicating that the complex system is becoming more stable, which is consistent with the calculation results of free energy above ( Figure S4).
In order to further analyze the local minimum area of the free energy distribution corresponding to simulation time, the conformational cluster analysis is performed with a cutoff of 1.1 Å for the PaLpxC-Cmpd # 290 and PaLpxC systems. The result is shown in Figure 8B,D. From Figure 8B,D, the number of clusters of the two systems is both two, and each cluster represents a relatively independent conformational mode. Statistical analysis shows that the distribution ratio of cluster 1 to cluster 2 of the PaLpxC protein is 25.2%:74.8%. After binding with Cmpd # 290, the corresponding distribution ratio changes to 21.1%:78.9%. The results indicate that the conformational transition space decreases in the complex system, and the association with Cmpd # 290 makes the system more stable. Figure S5 shows the slow-motion modes of PaLpxC (A) and PaLpxC-Cmpd # 290 (B). It can be seen from Figure S5 that after binding with the inhibitor Cmpd # 290, the local movement of the PaLpxC protein produce obvious changes. For example, the upper end of the α helices located in R1 (G263-A276) has a tendency to approach the inhibitor directly, as well as the β strand in the R2 (V137-A144) region. In addition, the rotational motion of the whole

Conformational Change of the Systems
FEL is a kind of probability distribution function of free energy obtained by principal component analysis and statistics based on MD equilibrium trajectories. Figure 8A,C shows the free energy distribution of PaLpxC-Cmpd # 290 and the PaLpxC systems at 300 K, respectively. The depth of the color in the figure represents the free energy level of the corresponding regions. According to Figure 8A, there are two lower free energy independent regions in the complex PaLpxC-Cmpd # 290, with M1 at the upper left and M2 at the upper right. The minimum free energy region mainly corresponds to 20-35 ns, 50-65 ns and 75-120 ns. There are also two independent lower free energy regions in the PaLpxC protein corresponding to M1 at the lower right and M2 at the upper left during 65-120 ns. Based on the distribution of the two reaction coordinates (i.e., PC1 and PC2), it is found that, after binding with the inhibitor Cmpd # 290, the low free energy region of the PaLpxC does not significantly increase or decrease, but the conformational sampling space reduces. Specifically, PC1 reduces from the original 0.70 to −0.70 to 0 to −0.60, indicating that the complex system is becoming more stable, which is consistent with the calculation results of free energy above ( Figure S4).
In order to further analyze the local minimum area of the free energy distribution corresponding to simulation time, the conformational cluster analysis is performed with a cutoff of 1.1 Å for the PaLpxC-Cmpd # 290 and PaLpxC systems. The result is shown in Figure 8B,D. From Figure 8B,D, the number of clusters of the two systems is both two, and each cluster represents a relatively independent conformational mode. Statistical analysis shows that the distribution ratio of cluster 1 to cluster 2 of the PaLpxC protein is 25.2%:74.8%. After binding with Cmpd # 290, the corresponding distribution ratio changes to 21.1%:78.9%. The results indicate that the conformational transition space decreases in the complex system, and the association with Cmpd # 290 makes the system more stable. Figure S5 shows the slow-motion modes of PaLpxC (A) and PaLpxC-Cmpd # 290 (B). It can be seen from Figure S5 that after binding with the inhibitor Cmpd # 290, the local movement of the PaLpxC protein produce obvious changes. For example, the upper end of the α helices located in R1 (G263-A276) has a tendency to approach the inhibitor directly, as well as the β strand in the R2 (V137-A144) region. In addition, the rotational motion of the whole protein becomes slightly irregular and divergent after the association with inhibitors. It is speculated that the partial loss of the whole rotational motion and the contraction of R1/R2 local conformation both confirm the results of the PCA analyses. Interestingly, the key residues G263 and H264 obtained from energy decomposition calculation (see Figure 7 and Table S6) are located in the R1 region; nevertheless, the relationship between the motion mode change of PaLpxC and its binding with inhibitors remains to be further explored. protein becomes slightly irregular and divergent after the association with inhibitors. It is speculated that the partial loss of the whole rotational motion and the contraction of R1/R2 local conformation both confirm the results of the PCA analyses. Interestingly, the key residues G263 and H264 obtained from energy decomposition calculation (see Figure 7 and Table S6) are located in the R1 region; nevertheless, the relationship between the motion mode change of PaLpxC and its binding with inhibitors remains to be further explored.

Water Mediated H-Bonds
Water is an important reactant of the hydrolysis reaction catalyzed by PaLpxC. In order to further explore the possible inhibition mechanism of Cmpd # 290, the water molecules within 4 Å of the center of the active site of Zn 2+ were sampled and analyzed before and after the association of PaLpxC with Cmpd # 290. The sampling interval is 400 ps and the sample capacity is about 300. Figure 9A shows the changes in the number of water molecules near the active site of PaLpxC in the absence of inhibitors. The result indicates that in the process of simulation, the number of water molecules within 4 Å around Zn 2+ numbered between one and five. The probability of occurrence of one, two, three, four, and five water molecules are about 6.98%, 41.20%, 36.88%, 13.62%, and 1.33%, respectively. The same analysis of the PaLpxC-Cmpd # 290 complex system shows that with the presence of the inhibitor molecule, in the same range of the active site, there is no water molecule. This indicates that Cmpd # 290 can effectively prevent the water molecules from approaching the reaction center.

Water Mediated H-Bonds
Water is an important reactant of the hydrolysis reaction catalyzed by PaLpxC. In order to further explore the possible inhibition mechanism of Cmpd # 290, the water molecules within 4 Å of the center of the active site of Zn 2+ were sampled and analyzed before and after the association of PaLpxC with Cmpd # 290. The sampling interval is 400 ps and the sample capacity is about 300. Figure 9A shows the changes in the number of water molecules near the active site of PaLpxC in the absence of inhibitors. The result indicates that in the process of simulation, the number of water molecules within 4 Å around Zn 2+ numbered between one and five. The probability of occurrence of one, two, three, four, and five water molecules are about 6.98%, 41.20%, 36.88%, 13.62%, and 1.33%, respectively. The same analysis of the PaLpxC-Cmpd # 290 complex system shows that with the presence of the inhibitor molecule, in the same range of the active site, there is no water molecule. This indicates that Cmpd # 290 can effectively prevent the water molecules from approaching the reaction center.
With the simulation trajectory of the PaLpxC system monitored by Visual Molecular Dynamics (VMD) package (University of Illinois at Urbana-Champaign, Urbana-Champaign, IL, USA), it is found that the variation of water molecules near the active site is obvious. Water molecules are present about 2.27 Å away from Zn 2+ on average and tetracoordinate to Zn 2+ . It is consistent with the extended X-ray absorption fine structure (EXAFS) experimental results reported by McClure et al. [34]. In order to further discuss the binding characteristics of the water molecules in this area, we focused on the binding mode of the water molecules shown in Figure 9B. According to Figure 9B, in addition to coordinating with Zn 2+ , the water molecule also forms a stable H-bond with the carboxyl oxygen anion of E77, and a weak H-bond with H264. Previous studies [30,31,35] show that E77 serves as the key residue for the attack and deprotonation of the water molecule, leading to the hydrolysis reaction. By calculation, the distance between the carboxyl oxygen anion of E77 and the hydrogen atom of water molecules is about 1.73 Å, which is beneficial for the hydroxyl ion's attack on the proton. The H-bond between H264 and the water molecule is another important factor stabilizing the hydroxyl ions, and assisting the nucleophilic attack on the acetyl of the natural substrate [35,36]. In addition, it is found that there is a certain interaction between the water molecule and M62 of the PaLpxC system. In some cases, an H-bond was also formed between the water molecule and T190. In general, the stability of the water molecule in this region mainly depends on the localization of multiple H-bonds near the active site of the PaLpxC system. The stable existence of these H-bonds is related to the inherent low flexibility of the amino acids, which is mentioned in the RMSF analysis above. With the simulation trajectory of the PaLpxC system monitored by Visual Molecular Dynamics (VMD) package (University of Illinois at Urbana-Champaign, Urbana-Champaign, IL, USA), it is found that the variation of water molecules near the active site is obvious. Water molecules are present about 2.27 Å away from Zn 2+ on average and tetracoordinate to Zn 2+ . It is consistent with the extended X-ray absorption fine structure (EXAFS) experimental results reported by McClure et al. [34]. In order to further discuss the binding characteristics of the water molecules in this area, we focused on the binding mode of the water molecules shown in Figure 9B. According to Figure 9B, in addition to coordinating with Zn 2+ , the water molecule also forms a stable H-bond with the carboxyl oxygen anion of E77, and a weak H-bond with H264. Previous studies [30,31,35] show that E77 serves as the key residue for the attack and deprotonation of the water molecule, leading to the hydrolysis reaction. By calculation, the distance between the carboxyl oxygen anion of E77 and the hydrogen atom of water molecules is about 1.73 Å, which is beneficial for the hydroxyl ion's attack on the proton. The H-bond between H264 and the water molecule is another important factor stabilizing the hydroxyl ions, and assisting the nucleophilic attack on the acetyl of the natural substrate [35,36]. In addition, it is found that there is a certain interaction between the water molecule and M62 of the PaLpxC system. In some cases, an H-bond was also formed between the water molecule and T190. In general, the stability of the water molecule in this region mainly depends on the localization of multiple H-bonds near the active site of the PaLpxC system. The stable existence of these H-bonds is related to the inherent low flexibility of the amino acids, which is mentioned in the RMSF analysis above.

Possible Inhibitory Mechanism of Cmpd # 290
Based on the analyses of the molecular recognition between PaLpxC and Cmpd # 290, as well as the conformational change of the complex and the water molecules near the active site of PaLpxC, the following possible inhibition mechanism of Cmpd # 290 was proposed. Firstly, Cmpd # 290 is driven by nonpolar interactions, depending on the H-bond interaction with M62, T190, K238, and H264 in the active pocket, and influences the localization of water molecules at the catalytic center. After that, Cmpd # 290 chelates to Zn 2+ with its hydroxamate moiety and forms a pentacoordination configuration, which prevents the water molecule from combining with it. Meanwhile, the hydrophobic terminus of Cmpd # 290 occupies the hydrophobic channel of PaLpxC via mixed hydrophobic interaction with F191 and I197, blocking the entry of other water molecules and natural substrates. Finally, the association with Cmpd # 290 also reduces the conformational transition space of PaLpxC, which makes the PaLpxC molecule more stable, ultimately resulting in inhibition (see Figure 10). In Figure 10, the crystal structure of the AaLpxC-substrate (green) was superimposed well onto the structure of the complex PaLpxC-Cmpd # 290 (cyan), due to high conservation. In the Figure 9. Variation of the number of water molecules within a 4 Å range near the active zinc ion of the PaLpxC system versus simulation time (A); and the binding mode of the key water molecule (B). The unit of distance between the two atoms exhibited in Figure B is Angstrom.

Possible Inhibitory Mechanism of Cmpd # 290
Based on the analyses of the molecular recognition between PaLpxC and Cmpd # 290, as well as the conformational change of the complex and the water molecules near the active site of PaLpxC, the following possible inhibition mechanism of Cmpd # 290 was proposed. Firstly, Cmpd # 290 is driven by nonpolar interactions, depending on the H-bond interaction with M62, T190, K238, and H264 in the active pocket, and influences the localization of water molecules at the catalytic center. After that, Cmpd # 290 chelates to Zn 2+ with its hydroxamate moiety and forms a pentacoordination configuration, which prevents the water molecule from combining with it. Meanwhile, the hydrophobic terminus of Cmpd # 290 occupies the hydrophobic channel of PaLpxC via mixed hydrophobic interaction with F191 and I197, blocking the entry of other water molecules and natural substrates. Finally, the association with Cmpd # 290 also reduces the conformational transition space of PaLpxC, which makes the PaLpxC molecule more stable, ultimately resulting in inhibition (see Figure 10). In Figure 10, the crystal structure of the AaLpxC-substrate (green) was superimposed well onto the structure of the complex PaLpxC-Cmpd # 290 (cyan), due to high conservation. In the binding pocket, there is a large volume clash caused by Cmpd # 290 (blue) and the natural substrate (red) and the hydroxamate moiety chelates to Zn 2+ , thus occupying the coordination position of the water molecule (brown), which also leads to a clash at the hydrophilic terminus of Cmpd # 290. binding pocket, there is a large volume clash caused by Cmpd # 290 (blue) and the natural substrate (red) and the hydroxamate moiety chelates to Zn 2+ , thus occupying the coordination position of the water molecule (brown), which also leads to a clash at the hydrophilic terminus of Cmpd # 290. In summary, the possible structure-activity relationships of pyridone methylsulfone hydroxamate PaLpxC were proposed through 3D-QSAR modeling. Then, the binding mode of the complex PaLpxC-Cmpd # 290 and potential inhibitory mechanism of inhibitors both were explored theoretically. This work will provide a theoretical basis for the molecular design of LpxC-targeting antibiotics. Nevertheless, the essential biological experiment is decisive. We plan to gradually introduce relevant experiments in the follow-up work to verify these modeling results.

Molecular Docking
ChemBio3D Ultra 12.0 (Cambridge Soft, Cambridge, MA, USA) was used to construct the molecular structures of the PMH LpxC inhibitors except for Cmpd # 290, and energy optimization of these structures was performed under the MM2 molecular force field with the convergence condition of root mean square (RMS) less than 0.0001 kcal·mol −1 ·Å −1 . After optimization, the structures were imported into SYBYL-X 1.3 [37] and molecular charges were calculated by the Gasteiger-Hückel method, and then another optimization under the Tripos force field was performed with the energy convergence condition of RMS less than 0.05 kcal·mol −1 ·Å −1 to generate the conformers for molecular docking. The AutoDock4.2 software package [38] was used to perform molecular docking, with sampling based on the Lamarckian genetic algorithm (LGA), total and local search results of energy combined, and the semi-empirical function of the binding free energy used to rank the intermolecular energy scores between the receptor and ligand. The docking force field used in this work was an exclusive force field developed for zinc [Zn(II) or Zn 2+ ] containing systems (i.e., AutoDock4Zn force field). In this force field, a pseudo-atom was introduced to guide the ligand coordinating with zinc, and the geometrical characteristics and bond strength of Zn(II) or Zn 2+ were considered sufficiently, which greatly improved the accuracy of molecular docking for zinc-containing molecular systems [39]. In the current work, we determined the box center according to the Zn 2+ coordinates of PaLpxC (PDB entry: 3UHM), and built a 40 Å × 90 Å × 60 Å rectangular box with a grid space of 0.375 Å. 128 conformations were collected in each docking, and the conformation of the lowest energy in the largest conformation cluster was defined as the near native conformation for the subsequent 3D-QSAR analyses. In summary, the possible structure-activity relationships of pyridone methylsulfone hydroxamate PaLpxC were proposed through 3D-QSAR modeling. Then, the binding mode of the complex PaLpxC-Cmpd # 290 and potential inhibitory mechanism of inhibitors both were explored theoretically. This work will provide a theoretical basis for the molecular design of LpxC-targeting antibiotics. Nevertheless, the essential biological experiment is decisive. We plan to gradually introduce relevant experiments in the follow-up work to verify these modeling results.

Molecular Docking
ChemBio3D Ultra 12.0 (Cambridge Soft, Cambridge, MA, USA) was used to construct the molecular structures of the PMH LpxC inhibitors except for Cmpd # 290, and energy optimization of these structures was performed under the MM2 molecular force field with the convergence condition of root mean square (RMS) less than 0.0001 kcal·mol −1 ·Å −1 . After optimization, the structures were imported into SYBYL-X 1.3 [37] and molecular charges were calculated by the Gasteiger-Hückel method, and then another optimization under the Tripos force field was performed with the energy convergence condition of RMS less than 0.05 kcal·mol −1 ·Å −1 to generate the conformers for molecular docking. The AutoDock4.2 software package [38] was used to perform molecular docking, with sampling based on the Lamarckian genetic algorithm (LGA), total and local search results of energy combined, and the semi-empirical function of the binding free energy used to rank the intermolecular energy scores between the receptor and ligand. The docking force field used in this work was an exclusive force field developed for zinc [Zn(II) or Zn 2+ ] containing systems (i.e., AutoDock4 Zn force field). In this force field, a pseudo-atom was introduced to guide the ligand coordinating with zinc, and the geometrical characteristics and bond strength of Zn(II) or Zn 2+ were considered sufficiently, which greatly improved the accuracy of molecular docking for zinc-containing molecular systems [39]. In the current work, we determined the box center according to the Zn 2+ coordinates of PaLpxC (PDB entry: 3UHM), and built a 40 Å × 90 Å × 60 Å rectangular box with a grid space of 0.375 Å. 128 conformations were collected in each docking, and the conformation of the lowest energy in the largest conformation cluster was defined as the near native conformation for the subsequent 3D-QSAR analyses.

3D-QSAR
The PMH compounds were randomly divided into a training set (31 molecules in total) for constructing the 3D-QSAR models and a test set (7 molecules in total) for testing the predictive ability of the biological activity of the models. The IC 50 value of the experiment was converted to a negative logarithm pIC 50 form. Before constructing the 3D-QSAR models, it was necessary to align the compound molecules. The quality of molecular alignment has a great influence on the reliability of the 3D-QSAR models. In molecular alignment, ensuring all molecules using the construction of 3D-QSAR models with the largest similar spatial orientation of molecular fields was helpful to improve the similarity of the steric and electrostatic fields computed. Considering that the 3D-QSAR models assume that all molecules are in the same active binding mode at the same site, the crystallographic structure of Cmpd # 290 was selected as the template molecule, and molecular alignment was based on the pyridone methylsulfone hydroxamate commom substructure. 3D-QSAR analyses were completed in SYBYL-X 1.3 [37]. To set up the CoMFA model, the PMH molecules with alignment were placed in a spatial grid, a sp 3 hybridized C + probe walks in the space around the molecules, the interaction between the probe and molecules was calculated, and the energy values of the interaction with different space coordinates were recorded to obtain the steric and electrostatic field data; CoMSIA was similar to CoMFA, but in addition to the analyses of the steric field (S) and electrostatic field (E), the hydrophobic field (H), H-bond acceptor (A) field, and H-bond donor (D) field were also calculated.
The partial least squares (PLS) method was used to analyze the PMH molecules in training set. Firstly, the leave-one-out (LOO) method was used for cross validation (parameters are default) to obtain the optimal numbers of component (ONC) and the determination coefficient q 2 . Secondly, according to the ONC, the 3D-QSAR models were established by non-cross validation, and the correlation coefficient r 2 , standard error for estimate E s , root mean squared error (RMSE) and the F-test values F were obtained. Parameters based on the PLS analyses above can be used to evaluate the predictive ability and stability of the models, and predict the biological activity of the molecules in the test set. The predicted correlation coefficient for the test set r 2 p was calculated by Equation (2): In Equation (2), SD represents the sum of the squared deviations of the experimental biological activity of the test set and the average molecular biological activity of the training set, and PRESS indicates the sum of the squared errors of the predictive biological activity values of the test set and their biological activity as determined by the experiment.

MESP and Mulliken Charge Analyses
Gaussian 09 (Gaussian Inc., Wallingford, CT, USA) was used to calculate the MESP and Mulliken charges based on DFT. The molecular conformation of the ligand used in the calculation was generated by energy optimization, and the method and criterion of energy optimization were consistent with those of the previous molecular docking. The B3LYP/6-31G (d,p) basis set [40] was applied to optimize the geometry of the PMH LpxC inhibitors. In order to simulate the physiological environment as much as possible, all the energy was calculated in the solvent environment in the CPCM model [41]. The calculation of the MESP and Mulliken atomic charges was completed at the same theoretical level. The MESP iso-energy surface was superimposed by the color coded iso-potential surface onto the iso-electron density surface (0.001 e au −3 ). Different colors are used to represent the distribution of the MESP on the surface of the molecule. For example, dark red represents the region of the most negative electrostatic potential and dark blue represents the region of the most positive electrostatic potential.

Molecular Dynamics Simulation
The structural models of PaLpxC and complex PaLpxC-Cmpd # 290 were further subjected to MD simulations in the solution system for 120 ns, respectively, using the AMBER16 software package [42] with the AMBER ff14SB force field. The force field parameters were fitted based on the experimental results [43], and the simulated temperature was 300 K. The solvent was added with the TIP3P water model [44], and formed a 10 Å thick layer of water molecules on the surface of the systems. Using MCPB.py [45] to call Gaussian 09, according to the B3LYP method, the Zn 2+ bond model and the force field parameters were constructed by using 6-31G* basis set. Energy minimization before simulation was completed by two steps. Firstly, with the constraint (constraint force is 500 kcal·mol −1 ·Å −2 ) of the solute, 5000 steps of the steepest descent method was performed and followed by another 5000 steps of the conjugate gradient method; then unrestrained minimizations of the whole system were executed by the same stages above with the convergence condition that the energy gradient is less than 0.01 kcal·mol −1 ·Å −2 .
MD simulation consists of two steps: constrained MD simulation was performed for 1 ns with the constraint force constant at 10 kcal·mol −1 ·Å −2 and the temperature of the system gradually increased from 0 K to 300 K. Then free isothermal MD simulation of the systems was performed for 119 ns. The conformational changes in the whole process of simulations were tracked and observed by VMD, and using the SHAKE algorithm [46] to restrain the bond length, the radius of non-bond was 0.1 Å, and the time integration step was set to 2 fs. 12,768 and 11,162 water molecules were added to the PaLpxC system and the complex PaLpxC-Cmpd # 290 system. To balance the total charge of the systems, 7 and 6 Na + were also added, and the total numbers of atoms in the PaLpxC and PaLpxC-Cmpd # 290 systems were 38,310 and 38,211, respectively.

Energy Decomposition Analysis
Energy decomposition analysis was performed by the molecular mechanics/generalized Born surface area (MM/GBSA) method [47]. The basic idea of this method is that: (1) the energy contribution of each residue can be divided into: the internal energy in a vacuum; the polar solvation energy calculated by the Generalized Born model [48]; and the nonpolar solvation energy calculated by LCPO algorithm [49]; (2) energy is decomposed onto the atoms of the backbone and side chain of each residue. Here, the internal energy in a vacuum can also be divided into two parts: polar electrostatic interaction and nonpolar van der Waals interaction. According to the LCPO algorithm, there is a positive correlation between nonpolar solvation energy and solvent accessible surface area (SASA). Thus, the contribution of the main residues in PaLpxC to the binding of Cmpd # 290 can be observed by MM/GBSA energy decomposition.

Free Energy Landscape Analysis
FEL [50] can promote our understanding of molecular recognition, folding, and aggregation by studying the free-energy minima and the free-energy barriers. The free-energy minima represent the conformational ensemble in stable states which are accessible to a biomolecule under physiological conditions and the free-energy barriers denote the transient states connecting them. The FEL can be constructed based on PCA [51]. The corresponding expression is where X stands for the PCs and thus P(X) is the probability distribution of the molecular system along the PCs, and k B and T are the Boltzmann constant and absolute temperature, respectively.

Conformational Cluster Analysis
On the basis of the conformations obtained by MD simulation, conformational cluster analysis was performed using MMTSB tools [52]. The specific analysis process is done according to the view of Daura et al. [53]: calculating the RMSD value between the conformations one by one, as a basis for the establishment of RMSD matrix (N × N, N stands for the numbers of conformations). After setting a RMSD threshold manually, and comparing any two RMSD values of conformations, if the difference is less than the threshold value, it is classified as a cluster. The lowest energy conformation in each cluster is used as the representative conformation. Here, the matrix was built by the RMSD of the whole system, and the situation of conformational cluster of overall systematic conformations is discussed.

Conclusions
In the current work, theoretical molecular simulation studies on a broad-spectrum LpxC inhibitor, a type of pyridone methylsulfone hydroxamate compound, were performed to investigate the relationship between its chemical structure and bioactivity. CoMFA and CoMSIA were used to study the 3D-QSAR of PMH LpxC, and QSAR models with good predictive ability were obtained. According to the results of the contour map, the introduction of bulky groups and electron withdrawing groups in PMH LpxC inhibitors can effectively improve the inhibitory activity of the inhibitors. Based on the analyses of the molecular electrostatic potential of the PMH LpxC inhibitors, the effect of substituents on the inhibitory activity was explained at the molecular level.
Comparative MD simulations of the PaLpxC and complex PaLpxC-Cmpd # 290 were carried out for 120 ns. The results of binding free energy calculations based on MD trajectories show that the binding free energy calculation results agree well with the experimental values, and the van der Waals interaction is the main driving force for the formation of the complex. The results of H-bond and energy decomposition show that Cmpd # 290 mainly interacts with M62, T190, F191, H264, I197, and K238 of PaLpxC.
Finally, through the simulation of trajectory sampling analyses, the binding mode of key water molecules in PaLpxC was investigated. Without any constraint, the water molecules can freely access the active pocket of PaLpxC; and with the assistance of E77, H264, T190, and M62, it can coordinate and form a tetrahedral configuration with Zn 2+ (the average length of O-Zn is around 2.27 Å), to provide the reagent for the hydrolysis reaction constantly. In the complex structure, Cmpd # 290 chelates Zn 2+ by hydroxamate moiety, which causes the destruction of the binding mode of water molecules. It further blocks the entry of the reaction substrate by enclosing the pocket through its hydrophobic chain, so as to exert an inhibitory effect. This work provides some theoretical guidance for the design of drug molecules targeting PaLpxC.