VIRTUAL SCREENING OF FDA-APPROVED DRUGS BY MOLECULAR DOCKING AND DYNAMICS SIMULATION TO RECOGNIZE POTENTIAL INHIBITORS AGAINST MYCOBACTERIUM TUBERCULOSIS ENOYL-ACYL CARRIER PROTEIN REDUCTASE ENZYME

Objective: This in silico study is aimed at identification of new possible inhibitors against Mycobacterium tuberculosis InhA enzyme by screening a library of FDA-approved drugs. Methods: In this in silico study, a library of FDA-approved drugs was screened by molecular docking against the monomer of enoyl-acyl carrier protein reductase to recognize potential inhibitors. Then, those best drugs with minimum docking energy were subjected to molecular dynamics simulation. Results: Out of the top ten docking hits, only revefenacin was able to maintain the closet proximity to InhA enzyme binding pocket during the two rounds of dynamics simulation. Analysis of molecular dynamics (MD) simulation data indicated that the antimuscarinic drug revefenacin has a ligand movement Root-Mean-Square Deviation (RMSD) that didn’t exceed 4 Angstrom. Also, in this MD study, revefenacin has a superior binding energy of -35.59 Kcal/mol as compared to -13.88 Kcal/mol for the other hit ergotamine. These favorable MD simulation records for revefenacin can be explained by its ability to continuously interact with enzyme binding pocket by two hydrogen bonds. Conclusion: We report that the antimuscarinic drug revefenacin may have the potential to inhibit the enoyl-acyl carrier protein reductase for Mycobacterium tuberculosis . However, these preliminary results must be further evaluated by in vitro and in vivo studies.


INTRODUCTION
Tuberculosis (TB) is one of the earliest recognized infective diseases in human history that is caused by the aerobic bacillus bacterium called Mycobacterium tuberculosis [1].This acid-fast pathogen can be transmitted by the inhalation of respiratory droplets that usually generated by the coughing of active TB patients [2].Once the bacilli of Mycobacterium tuberculosis reach the respiratory alveoli, it will then get phagocytized by the resident alveolar macrophages.Some of these engulfed bacilli will continue to multiply inside macrophages, leading to the lysis of these innate immune cells.If the alveolar macrophages failed to eliminate Mycobacterium tuberculosis infection within one to two weeks, then local inflammation will recruit monocytes to the site of infection with subsequent establishment of specific T-cells mediated immune response.When immune cells are unable to eradicate infection by Mycobacterium tuberculosis, then the best strategy to prevent further spread of infection is by the formation of granuloma around mycobacterial bacilli [3].It is estimated that in about 10% of Mycobacterium tuberculosis infection cases, immune cells will succuss to eliminate infection.However, in most infection cases, some mycobacterial bacilli will escape immune response and enter a latent and nonreplicating stage in lesion areas of the lung, leading to a latent tuberculosis infection (LTBI).These dormant bacilli can regain their ability to multiply and spread when immune system is disrupted, leading to active tuberculosis infection [4].According to the estimates of World Health Organization (WHO), about one-quarter of world population is latently infected with Mycobacterium tuberculosis and 5-10 % of them will later develop active TB (https://www.who.int/news-room/fact-sheets/detail/tuberculosis)[5,6].In active TB cases, the immune cells will fail to contain mycobacterial infection and bacilli can spread through blood circulation and lymphatic channels into other distal parts of the body.This can result in a serious damage in lung and other parts of the body like the brain, spinal cord and bones.Usually, individuals with latent TB are asymptomatic and can't transmit infection.On the contrary, actively infected patients can spread the infection to others and experience clinical symptoms like continuous coughing, bloody sputum, fatigue, fever and chest pain [7].For the effective management of TB cases, rapid and reliable diagnostic methods must be utilized.The frequently used methods for diagnosis of TB cases are clinical diagnosis, chest X-ray, microscopic identification, bacterial culture and immunologic tests like skin test [8].Also, early and effective treatment is important in the management of TB to avoid the emergence of drug-resistant strains of Mycobacterium tuberculosis.The standard treatment regimen recommended by WHO for active pulmonary TB involves the administration of four antibiotics in combination for two months and these are: isoniazid (INH), rifampicin, ethambutol and pyrazinamide.After that, a combination of only isoniazid and rifampicin is used for another four months [9].Unfortunately, this long and complex TB treatment regimen has led to issues like patient non-compliance with subsequent treatment failure [10].The failure of treatment can lead to TB relapse and emergence of either multidrug-resistant tuberculosis (MDR-TB) or extensively drug-resistant tuberculosis (XDR-TB) [11].Moreover, the only available tuberculosis vaccine is the Bacillus of Calmette and Guérin (BCG) vaccine, whose efficacy is controversial and can range between 0-80 % according to clinical trials [12].Therefore, it is relevant to introduce new therapeutic agents that can reduce the current burden of public health threat imposed by TB.In this direction, one of the potential molecular targets to design new anti-tuberculosis agents is the mycolic acids (MAs) synthesis pathway.The mycolic acids are long chain fatty acids that are essential for the buildup and integrity of mycobacterium cell wall.These fatty acids can act as a permeability barrier and contribute to bacterial resistance against both host defense mechanisms and different antibiotics [13,14].In Mycobacterium tuberculosis, one of the key enzymes for the reduction stage of fatty acids production and biosynthesis of mycolic acids is the enoyl-acyl carrier protein (ACP) reductase (InhA).The InhA is an NADH-dependent reductase and is a well-known target for isoniazid (INH) [15].As such, we have chosen in this in silico study the enoyl-ACP reductase of Mycobacterium tuberculosis as a molecular target to screen a library of FDA-approved drugs.This virtual screening study is aimed at the identification of new potential inhibitors against Mycobacterium tuberculosis InhA enzyme.

Setting up a virtual screening study plan
The main steps and options applied in this in silico study are outlined in (fig.1).As noted in this figure, both molecular docking and two rounds of dynamics simulation were employed to virtually screen a library of FDA approved drugs against Mycobacterium tuberculosis InhA enzyme.

Molecular docking
For this computational study, we have used the online drug discovery web server named DrugRep to virtually screen a library of 2,315 FDA-approved drugs against chain A of enoyl-ACP reductase for Mycobacterium tuberculosis [16].Initially, chain A only of the enoyl-ACP reductase (PDB: 1BVR) was extracted by using UCSF Chimera version 1.15 [17,18].Then, this extracted monomer of InhA enzyme was uploaded to the DrugRep server to perform a structure-based virtual screening.It is worthwhile to mention that the DrugRep website utilizes both AutoDockTools version 1.5.6 and AutoDock Vina version 1.1.2to prepare the uploaded target and perform the docking step, respectively [19,20].In this study, the used docking coordinates were X: 17.0, Y: 14.0, Z: 8.0 while the dimensions of grid box were 22*22*22 Angstrom.Moreover, the accuracy of this docking protocol was assessed by redocking the co-crystalized NAD + into chain A of enoyl-ACP reductase.Then, PacDOCK web server was used to compare the conformation of co-crystalized and docked NAD + by calculating the root mean square deviation (RMSD) [21].Finally, we have selected the best 10 drugs with least docking energy for more assessment by molecular dynamics (MD) simulation.Additionally, the orientation of least binding energy pose for each drug-enzyme docking complex was visualized by using PyMOL version 2.4.1 (https://pymol.org/2/)and protein-ligand interaction profiler (PLIP) [22].It should be noted that during selecting the top 10 hits, we have removed anticancer, antihypertensive and psychoactive drugs from the results of this virtual screening as the use of these agents can be associated with unacceptable adverse effects.Thus, the repurpose of these drugs may be clinically inappropriate for tuberculosis patients.

Molecular dynamics (MD) simulation study
After molecular docking, YASARA Dynamics version 20.12.24 was utilized to carry out two rounds of molecular dynamics simulation for 20 and 50 nanoseconds intervals [23].In the first round of MD simulation, the drug-enzyme docking complex with least energy binding pose for each of the top 10 hits was submitted to undergo 20 nanoseconds simulation.After that, only those drugs with maximum proximity RMSD to enoyl-ACP reductase binding pocket of less than 4 Angstrom were submitted to the second round of 50 nanoseconds simulation.Again, the ligand maximum proximity RMSD to the enzyme active site was calculated to validate results of MD simulation second round.In addition, the binding energy of molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) was computed for each drug by using AMBER14 force field [24].For this MD study, the detailed procedure and applied options are identical to what we have used in our previously published studies [25][26][27].Concisely, NaCl was employed in this simulation with a concentration of 0.9%.Also, an excess concentration of either Na + or Cl ˉ ions was applied to the simulation environment to ensure the neutralization of the drugenzyme complex.This MD simulation was carried out with the help of the following force fields: AM1BCC and GAFF2 for the ligand, AMBER14 for the solute, TIP3P for water [24,28,29].Moreover, any probability of clashes in this simulation was prohibited through the reduction of the steepest descent and simulated annealing.

RESULTS AND DISCUSSION
Before starting the virtual screening study, the precision of the followed docking protocol was assessed first by using the redocking method.In this method, the co-crystalized NAD + was removed from chain A of InhA enzyme and docked again into the same binding pocket by using similar procedure and options to those applied in the virtual screening project.Then, the native conformation of NAD + was aligned with the docked conformation and RMSD value was calculated to estimate the variation degree between these two conformations.It is well-known that a lower RMSD value for conformations variation reflects good accuracy of docking protocol.According to literatures review, a precise docking protocol is usually coincided with a conformations difference RMSD between 1.5 and 2.0 Angstrom [30].In this study, the alignment of NAD + native conformation with the docked one can be observed in (fig.2).As seen in this figure, the computed RMSD value for conformations difference was 1.58 Angstrom.Thus, the applied docking protocol for this study is deemed to be accurate.It should be pointed out that the docking energy of NAD + into chain A of enoyl-ACP reductase was -9.4 Kcal/mol.An overview for the results of molecular docking and dynamics simulation study for the best 10 drugs can be observed in (table 1).Based on findings in this table, the docking energy for these drugs is ranging between -12.1 and-9.6Kcal/mol.interestingly, the first three hits in (table 1) with least binding energy are anti-migraine medications.Also, two of the hits listed in (table 1) are penicillin antibiotics and these are carindacillin and sultamicillin.According to several published studies, many of the listed drugs in (table 1) may have activity against viruses like SARS-CoV-2 and West Nile virus [31][32][33][34][35][36].Moreover, based on an in vitro study, the anti-asthma drug zafirlukast may have antituberculosis activity.In this study, zafirlukast was able to inhibit the growth of Mycobacterium tuberculosis by causing a dysregulation in the expression of several genes inside bacteria [37].However, the maximum ligand movement RMSD in the first round of MD simulation in (table 1) indicated that only ergotamine and revefenacin were able to stay close to binding site with RMSD value less than 4 Angstrom throughout 20 nanoseconds interval.It is well-known that a low ligand movement RMSD means a close proximity of the drug to target site and, thus stronger binding [38].Consequently, the second round of MD simulation for 50 nanoseconds were executed only for ergotamine and revefenacin.During this later simulation, only the antimuscarinic drug revefenacin was able to maintain a ligand movement RMSD value that didn't exceed 4 Angstrom as seen in the last column in (table 1).A detailed plot of ligand movement RMSD throughout 50 nanoseconds interval of MD simulation can be viewed in (fig.3) for both ergotamine and revefenacin.Based on this plot, it is expected that revefenacin was able to stay closer than ergotamine to the binding site of InhA enzyme monomer.As a result, revefenacin may exhibit a stronger binding potential to enzyme active site because it has a lower ligand movement RMSD values during MD simulation.Moreover, the computed average MM-PBSA binding energy during 50 nanoseconds simulation was -13.88 and-35.59Kcal/mol for ergotamine and revefenacin, respectively.

Fig. 3: A plot of ligand movement RMSD versus molecular dynamics simulation interval for both ergotamine and reverencing
This difference in ligand movement RMSD and MM-PBSA binding energy during MD simulation between ergotamine and revefenacin can be explained through examination of types of interactions between drug and residues of enzyme binding pocket.The visualization of docking complexes in (fig.4), indicated that both ergotamine and revefenacin can form two hydrogen bonds with residues of enzyme binding site.
However, only revefenacin was able to keep these two hydrogen bonds with enzyme active site throughout 50 nanoseconds of MD simulation as can be noted in (fig.5).Hence, during MD simulation, the closer proximity and favorable binding energy for revefenacin can be explained by the ability of this antimuscarinic agent to maintain these two hydrogen bonds with InhA enzyme.

Fig. 1 :
Fig. 1: An outline for the main steps and options of the applied virtual screening study

Fig. 4 :Fig. 5 :
Fig. 4: A graphical illustration of docking complex for (A) ergotamine and (B) revefenacin.Hydrogen bond is illustrated as blue continuous line while hydrophobic bond is represented by the dashed line