Binding pattern analysis and structural insight into the inhibition mechanism of Sterol 24-C methyltransferase by docking and molecular dynamics approach

Sterol 24-C methyltransferase (SMT) plays a major role during the production of steroids, especially in the biosynthesis of ergosterol, which is the major membrane sterol in leishmania parasite, and the etiological basis of leishmaniasis. Mechanism-based inactivators bind irreversibly to SMT and interfere with its activity to provide leads for the design of antileishmanial inhibitors. In this study, computational methods are used for studying enzyme–inhibitor interactions. fifty-seven mechanism-based inactivators are docked using 3 docking/scoring approaches (FRED, GoldScore, and ChemScore). A consensus is generated from the results of different scoring functions which are also validated with already reported experimental values. The most active compound thus obtained is subjected to molecular dynamics simulation of length 20 ns. Stability of simulation is analyzed through root-mean-square deviation, beta factor (B-factor), and radius of gyration (Rg). Hydrogen bonds and their involvement in the structural stability of the enzyme are evaluated through radial distribution function. Newly developed application of axial frequency distribution that determines three-particle correlation on frequency distributions before and after simulation has provided a clear evidence for the movement of the inhibitor into active pocket of the enzyme. Results yielded strong interaction between enzyme and the inhibitor throughout the simulation. Binding of the inhibitor with enzyme has stabilized the enzyme structure; thus, the inhibitor has the potential to become a lead compound.


Introduction
Ergosterol is the major membrane sterol in leishmania parasites, which is introduced by the enzyme sterol methyltransferase (SMT). Lead discovery is currently an important area requiring diligent efforts for exploring the novel antileishmanial drugs (Siqueira-Neto et al., 2010). Inhibition of ergosterol biosynthesis provides safe and efficient therapy for leishmania (Nes et al., 2009). Both the structure and function of SMT are ideal for the design of species-specific inhibitors which have the advantage of not being harmful to mammalian physiology. Thus, it has been treated as a potent drug target in leishmaniasis and other diseases caused by fungal agents (Lorente et al., 2004). Pentavalent antimonials, polyene antibiotic agents, antiretroviral drugs, immunomodulators, and combined therapy comprise the commonly used treatment against leishmaniasis (Monzote, 2009). Polyene antibiotic agents and ergosterol biosynthesis inhibitors (EBIs) such as itraconazole interfere with ergosterol homeostasis in independent ways. Amphotericin B, a representative drug from the class of polyene antibiotic agents, interferes with membrane permeability by binding to ergosterol in the lipid leaflet, whereas the EBIs produce a significant loss of ergosterol by targeting the enzymes that interrupt carbon flux. SMT, a 353-amino acid enzyme, is a member of one-carbon group transferring methyltransferase family and employs glutathione as a cofactor (Ator, Schmidt, Adams, & Dolle, 1989;Davies, Kaye, Croft, & Sundar, 2003). The active site is a deep hydrophobic groove where the inhibitors with hydrophobic substituents would bind (Azam, Abro, Raza, & Saroosh, 2014). A series of inhibitors against SMT have been started to evolve by the end of 1980s (Nes, 2000;Nes, Xu, & Parish, 1989). Mechanism-based inactivators bind irreversibly to SMT which form the basis of this study. Sterol analogs have been studied in vitro and in vivo to interfere with SMT activity irreversibly which provide leads for the design of antileishmanial inhibitors. Excessive experimental work is being carried out on this target and its inhibitors (Ator et al., 1992;Guo et al., 1997;Mangla & Nes, 2000;Nes, He, & Mangla, 1998;Nes et al., 1999;Rahman & Pascal, 1990;Zhou & Nes, 2003;Zhou, Song, Liu, Miller, & Nes, 2004). Many inhibitors designed have ammonium or sulfonium function incorporated in their side chain which is mainly responsible for the inhibitory activity. The sterol molecule is predicted to bind to the donor site through a combination of polar and hydrophobic interactions. They mimic the reactive intermediates in the mechanism of methylation. After binding with SMT, the inhibitors undergo C-methylation but no other enzyme can catalyze the reaction . Although several drugs are in use for treating leishmaniasis, there is a need for new antileishmanial agents as already present drugs have limitations such as high cost, toxicity, and drug resistance. In silico, approaches are extensively used for the design of novel enzyme inhibitors. In terms of speed and accuracy, computational techniques that are most important in modeling protein-ligand binding are molecular docking, molecular mechanics, and molecular dynamics (Gilson & Zhou, 2007).
Herein, we used two different docking tools FRED and GOLD to show the interactions between protein and ligands so as to contribute to the drug design against SMT. This will also provide an estimate of the binding affinities of inhibitors before they are subjected to more advanced steps of synthesis and biological elucidation. We compare docking results from two different tools and investigate to what extent and how adequately do these tools explain the biological results. Besides docking, another more detailed and much valued tool in drug design is molecular dynamics simulation. In the whole process of in silico drug discovery, molecular dynamics provide a spectrum of solutions for different problems that are encountered. It is employed for the refinement of enzyme-inhibitor complexes and the calculation of binding energies (Alonso, Bliznyuk, & Gready, 2006). The best-docked pose of the most active inhibitor is investigated using molecular dynamics simulation in order to understand the dynamic pattern and conformational variability with respect to time in nanosecond scales. The simulation is analyzed with respect to statistical parameters such as root-mean-square deviation (RMSD), B-factor, radius of gyration (Rg), and radial distribution function (RDF). Hydrogen bonds distance is investigated using the statistical parameter RDF, and for further insight into coordination geometry of the enzyme-inhibitor complex, we utilized a novel analysis approach of axial frequency distribution (AFD). Movement of inhibitor was identified inside the binding pocket of enzyme during simulation using AFD. Furthermore, we attempted to understand the structure of enzyme-inhibitor complex in perspective of the network of hydrogen bonds which has great association with the stability of the complex.

Materials and methods
The protein utilized in the docking study was the product of homology modeling. Fifty-seven potential sterol analogs including the relevant substrates and two standard drugs were selected after thorough literature search (Mangla & Nes, 2000;Siqueira-Neto et al., 2010).
Structures were drawn using chemical structure drawing package, ChemOffice 2004. The energy minimization was carried out using MM2 force field. Ten thousand iterations were given with a convergence criterion of 1.000 Kcal/atom/ps. Binding regions were proposed using binding site predictors; also, a careful interpretation of multiple sequence alignment data of the SMT enzyme with its homologs was performed in a previous study by Azam et al. (2014). GOLD and FRED with genetic algorithm and exhaustive search respectively were utilized to dock inhibitors in the active site of SMT. Speed in docking is important in a well-performed study of a set of compounds. Default parameters were used in all docking experiments. The average time to dock a ligand was approximately 1-2 min by FRED and 2-5 min by GOLD.

GOLD
While docking with GOLD, for each inhibitor, ten was the highest number of genetic algorithm runs. The default parameters for genetic algorithm were designated. For van der Waals and hydrogen bond distances, cutoff values of 2.5 and 4.0 Å, respectively, were taken as used by default. The scoring functions taken into account in the study were GoldScore and ChemScore (Jones, Willett, Glen, Leach, & Taylor, 1997;Verdonk et al., 2004). The GoldScore function is made up of four components: where S hb-ext is the enzyme-inhibitor hydrogen bond score and S vdw-ext is the enzyme-inhibitor van der Waals score. S hb-int is the contribution to the function due to intramolecular hydrogen bonds in the ligand; S int is the scoring function contribution due to intramolecular strain in the ligand. whereas the free energy of binding an inhibitor to an enzyme as estimated by ChemScore function is given as follows: where S hbond , S metal , and S lipo are scores for hydrogen bonding, acceptor-metal, and lipophilic interactions, respectively. H rot is a contribution to the overall score that symbolizes the loss of conformational entropy that occurs when inhibitors bind to the enzyme. A clash penalty acting upon close contacts and internal torsion terms for poor internal conformations are added to make up the final ChemScore value that includes covalent and constraint scores also:

FRED
The first step in docking with FRED was the shape fitting process where inhibitors were taken as input and were tested where they can potentially be placed. The grid box was centered at dimensions 20.33 Å × 20.33 Å × 20.00 Å, and the total volume of the box was 8268 Å 3 . The next step was the optimization. For this purpose, Chemgauss4 was employed and it is a scoring function with Gaussian type fitting. The potentials between chemically matched positions around the docked pose of an inhibitor are used by Gaussian function. These matched positions are paired with the neighboring particular groups in the enzyme. Commonly, the interactions are hydrogen bonds; thus, an acceptable hydrogen bond score is assigned when a lone pair position of a molecule is overlapped by a polar hydrogen position on another molecule. Other types of interactions that can be estimated and scored by Chemgauss scoring function are acceptor, donors, steric, lone pairs, metals, coordinating groups, and chelator coordinating groups (Kitchen, Decornez, Furr, & Bajorath, 2004).

Molecular dynamics simulation protocol
MD simulation of SMT with most active compound was carried out using the AMBER software (Case et al., 2008). Preprocessing of the inhibitor involved the application of antechamber program. For ligands, the general AMBER force field (GAFF) was chosen and for enzyme the ff99SB force field (Duan et al., 2003;Lee & Duan, 2004;Wang, Wolf, Caldwell, Kollman, & Case, 2004). To record the topology of the enzyme and inhibitors, the LEaP module was employed (Case et al., 2008). The system was minimized before solvation with 750 steps of both steepest descent and conjugate gradient. The system was placed in TIP3P tetrahedral box of water molecules with dimension of 65Å × 70Å × 70Å (Jorgensen, Chandrasekhar, Madura, Impey, & Klein, 1983). Histidine residues were present as neutral species with hydrogen at epsilon position. In the final form, these were named as HIE 8,HIE 25,HIE 40,HIE 79,HIE 117,HIE 156,HIE 161,and HIE 244. SMT docked with an inhibitor consisted of 4605 atoms surrounded with 10,636 water molecules. Minimization of the solvated enzyme system was performed before undertaking MD simulations. 2500 cycles of energy minimization (steepest descent algorithm of 1500 steps was followed by 1000 steps of conjugate gradient) were conducted through SANDER module to relieve unfavorable clashes in the docked complex of SMT. It was followed by heating with initial period of 100 picosecond starting gradually from 0 K to a temperature of 300 K and pressure of 1 atm. Heating sets the system and assigns the initial coordinates and velocities. After heating, an initial round of 100 ps of equilibration at a constant temperature of 300 K is required before the production phase of the simulation. Before launching the production run for the system, we need to equilibrate the system, which is normally achieved by an exchange between potential and kinetic energies. This exchange then tends to stabilize the system which is apparent when kinetic and potential energies, and the total energy seems to converge with oscillations about the mean value. This behavior is achieved from a microcanonical (NVE) ensemble (fixed total energy, volume, and particle number). So our system was equilibrated on NVE ensemble before the production was started. Production run of 20 ns followed the equilibration for SMT docked complex, to bring together the correct results, statistically and numerically. The SHAKE algorithm was used for constraining the bond lengths that involved hydrogen bonds (Ryckaert, Ciccotti, & Berendsen, 1977). Simulation box was kept under periodic boundary conditions, and the Berendsen coupling integration algorithm was employed to keep the temperature constant; a non-bonded cutoff of 8.0 Å was applied. To explain long-range electrostatic interactions, Ewald summation method was used for performing MD simulations (Darden, York, & Pedersen, 1993). Files were saved after every .5 ps for analyzing the complex structure. PTRAJ module of the AMBER saved these files.

Structural evaluation
The preliminary analysis to obtain fundamental structural attributes is RDF, which is plotted between O and N of enzyme with H of inhibitor. The RDF, represented as g (r), shows the likelihood of finding an atom or a molecule in the distance r from another atom or molecule. If we find the frequency of two molecules at a distance of r, ranging from 0 to ∞, RDF can be calculated. It is a prominent analytical tool in order to describe the structural features of a system. Secondary structures of large polypeptide chains are routinely determined using this valuable tool (Donohue, 1954). Docking study highlighted the enzyme residues Asp 95, Ile 112, Glu 113, and Ala 114 which were visible in the course of the inhibitor binding process. These residues forming hydrogen bonds were further analyzed using RDF. We investigated the RDF between SMT and its inhibitor hydrogen atom to see the influence of ligand on enzyme activity. While analyzing the association of an inhibitor with the enzyme, RDF provides the estimate average distance between two atom masks on the basis of an orientation vector. It does not calculate the coordination geometry of the complex as it averages the density around all directions of the initial atom masks. To reveal this precise information relevant to the structure of the enzymeinhibitor complex, another parameter was used. This parameter utilizes the atomic molecular distribution of the atomic mask depending on the relative coordinate axis, which can give the insight into coordination geometry of the enzyme-inhibitor complex and further establish their stability. Initially, orientation vector of 5 Å was utilized to scrutinize non-bonding atoms. Figure 1 represents the initial scrutiny of the non-bonding atoms. Only those distances are taken which are within the boundary. This is represented by sphere and line model where the blue sphere at the center represents the protein atom showing the interaction during the course of simulation. The white spheres represent the 3D distribution of ligand atom at different time steps of simulation. The lines represent the x-, y-, and z-axis. Protein atom is our point of reference, and the coordinates of ligand atoms are mapped on to XY plane. Thus, a distribution is obtained on the basis of their relative position to X-and Y-axis coordinates. This distribution is plotted in the form of a 3D histogram where the axial distribution of the atoms is represented and will be referred in this work as AFD. AFD could be represented as follows: i and j represent the relative coordinate on X-and Y-axis.
The k and l values represent the cutoff value on X-and Y-axis relative to the protein atom, respectively and m i, j counts the number of observations that fall in the coordinate (i, j).
In principle, the information about particles positions is available in the simulation trajectory. In case of RDF which is a two-particle correlation, one can observe local displacement of atom from the aspect of a protein residue. One can also expand from two-particle to three-particle correlations in order to have detailed observation of the system. In this manner, local reorganizational attributes can be measured, along with local displacement, thus allowing more information to be obtained for structures. A structural network comprising of hydrogen bonds is deduced from RDFs of O-H and N-H. This structural framework is influenced by a potent inhibitor when it binds in the active site of the enzyme structure. The enzyme-inhibitor pair distribution functions are not significantly responsive to local structure reorganizations events, as the majority of the contributions result from the movement of an inhibitor in the active pocket of the enzyme to attain an energetically stable state as the simulation advances. AFD is sensitive to such local movements and hence determines the structural change.

Assessment of the docking results
Molecular docking studies were performed for the calculation of enzyme-inhibitor interactions, which is an efficient technique to predict the potential binding site(s) on the whole enzymatic target (Azam, Uddin, & Syed, 2012). For the estimation of the binding affinities of certain inhibitors, molecular docking is widely used. It is also employed for the prediction of the relationship between activity data and the docking scores (Saleem, Azam, & Zarina, 2012). Predicting binding affinities using scoring function has proved to be more challenging as compared to the prediction of their binding modes. This is where different docking tools and routines perform differently. Fifty-five mechanism-based inactivators and two FDA-approved drugs were used in this study. Different scoring functions such as GoldScore, ChemScore, and Chemgauss4 score were employed for scoring. Table 1 lists the top-ranked complexes from this study along with the different scoring functions' values. Upon visual inspection of the binding modes and key enzyme-inhibitor interactions, the topscoring inhibitors docked into the binding pocket of SMT appeared to be the same (as shown in Figure 2), but difference in the scoring function resulted in the difference of the docking scores. Hydrophobic interactions are formed in majority among all the enzyme-inhibitor interactions. Therefore, an inhibitor with hydrophobic substituent would be better as it would project into the hydrophobic groove and thus effectively inhibit the enzyme activity. It is clearly understood that the placement of each mechanism-based inactivator in the binding pocket of the enzyme is enhanced by hydrophobic substituents. Here, a stable conformation for each inhibitor is considered to be maintained by the increase of the surface complementarity. The presence of a small number of hydrogen bonds further stabilizes the interaction with the receptor to its solvent-exposed binding site. These results are consistent with the experimental findings (Nes et al., 2009) as the IC 50 can only be justified by the interaction of ligand molecule with active site or the allosteric site which hampers the activity of the protein or enzyme. In-depth analysis of the top-scoring compounds is performed. GOLD scoring functions ranked the other compounds higher as compared to the standard drugs (amphotericin B and itraconazole). In the case of FRED, majority of the compounds were ranked lower than the standard drugs. A better fit in the binding site is  shown by more negative value. To elucidate the binding mode of mechanism-based inactivators in the binding pocket of SMT, the docking pose of active compounds (19 in GOLD results with GoldScore -64.4) and (26 in FRED results with binding affinity value −11.3) has been discussed in the forthcoming section.
The best docking conformation of molecule 19 (24 (R,S),25-epiminolanosterol), which seems to be the most active compound, is depicted in Figure 3(a). The cutoff value for hydrogen bond is usually considered in range of 1.5-2.5 Å. The range represents the distance between the hydrogen joined to highly electronegative atom and another highly electronegative atom. The distance usually taken is from one electronegative element to another, so the hydrogen van der Waal's radius is added to the cutoff and the cutoff boundary is set to 4 Å. If the distance is less than 2.5 Å, they are considered as strong; if it is between 2.5 and 3.2 Å, it is considered as moderate; if it is between 3.2 and 4 Å, it is classified as weak (Jeffrey, 1997). So, according to this, the hydrogen bonds observed from the docking are further discussed in the following lines. The -OH group of compound 19 forms a weak hydrogen bond with both backbone oxygen and nitrogen of Glu 113 (3.42 and 3.90 Å), while the distance of -OH from backbone carbonyl oxygen of Ile 112 is 2.01 Å. The backbone -NH group of Ala 114 also forms weak hydrogen bond with -OH of the same compound with the distance of 3.38 Å. In this study, all the mechanism-based inactivators indicate isoleucine 112 and glutamate 113 residues in their vicinity which gives a clear idea that a potent SMT inhibitor must be interacting with these isoleucine and glutamate amino acids.
Prediction of the correct binding of a ligand into its active site was evaluated by comparing the binding affinity values. FRED performs docking faster as compared to GOLD. It is efficient in terms of speed; however, results are different from those found in the literature for the same protein in different organism, while GOLD results are more favorable and coherent with the in vitro experiments conducted against opportunistic pathogens (Nes et al., 2009). The enzyme-inhibitor complex with minimum binding affinity (compound 26, that is, 25-azalanosterol) value is depicted in Figure 3(b). Majority of the residues are polar and include Tyr 1, Gly 4, Gln 5, His 8, Gly 45, Gly47, Asn 67, Asn 68, Asn 69, Gln 72, His 117, and Tyr 124. Glu 113 is a single acidic residue which is very important in the binding site as it acts as a side-chain donor and involved in hydrogen bonding with the inhibitor. Phe 7, Val 44, Phe 96, Ile 112, Ala 114, and Ala 118 are the greasy or hydrophobic residues. Majority of the interactions formed are hydrophobic, and a few hydrogen bonds are also observed. Oxygen of both Glu 113 and Ile 112 forms hydrogen bonds with inhibitor with distances 3.20 Å and 2.86 Å, respectively. All the other compounds show similar interactions with the active site amino acid residues.
A correlation between available IC 50 data (Ator et al., 1992;Guo et al., 1997;Mangla & Nes, 2000;Nes, et al., 1997Nes, et al., , 1998Nes, et al., , 1999Rahman & Pascal, 1990;Zhou & Nes, 2003;Zhou, Song, Kanagasabai, et al., 2004; and the respective docking score or binding affinity values is determined. The accuracy of a docking protocol is determined by the assessment of the scoring function that how well it assigned scores to certain docked pose. The best pose of the inhibitors was selected from the 10 GOLD-generated poses on the basis of two gauges, inhibitor fit in the active site and fitness function scores. The fitness function scores are given in supplementary information. From the poses as viewed in VMD, it is clear how accurately GOLD predicts the correct binding conformation of ligands, GoldScore showed better predicting ability. GoldScore has compatibility problem when calculating the large ligand molecule interaction as the large molecules rely primarily on hydrophobic interactions. GoldScore has decreased sensitivity for hydrophobic interactions (Perola, Walters, & Charifson, 2006). The term S(int) in GoldScore equation represents the intramolecular strain in the ligand which can have negative value if the ligand is large as compared to the pocket or other ligands. Amphotericin B is about 923 Da in size which is comparatively higher than other ligands that is the reason value for S(int) is approximately −9.4 which removes it from the top-ranked ligands. Docking results reported by Azam et al. from AutoDock Vina predicted binding affinity of −11.1 Kcal/mol. Large molecules are usually an exception when ranking their docked conformation on the basis of GoldScore. Ergosterol was docked to find the cutoff value to determine the mechanistic inhibitors.
The compounds with low IC 50 values correspond to high GoldScores, while inhibitors scored with Chem-Score were not found in any such correlation and the scatter appeared to be random for ChemScores. Gold-Score has performed better than ChemScore fitness function, and a qualitative agreement is found with the reported IC 50 values of mechanism-based inactivators. A linear regression model was used for establishing the correlation where the GoldScore has shown the higher scoring consistency. Results yielded for current study provided the following linear regression model for explaining the data: GoldScore ¼ 15:82 þ 6:59 pIC 50 (5) where pIC 50 ¼ Àlog IC 50 (6) R 2 , the calculated correlation coefficient, shows the strength of this model. The value of R 2 came out to be .50 which means that IC 50 values and GoldScore are correlated. Additionally, this also implies that variation in IC 50 can be described by 50% of the variation in Gold-Score. The potency of newly designed SMT inhibitors can successfully be predicted using the given linear regression model. GOLD provided a satisfactory correlation value when estimated using GoldScore. Hence, for further molecular dynamics study, docked pose from GOLD is used for observing dynamics of SMT enzyme.
In an earlier in vitro study by Nes and coworkers, the same inhibitor is reported to be the most potent against the same protein in different organism. It was proved to be the most active antifungal agent while inhibiting SMT activity in Cryptococcus neoformans, and in our results, the same compound is exhibiting best antileishmanial activity. Due to sequence similarity in both the proteins, their functions are alike, and hence, the inhibitors acting upon their active sites are also same. (Nes et al., 2009).

Simulation analysis
A 20 ns molecular dynamics simulation is carried out on enzyme-inhibitor complex for the assessment of its dynamic behavior and stability. The inhibitor is stabilized mainly by hydrogen bonds. Time-dependent behavior of hydrogen bonds is analyzed through molecular dynamics. A hydrogen bond is formed between OH of the ligand and Ile 112 at the beginning of the simulation (60 ps). However, as the simulation proceeds (around 760 ps), the ligand is internalized and the OH is able to form two new hydrogen bonds with Glu 113 and Ala 114, which are present in the catalytic site. It is worth mentioning that the formation of the hydrogen bond between the ligand and the Asp 95 is stabilized in the later steps of simulations. It is noteworthy that the ligand remains inside the binding site and once formed, the hydrogen bonds are present throughout the course of simulation. Furthermore, it is observed that the enzyme-inhibitor conformation does not involve crucial changes. Aforementioned findings provide essential information regarding interactive dynamics of enzyme-inhibitor complex. RMSD is the foremost and basic measure that is used to monitor structural changes in simulations of protein dynamics. The overall RMSD of SMT enzyme-inhibitor complex is shown in Figure 4(a). Fluctuations have been observed in RMSD values; the maximum value for RMSD noted is 3.71 Å. The RMSD stays around 3.0 Å which is its average value also. It has been suggested by Azam et al. (2014) in a previous study that the ligand is stabilized in the active site of SMT primarily by hydrogen bonding which is observed in the molecular docking. The RMSD plot has moderate mean value of 3.0 Å depicting complex stability. The enzyme-inhibitor complex shows stability due to strong binding which is possible when we have high hydrogen bonds or hydrophobic interaction. In this system, we have four hydrogen bonds which are attributed to complex stability.
B-factors have been exploited to predict flexibility in structure, assess their thermal stability, and predict disordered regions in protein dynamics (Kuzmanic & Zagrovic, 2010). During any structural analysis, B-factors are very important to scrutinize: For an atom, value for B-factor that is below 30 Å 2 points to the strength of its atomic position, but a value greater than 60 Å 2 shows that the position of an atom may be disordered. B-factor plot (depicted in Figure 4(b)) shows that half of the residues lay below 30 Å 2 and the other half above the 60 Å 2 . This indicates that the structure of the protein overall is stable, and at the same time, some of the residues are disordered. The structure of the protein seems to be highly stable where the active site is present and the inhibitor is bound as the B-factor values are below 30 Å 2 . This is validating the fact that the active site is the region which is evolutionarily highly conserved and thus stable as well. Ala 288 (579 Å 2 ) has the highest value for B-factor of all the 288 residues of SMT.
Next was the Rg, which gives an understanding of the equilibrium conformation of the system. It reflects the compactness of a protein. Low Rg suggests rigid packing, while higher radii suggest the packing flexibility (Lobanov, Bogatyreva, & Galzitskaya, 2008). The enzyme-inhibitor complex reaches structural equilibrium at 20.13 Å according to the Rg calculations (see Figure 4(c)). This confirms the results of RMSD and B-factor reported herein. Overall, the structure bound with the ligand gains stability. To date, there is no experimental data available for the Rg values for SMT; hence, experiments may be conducted to authenticate the value predicted through simulation.

Radial distribution function
The RDF between N and O atoms of SMT with bound inhibitor is discussed in detail. The sharpest RDF peak for the ligand with the Ile112 interaction appears at distance of 1.774 Å with a probability value of 2.207 (see Figure 5(a)). The ligand H1 atom is 2.203 times more likely to be found at a separation distance of 1.774 Å from carbonyl oxygen of Ile 112. The highest peak in RDF graph for ligand and Glu 113 occurs at 3.527 Å and has a g(r) value of .457 indicating a .457 times greater probability of finding the ligand atom H1 and the N atom of backbone of residue Glu 113 at a distance of 3.527 Å (see Figure 5(b)). The sharpest peak for the ligand and Ala 114 was found at 4.029 Å in the graph having a g(r) of .285 ( Figure 5(c)), and the sharpest peak for ligand and Asp 95 appears at 2.124 Å having a g(r) value of .210 (see Figure 5(d)). The ligand and Ile 112 have the maximum g(r) value amid all of the aforementioned interactions. The comparatively higher values of g (r) inferred at the end of simulation signify the confidence at the distance between atoms of the inhibitor and enzyme residue. It is noticed that the presence of inhibitor enhanced the peaks of RDF for O and N of enzyme with H of inhibitor. This indicates the strengthening of enzyme-inhibitor complex structure which can be attributed to the shortening of the hydrogen bond angle and the distance between enzyme O and N atoms with H atom of the inhibitor. The results obtained from the RDF graphs depict a movement of the inhibitor structure. It appeared that a new hydrogen bond is formed toward the end of simulation run as mentioned before. For the definition of hydrogen bond between donor and acceptor, a cutoff radius of 3.5 Å was taken. The distribution of hydrogen bond distances between enzyme and its inhibitor in the complex was inspected, and it was found that at the peaks of the graphs, there is a significant increase in the probability distribution. It is clearly observed that the addition of a mechanism-based inactivator of an enzyme forces the shorter angles of the hydrogen bond distances; thereby, the ratio of strong hydrogen bonds is significantly increased by the inhibitor.

Axial frequency distribution
In order to have deeper insight, we investigated the hydrogen bonds through AFD and rationalized the formation of new hydrogen bonds observed at bonds with Glu 113 and Ala 114 during the simulation. As mentioned earlier, statistical parameter such as RMSD and RDF can give an insight into the structural stability of protein-ligand complex. Structure stability depicted cannot account for the formation of new hydrogen bonds. Appearance of hydrogen bonds between ligand and protein during the course of simulation is due to the movement of ligand within the active site cavity. Results yielded from AFD employed on the simulation trajectories agree well with the fact that inhibitor has gone through change in its position. No net displacement of the inhibitor with respect to the enzyme is shown by the RDF, but one new hydrogen bond appears at the end of simulation. The PTRAJ results could not explain the anomaly; however, AFD assists in the characterization of the aforesaid behavior. Keeping in view this fact, AFD has been defined and explored for the first time in this study. Figure 6 shows the RDF and AFD graphs taken at different intervals, where graph is depicting a clear movement of the ligand with the progression of simulation while RDF suggests that the atoms are more or less vibrating around the same position. All four hydrogen bonds that were analyzed through RDF are subjected to analysis through AFD. The results are determined at the beginning and after 20 ns of the simulation for comparison of hydrogen bonds. The asymmetrical orientation difference observed in hydrogen bond between residue 112, 113, and 114 with respect to first hydrogen atom of the ligand is due to the movement of backbone and the rotation of ligand during the course of simulation. In Figure 7, AFD graph of ligand H1 atom with respect to carbonyl O of Ile 112 depicts a clear movement of ligand atoms with respect to protein atoms. Conservation of distance observed from RDF is also visible in AFD, but the ligand atoms have taken a turn of 180°with reference to the starting position of simulation. Relative distance and surface area of distribution of H1 with respect to O atom at the start of simulation is more or less equal to the end of simulation. Similar pattern is observed in AFD graph between ligand H1 atoms with respect to backbone N of Glu 113 (Figure 8), with some deviation of distance between the atoms at the end of simulation which is coherent with RDF graph in Figure 5(b). As it is observed in RDF, there are two peaks of the graph at the end of simulation and these are also visible in the AFD graph. But the atom distribution of ligand had taken a turn of nearly 180°depicting a rotation in ligand molecule. This twist in ligand is further established, is  visible from the AFD graph of ligand H1 atom with respect to backbone N of Ala 114, and is also observed in Figure 9 with slight reduction in angle of rotation with reference to the starting point of simulation. In all of the above mentioned AFD graphs, the sharpness in the peaks of the graph with respect to XY-axis has slightly increased at the end of simulation which could depict the stabilization of hydrogen bonds at the end of simulation. Relatively higher fluctuation in the hydrogen bond between residue 95 and second hydrogen atom suggests that the hydrogen bond is less stable as compared to the other three hydrogen bonds formed with first hydrogen atom ( Figure 10). The rotation as shown in Figure 11 played its role in bringing the second hydrogen atom of ligand in vicinity of the protein residue Asp 95 which resulted in the formation of a hydrogen bond. Due to movement of the ligand and variation in interaction pattern, one can induce the behavioral evidence for a new residue involving in hydrogen bond formation from the interaction diagram taken at both the start and end of simulation ( Figure 12). This behavioral change in ligand is easily observed through AFD; thus, it gives an edge over the representation of atom distribution on the basis of a single vector. It is observed that RDF has limitations in determining the relative movement and rotation of atoms. AFD proves to be more explanatory regarding aforementioned dilemma, as it shown to be changing the numeric values with respect to rotation in the molecule. The main factor which helps in understanding the underlying mechanism is that it utilizes two vectors instead of one in plotting the distribution. RDF is useful in determining the net distance of two atoms, but AFD can additionally provide the detailed description of the geometric distribution of the atoms of ligand with respect to protein.

Conclusion
The data we had on SMT are based on the various bioactivity analyses and IC 50 values of the compounds against SMT. No crystal structure has been reported yet to establish any molecular mechanistic explanation of the Figure 11. Rotation of docked ligand toward the protein residue during the simulation. Red represents the docked conformation of ligand and protein which is also the initial point of simulation. Blue represents the conformation adopted by protein and ligand after 20 ns. inhibition mechanism. The study is designed to elaborate the molecular approach adopted by the ligands in inhibiting the SMT of leishmania infantum. In this study, we investigated the binding modes of inhibitors docked into the binding site of SMT enzyme along with the detailed analysis of the most potent complex by taking into account various parameters mainly the hydrogen bonding. As discovery of novel inhibitor templates remains a clear need, the best-docked enzyme-inhibitor complex is subjected to molecular dynamics to further our understanding of the system. Simulation is performed to check the dynamicity of ligand moiety in the presence of its designated target protein which is SMT. RDF and a novel analysis tool AFD were used to have a deeper understanding of enzyme stabilization once it is bound with the inhibitor. The increase in the final peaks of RDFs in the presence of an inhibitor as well as hydrogen bonding present in the complex indicated the strengthening of the enzyme-inhibitor complex. AFD clearly depicted the movement of the ligand indicating its strength over RDF that has its limitations for this sort of explanation. AFD, is a novel in-house analysis tool for structural analysis of the system. We are currently working on the performance evaluation of AFD and we plan to report the results in near future. Only 20 ns simulation is considered because ligand rotation is observed in first 15 ns. Later on, the ligand had no considerable movement with respect to protein. Data obtained from AFD show the stability of ligand after 15 ns as the graph depicts the relative motion of ligand atoms to be more or less within a certain limit. So the dynamicity of ligand atoms with respect to protein is depicted, and therefore, the study contained just 20 ns. This is also concurred from the current study that after applying sophisticated and sensitive tool like AFD, it is hydrogen bonding playing crucial role in the whole binding mechanism. Hydrogen bonding not only has attracted the ligand to bind with the receptor but also plays an important role in maintaining the underlying dynamics of the system. Both receptor and ligand are dwelling in a dynamic environment due to efficient hydrogen bonding pattern with respect to time which is indicated only by molecular dynamics simulation.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.10.1080/07391102.2014.1002423.