key: cord-0023575-5465z96e
authors: Almihyawi, Raed A. H.; Al-Hasani, Halah M. H.; Jassim, Tabarak Sabah; Muhseen, Ziyad Tariq; Zhang, Sitong; Chen, Guang
title: Molecular Insights into Binding Mode and Interactions of Structure-Based Virtually Screened Inhibitors for Pseudomonas aeruginosa Multiple Virulence Factor Regulator (MvfR)
date: 2021-11-11
journal: Molecules
DOI: 10.3390/molecules26226811
sha: b1c16be28033ce5b1018594524c3d59cd5daf37e
doc_id: 23575
cord_uid: 5465z96e

Multi-drug resistance (MDR) bacterial pathogens pose a threat to global health and warrant the discovery of new therapeutic molecules, particularly those that can neutralize their virulence and stop the evolution of new resistant mechanisms. The superbug nosocomial pathogen, Pseudomonas aeruginosa, uses a multiple virulence factor regulator (MvfR) to regulate the expression of multiple virulence proteins during acute and persistent infections. The present study targeted MvfR with the intention of designing novel anti-virulent compounds, which will function in two ways: first, they will block the virulence and pathogenesis P. aeruginosa by disrupting the quorum-sensing network of the bacteria, and second, they will stop the evolution of new resistant mechanisms. A structure-based virtual screening (SBVS) method was used to screen druglike compounds from the Asinex antibacterial library (~5968 molecules) and the comprehensive marine natural products database (CMNPD) (~32 thousand compounds), against the ligand-binding domain (LBD) of MvfR, to identify molecules that show high binding potential for the relevant pocket. In this way, two compounds were identified: Top-1 (4-((carbamoyloxy)methyl)-10,10-dihydroxy-2,6-diiminiodecahydropyrrolo[1,2-c]purin-9-yl sulfate) and Top-2 (10,10-dihydroxy-2,6-diiminio-4-(((sulfonatocarbamoyl)oxy)methyl)decahydropyrrolo[1,2-c]purin-9-yl sulfate), in contrast to the co-crystallized M64 control. Both of the screened leads were found to show deep pocket binding and interactions with several key residues through a network of hydrophobic and hydrophilic interactions. The docking results were validated by a long run of 200 ns of molecular dynamics simulation and MM-PB/GBSA binding free energies. All of these analyses confirmed the presence of strong complex formation and rigorous intermolecular interactions. An additional analysis of normal mode entropy and a WaterSwap assay were also performed to complement the aforementioned studies. Lastly, the compounds were found to show an acceptable range of pharmacokinetic properties, making both compounds potential candidates for further experimental studies to decipher their real biological potency.

Infectious diseases are a major reason for human disorders, particularly in low income countries [1, 2] . Infectious diseases have been the top cause of deaths around the globe for a long time and have high economic costs [3, 4] . Multi-drug-resistant bacterial species emerged as a serious threat to public health and are classified by the World Health Organization (WHO) as one of the top 10 health problems that humanity is currently facing [5] [6] [7] . Antibiotic resistance, in particular, is of great concern in six highly virulent bacterial species (Enterococcus faecium, Staphylococcus aureus, Klebsiella pneumoniae, Acinetobacter baumanii, Pseudomonas aeruginosa, and Enterobacter spp.) (commonly referred to as ESKAPE pathogens) [6, 8] . The design of new drugs against the mentioned antibioticresistant bacterial pathogens involves a constant search and the unveiling of new chemically diverse molecules to tackle ESKAPE pathogens requires more time [9] .

Gram-negative bacilli of the genus Pseudomonas are found in freshwater, soil, and marine environments [10] . P. aeruginosa is a frequent causative pathogen of nosocomial infections such as bacteremia, pneumonia, urinary tract infections, meningitis, and damaged mucous membranes or skin, the latter enabling pathogens to enter the blood circulation and cause septicemia [11] . The infectious bacteremia caused by P. aeruginosa has a higher mortality rate than other species of Pseudomonas due to its higher resistance spectrum against many of the antibiotics [12] . It is a ubiquitous pathogen that has the natural ability to thrive in moist environments and show resistance to many antiseptics and antibiotics, and thus, is commonly found in hospital intensive care units [13] . The resistance is multifactorial and is mediated by porins, penicillin-binding proteins, efflux pumps, chromosomal β-lactamases, and aminoglycoside-modifying enzymes, all of which contribute to resistance against antibiotics that are commonly used for treating P. aeruginosa infections [14] . The multi-drug resistance in this pathogen has made it critical to come up with new antimicrobial drugs.

P. aeruginosa survives the action of antibiotics through the formation of dormant cells known as antibiotic-tolerant/persister (AT/P) cells [15] . In these cells, the metabolic state is suppressed, enabling tolerance to lethal antibiotic concentrations. It was demonstrated that multiple virulence factor regulator (MvfR) plays a key role in the formation of AT/P cells and the regulation of different virulence functions in P. aeruginosa [16] . In order to block the function and to design anti-virulent drugs, the current study uses different applications of computer aided drug design (CAAD) [17] . Computational approaches are of significant importance in the process of drug discovery and development [18] [19] [20] . The search for specific and selective novel drug targets against bacterial pathogens is an important step in the design of new drug molecules to fight bacterial infections. This in silico study aims to identify potential inhibitory molecules against P. aeruginosa that can be developed as drugs.

The objective is to screen high-affinity binders from antibacterial and natural databases. Virtual screening was performed to prioritize the best-docked molecule for the MvfR, followed by a biophysical analysis of molecular dynamics simulation and binding free energies to validate the docking predictions. The findings of this study can help in the identification of novel leads against nosocomial P. aeruginosa infections.

Initially, the crystal structure of P. aeruginosa MvfR was retrieved from the protein data bank (PDB) using the PDB ID of 6B8A [21] . The MvfR crystal structure was of 2.65 Å resolution, and had an R-Value Free score of 0.251 and an R-Value Work score of 0.216 [16] . The enzyme was visualized in UCSF Chimera version 1.15 [22] , and was analyzed to prepare it for the molecular docking study. The water molecules and associated co-crystallized ligand (M64 compound) were deleted from the protein structure. The structure then entered the energy minimization phase of 2000 steps: 1000 steps of the steepest descent algorithm (to ease highly unfavorable clashes) and 1000 steps of the conjugate gradient algorithm (a slower algorithm that is effective at reading the energy minimum). The said algorithms were run at a default step size of 0.02 Å. AMBER ff14SB [23, 24] was used to assign charges to the protein residues.

In order to discover novel chemical entities of good safety and strong antibacterial activity, valuable starting leads are indeed of significance. To obtain these, several antibacterial libraries were utilized for use in structure-based virtual screening. The libraries employed in this study included the Asinex antibacterial library (~5968 molecules) (http://www.asinex.com/?page_id=14/, accessed on 7 August 2021) and the comprehensive marine natural products database (CMNPD) (~32 thousand compounds) [25] . Both libraries were imported to FAFDrugs4 (https://mobyle.rpbs.univ-paris-diderot.fr/cgibin/portal.py?form=FAF-Drugs3#forms::FAF-Drugs4, accessed on 8 August 2021) [26] to be filtered based on Lipinski's rule of five [27, 28] . The filtered compounds were then transferred to PyRx v0.8 [29] and converted into .pdbqt, and their energy was minimized via an MM2 force field [30] .

After preparing both the MvfR receptor and the drug libraries, SBVS was performed, targeting the ligand-binding domain of MvfR, which was collectively formed by two subdomains connected through antiparallel β-sheets. The ligand-binding domain was hydrophobic and comprised active residues (Thr265 (X: 21.693 Å, Y: −20.196 Å and Z: 5.844 Å) and Leu189 (X: 15.868 Å, Y: −31.230 Å and Z 0.342 Å)) that were reported to be in regular contact with the M64 co-crystallized compound [16] . Virtual screening of the libraries was achieved using the AutoDock Vina program [31] and GOLD 5.2 [32] , where the grid box was centered at the above-mentioned residues with dimensions along the XYZ axes of 25 Å. To be certain about the docking protocol, the co-crystalized ligand was extracted and docked to the MvfR blindly. After confirmation of the docking method, the ligand libraries were screened against the targeted pocket of the MvfR. The number of poses generated for each compound was tuned to 100; these were clustered, and the ones with the lowest binding energy scores and the greatest numbers of hydrophobic and hydrophilic interactions were selected for complex formation. In total, 3 complexes were selected, including one control (M64), for further analysis.

Molecular dynamics simulations were performed to evaluate the binding mode of the leads and the control. Docking results are usually not satisfactory and post-moleculardocking analyses, including molecular dynamics simulation and binding free energies, are widely applied to validate docking predictions [33] [34] [35] [36] . The AMBER20 simulation package [37] was used to perform all atom simulations. The Antechamber program [38] was employed to recognize the atom type and bond type, to find the missing force field parameters and provide similar substitutes, and to generate the topology files. This module was considered to automatically generate drug molecules and protein input parameters for simulation. Further, parametrization of the compounds and the MvfR was conducted using AMBER GAFF [39] and the ff14SB force field [23] , respectively. After preparation, each complex was placed in a TIP3P water box of 12 Å dimensions (to ensure that the box size was sufficient to allow proper complex dynamics and that the opposite parts of the complex from the adjacent cells did not see each other), which was then treated with a suitable number of Na + counter ions (9 in total) to obtain a neutral charge system ( Figure 1 ). To prepare the systems for the production run, the complexes were subjected to different phases of energy minimization, as follows: energy minimization for a total of 3000 steps (hydrogen atom energy minimization, water box energy minimization, complex energy minimization, and non-heavy atom minimization) with different set restraints as used previously [40] . The number of cycles of each energy minimization step was adjusted according to the studied systems. Further, an NVT ensemble was used to gradually heat the systems from 0 to 300 K at a time interval of 2 femtoseconds for 20 ps, keeping a constraint of 5 kcal/mol. During this phase, the temperature was kept constant through the use of Langevin dynamics [41] . The SHAKE algorithm [42] was used to constrain the geometry of bonds that involved hydrogen atoms. The systems were then equilibrated using a two-femtosecond time step. After this, each system was treated at constant pressure and temperature for 1 ns, at 1 bar pressure and 300K temperature. Another round of equilibration was performed for 1 ns. The production run was performed for 200 ns. The CPPTRAJ module [43] was applied to perform an analysis of the simulation trajectories in order to evaluate the stability of the system structures. 

The radial distribution function (RDF) is denoted by g(r) and was applied after molecular dynamics simulation to illustrate the variation in the density of interatomic interactions with respect to time [44] . RDF was performed through the CPPTRAJ module of AMBER considering only hydrogen bonds between the compounds and MvfR residues. The hydrogen bond interactions were examined using an in-house Visual Molecular Dynamics (VMD) Perl script [45] .

Complexes were then subjected to the AMBER MMPBSA.py method to calculate the net binding free energies of the complex, enzymes and compounds [46] . The net system binding free energy was considered by subtracting the combined binding energy of the enzyme and compound from the complex binding energy. From simulation trajectories, 1000 frames were picked at a regular time interval and investigated for molecular mechanical energies and solvation free energies. The binding free energy estimation was conducted through two methods: MM-PBSA and MM-GBSA [47, 48] . Statistically, both methods estimate binding free energy as, ∆G net binding free energy = ∆G binding free energy of complex − ∆G receptor + ∆G ligand

The AMBER NMODE module was employed to compute the contribution of entropy to the net binding MM-PB/GBASA energy of the complexes [49] . Only 10 frames of the trajectories were analyzed.

Further confirmation on the intermolecular stability of the MvfR-compound complexes was achieved by estimating the absolute binding free energies using WaterSwap from the Sire Package [50, 51] . WaterSwap works on the idea of swapping ligand dimensions at the active pocket of the enzyme with water molecules of equal volume and size from the explicit environment. One thousand iterations were completed for each system, which is reported to be enough to obtain converged binding energy values. The absolute binding free energy was calculated using four efficient methods: thermodynamic integration (TI), free energy perturbation (FEP), quadrature-based integration of TI, and Bennett's acceptance ratio (BAR) method. The stability of complexes was considered to be high if they had net energy values <1 kcal/mol [52] .

The SwissADME [53] and pkCSM [54] online software applications were used to predict the ADMET properties of the lead molecules.

SBVS screening studies were performed at the active pocket of the MvfR enzyme using two highly reliable docking software packages: AutoDock Vina and GOLD. The results of both virtual screenings were then sorted on the basis of the docking function, and top two leads (in comparison with the M64 control) were selected. Top-10 hits that were virtually screened and had higher binding affinity binders for the targeted MvfR are tabulated in Table 1 . The top two hit compounds: (4-((carbamoyloxy)methyl)-10,10dihydroxy-2,6-diiminiodecahydropyrrolo[1,2-c]purin-9-yl sulfate) and 10,10-dihydroxy-2,6diiminio-4-(((sulfonatocarbamoyl)oxy)methyl)decahydropyrrolo[1,2-c]purin-9-yl sulfate were consistently observed to have excellent binding; therefore, they were considered for additional computational analysis. Selection of these compounds in complex with MvfR was conducted considering their structural stability in molecular dynamics simulations and MM-GBSA analysis. The 50-ns MD simulation validated the stability of the docked conformation of the lead compound with MvfR, and the structural deviations of the Cα atoms were plotted against time as the RMSD (presented in Section 3.3). Contrary to the control M64, the screened lead-MvfR complexes were stable and followed a somewhat similar RMSD trend within 2 Å. Small numbers of frames from the molecular dynamics simulations were then selected to estimate the MM-GBSA binding free energy for further validation of the binding strength of the lead molecules for MvfR. The estimated net MM-GBSA binding free energies of Top-1, Top-2 and control M64 were −24.15 kcal/mol, −45.47 kcal/mol and −68.89 kcal/mol, respectively. 

The AutoDock binding free energy score and GOLD fitness score of the control were −8.14 kcal/mol and 55.14, respectively. The majority of the interactions produced by M64 with MvfR residues were hydrophobic and were similar to those reported previously. In our study, we found that the M64 formed a closed-distance hydrogen bond with Gln102 through the central chemical moiety oxygen atom. This finding is in line with the cocrystallized structure, which is a strong indicator of the soundness of the docking methodology applied herein. The rest of the M64 structure, including Ala10, Ile57, Thr74, Hie92, Ser93, Ile94, Asn114, Arg117, Pro118, Phe129, Trp142, Ala145, Pro146, Leu162, Ser163, and Thr173, was strongly entangled by van der Waals residues. Along with this, several pi-alkyl, pi-sigma and pi-pi stack interactions could be noticed between M64 and MvfR ( Figure 2A ). M64 is a highly competitive antagonist and has an enhanced affinity for MvfR than natural binding substrates. The Top-1 lead, in contrast to M64, produced a higher number of hydrogen bonds than van der Waals and other chemical contacts. The AutoDock binding energy of the compound was −9.18 kcal/mol and its GOLD fitness value was 61.4. The 10,10-dihydroxy-2,6-diiminiodecahydropyrrolo[1,2-c]purin-9-yl 52.2 −6.78

The AutoDock binding free energy score and GOLD fitness score of the control were −8.14 kcal/mol and 55.14, respectively. The majority of the interactions produced by M64 with MvfR residues were hydrophobic and were similar to those reported previously. In our study, we found that the M64 formed a closed-distance hydrogen bond with Gln102 through the central chemical moiety oxygen atom. This finding is in line with the co-crystallized structure, which is a strong indicator of the soundness of the docking methodology applied herein. The rest of the M64 structure, including Ala10, Ile57, Thr74, Hie92, Ser93, Ile94, Asn114, Arg117, Pro118, Phe129, Trp142, Ala145, Pro146, Leu162, Ser163, and Thr173, was strongly entangled by van der Waals residues. Along with this, several pi-alkyl, pi-sigma and pi-pi stack interactions could be noticed between M64 and MvfR ( Figure 2A ). M64 is a highly competitive antagonist and has an enhanced affinity for MvfR than natural binding substrates. The Top-1 lead, in contrast to M64, produced a higher number of hydrogen bonds than van der Waals and other chemical contacts. The AutoDock binding energy of the compound was −9.18 kcal/mol and its GOLD fitness value was 61.4. The 10,10-dihydroxy-2,6-diiminiodecahydropyrrolo[1,2-c]purin-9-yl sulfate region, in particular, dominated the bindings and favored the compounds that showed hydrogen bonding with active pocket residues such as Gln102, Ser104, Asn114, and Arg117. The opposite carbamic acid chemical moiety of the compound favored the formation of hydrogen bonds with Asp172. The van der Waals interactions of this compound involved the residues Ile57, Ile103, Leu116, Pro118, Val119, Trp142, Gly143, Ile144, Ile171, and Thr173 ( Figure 2B ). The Top-2 compound had a GOLD score of 59.2 and a binding energy value of -9.0 kcal/mol. Similarly to Top-1, Top-2 dominated its interactions with the MvfR through hydrogen bonds. The compound formed hydrogen bonds with Gln102, Leu105. Asn114, Leu115, Arg117, and Val119. The rest of the interactions can be seen in Figure 2C . Overall, the docking study predicted that the lead compounds and the control would favor binding to the same ligand-binding domain and produce interaction networks of the same kind. 

The initial 50 ns of the molecular dynamics simulation was extended to 200 ns to obtain confidence values of the complexes' stability and the strength of leads' interactions with the MvfR. Different analyses conducted on the simulation trajectories of the leads and control trajectories are presented in Figure 3 . Molecular dynamic simulation is a highly useful technique for determining the time-dependent stability of the ligand-receptor interactions and docked conformation. Different statistical analyses based on Cα atoms were performed, such as the radius of gyration (Rg) [55] , the root mean square deviation (RMSD) [56] and root mean square fluctuations (RMSF) [57] . Analyses were started by calculating the RMSD for all three complexes by superimposing all frames of the molecular dynamics simulation over the initial docked conformation, and a plot was generated for the entire simulation time using the XMGRACE software [58] . From Figure 3A , consistent stability of the systems could be witnessed and the systems' RMSD were close to the 2 Å mark throughout the simulation time. This indicates that the leads were enjoying the affinity inside the active pocket of the MvfR, similarly to that of the co-crystalized M64, and were strongly attached to the active side residues via hydrophobic and hydrophilic interactions. The different interactions that allowed stable binding of the compounds at the enzyme's active pocket are mentioned in Table 2 and the Radial Distribution Function section. For compound 1, Leu71, Tyr73, Arg197, and Leu200 were reported to be in consistent contact with the ligands, while Ser104, Leu115, Arg117, Ser163, Gln190, and Ile194 were reported for compound 2. Due to the stable binding conformation of the compounds during the simulation time, the interacting residue pattern remain mostly the same and did not experiences any major shifts. Next, RMSF analysis was performed on the Cα atoms, visualizations of which are shown in Figure 3B . All the systems' residues revealed highly stable RMSF values and were within the range of <2 Å with the exception of a few 3 Å spikes. It was observed that the loop regions of MvfR represented high fluctuations in the presence of compounds/control; however, they were still in an acceptable range and allowed proper accommodation of the leads/control inside the pocket for enhanced docked stability. Lastly, the MvfR's compactness and structure equilibrium were tested in the presence of the leads/control by means of Rg analysis of the Cα atoms. As can be noticed in Figure 3C , all three systems were in good equilibrium (<19 Å). The lead complexes, in particular, were more stable than the control complex. In summary, the simulation analysis confirmed the stability of the systems and identified the compounds as suitable candidates for further experiments in order to unravel their real affinity for the MvfR. 

Acceptor 

In order to find out the occupancy of hydrogen bonding between the ligand molecules and the protein, we performed hydrogen bond analysis [59] . These interactions determined the intermolecular specificity and were important to stabilize the ligand-protein complexes. Table 2 . The control, Top-1 and Top-2 leads were found to form 12, 68, and 28 hydrogen bonds, respectively, with the MvfR. Several key residues already predicted by molecular docking studies were unveiled to play crucial roles in ligand binding throughout the length of the simulation time. Several previous studies reported the importance of hydrogen bonds while designing new drug molecules against a given biological target [60, 61] . For example, Khalid et al. [24] demonstrated several key residues of soluble guanylate cyclase H-NOX domain with ligand molecules. 

Furthermore, the RDF analysis was conducted using strong intermolecular interactions between MvfR and the compounds to understand the intensity of interactions versus time. RDF has been frequently employed in studies to highlight the critical intermolecular interactions that are key in the recognition and binding of good affinity binders [9, 57, 62] . Several residues were filtered that favored continuous contacts with the compounds throughout the simulation time ( Figure 5 ). These interactions were plotted in terms of density versus distance. In the case of Top-1, residues such as Leu71, Tyr73, Arg197, and Leu200 were among the high-density interactions with MvfR, while, in the case of Top-2, Ser104, Leu115, Arg117, Ser163, Gln190, and Ile194 were among the high-density residues that were in consistent interactions. Interactions that remained constant after specific time periods are not provided while those of bond distance variations are plotted. compounds throughout the simulation time ( Figure 5 ). These interactions were plotted in terms of density versus distance. In the case of Top-1, residues such as Leu71, Tyr73, Arg197, and Leu200 were among the high-density interactions with MvfR, while, in the case of Top-2, Ser104, Leu115, Arg117, Ser163, Gln190, and Ile194 were among the highdensity residues that were in consistent interactions. Interactions that remained constant after specific time periods are not provided while those of bond distance variations are plotted. 

The estimation of binding free energy via the MM-PBSA and MM-GBSA gives reliable predictions about a compound's affinity for a given biological macromolecule [47, 48] . Both of the mentioned techniques are widely employed in drug design processes because of their low computational demands and the fact that the results can be easily correlated to experimental results [63] . Both of the methods have been regularly employed to validate molecular dynamics simulations and docking predictions [64, 65, 66] . Complete results of the binding free energies of the complexes are tabulated in Table 3 . For both the lead and control complexes, electrostatic and van der Waals energy were revealed to play a critical role in binding, as predicted by the docking findings. The control was found to show higher van der Waals domination in net interactions compared to the lead molecules, which exhibited higher electrostatic contributions as well equal contributions from van der Waals energy. The solvation energy in the case of the control was found to be more dominated by non-polar energy, whereas, for the leads, the polar solvation energy was three times more stabilized than the non-polar solvation energy. Overall, the net binding energies of the systems were very high, thus demonstrating the stability of the systems. The net binding free energies of the systems are in following order: control 

The estimation of binding free energy via the MM-PBSA and MM-GBSA gives reliable predictions about a compound's affinity for a given biological macromolecule [47, 48] . Both of the mentioned techniques are widely employed in drug design processes because of their low computational demands and the fact that the results can be easily correlated to experimental results [63] . Both of the methods have been regularly employed to validate molecular dynamics simulations and docking predictions [64] [65] [66] . Complete results of the binding free energies of the complexes are tabulated in Table 3 . For both the lead and control complexes, electrostatic and van der Waals energy were revealed to play a critical role in binding, as predicted by the docking findings. The control was found to show higher van der Waals domination in net interactions compared to the lead molecules, which exhibited higher electrostatic contributions as well equal contributions from van der Waals energy. The solvation energy in the case of the control was found to be more dominated by non-polar energy, whereas, for the leads, the polar solvation energy was three times more stabilized than the non-polar solvation energy. Overall, the net binding energies of the systems were very high, thus demonstrating the stability of the systems. 1 kcal/mol) ). For Top-1, electrostatic energy was found to be predominant in the docking studies as the hydrogen bond distances were very close, mostly within 3 Å. The hydrophobic interactions were also reported to equally favor the stable binding of compounds and played a key role, along with hydrophilic interactions, in holding the compound conformation at the active pocket. Both MM/GBSA and MM/PBSA agree on the significant electrostatic contribution and equally favorable binding of the van der Waals energy. In dynamics simulation analysis, the same findings were revealed, where both the hydrophilic and hydrophobic interactions (as shown in Table 2 and RDF analysis) remained critical for the stability of the RMSD plot. In case of Top-2, the van der Waals energy dominated over the electrostatic energy by a very low margin; the same was observed in the docking analysis. The van der Waals and other hydrophobic interactions pushed the more electronegative chemical moieties of the compound towards the inside of the pocket. This resulted in good interaction networks of both the electrostatic and van der Waal contacts. The binding conformation stabilities and binding interaction profiles of the compounds with the enzyme remained consistent in all of the analyses performed in this study, all of which classified the compounds as strong binders of MvfR. 

Further analysis was conducted to determine the key hotspot residues of MvfR that contributed significantly in terms of binding and holding the leads/control at the active pocket. Identification of hotspot residues was performed in many previous studies to report key interactions between ligands and residues that were vital in stabilizing the ligands at the docked site [57, 67] . The net MM-GBSA binding energies of the systems were decomposed into residues of the MvfR, and only the common residues that were critical in binding the ligands were shortlisted, as shown in Table 4 . Gln102, Asn114, Arg117 and Val199 were common in all complexes and were found to be major contributors to the ligand interactions. Gln102 was a key hydrogen bonding residue and was reported previously in hydrogen-bonding interactions with ligand leads. It was observed that the rest of the residues involved both hydrogen bonding as well as van der Waals interactions. 

To compensate for the missing approximation of binding entropy in MM-PBSA and MM-GBSA, the entropy calculation was implemented via normal mode in the AMBER package. As the calculation was very slow, only a limited number of frames were analyzed. The net entropy of the systems was in the following order: control (−8.89 kcal/mol), Top-1 (−10.10 kcal/mol) and Top-2 (−11.00 kcal/mol).

Although the MM-PBSA and MM-GBSA methods are very successful in determining free energies, they have several limitations; therefore, another validation method, Wa-terSwap, was applied in the study. The WaterSwap-based binding free energy values, calculated using different algorithms, are illustrated in Figure 6 . Both of the lead molecules were disclosed as better binders than control M64. As can be seen, the net WaterSwap energies calculated the using algorithms for all three systems differed by no more than 1 kcal/mol, which demonstrated highly converged systems. 

Unfavorable pharmacokinetics of compounds in the process of drug discovery can lead to drug failure, and thus, can increase the time and cost involved in the development of potent and safe drugs [53, 54] . For this purpose, pharmacokinetics predictions are important in the early stages of drug development using available in silico tools to enhance the chances of selecting the correct molecules for development. Medicinal chemistry focuses on drug absorption, and this was evaluated as the first step in these in silico studies. It was observed that both compounds were highly water soluble, as predicted by the ESOL, Ali and SILICOS-IT methods in the SWISSADME server. For this reason, the compounds are excellent candidates in terms of oral bioavailability. Further, the compounds had no Pan-assay interference compounds (PAINS) structure; thus, they targeted one specific biological target and had one desired effect [68] . From a synthetic chemistry perspective, the compounds had a good synthetic accessibility score of~5, meaning they will be easy to synthesize in future experimental analyses. The compounds also had high gastrointestinal absorption and did not act as substrates for the P-glycoprotein transporters. The transdermal deliveries of the compounds are also predicted to be very good, making them suitable for skin-related products. They had volume distribution values that indicated their low tissue distribution as compared to their distribution in the plasma. Likewise, they also had low fraction unbound values, which indicate that they could lower their serum protein binding affinities and could enhance their distribution efficiency through the cell membranes. The blood-brain barrier crossing abilities of drugs are important in terms of evaluating their side effects and toxicity, as well as the efficiency of their pharmacological action in the brain [69] . These compounds had poor blood-brain barrier penetration, and thus, they could not move through the central nervous system easily. Additionally, they did not inhibit the detoxification of cytochrome P450, and thus, were involved in the oxidation of xenobiotics to help in their removal. The renal and hepatic clearance of the compounds were projected to be~0.53 log mL/min/kg. This total clearance of compounds is an important factor in terms of evaluating their bioavailability and calculating the rate of dosage for their steady-state concentration. They were found to be AMES non-toxic based on their LD50 values during oral administration to rats, and were anticipated to demonstrate no sensitization of the skin and to not inhibit hERGI and hERGII, which can reduce the likelihood of QT syndrome development. Detailed pharmacokinetic data of both lead molecules are tabulated in Table 5 . 

Computer aided drug design is an integral part of modern drug discovery and has played a significant role in the development of drugs that are in clinical trials and in clinical use [17] . In this study, we identified two leads against MvfR of P. aeruginosa, which, as with the co-crystallized M64, had high binding affinity for the relevant enzyme. Additionally, the compounds met prominent druglike rules and had good pharmacokinetics and acceptable safety. The docking and subsequent molecular dynamics simulations revealed the formation of strong interactions by the compounds with active pocket residues and significant conformation stability. There were multiple van der Waal and hydrogen bond interactions of the compounds with the hotspot residues of the pocket, which resulted in increased intermolecular affinity. In particular, strong van der Waals interactions and hydrogen bonding was observed between all the screened molecules and MvfR residues (Gln102, Asn114, Arg117, Val119 and Asp172). These residues were present in all the crystal structures of the MvfR in P. aeruginosa and formed bonds with co-crystallized ligands [16] . These residues are conserved among the different MvfR and are considered to be key for enzyme functionality [16] . Thus, there are fewer chances for the enzyme to escape the compound's action. In addition, the compounds revealed favorable druglike and lead-like properties and were reported to have good pharmacokinetic profiles. This also increases the chances of the compounds reaching the market. Since these compounds showed promising results and are easily available from commercial sources, they can be used in further quick in vivo experiments to determine their real binding affinity and MvfR inhibition potential. Moreover, it is suggested that the anti-P. aeruginosa activity of the compounds be investigated to affirm that the molecules do not infer with the essential biological pathways of the pathogens and do not harm bacterial cells. This will ensure that the compounds will only de-weaponize the pathogen and will not mediate the evolution of new resistant mechanisms. To summarize, these screened compounds provide promising starting structures in the search for novel anti-virulent compounds against the superbug P. aeruginosa. Although very stringent criteria for docking score functions, as well as for the dynamics simulations and binding free energy methods, were applied for the compounds in this study, the shortlisting and selection of stable conformations can still be undertaken with confidence, and the selection of the ligand conformation can be enhanced via such approaches as using virtual screening rather than docking-based screening in long-length traditional molecular dynamics simulations [70, 71] , thoroughly analyzing the binding of free energies throughout the frames of simulations, using the hybrid QM/MM approach for further enhancing the prediction of stable conformation accuracy [72] , using the Selective Ligand-Induced Conformational Ensemble (SLICS) method [73] , using different analytical methods such as axial frequency distribution (AFD) [74] and the Binding Free Energy-Based Footprint Pharmacophore Model method [75] , etc. 

The data presented in this study are available within the article.

Access to effective antimicrobials: A worldwide challenge

Tahir ul Qamar, M. Integrated Core Proteomics, Subtractive Proteomics, and Immunoinformatics Investigation to Unveil a Potential Multi-Epitope Vaccine against Schistosomiasis. Vaccines 2021

The global burden of chronic diseases: Overcoming impediments to prevention and control

Development of a Novel Multi-Epitope Vaccine Against Crimean-Congo Hemorrhagic Fever Virus: An Integrated Reverse Vaccinology

Microbial resistance movements: An overview of global public health threats posed by antimicrobial resistance, and how best to counter. Front

Tahir ul Qamar, M. Pan-Vaccinomics Approach Towards a Universal Vaccine Candidate Against WHO Priority Pathogens to Address Growing Global Antibiotic Resistance

Designing multi-epitope vaccine against Staphylococcus aureus by employing subtractive proteomics, reverse vaccinology and immuno-informatics approaches

Mechanisms of antimicrobial resistance in ESKAPE pathogens

Comparative subtractive proteomics based ranking for antibiotic targets against the dirtiest superbug: Acinetobacter baumannii

Prediction of vaccine candidates against Pseudomonas aeruginosa: An integrated genomics and proteomics approach

Growing emergence of drug-resistant Pseudomonas aeruginosa and attenuation of its virulence using quorum sensing inhibitors: A critical review

Infections caused by resistant gram-negative bacteria: Epidemiology and management

Antibiotic resistance in Pseudomonas aeruginosa: Mechanisms and alternative therapeutic strategies

Tolerance and resistance of Pseudomonas aeruginosa biofilms to antimicrobial agents-How P. aeruginosa can escape antibiotics

Assessing Pseudomonas aeruginosa persister/antibiotic tolerant cells

Molecular insights into function and competitive inhibition of Pseudomonas aeruginosa multiple virulence factor regulator

Computer-aided drug design methods

Mutational Landscape of Pirin and Elucidation of the Impact of Most Detrimental Missense Variants That Accelerate the Breast Cancer Pathways: A

Discovery of anti-MERS-CoV small covalent inhibitors through pharmacophore modeling, covalent docking and molecular dynamics simulation

Ullah Saeed, H.F. Immuno-Informatics Analysis of Pakistan-Based HCV Subtype-3a for Chimeric Polypeptide Vaccine Design

Protein Data Bank (PDB): Database of threedimensional structural information of biological macromolecules

UCSF Chimera-A visualization system for exploratory research and analysis

The FF14SB force field

Comparative Studies of the Dynamics Effects of BAY60-2770 and BAY58-2667 Binding with Human and Bacterial H-NOX Domains

CMNPD: A comprehensive marine natural products database towards facilitating drug discovery from the ocean

Free ADME-tox filtering computations for chemical biology and early stages drug discovery

Lead-and drug-like compounds: The rule-of-five revolution

Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings

Small-molecule library screening by docking with PyRx

Merck Molecular Force Field

Using AutoDock 4 and AutoDock Vina with AutoDockTools: A Tutorial; The Scripps Research Institute, Molecular Graphics Laboratory

Improved protein-ligand docking using GOLD

Discovery of selective inhibitors for cyclic AMP response element-binding protein: A combined ligand and structure-based resources pipeline

Discovery of human coronaviruses pan-papain-like protease inhibitors using computational approaches

Investigating the molecular mechanism of staphylococcal DNA gyrase inhibitors: A combined ligand-based and structure-based resources pipeline

Promising terpenes as SARS-CoV-2 spike receptorbinding domain (RBD) attachment inhibitors to the human ACE2 receptor: Integrated computational approach

Antechamber: An accessory software package for molecular mechanical calculations

Development and testing of a general amber force field

Imtiaz-Ud-Din

A one-pot multicomponent facile synthesis of dihydropyrimidin-2(1: H)-thione derivatives using triphenylgermane as a catalyst and its binding pattern validation

Constant pressure molecular dynamics simulation: The Langevin piston method

A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations

Software for processing and analysis of molecular dynamics trajectory data

Radial Distribution Functions of Some Structures of the Polypeptide Chain

VMD: Visual molecular dynamics

py: An efficient program for end-state free energy calculations

The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities

Assessing the Performance of the MM_PBSA and MM_GBSA Methods. 1. The Accuracy of Binding Free Energy Calculations Based on Molecular Dynamics Simulations

The normal-mode entropy in the MM/GBSA method: Effect of system truncation, buffer region, and dielectric constant

Rapid decomposition and visualisation of protein-ligand binding free energies by residue and by water

Computational prediction of drug solubility in water-based systems: Qualitative and quantitative approaches used in the current drug discovery and development setting

Screening pipeline for Flavivirus based inhibitors for Zika virus NS1

SwissADME: A free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules

Predicting small-molecule pharmacokinetic and toxicity properties using graph-based signatures

Radius of gyration as an indicator of protein structure compactness

Significance of Root-Mean-Square Deviation in Comparing Three-Dimensional Structures of Globular Proteins

Binding mode analysis, dynamic simulation and binding free energy calculations of the MurF ligase from Acinetobacter baumannii

The role of hydrogen-bonds in drug binding

Hydrogen bond analysis of the EGFR-ErbB3 heterodimer related to non-small cell lung cancer and drug resistance

Toward novel inhibitors against KdsB: A highly specific and selective broad-spectrum bacterial enzyme

Interaction mechanisms of a melatonergic inhibitor in the melatonin synthesis pathway

What is the current Value of MM/PBSA and MM/GBSA Methods in Drug Discovery

Structural insight into the binding pattern and interaction mechanism of chemotherapeutic agents with Sorcin by docking and molecular dynamic simulation

Synthesis and Identification of Novel Potential Molecules Against COVID-19 Main Protease Through Structure-Guided Virtual Screening Approach

Abrogation of SARS-CoV-2 interaction with host (NRP1) Neuropilin-1 receptor through high-affinity marine natural compounds to curtail the infectivity: A structural-dynamics data

Binding free energy based analysis of arsenic (+3 oxidation state) methyltransferase with Sadenosylmethionine

Growing PAINS in academic drug discovery

The blood-brain barrier: Clinical implications for drug delivery to the brain

Structure-based virtual screening to get new scaffold inhibitors of the Ser/Thr protein kinase PknB from mycobacterium tuberculosis

Structurally simple synthetic 1, 4-disubstituted piperidines with high selectivity for resistant Plasmodium falciparum

In silico methods and tools for drug discovery

Accelerated Structural Prediction of Flexible Protein-Ligand Complexes: The SLICE Method

AFD: An application for bi-molecular interaction using axial frequency distribution

Binding free energy-based footprint pharmacophore model to enhance virtual screening and drug discovery: A case on glycosidases as anti-influenza drug targets

The authors would like to acknowledge Jilin Agricultural University, Ministry of Education, Jilin, China, for providing the facilities for this research.

The authors declare no conflict of interest.