Estimation of Bond Free Energy with Gmx_Mmpbsa in Ndm-1 Complexes

We report the values obtained for the binding free energies (ΔG_bind) of the complexes: NDM-1-M25 (43.80 kcal/mol), NDM-1-M26 (12.71 kcal/mol), NDM-1-M35 (19.92 kcal/mol), NDM-1-M37 (2.46 kcal/mol) and the reference system NDM-1-Meropenem (-10.09 kcal/mol). These results are based on previous absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties, docking, and molecular dynamics calculations performed on a small library of potential NDM-1 enzyme inhibitors reported in the literature. Subsequently, the selected systems were evaluated by the Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) method using the gmx_MMPBSA software. The MMPBSA method has been widely adopted to estimate protein-ligand binding affinities due to its efficiency and high correlation with experimental data. The obtained results suggest that, unlike the reference complex, the studied ligands probably do not bind to the receptor under the simulation conditions, which is in significant agreement with experimental data.

Introduction

Bacterial resistance to β-lactam antibiotics represents a significant threat to global public health, in part due to the activity of the New Delhi metallo-β-lactamase (NDM-1) enzyme. This enzyme, present in bacteria such as Klebsiella pneumoniae and Escherichia coli, has the ability to hydrolyze the β-lactam core of antibiotics, resulting in their inactivation and drastically reducing the therapeutic options available to treat serious infections [1]. Due to its clinical impact and rapid spread in hospital settings, the development of new NDM-1-specific inhibitors has become a priority in the fight against antibiotic resistance [2]. Research on NDM-1 inhibitors not only seeks to restore the efficacy of antibacterial treatments but also to curb the spread of these multidrug-resistant pathogens [3].

Within the process of developing new inhibitors, the calculation of molecular docking stands out as an essential computational tool. This method allows for predicting how ligands interact with the active site of NDM-1, revealing possible key interactions and preferred binding modes [4]. By obtaining predictions about the affinities of compounds to the enzyme, docking facilitates the identification of promising candidates, thus optimizing the drug discovery process. In addition, this approach allows multiple conformational configurations to be explored, providing a broad view of the inhibitory potential of the molecules studied [5].

The importance of Molecular Dynamics (MD) simulations lies in their ability to provide a more detailed understanding of the structural stability and conformational flexibility of protein-ligand complexes over time [6]. Unlike docking, which provides a static view of the interaction, molecular dynamics allows one to assess how these interactions evolve under conditions that simulate the real biological environment. This is particularly useful in the case of NDM-1 since the stability of the inhibitor-enzyme complex is crucial for its efficacy. MD simulations, moreover, provide information on conformational changes of the protein, which can influence the design and optimization of more effective inhibitors.

Finally, the calculation of binding energies using the MMPBSA (Molecular Mechanics Poisson-Boltzmann Surface Area) method is an essential step to validate the results obtained in the previous simulations. This method allows for estimating the binding free energy in an aqueous environment, taking into account both electrostatic interactions and non-polar forces [7]. By applying MMPBSA to the complexes generated in molecular dynamics simulations, it is possible to more accurately predict the affinity of ligands for NDM-1, providing a closer correlation with experimental results. This approach not only complements computational findings but also guides the rational design of inhibitors, optimizing their therapeutic potential [8].

In the present investigation, a group of potential NDM-1 enzyme inhibitors, selected from a comprehensive literature review [9-14], were evaluated. These studies identified compounds derived from cyclic borates and chromones as potential inhibitors of the metalloenzyme, in addition to highlighting their promising values of oral bioavailability and low toxicity [9,12]. As a result, four molecules were selected (Figure 1), to which binding affinity energy values were calculated using the MMPBSA method, in complexes formed with NDM-1, and these data were compared with those obtained for the NDM-1-Meropenem reference system. The aim of this analysis was to determine whether the values obtained in the molecular docking and molecular dynamics studies corroborated with the free binding energies (ΔG_bind) between the ligands and the receptor. Finally, it was established that, although compounds M25, M26, M35, and M37 presented promising values in terms of docking energy and showed stable molecular dynamics comparable to the reference system, the MMPBSA calculations were conclusive in indicating that these inhibitors were unlikely to bind to the receptor under the conditions studied. This is because all binding energy values obtained were positive, in contrast to the total ΔG_bind value of the reference system, which was -10.09 kcal/mol, evidencing strong binding in that complex.

Figure 1: Structures at pH=7.5 of the 4 compounds with the best coupling energy values.

Materials and methods

Data obtained from previous research by the authors, including ADMET properties, molecular docking energies, and 10 ns molecular dynamics simulations, were used for a group of molecules with inhibitory potential against the NDM-1 enzyme (PDB code: 5ZGZ) (www.rcsb.org). These systems showed promising results both in terms of docking energies and RMSD, RMSF, and hydrogen bonding (H-BOND) calculations derived from molecular dynamics, compared to the reference NDM-1-Meropenem complex. To estimate the binding free energies of the systems under study, the Molecular Mechanics/Poisson-Boltzmann Surface (MM/PBSA) method was employed and implemented using the gmx_mmpbsa software [15]. The analysis was carried out using the files generated by the molecular dynamics of the complexes from GROMACS 2024.1 [16]. In the calculations of binding free energies, parameters such as the number of frames were considered, with a total of 300 frames analyzed at a temperature of 298.15K. The inter-frame interval was 1 frame, while the rest of the parameters were kept as specified in the mmpbsa.in file. There are several quantum mechanical methods to evaluate systems, however, for biochemical systems, the most efficient and reliable method is the MM/PBSA, especially in the estimation of binding free energies, as previously mentioned [17].

Results and discussion

Four systems (NDM-1-M25, NDM-1-M26, NDM-1-M35, and NDM-1-M37), which had previously shown promising results in terms of docking energies and stability in molecular dynamics, evaluated in terms of structural stability and flexibility, were analyzed. These systems were compared with the NDM-1-Meropenem reference system (Table 1). Docking energy values obtained by docking showed negative interactions on the kcal/mol scale, suggesting favorable binding of ligands to the active site of NDM-1. For example, the NDM-1-M25 complex presented a ΔG value of -10.6 kcal/mol, as did the NDM-1-M26 complex. However, these results are preliminary and do not reflect physiological conditions, as they do not include solvation effects or full protein flexibility.

To address these limitations, binding energy calculations were performed using the MMPBSA method, which decomposes the energetic contributions into components such as van der Waals energy (ΔVDWAALS), gas-phase electrostatic energy (ΔGGAS), solvation energy (ΔGSOLV), and total energy (ΔTOTAL). The results of these energies, presented in Table 1, indicate that the complexes formed between NDM-1 and potential inhibitors have a range of binding energies from 2.46 to 43.80 kcal/mol, suggesting that compounds M25, M26, M35, and M37 would not bind favorably to the receptor, since their binding energies are positive. In contrast, the NDM-1-Meropenem system presented a total binding energy of -10.09 kcal/mol, which reinforces the hypothesis that those systems whose binding free energy values do not approach this negative value would not be good candidates as potential inhibitors of the NDM-1 metalloenzyme.

The results show that van der Waals energies (ΔVDWAALS) and gas-phase electrostatic energy (ΔGGAS) are key components influencing ligand binding affinity. In particular, van der Waals interactions have been shown to be favorable in all complexes analyzed, with negative values ranging from -20.13 to -31.89 kcal/mol, suggesting relative stability in binding to the carbapenemase active site. For the gas-phase electrostatic energy, negative values were obtained for all the systems under study, which is also favorable, except for the NDM-1-M35 complex, which presented a value of 23.7 kcal/mol. However, the values for the solvation energy (ΔGSOLV) were positive, with the exception of the NDM-1-M35 system, which had a value of -3.77 kcal/mol, resulting in unfavorable total binding energies compared to the reference system.

Plots of the energetic variations (ΔVDWAALS, ΔGGAS, ΔGSOLV, and ΔTOTAL) for the complexes under study are presented in Figure 2. For the NDM-1-M25 complex, an energy fluctuation ranging from approximately 22 to 59 kcal/mol is observed throughout the simulation. However, since the energy range is positive, this indicates that the inhibitor will not bind favorably to the receptor, suggesting instability in this system. For reference, the plots corresponding to the NDM-1-Meropenem system, composed of the commercial antibiotic Meropenem and the NDM-1 enzyme, show a total energy fluctuation between -20 and -0.2 kcal/mol during the simulation, stabilizing around -10 kcal/mol. In the NDM-1-M26 and NDM-1-M35 systems, a behavior similar to that of NDM-1-M25 is evident, with energetic variations in positive ranges throughout the simulation time, which also indicates unfavorable binding. In the case of NDM-1-M37, it is important to evaluate in greater detail the energetic fluctuations and compare them with the reference system plots, since they present negative values reaching almost -9 kcal/mol in frames 220-222. However, the ΔTOTAL value stands at 2.46 kcal/mol, which is thermodynamically unfavorable and indicates the lack of binding to the receptor. These results suggest that M37 could be considered as a starting compound for the development of new structures. This could improve the binding efficiency of potential new inhibitors of the NDM-1 metalloproteinase by making specific structural modifications that optimize van der Waals interactions, hydrogen bridging, and hydrophobic interactions at the enzyme active site. These modifications have been shown to improve various drug-likeness properties, as reported in previous research [18-20]. Modifications include the addition of electron-donating groups, such as hydroxyl, alkyl groups, amines, or aryl groups (in certain cases of conjugated systems), in proximity to the polar regions of the active site, which could intensify hydrogen bond formation. In addition, hydrophobic adjustments that improve compatibility with non-polar regions of the active site could take advantage of pi-type interactions or hydrophobic contacts, thus strengthening the bond [21]. On the other hand, higher structural rigidity of ligands could reduce the entropy of excited states, favoring the stability of the protein-ligand complex. It has been shown that higher ligand stiffness correlates with a decrease in entropy, resulting in more stable binding to the target protein. Also, a detailed analysis of how these modifications impact specific energetic terms, such as solvation energy and van der Waals interactions calculated using MMPBSA, would provide a deeper understanding of their effect on binding affinity [7,8].

Figure 2: Variation of the energies (kcal/mol) of gas, solvent, total, and formation of all interactions occurring throughout the simulation of systems A (NDM-1-25), B (NDM-1-26), C (NDM-1-35), D (NDM-1-37) and E (NDM-1-Meropenem).

On the other hand, the accuracy of the MMPBSA method is subject to several limitations that it is essential to discuss in order to improve the understanding of the obtained results. First, the use of the AMBER force field in the calculations may influence the accuracy of the calculated energies, since the choice of the force field is crucial for evaluating protein-ligand interactions, especially electrostatic and van der Waals interactions [7]. Moreover, the length of the simulation, which in this study was 10 ns, may not be sufficient to capture long-term conformational events, which could lead to underestimation or overestimation of the stability of the complex. Extending the duration of the simulations, e.g., to 100 ns, could provide a clearer picture [22].

It is also important to note that the solvation treatment, using the implicit Poisson-Boltzmann method, might not adequately reflect solvent interactions in complex systems; in this context, the inclusion of explicit solvation models could improve the accuracy of the calculations [8]. Addressing these limitations will allow us to justify the results and propose strategies to minimize biases in future research.

From the study, it can be considered that the results of the NDM-1-Meropenem system can be used as a reference, and any significant deviation in terms of stability and free energy in the other systems suggests a lower inhibitory potential. This was clearly evidenced in the systems consisting of NDM-1 and the ligands M25, M26, M35, and M37.

For future research, it is recommended that empirical validation be incorporated to significantly complement the computational calculations performed in this study. This would involve the synthesis of the proposed ligands and experimental evaluation of their affinities for NDM-1 by enzymatic binding or inhibition assays. In addition, it would be beneficial to develop inhibitors from the most promising ligands, based on the results obtained from molecular dynamics (MD) and MMPBSA simulations. It is also suggested that a more comprehensive analysis including a larger number of structural analogues of NDM-1 inhibitors be carried out, which would provide a more robust understanding of the binding properties and facilitate the identification of patterns in ligand efficacy, thus guiding the design of new and more effective inhibitors. Finally, experimental evaluation of the observed interactions using techniques such as X-ray crystallography or spectroscopy (NMR) is recommended, which would allow corroboration and validation of the established theoretical models.

Conclusion

The study of NDM-1 enzyme inhibitors is fundamental in the fight against bacterial resistance, highlighting the importance of the use of computational tools such as molecular docking, molecular dynamics simulations, and MMPBSA calculations which provide a solid foundation for effective drug development. Moreover, these approaches not only improve our understanding of key molecular interactions but also accelerate the process of discovering new treatments, offering an innovative response to one of the major global public health crises. On the other hand, results obtained using the MMPBSA method indicate that, although ligands M25, M26, M35, and M37 showed relatively stable molecular dynamics during 10 ns simulations, their interaction energies fluctuated within positive ranges, suggesting instability and a lack of effective binding to the NDM-1 receptor under the simulation conditions used. In fact, all ligands presented positive binding free energy values, which contrasts with the NDM-1-Meropenem reference system, which showed a negative value of -10.09 kcal/mol, indicating a strong binding affinity. However, although the compounds studied did not show favorable binding under the present conditions, the M37 ligand exhibited energetic fluctuations reaching almost -9 kcal/mol at certain times during the simulation, suggesting that it could be a good starting point for the development and optimization of more effective inhibitors. Overall, these results highlight the potential of computational simulations to guide the rational design of NDM-1 inhibitors, despite the challenges that remain.

Authors' contributions

Eduvan Valencia: Writing – original draft, Formal analysis, Methodology, Conceptualization Data curation.

Mauricio Galvis Patiño: Supervision, Conceptualization.

1. Haque M, Sartelli M, McKimm J, Abu Bakar M. Health care-associated infections - an overview. Infect Drug Resist. 2018;11:2321-2333. Available from: https://doi.org/10.2147/IDR.S177247
2. Yang SK, Kang JS, Oelschlaeger P, Yang KW. Azolylthioacetamide: a highly promising scaffold for the development of metallo-β-lactamase inhibitors. ACS Med Chem Lett. 2015;6:455-460. Available from: https://doi.org/10.1021/ml500534c
3. Salih T, Ali PG. Rational design of novel compounds to serve as potential NDM-1 inhibitors using molecular docking, molecular dynamics simulation, and physicochemical studies. Mol Simul. 2023;49:1373-1387. Available from: https://doi.org/10.1080/08927022.2023.2232468
4. Morris GM, Huey R, Lindstrom W, Sanner MF, Belew RK, Goodsell DS, et al. AutoDock4 and AutoDockTools4: automated docking with selective receptor flexibility. J Comput Chem. 2009;30(16):2785-2791. Available from: https://doi.org/10.1002/jcc.21256
5. Etruri B, Jardini A, Maciel R. In silico evaluation of most used drugs on the treatment of slight and mild COVID-19 cases in Brazil. Int J Dev Res. 2021;11:48146-48150. Available from: https://www.journalijdr.com/silico-evaluation-most-used-drugs-treatment-slight-and-mild-covid-19-cases-brazil
6. Kar B, Kundu CN, Pati S, Bhattacharya D. Discovery of phyto-compounds as novel inhibitors against NDM-1 and VIM-1 protein through virtual screening and molecular modelling. J Biomol Struct Dyn. 2023;41:1267-1280. Available from: https://doi.org/10.1080/07391102.2021.2019125
7. Wang C, Nguyen PH, Pham K, Huynh D, Le TBN, Wang H, et al. Calculating protein-ligand binding affinities with MMPBSA: method and error analysis. J Comput Chem. 2016;37(27):2436-2446. Available from: https://doi.org/10.1002/jcc.24467
8. Genheden S, Ryde U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin Drug Discov. 2015;10(5):449-461. Available from: https://doi.org/10.1517/17460441.2015.1032936
9. Wang T, Xu K, Zhao L, Tong R, Xiong L, Shi J. Recent research and development of NDM-1 inhibitors. Eur J Med Chem. 2021;223:113667. Available from: https://doi.org/10.1016/j.ejmech.2021.113667
10. Liu S, Jing L, Yu ZJ, Wu C, Zheng Y, Zhang E, et al. ((S)-3-Mercapto-2-methylpropanamido) acetic acid derivatives as metallo-β-lactamase inhibitors: synthesis, kinetic and crystallographic studies. Eur J Med Chem. 2018;145:649-660. Available from: https://doi.org/10.1016/j.ejmech.2018.01.032
11. Skagseth S, Akhter S, Paulsen MH, Muhammad Z, Lauksund S, Samuelsen Ø, et al. Metallo-β-lactamase inhibitors by bioisosteric replacement: preparation, activity and binding. Eur J Med Chem. 2017;135:159-173. Available from: https://doi.org/10.1016/j.ejmech.2017.04.035
12. Christopeit T, Albert A, Leiros H-KS. Discovery of a novel covalent non-β-lactam inhibitor of the metallo-β-lactamase NDM-1. Bioorg Med Chem. 2016;24(13):2947-2953. Available from: https://doi.org/10.1016/j.bmc.2016.04.064
13. Zhang YL, Yang KW, Zhou YJ, LaCuran AE, Oelschlaeger P, Crowder MW. Diaryl-substituted azolylthioacetamides: inhibitor discovery of New Delhi metallo-β-lactamase-1 (NDM-1). ChemMedChem. 2014;9(11):2445-2448. Available from: https://doi.org/10.1002/cmdc.201402249
14. Shen B, Yu Y, Chen H, Cao X, Lao X, Fang Y, et al. Inhibitor discovery of full-length New Delhi metallo-β-lactamase-1 (NDM-1). PLoS One. 2013;8(5). Available from: https://doi.org/10.1371/journal.pone.0062955
15. Valdés-Tresanco MS, Valdés-Tresanco ME, Valiente PA, Moreno E. gmx_MMPBSA: a new tool to perform end-state free energy calculations with GROMACS. J Chem Theory Comput. 2021;17(10):6281–91. Available from: https://doi.org/10.1021/acs.jctc.1c00645
16. Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, Berendsen HJC. GROMACS: fast, flexible, and free. J Comput Chem. 2005;26:1701–18. Available from: https://doi.org/10.1002/jcc.20291
17. da Silva Moreira LP. Protein using molecular dynamics [dissertation]. Coimbra: Universidade de Coimbra; 2024.
18. Cuervo Prado PA. Rational design of new purine α, γ-diketo acids and their bioisosteric derivatives with potential inhibitory activity of the HIV integrase enzyme [dissertation]; 2014.
19. Valencia Cristancho ED. Rational design of spirothiazolidonic, pyrazolothiazolidonic and pyrazolo β-lactam compounds with potential antimicrobial activity [dissertation]; 2014.
20. Srinivas A, Nagaraj A, Reddy CS. Synthesis of some novel methylene-bis-pyrimidinyl-spiro-4-thiazolidinones as biologically potent agents. J Heterocycl Chem. 2008;45(4):1121-1125. Available from: https://doi.org/10.1002/jhet.5570450428
21. Young T, Abel R, Kim B, Berne BJ, Friesner RA. Motifs for molecular recognition exploiting hydrophobic enclosure in protein–ligand binding. Proc Natl Acad Sci USA. 2007;104(3):808-813. Available from: https://doi.org/10.1073/pnas.0610202104
22. Jorgensen WL, Tirado-Rives J. Molecular modeling of organic and biomolecular systems using BOSS and MCPRO. J Comput Chem. 2005;26(16):1689-1700. Available from: https://doi.org/10.1002/jcc.20297