In silico screening of W. salutaris bioactive constituents reveal betulinic acid, ursolic acid acetate and eucosterol as potential M. tuberculosis antagonists targeting FabF


 Tuberculosis (TB) remains a long-standing burdening disease to control worldwide. The lengthy current TB treatment, which boasts with unbearable adverse effects, and frequent emergence of drug resistant strains of M. tuberculosis lays an increasing burden. This behests urgent discovery and development of alternative novel medicine to alleviate TB. In this report, in silico methods were applied to examine the propensity of W. salutaris active compounds as potential inhibitors of M. tuberculosis fatty acid biosynthesis protein (FabF). Thirteen compounds were virtually screened against FabF and subjected to molecular dynamics simulations and post-dynamics analyses to examine their inhibitory potential. Betulinic acid, ursolic acid and ursolic acid acetate had the best binding energies and hence the best inhibitory potential against FabF and desirable cytotoxicity profile. These compounds bind and interact with FabF active site residues to exert their inhibitory potential. Findings in this preliminary report warrant further experimental validation towards the development of these compounds as potential drugs targeting FabF in the treatment of tuberculosis.


Introduction
Mycobacterium tuberculosis (Mtb) is a Tuberculosis (TB) causing agent which, if not properly diagnosed, can be fatal. TB has been reported as one of ten greatest killer diseases globally which affected about 10 million people in 2017 (Who, 2019). TB is also considered a leading opportunistic disease and cause of death in HIV/AIDS patients as compromised immune system in HIV infected individuals enhances the growth of Mtb (Trinh et al., 2015).
Mtb invades the human cells through encapsulation with a cell wall made up of long-chain fatty acids known as mycolic acids (Takayama et al., 2005). Encapsulation by mycolic acids effectively allows the bacterium to grow inside the microphage, hidden from the external environment. Mycolic acids are crucial for growth and survival of Mtb, therefore providing a potential underlying mechanism for Mtb pathogenesis (Gajiwala et al., 2009). As a result, Mtb fatty acid biosynthesis pathway provides a great focal point for the discovery of anti-tubercular drugs.
Unlike mammals with only type I fatty acid synthesis, mycobacteria have both type I and II systems (Choi et al., 2000). Type I synthesis involves a single multifunctional enzyme that carries out all the reactions of chain elongation. The type II pathway has each enzymatic activity found on a separate protein. The amino acid sequences and the active sites of type I and II fatty acid biosynthesis enzymes are completely different from each other (Gajiwala et al., 2009). This allows a possibility to design speci c and potent inhibitors of enzymes from the type II pathway, with little or no effect on type I fatty acid biosynthesis enzyme. β-ketoacyl-acyl carrier protein synthase III (FabH) and β-ketoacyl-acyl carrier protein synthase II (FabF) are the proteins of fatty acid synthesis type II pathway (Daffé et al., 2017). FabH catalyses the initial cycle that involves condensation of malonyl-ACP and Acetyl-CoA, while FabF catalyses subsequent cycles of elongation that result to mycolic acids (Sridharan et al., 2008).
Currently, TB is treated with a 6 months multi-drug regimen starting with rifampicin (RIF), isoniazid (INH), pyrazinamide (PZA) and ethambutol (EMB) taken for 2 months; followed by RIF and INH for 4 months (Jnawali et al., 2013). These drugs mainly inhibit the biosynthesis of essential bacterial biomolecules including mycolic acids. However, inappropriate use of anti-TB drugs and failure to complete the full course of treatment have been reported to encourage the development of drug-resistant Mtb strains (Jnawali et al., 2013). Frequent emergence of drug-resistant strains compromises the already di cult anti-tubercular treatment (Chollet et al., 2016). Thiolactomycin (TLM) is an FDA approved antibiotic known to inhibit type II fatty acid biosynthesis and has shown great activity against a broad spectrum of pathogens  both in vitro and in vivo (Douglas et al., 2002). The major problem with these synthetic drugs is that their manufacture is cost-effective and time consuming, and they are associated with drastic adverse effects on human health, thus they do not reach and pass preclinical stages. Hence, there is a great need for alternative sources of medicine to provide safe and affordable drugs.
These compounds have been reported to exhibit a wide range of biological activities including antioxidant (Soyingbe et al., 2018), antidiabetic, antimalarial, antimicrobial, antifungal (Soyingbe et al., 2018 and antibacterial properties (Maroyi, 2013). Among these compounds, presence of hydroxyl group in the A or B ring combined with the presence of either an acid, propenyl, hydroxymethyl or 2hydroxypropyl in the E ring of the triterpene was reported to be the potential structures that can result in anti-tubercular activity (Wachter et al., 1999). Although anti-bacterial activity of these compounds have been reported in some species such as Staphylococcus aureus, Escherichia coli, Klebsiella pneumoniae and Bacillus subtilis, but to the extent of our knowledge there has not been any anti-tubercular properties reported on any of these compounds. Hence, channelling our interest in nding the active ingredient responsible for W. salutaris`s previously reported anti-mycobacterial activity. This study aims to explore the capacity of twelve W. salutaris compounds as potential inhibitors of Mtb FabF through the application of robust computational methods.

Processing of structures for molecular docking
The crystal structure of FabF exists as a homodimer with chain A and B. To reduce the computation cost, only chain A was considered for this study. Hydrogen atoms were removed from all ligands prior to energy minimization using MMFF94 force eld in Avogadro (Hanwell et al., 2012) software. UCSF Chimera (Pettersen et al., 2004) was used to separate TLM from FabF and to add hydrogen atoms to the protein in preparation for subsequent molecular docking.

Molecular docking
AutoDock Vina (Trott et al., 2010) was applied to dock all ligands into FabF. Gasteiger (Steffen et al., 2010) partial charges were allocated during docking. AutoDock graphical user interface issued by MGL tools was utilised to outline the AutoDock atom types. The active site docking was carried out in which the grid box was drawn encapsulating all amino acids that participate in the binding of a ligand to FabF crystal structure. The grid box was analysed with the following grid parameters: x = 74 Å, y = 68 Å, and z = 60 Å for the dimension, while for the centre grid was x = -7.321 Å, y = 44.488 Å, and z = -0.282 Å with the space being 0.431 and exhaustiveness = 8. The Lamarckian genetic algorithm (Morris et al., 1998) was used to generate docked conformations in accordance with their docking energy in a descending order. The resultant complexes were saved using UCSF Chimera software.

Settings for molecular dynamics simulation
The top scoring complexes of each compound from molecular docking were prepared for molecular dynamics (MD) simulations. Avogadro and UCSF Chimera software were used to prepare both the ligand and the receptor. Preparations included separating the ligand from the receptor, adding hydrogen atoms to the ligand while deleting them from the receptor and correcting the bonds of the ligand where necessary.

MD simulations
Molecular dynamic (MD) simulation of all complexes were carried out using the Amber 14 software (Salomon-Ferrer et al., 2013). The Antechamber module was used to generate atomic partial charges for the all compounds utilising Amber force eld (GAFF) and restrained electrostatic potential parameters (Wang et al., 2004). To describe the protein systems, Amber force eld ff14SB was applied (Perez et al., 2015;Tan et al., 2015). Herein, hydrogen atoms and counter ions were added to the protein, using LEAP module implemented in Amber14. Before simulations were started, all systems were enclosed in a TIP3P water box with a located 10 Å. Hence, the cubic periodic boundary conditions were implemented in all the systems. Particle-mesh Ewald approach (Salomon-Ferrer et al., 2013) was applied to treat a long-range electrostatic interactions with a cut-off distance of 12 Å. The initial energy minimization step of all systems were applied with a restraint potential. Unrestrained conjugated gradient minimization was carried out for all system using Sander module of Amber14. Canonical ensemble (NVT) MD simulations were performed for 50 ps; where the system was gradually heated from 0 to 300 K with harmonic restraints of 5 kcal/mol Å −2 for solute atoms, with the Langevin thermostat (Johnson et al., 2012) at 1 ps random collision frequency. The energy minimization, heating and equilibration steps were performed using the standard parameters and a standard MD production run for 100 ns was carried out in all systems. After MD simulations, trajectories of all systems were saved and analysed using the CPPTRAJ and PTRAJ modules of the Amber14 suite.

Binding energy calculations
The binding free energy calculations of all systems were performed using Molecular Mechanics/Generalised-Born Surface Area (MM/PBSA), an important method that gives detailed information on the binding a nity between the ligand and receptor (Genheden and Ryde, 2015). Hence, MM/PBSA approach was carried out in the current study. The approach is based on receptor-ligand complex, after the MD simulation. For a 100 ns trajectory, 1000 snapshots were evaluated during this calculation. The following set of equations describes the binding free energy calculation: (1) From the equation above, ΔE ele and ΔE vdW are electrostatic and van der Waals interactions of the ligands with the proteins in gas phase, respectively. ΔG pol is the polar solvation free energy calculated using MM-PBSA program. The ΔG nonpolar signi es nonpolar solvation free energy. The entropy contributions (− TΔS) of the binding free energies were calculated using normal mode calculation (Genheden and Ryde 2012) in the translational, rotational, and vibrational entropy components by conformational snapshot derived from MD simulations.

Per-residue energy decomposition analysis
The contribution of each residue to the total binding free energy pro le was carried out for all amino acids that interact with the compound, using the MM/PSA method in Amber 14 software.

Post-MD analyses
Generated MD trajectories were analysed using robust in silico analytical methods including Root Mean Square Deviations (RMSD), Root Mean Square Fluctuations (RMSF), hydrogen-bond network pro le, Dynamic Cross Correlation (DCC) and Principal Component Analysis (PCA) using CPPTRAJ module in Amber14. Visualisation of trajectories was performed using UCSF Chimera. The results were analysed and plots were generated using Origin software (Seifert, 2014).

Principal Component Analysis
An effective method of explaining the dynamic nature of proteins is through the essential dynamics analysis method, also known as PCA (Desdouits et al., 2015;Amadei, Linssen, & Berendsen, 1993). In the current study, PCA was performed on C-α atoms on 1000 snapshots at 100 ps time interval utilising an inhouse script. The rst two principal components (PC1 and PC2) were calculated and matrices were carried out using Cartesian coordinates of Cα atoms. PC1 and PC2 correspond to rst two eigenvectors of covariance matrices. Origin software was utilised to construct PCA plots.

Toxicity analysis
Selected ligands from MD simulation were subjected to toxicity analysis using the toxicity prediction tool ProTox-II server (Banerjee et al., 2018). This method predicts lethal doses of compounds in LD 50 in mg/kg of test subject and assign LD 50 values according to the toxicity class. This server has been previously used successfully to assess the toxicity pro les of some anti-tubercular compounds (Cinaroglu and Timucin, 2019).

Stability of simulated complexes
Based on the calculated binding free energies of all the simulated compounds, only BA, EUC and UAA had the best binding a nities relative to TLM and these were selected for further post-MD analyses. To determine the structural stability of each receptor-ligand complex, RMSD of the above-mentioned systems was calculated and presented in Fig. 2.
The average RMSD was 1.4 Å and 1.15 Å for FabF-TLM and FabF-BA respectively; and 1.0 Å for both FabF-EUC and FabF-UAA. These gures account for system stability since a standard parameter that signify a stable system is a maximum RMSD of 2 Å and below (Carugo, 2001). Therefore, we can conclude that in the duration of 100 ns simulation, all the systems formed an acceptable stable conformations.

Protein exibility during MD simulation
The RMSF of amino acids C-alpha atoms of the ligand-receptor complexes provide structural insights into the residual exibility of different regions of the studied complexes. Amino acids determine the conformational features of the protein. However, conformational changes in proteins are brought by chemical or mechanical alterations in amino acid sequence and ligand binding. Studies have shown that conformational changes occur as a result of ligand-induced motion during ligand binding to the receptor. The RMSF of the simulated complexes is presented in Fig. 3.
In all systems, the region of residues 116-150 shows high uctuation denoting high residue mobility in this region, ranging from 4.3-8.8 Å. Although binding a nity of FabF-TLM (-15.92 Kcal/mol) is lower than that of FabF-BA (-31.17 Kcal/mol), FabF-EUC (-26.27 Kcal/mol) and FabF-UAA (-29.52 Kcal/mol), FabF-TLM complex showed lower mobility in residues 211-229. However, binding of TLM to Fab-F seemed to have favoured an increased uctuation all other residues when compared FabF-BA, FabF-EUC and FabF-UAA. From these analyses, it can therefore be inferred that the binding a nity of FabF-TLM is associated with decreased residue stability. A decrease uctuation was also observed in residues 100-115 (which include residues Gly 112 and Leu 113 which participate in ligand binding) in FabF-BA, FabF-EUC and FabF-UAA compared to FabF-TLM complex. This decrease may be associated with the presence of numerous cyclohexane rings in these compounds which limit the capacity of residues in the vicinity to uctuate, hence resulting to the rigidity of FabF complexed with natural compounds. However, the binding of BA, EUC and UAA rigidi ed the active site residues of FabF, suggesting that natural compounds form stable interactions with FabF relative to TLM.
3.3 Binding free energy analysis MM/PBSA method is a popular approach used to calculate the binding free energy of ligands to biological macromolecules (Hou et al., 2011). Table 1 presents the total binding free energy of twelve natural compounds and TLM calculated from 100 ns MD trajectories. FabF-UAA (-29.52 Kcal/mol) compared to FabF-TLM (-15.92 Kcal/mol) were considered for further analysis. The difference in binding free energy between systems indicate that the force of interactions contributes higher energy in the binding of natural compounds to FabF as compared to the TLM. Electrostatic contribution to the total binding free energy was higher in FabF-BA (-4.4420 kcal/mol) and FabF-EUC (-10. 3945) compared to FabF-UAA (-2.3505 Kcal/mol) and FabF-TLM (-2.1481 Kcal/mol). However, FabF-BA, FabF-EUC and FabF-UAA systems exhibited a relatively higher van der Waal energy (-38.0026, -32.7515, -36.5090 Kcal/mol respectively) contribution to the total binding free energy compared to the FabF-TLM system with van der Walls energy contribution of -22.0500 kcal/mol. Based on these results, vdW forces are potentially important binding forces between studied ligands and FabF.
Hence, vdW contribution is comparatively higher than the electrostatic contribution as displayed in Table 1.

Per-residue energy decomposition analysis
Assessing energy contribution of individual active site residues to the total binding free energy was achieved by decomposing binding pro le into the contributions of each active site amino acids in all systems. In addition, this also provided a detailed understanding of the protein dynamics on a degree of different binding forces. Figure 4 shows that major energy contributors in FabF-BA system were Phe207, Leu113 and Phe403, with − 1.

Principal Component Analysis (PCA)
To further understand conformational motions for all sysrems, PCA calculation was carried out. The clustering system of PCA was preferred due to its ability to describe various conformational states during the simulation period by grouping molecular structure into a subset based on their conformational similarities (Wolf and Kirschner, 2013). The conformational changes of all systems were projected along the rst two principal components (PC1 vs PC2) to give a better understanding of conformations. The scatter plot (Fig. 5) shows the dominant changes in motion across PC1 vs PC2 of natural compounds and TLM bound to Mtb FabF.
A recognisable separation of motion was observed in the FabF-UAA system showing a greater correlated motion along PC1 and PC2 components, while other systems showed a signi cantly smaller or interrelated movement. It is obvious from the PCA plot that FabF-UAA system appeared to occupy less space, meaning that the binding of UAA to the residues in the active site of FabF induced the residue dynamics that resulted into conformational rigidity, thereby enhancing ligand residue interaction in this system.

Dynamic Cross-Correlation
In order to further analyse the internal dynamics of in the binding of natural compounds and TLM to FabF and to assess the existence of associated motions, DCCM analysis was carried out on the position of the Cα atoms using the 100 ns MD trajectories. The highly positive correlated movement of speci c residues are represented by yellow to red colored regions, whereas highly negative anti-correlated movement of speci c residues is represented by blue to black colored regions. The diagonal part shows obvious correlated movements, while other regions rarely exhibit such highly correlated movements. Figure 6 shows the DCC matrixes of Cα atoms uctuations of the four systems.
Overall DCCM analyses showed that the binding of BA, EUC, UAA and TLM alters the structural conformation of FabF as evidenced by changes in associated movements and dynamics. Positive correlated regions in plot A are more prominent in marked regions between residues 375-400 related to plot B, C and D. In FabF these residues are seen located away from the active site. Upon natural compound binding to FabF correlation increases relative to FabF-TLM system suggesting that natural compounds binding to FabF may have resulted in conformational changes in the protein. In relation to plot A, plot B, C and D exhibit highly positive correlated motion in residues between 200-250 whereas slightly correlated motion are seen between 250-275 region. These prominent correlated regions reside a majority of hydrophobic active cite residues hence binding of BA, EUC and UAA to FabF may have induced residue dynamics that result to more to more correlated motion among the active site residues.
These results are in agreement with observed RMSF suggesting that BA, EUC and UAA are experimentally more potent than TLM since their binding to FabF exhibit higher binding a nity to the active site residues.

Hydrogen bond network pro le
Throughout biological systems, hydrogen bonds play an important role in maintain the structural integrity of proteins, protein-ligand interactions and catalysis. This is because of their special and essential existence (Chen et al., 2016). In order to investigate the impact of BA, EUC, UAA and TLM binding on FabF, evaluation of hydrogen bond distances between amino acid residues interacting with these natural compounds and TLM in the active site was monitored for 100-ns simulations (Table 2a-2d).    (Table 2a) system showed an increase in hydrogen bonds distance between residue atoms interacting with TLM. The increase is about 0.007 Å and 0.003 Å greater than that of FabF-EUC and FabF-UAA system respectively. These differences may be responsible for relative high potency exhibited by these natural compounds in inhibiting FabF compared to TLM. Frequent occurrence of shorter hydrogen bond distance in FabF-UA (2.8 Å) and FabF-EUC (2.86 Å) relative to FabF-TLM (2.9 Å) may also have resulted to the high potency of natural compounds over TLM (Fig. 7). Also the presence of multiple cyclohexane rings in natural compounds may have in uenced ligand orientation in the active site region and decreased the average angle formation, thereby increasing ligand-residue interaction.

Predicted toxicity analysis
Some chemical compounds are promiscuous in nature, meaning that they may interact with biological receptors different from experimentally recognised receptors. Compound promiscuity re ects the molecular basis of the pharmacological effects therefore, in-depth assessment of toxicity of compounds is essential to establish the degree of promiscuity among compounds at different level of drug research (Hu et al., 2014). Toxicity of a compound de nes a quality of it being poisonous and estimating the toxicity and biological activity of a compound provides a detailed knowledge of other properties of a drug that likely have not been anticipated (Ndagi et al., 2014). The ProTox-II server toxicity prediction tool has been previously used to e ciently assess toxicity pro les of some anti-tubercular compounds (Cinaroglu and Timucin, 2019). It assigns toxicity classes of compounds ranging from 1-6: One being the most lethal and six being the non-toxic class. In this study, this sever was used to investigate the toxicity pro le of TLM, BA, EUC and UAA (Table 3). The toxicity analyses showed that the investigated natural compounds are generally less toxic, with BA and UAA being less toxic compared to TLM. This was derived from their toxicity class obtained from the ProTox-II server toxicity prediction tool. LD50 is the dose at which 50% of the test subject die upon exposure to the compound. The LD50 are associated with toxicity class of the compound. These results indicated the potential of these ligands as plausible TB drug leads. These ndings are in agreement with who reported that extracts of W. salutaris are not toxic to mammalian cell line and possess great antimycobacterial activity (Madikane et al., 2007, Oluwagbemigia et al., 2018.

Conclusion
The discovery of compounds from natural sources has given hope to the development and designing of potent M. tuberculosis FabF inhibitors for the treatment of tuberculosis. Understanding conformational features of FabF induced by natural compounds, which could open alternative avenues in the treatment of TB is unexploited. In this report, computational approaches have been employed to provide an in-depth understanding of the in uence of W. salutaris bioactive compounds. BA, EUC and UAA had the highest binding energies amongst other compounds compared to TLM, meaning they have a higher binding a nity to FabF and consequently a higher inhibitory potential. In addition, BA and UAA further boasts with another attractive character of a relatively low toxicity pro le compared to TLM. The qualities of these three compounds warrant further experimental validation as they bear promising inhibitory properties against M. tuberculosis FabF. This report serves as a starting point to the development of these compounds as drugs to potentially replace TLM as better drugs with reduced toxicity pro les in the treatment of TB.

Declarations
Funding: The authors acknowledge the University of KwaZulu-Natal and the National Research Foundation (NRF) for nancial assistance Con icts of interest/Competing interests: The authors declare no con ict of interest