Virtual Screening of Natural Products , Molecular Docking and Dynamics Simulations on M . tuberculosis S-adenosyl-L-homocysteine Hydrolase ABDUL-RASHID

The activated methyl cycle of Mycobacterium tuberculosis (Mtb)is responsible for the regeneration of S-adenosyl methionine (SAM) from S-adenosyl-L-homocysteine (SAH). Inhibition of the key enzymes in this transformation may lead to accumulation of SAH and depletion of SAM in the Mtb cell. This has detrimental effects onthe bacterium’s cellular processes. Virtual screening of natural products from the Philippines and those in Ambinter database against S-adenosyl-Lhomocysteine hydrolase (SAHH) yielded the tautomer of the molecule, methyl 4-({2-[(4-hydroxy2-oxo-1,2-dihydro-3-quinolinyl)carbonyl]hydrazino}sulfonyl)phenylcarbamate, which displays better binding energy (-307.64 kcal/mol) than the substrate, SAH (-270.601 kcal/mol). Molecular dynamics simulations at body temperature indicated that the hit-SAHH complex is more stable than the enzyme-substrate complex.


INTRODUCTION
Tuberculosis (TB), one of the world's most deadly infectious diseases, is caused by the bacterium called Mycobacterium tuberculosis(Mtb).According to the World Health Organization (WHO), Mtb is ranked second to HIV/AIDS as the deadliest single infectious agent worldwide.In 2013, approximately 9.0 million people developed TB, which caused 1.5 million deaths in the same year. 1 The current treatment for TB consists of the administration of four first-line drugs -isoniazid, rifampicin, ethambutol, pyrazinamide -for two months, followed by 4 months of treatment with isoniazid and rifampicin. 1 This treatment however, is too time-consuming and is not usually compatible with medication for treatment of other conditions like HIV.Improper and immature use and monitoring of the drug regimen have led to the emergence of multidrug resistant, extensively drug resistant, and totally drug resistant TB strains [2][3][4] .
Part of discovering new drugs to combat TB is by targeting novel enzymes.One of the new attractive drug targets in Mtb is S-adenosyl-Lhomocysteine hydrolase (SAHH), an enzyme in the activated methyl cycle.The activated methyl cycle is responsible for the regeneration of S-adenosyl methionine (SAM) from S-adenosyl-Lhomocysteine (SAH) 5 .Compounds that hamper the activated methyl cycle cause the accumulation of SAH and depletionof SAM.SAM, the end product of the cycle, is a donor of active methyl groups in several essential cellular reactions and also a cofactor of certain enzymes 6 .Particular ratio of SAM to SAH in bacterial cell should be maintained for survival, and perturbation of this ratio level leads to growth arrest. 7The processes that require methyl groups from SAM include DNA methylation, polyamine and N-acylhomoserine lactone biosyntheses, biotin formation, 8 and ribosomal RNA methylation.9,10 In M. tuberculosis, the SAMdependent methyltransferases include Hma 11 , which catalyzes keto and methoxy-mycolic acids, and Rv2952. 12Another SAM dependent methyltransferase in Mtb encoded by umaA catalyzes the conversion of the direct methylation ofphospholipid-linked esterified oleic acid to essential fatty acid, 10-methylstearic acid(or Tuberculostearic acid, TSA). 13nce in vitro and in vivo testing of binding affinities of a vast collection of compounds to TB protein targets are time consuming and expensive, computer-aided drug discovery has been used as a more practical and efficient alternative for speedy identification of potential leads.This study has been focused on SAHH as thedrug target, based on which a pharmacophore was developed and used to screen a database of natural products.The interaction of the natural products with SAHH was modeled using molecular docking and molecular dynamics techniques.

EXPERIMENTAL
The visualization andoptimization of protein and ligand structures, virtual screening, molecular docking,calculation of binding energies, and molecular dynamics were performed by the use of BIOVIA's Discovery Studio® (DS) package. 14

Protein Preparation and Minimization
The crystal data for S-adenosyl-Lhomocysteine hydrolase were retrieved from Research Collaboration for Structural Bioinformatics (RCSB) Protein Data Bank (PDB) (http:// www.rcsb.org/).For M. tuberculosis, the protein (PDB ID: 3DHY) bound to ethylthioadenosine, a close analog of the substrate was used. 6The protein structure was preparedusing the Prepare Protein protocol of DS.The output file obtained from preparation was used as the input filein structure optimization using theMinimization protocol of DS.The resulting prepared and minimized proteinstructure was used for all other subsequent protocols.

Site Sphere Definition
The site sphere locates and defines the binding site based onthe position of the bound ligand on the original protein crystal structure.Since the protein in thisstudy was homotetrameric, only onesite sphere was generated.Thus,only one ligand wasselected and a site sphere was generated using the Define Sphere protocol.The coordinates of the site sphere were used to makeidentical copies on the prepared and minimized structure.

Pharmacophore Generation
Pharmacophore features based on the active site of Mtb SAHH weregenerated using the Interaction Generation protocol.The generatedpharmacophore features were clustered together using the Cluster CurrentFeatures and Keep Clusters Only tools based on their type as acceptor,donor, and hydrophobe.

Ligand Preparation and Building of 3D Databases
The molecules docked were a compilation of 1029Philippine natural products and the Ambinter Database (www.ambinter.com).Over 75,000 structures wereprepared using the Prepare Ligands protocol.The prepared ligands weredivided into smaller groups.Each group was built into respective 3Ddatabases using the Build 3D Database protocol.

Rigid and Flexible Fitting
The virtual screening of compounds was composed of two parts: rigidand flexible fitting.After rigid fitting, all ligands with Fit Values of 2.9 or higher were forwarded to screening by flexible fitting method.All ligands with Fit Value of 3.2 or higher were then docked to the target enzyme.

Molecular Docking and Calculation of Binding Energies
The CDOCKER protocol was used for docking each hit to SAHH.Calculate Binding Energies protocol was used for computing binding affinity of each complex.The binding energy value of the SAH-SAHH complex was used as baseline.

Molecular Dynamics
Standard Dynamics Cascade protocolwas used for the moleculardynamics simulations.Its output contains different conformations of theprotein at different temperatures.Three dynamics simulations weredone: 1.) Mtb SAHH alone, 2.) Mtb SAHH bound to substrate SAH, and 3.)Mtb SAHH bound to the top hit.The potential energies of their conformationsat the body temperature (310 K) were compared.Ligand interaction diagramswere also generated for analysis of interactions.

Preparation of the Target Protein
Out of five crystal structures of Mtb SAHH available in PDB, the one bound to ethylthioadenosine was chosen because of the close resemblance of this ligand to thesubstrate.Like SAH, ethylthioadenosine also possesses a ribose ring bondedto an adenine through an Nglycosidic bond.The only difference in theirstructure is that in ethylthioadenosine, the 5' of the ribose ring is attached toa thioethyl group rather than homocysteine.The target protein was prepared byinserting missing main-and side-chain atoms, optimizing the conformation onthe area where the missing atoms were added, and removing the watermolecules.The potential energy was subsequently minimizedthrough geometry optimization of the protein structure.The potential energyof the minimized structure of Mtb SAHH was found to be -135,696 kcal/mol.Moreover, the deviation of the prepared protein structure (Figure 1) from the original crystal structure was evaluated.
The two structures were overlaid using a carbon stick diagram (Figure 2)for ocular comparison.The overlay diagram shows very little deviation intheir structures (RMSD = 0.762Å for the main chain atoms).SAHH is a tetramer with four identical subunits.Each subunit contains onebinding site and the residues on the binding site of onesubunit is identical to the residues on the other subunits.Experimentally, the adenine ring of SAH was surroundedby Leu68, Thr71, Gln73, Leu410, Met421, and Phe425 residues. 6The twohydroxyl groups of the ribose ring interacted with ionizablegroups of Asp156, Glu218, Lys258, and Asp252.Anotherinteraction was found between the carboxylate group and His363.
A ligandinteraction diagram for ethylthioadenosine in complex with SAHH also showed that the key amino acids at the active site also interacted with ethylthioadenosine.The only aminoacid that was not seen in the diagram is His363.This observation is justifiedbecause His363 is supposed to interact with the carboxylate group of SAH,which is absent in ethylthioadenosine.
To define and mark the binding site on the structure, a site spherearound the bound ligand was generated (Figures 3).

Pharmacophore Generation and Virtual Screening against Mtb SAHH
A pharmacophorewas generated based on the active site of Mtb SAHH, by analyzing the active sitefor donors, acceptors, and hydrophobes.A total of 1,775 features weregeneratedand subsequently clustered down to only 17 features consisting of 6 donors, 7 acceptors, and 4hydrophobes (Figure 4).

Fig. 3: The site sphere on Mtb SAHH viewed closely Ligand Preparation and Database Construction
The molecules that were screened came from Philippine natural sources and the Ambinter Natural Product Database.These two databasescontain a total of 75,702 molecules.Prior to any computational protocol, theywere first prepared by enumerating their isomers, tautomers, and ionizationstates.3D conformations of each were then generated.After preparation, themolecules were consolidated as databases ready for screening.

Fig. 4: Pharmacophore features of the Mtb SAHH site sphere
The 17-feature pharmacophore was used to screen the databases of natural products against SAHH.The firstscreening process was done through rigid fitting.The arbitrary fit value cut-off was setto 2.9.All molecules that did not reach the cut-off were eliminated.Out of theinitial 75,702 molecules that were screened, a total of 3,746 reached the cutoffand were passed on to the next step.All of these molecules werere-screened through flexible fitting method.With thearbitrary cut-off set at 3.2, only 1,509 molecules passed through the screen and were set for the subsequent molecular docking procedure.

Molecular Docking and Binding Energy Calculations
The CDOCKER docking protocol works by first generating several conformations of the ligand bycontinuously adding thermal energy.This is a molecular dynamics method thattreat the ligand as a flexible structure. 5Each conformation isthen directed into an area in the site sphere where it is randomly rotatedseveral times.Each resulting orientation is saved as one pose.The program was set such that only the top ten poses in terms of binding energy are reported.The 1,509 molecules that passed both rigid and flexiblefitting were docked onto the binding site of Mtb SAHH.In addition, the preparedSAH was also docked onto the protein for reference.
The binding energy of the substrate, SAH, to the protein was found to be -270.60kcal/mol.Analysis of its ligand interaction diagram (Figure 5),revealed that all amino acid residues found in experiment werealso present in the diagram.Unlike ethylthioadenosine, the carboxylate group of SAHwas shown to exhibit charged interactions with His363.There are a total of four observable H-bonds andone charged interaction.His416 behaves as an H-bond donor to N-7 ofadenine and an H-bond acceptor to its amino group.Side chains of Thr220and Asp252 also behave as H-bond acceptors to the alpha amino group and aribose hydroxyl group, respectively.and Lys 248.There is also an interaction betweenthe pi electrons of the aromatic ring with the positively charged imidazole ringof His69.

Fig. 5: Ligand interaction diagram of SAH with Mtb SAHH
Outof the 1,509 molecules that were docked, only one was found to havebetter binding affinity to the protein than SAH.This molecule, a tautomer of methyl 4-({2-[(4-hydroxy-2-oxo-1,2dihydro-3-quinolinyl) carbonyl] hydrazino} sulfonyl) phenylcarbamate (Figure 6), has a binding energy of -307.64kcal/mol.Figure 7 shows thatall amino acidsreported to interact with SAH, save Leu68, were also found to interact withthe top hit.It is observed that its sulfone group plays an important role in itsbinding activity.The two oxygen atoms both participate in H-bondinteractions with Thr219

Molecular Dynamics
Molecular dynamics is a computational method in which the motions of amolecule brought about by a change in the environment arepredicted and simulated.In this study, the standardized dynamicsprotocol of DS was used.An input protein structure was minimizedthrough two intensive minimization protocols to bring its energy to theminimum.Thermal energy was continually added to the minimized structureuntil it reached a target temperature, i.e. body temperature.The structure was then equilibrated suchthat the energy applied was distributed over the whole structure.In this way,the most stable conformation at a specified temperature was generated.Theprocess was done repeatedly until 10 conformations of the protein at differenttemperatures ranging from 298 K to 310 K were generated.The potentialenergy of each conformation was also calculated.
The molecular dynamics simulation done on Mtb SAHH at 310K resulted in a slight deviation from the minimizedstructure.The RMSD calculations showed thatthere is a deviation of 1.066 Åon the main chain atoms and 1.024Å on the alphacarbons.Unlike molecular docking, molecular dynamics simulation does not "dock" severalconformations of a molecule to a rigid protein structure.Both the protein and theligand are treated as flexible structures.Basically, the program determines thelowest energy conformation of the complex at different temperatures.Therefore, two features of molecular dynamics explain why it is moreaccurate than molecular docking: 1) it treats the whole complex, both proteinand ligand, as a flexible structure, and 2) it simulates binding at severaltemperatures.In this study, the complexes of the protein with both the tophit and the substrate were simulated separately.Their conformations at 310 Kwere evaluated and their potential energies werecompared.The potential energy of the proteinsubstrate complex was -118,996 kcal/mol while that of the protein-top hit complex was -119,099kcal/ mol.This corroborates the results of the molecular docking procedure, the potential energy of the protein with the top hit being morenegative than that of the protein bound to the substrate.In otherwords, the complex formed by the protein and the top hit is more stable thanthe protein-substrate complex.
The ligand interaction diagram of the SAHH-SAH complex at 310 K was compared with that obtained from molecular docking.There are someinteractions that were consistent with the diagram generated from thedocking simulation.These interactions are the H-bond interactions betweenAsp252 and a hydroxyl group of the ribose, between Thr220 and the alphaamino group, and the charged interaction between His363 and thecarboxylate group.However, additional pi-pi interactionsbetween the aromatic adenine rings and the imidazole rings of His69 andHis416 were observed in dynamics simulation.The ligand interaction diagram of the SAHH-top hit complex was also compared to that generated from the docking simulation.Allnotable interactions observed in the docking simulation were also present inthe dynamics simulation.However, there are additional interactions that alsowarrant attention.His416 participated in two interactions: itacts as an H-bond donor to one of the nitrogens in the hydrazine group, andexhibits a charged interaction with the oxygen of the enolate.It should alsobe noted that unlike in the docking simulation, His69 exhibitscation-pi interactions with all aromatic rings of the ligand.Lastly, there aretwo additional amino acids that exhibit pi interactions with the ligand'simidazole ring, Phe364 and His363.

CONCLUSION
Pharmacophore-based virtual screening of over 75,000 natural products from Philippine sources and Ambinter database against Mtb SAHH was carried out.Both molecular docking and dynamics simulations established the top hit, a tautomer of methyl 4-({2-[(4-hydroxy-2-oxo-1,2dihydro-quinolinyl) carbonyl] hydrazino} sulfonyl) phenylcarbamate as a good inhibitor of SAHH and a potential lead compound against TB.