Deciphering the competitive inhibition of dihydropteroate synthase by 8 marcaptoguanine analogs: enhanced potency in phenylsulfonyl fragments

Abstract The emergence of sulfa-drug resistance and reduced efficacy of pterin-based analogs towards Dihydropteroate synthase (DHPS) inhibition dictate a pressing need of developing novel antimicrobial agents for immune-compromised patients. Recently, a series of 8-Marcaptoguanin (8-MG) derivatives synthesized for 6-Hydroxymethyl-7,8-dihydropterin pyrophosphokinase (experimental KD ∼ 100–.0.36) showed remarkable homology with the pteroic-acid and serve as a template for product antagonism in DHPS. The present work integrates ligand-based drug discovery techniques with structure-based docking, enhanced MD simulation, and MM/PBSA techniques to demonstrate the essential features of 8-MG analogs which make it a potent inhibitor for DHPS. The delicate balance in hydrophilic, hydrophobic substitutions on the 8-MG core is the crucial signature for DHPS inhibition. It is found that the dynamic interactions of active compounds are mainly dominated by consistent hydrogen bonding network with Asp 96, Asn 115, Asp 185, Ser 222, Arg 255 and π-π stacking, π-cation interactions with Phe 190, Lys 221. Further, two new 8-MG compounds containing N-phenylacetamide (compound S1, ΔGbind-eff = –62.03 kJ/mol) and phenylsulfonyl (compound S3, ΔGbind-eff = −71.29 kJ/mol) fragments were found to be the most potent inhibitor of DHPS, which stabilize the flexible pABA binding loop, thereby increasing their binding affinity. MM/PBSA calculation shows electrostatic energy contribution to be the principal component in stabilizing the inhibitors in the binding pocket. This fact is further confirmed by the higher energy barrier obtained in umbrella sampling for this class of inhibitors. Graphical Abstract Communicated by Ramaswamy H. Sarma


Introduction
Dihydropteroate Synthase (DHPS) is a pivotal enzyme in the bacterial folate biosynthetic pathway which catalyzes the production of folate intermediate, 7,8-dihydropteroate from the addition of 6-hydroxymethyl-7,8-dihydropterin-pyrophosphate (DHPP) and p-aminobenzoic acid (pABA) (Matherly, 2001). Folate is known to be essential in all living organisms for the synthesis of purines, pyrimidines, and amino acids (Bermingham & Derrick, 2002;Bourne, 2014;Swarbrick et al., 2008). However, the assimilation mechanism of folate is different in Human and pathogenic microorganisms. The higher eukaryotes consume folate from dietary sources, while plants, lower eukaryotes, and microorganisms synthesize folate by de novo pathway (Matherly & Goldman, 2003;Revuelta et al., 2018;Whetstine et al., 2002). Few enzymes in the de novo pathway, particularly DHPS are not human homologs, therefore makes them a promising target and topic of intense research for developing antimicrobial agents to reduce the mortality of HIV/AIDS and COVID-19 patients (Arooj et al., 2013;Shaw et al., 2014;Suthar et al., 2015;Tari, 2012;Zhang et al., 2021).
DHPS is reported to be encoded by folP gene (S anchez-Osuna et al., 2019) and the entire polypeptide chain is folded into the most commonly found TIM barrel-like structure (b/a) 8 in which eight a-helices bundles around the inner cylinder made of eight parallel b-strands (Babaoglu et al., 2004;Baca et al., 2000;Morgan et al., 2011;Yogavel et al., 2018). The N-terminal pole of the a/b barrel is interconnected by rigid loops while the opposite pole consists of extended and flexible loops (Achari et al., 1997), crucial for ligand recognition and binding ( Figure 1B,C). The DHPS catalytic cavity is consists of three distinct regions, such as: (a) pterin pocket, (b) pyrophosphate pocket, and (c) pABA loop. It has been reported that pABA binding is entirely dependent on the binding of DHPP at the deep cleft of the binding pocket which further induces the opening of the flexible loops located at the solvent-exposed part of the DHPS active site . The bi-substrate pocket of DHPS is therefore the target of the long-standing sulfonamide or sulfanilamide class of antibiotics to block the folate pathway by producing dead-end conjugates (Capasso & Supuran, 2019;Chakraborty et al., 2013). However, the point mutations at the flexible pABA binding loops of DHPS  in various pathogenic microorganisms, hiders the utilization of sulfonamide mediated antimicrobial therapy (Griffith et al., 2018;Ho & Juurlink, 2011;Sk€ old, 2000). Additionally, the use of sulfa-drugs against bacterial infections have been limited due to rigorous immunological reactions (i.e. sulfa allergy) and toxicity that may cause breathing problem, loss of appetite, vomiting, nausea, fever, etc (Mondal et al., 2015). As a consequence, the attention shifted towards targeting the structurally rigid pterin binding pocket to bypass sulfa drug resistance and its side effects. Although, there are many reports regarding the synthesis, structural and computational studies of pterin analogues on successful inhibition of DHPS (Azzam et al., 2020;Babaoglu et al., 2004;Hevener et al., 2010), the major bottleneck of the these inhibitors are lower solubility and higher selectivity towards the pterin binding pocket .
Therefore, a big thrust of the drug-discovery community has been redesigning DHPS inhibitors, which are capable of circumventing sulfa-drug resistance as well as have higher solubility. Accordingly, structure-based drug-discovery schemes (Dennis et al., 2014(Dennis et al., , 2016Hammoudeh et al., 2013;Zhao et al., 2016) have been applied in designing DHPS inhibitors, which not only bind the structurally rigid pterin binding pocket but also occupy the triosephosphate and flexible pABA binding pocket. Recent in silico efforts by Chakraborty and co-workers (Das et al., 2019) have also presented a free energy basis of FDA-approved sulfa drugs followed by the interpretation of crucial mutations towards sulfa drug resistance. The results obtained from this study also provided the idea to use product antagonists for retaining the inhibitory activity in face of some adverse point mutations at flexible loop regions. Additionally, a series of 8-Marcaptoguanine derivatives against 6-Hydroxymethyl-7,8dihydropterin pyrophosphokinase (HPPK), reported by Dennis and co-workers offer lower-order K D ($100-0.36 mM)value towards the enzymes of folate biosynthetic pathway (Dennis et al., 2016(Dennis et al., , 2018 and shown excellent homology with the natural product of DHPS. However, the key interaction of the compounds with the catalytic pocket and associated conformational dynamics of the receptor remains elusive and poses a fundamental challenge in the development of DHPS inhibitors. In this context, the present article addresses a pertinent query about the key factors that stabilize the 8-MG analogs at the simultaneous Pterin, pyrophosphate, pABA binding pockets of DHPS and increases the binding affinity. Here, we applied a robust in-silico approach to discover the potential hits for DHPS inhibition which have the ability to circumvent the sulfa drug resistance ( Figure 1D) as well as show higher affinity towards DHPS catalytic pocket. Specifically, we combined ligand-based drug discovery techniques such as Pharmacophore based virtual screening (Pal et al., 2019;, 3 D-QSAR (Ajmani et al., 2006;Peng et al., 2019), Density functional theory (DFT) (Andrade-Ochoa et al., 2018) to determine crucial functionalgroup substitution responsible for the potency of 8-MG compounds as well as screen large drug library. Next, structure-based molecular docking (Shahzad et al., 2020), all-atom MD simulation (Das & Chakraborty, 2020;Koneru et al., 2019), MM/PBSA (Genheden & Ryde, 2015;Kumari et al., 2014) studies were applied on top screened lead compounds to understand the conformational dynamics and free energetics of potent 8-MG analogs at the DHPS catalytic pocket. Furthermore, the enhanced sampling-based free energy simulation technique, namely umbrella sampling simulations (K€ astner, 2011;Sun et al., 2015;Zhou et al., 2015) were employed to cross-validate the affinity or resilience of DHPS inhibitors against sulfa-resistant mutations and dissociation from the catalytic cavity. Results obtained from this study reveal insights into the potency and future development of 8-MG derivatives as an antibacterial agent.

Data-set preparation
In order to develop a common pharmacophore and 3 D-QSAR model, the primary step was ligand preparation. A set of sixty-two 8-mercaptoguanine compounds (Table S2) having an inhibitory effect against dihydropteroate synthase (DHPS) were retrieved from literature (Dennis et al., 2018) and used in the present study. All compounds used for the modeling study were reported to share the same assay procedure (Seabrook & Newman, 2013). The biological activity data (K D ) of those compounds were found to span over four orders of magnitude and ranges from 420 to 0.39 lM. The biological activities of all the compounds were converted into pK D values (pK D ¼ 6-log 10 (K D )) for ease in modeling studies. The three-dimensional structures of the compounds were constructed in Maestro (Schr€ odinger Release 2017-2: Maestro, Schr€ odinger, LLC, New York, NY, 2017) builder panel. Next, all the 3 D structures of the inhibitors were imported to the Ligprep module (Schr€ odinger Release 2017-2: LigPrep, Schr€ odinger, LLC, New York, NY, 2017) in order to optimize the geometry and generate possible ionization state present at the pH of 6.9. OPLS_2005 force filed (Jorgensen et al., 1996;Kaminski et al., 2001) was employed to minimize the energy of each compound using a root mean square deviation (RMSD) cut-off of 0.01 Å. The resulting structures were used for further modeling studies.

Pharmacophore mapping
Pharmacophore is a part of molecular structure that is necessary for recognizing the ligand by a biological macromolecule. Phase (Dixon et al., 2006) was employed to find out the common pharmacophoric features of all chosen inhibitors of DHPS. Phase provides a set of five pharmacophoric features, hydrogen bond donor (D), hydrogen bond acceptor (A), hydrophobic group (H), negatively ionizable (N), positively ionizable (P), and an aromatic ring (R). Inhibitors having K D value greater than 5.3 was considered as active and the threshold of K D for inactive molecule was 4.2. The rest of the compounds were considered moderately active. This yields ten active molecules and ten inactive molecules which were used in pharmacophore model generation and corresponding scoring of that hypothesis. During hypothesis settings 'number of pharmacophoric features' was kept at five and a maximum of thirty hypotheses was generated. All the hypotheses were scored under the Phase Hypo Scoring method using default parameters for vector, site, volume, and energy terms. The best-scored hypothesis was chosen to understand the molecular fragment responsible for crucial non-bonding interactions.
To assess the predictive power of the model the bestscored hypothesis was tested against a set of decoy compounds (Kirchmair et al., 2008). Decoys are organic molecules that are assumed to be inactive against a pharmacological target. Due to the unavailability of decoys for DHPS, Decoy Finder (Cereto-Massagu e et al., 2012) was employed to generate target-specific decoy compounds for assessing the performance of generated hypothesis. The inputs for Decoy Finder include a set of active compounds for a biological macromolecule, known as 'queries' and a set of molecules from which the decoys will be selected. Next, the highly scored hypothesis was tested against 252 decoy compounds obtained from Decoy-Finder. A database containing 81,549 druglike molecules was generated by Phase-Database builder utility based on the hits collected from ZINC, PubChem and Schrodinger database. Finally, the validated hypothesis was subjected to database screening in order to screen new 8MG-compound with varied substitution.

Atom-based 3 D-QSAR modelling
The quantitative structure-activity relationship (QSAR) methods have been applied widely in computer-aided drugdesign approaches (Hoekman, 1996;Puzyn et al., 2010). It is believed that the classification of the training set and test set is the key to build a predictive 3 D-QSAR model Leonard & Roy, 2006). Many classification techniques, such as: (1) Random selection, (2) Selection based on biological activity, (3) Kennard-stone algorithm, and (4) kmeans clustering are reported in the literature. However, the first two methods are reported to fail in building QSAR models with optimal predictability. To the best of our knowledge, the nonhierarchical k-means clustering with partial least square (PLS) technique is highly reliable to build a statistically significant 3 D-QSAR model Martinez et al., 2017). The k-means clustering technique was implemented on the factor scores of the original variable matrix corresponds to each 8MG-compound considered in this study. The variable matrix consists of response variables (pK i ), predictor variable (molar refractivity (MR), ADME properties, reactivity descriptors), and topological descriptor variables predicted by PaDEL software (Yap, 2011). The details of the descriptors used in this study are given in Supplementary Section. Factor analysis on this variable matrix indicates that three factors can explain the maximum variance of the dataset. Therefore, k-means clustering on the three-factor scores of all sixty-two 8-MG derivatives was carried out to group molecules with maximum similarities. The ratio between the entire training set and test set compounds was 3:1 which is known to be standard for 3 D-QSAR modeling (Golbraikh & Tropsha, 2000;Umamatheswari et al., 2010). Atom-based ligand alignment tool of phase was employed to align all selected inhibitors prior to 3 D-QSAR calculation (Bemis & Murcko, 1996). The highest number of PLS in the 3 D-QSAR model is generally kept at N/5 which is $9 in the present case (N is the number of molecules in the training set). It is known that a higher number of PLS may cause statistical overfitting of the data (Pola nski et al., 2002;Tropsha et al., 2011). Therefore, we use six PLS factors as an incremental increase in the statistical significance to improve the predictability. The contour map obtained from the generated 3 D-QSAR model was analyzed to interpret the effect of threedimensional arrangements of structural features on DHPS inhibition. The quality of the 3 D-QSAR model was evaluated by predicting the activity of external test set compounds which were not considered during model development Martin et al., 2012).

Protein preparation
The co-crystal structure of DHPS with 8-MG analog (PDB ID: 5U10) was retrieved from Protein Data Bank (http://www. rcsb.org) and further processed with protein preparation wizard (Sastry et al., 2013) of Schr€ odinger software. The polar hydrogens were added to the protein structure at the pH of 6.9 to mimic the experimental ionization state. The hydrogen bond orientations were subsequently optimized and minimized. The conserved water molecule in the catalytic pocket was preserved. Lastly, the protein heavy-atoms were minimized by the AMBER99SB force field with an RMSD cut-off value of 0.30 Å by Gromacs-2016.5. In the course of minimization, the protein-heavy atoms were constrained, and hydrogen torsion parameters were turned off, to permit free rotation of hydrogen atoms. The receptor was minimized in vacuum.

Receptor grid generation
The minimized 3 D structure of 5U10 was imported to the Maestro GUI to prepare grid box for docking studies using Glide/receptor-grid generation tool Halgren et al., 2004). The crystal structure of the ligand was selected to identify the centroid of the active site of 5U10. A 3 D-grid box was defined within a 20 Å radius around the cocrystal ligand structure keeping van der walls scaling factor 1.0 Å and partial charge cut-off of 0.25. Prior to dock all the 8-MG derivatives, the co-crystal ligand structure of 5U10 was removed and re-docked to assess the accuracy of the docking calculation.

XP Docking
After achieving proper RMSD value between the co-crystal and re-docked ligand the grid box was selected for docking all sixty-two 8-MG derivatives at the catalytic pocket of 5U10. The van-der-walls scaling factor was kept at 0.80 during docking calculation with a partial charge cut-off value of 0.15. During the docking calculations, the receptor atoms were kept rigid while the rotatable bonds of the inhibitors were allowed to move flexibly inside the binding pocket of 5U10. The docking of all the inhibitors was performed in extra precision mode . After docking of each inhibitor, the strength of protein-ligand interactions was determined using the empirical scoring function implemented in Glide.

DFT calculation setup
The electronic properties of drug candidates are known to play a crucial role in eliciting pharmacological responses in the human body (Matysiak, 2007). The most active, least active and top three screened compounds (collected from database screening) identified by molecular docking were subjected to DFT calculation to study the quantum chemical properties of the molecules such as MO, density, molecular electrostatic potential (MEP), and dipole moments, etc. Primarily the geometry of chosen ligands was optimized in Gaussian09 (Gaussian 09 Citation j Gaussian.Com, n.d.) program using the combination of gradient correction of the exchange functional by Becke (B3LYP) (Becke, 1993) and the correlation function of Lee, Yang, and Parr (Lee et al., 1988) with both 6-31 G ÃÃ (d,p), 6-311 G ÃÃ (d,p) basis sets. All the quantum chemical calculations were carried out in the aqueous phase by employing the polarizable conductor-continuum model (CPCM) (Takano & Houk, 2005) and in vacuum. In order to explain the chemical reactivity of ligand, the frontier orbitals such as highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) were calculated. The location of HOMO indicates the electron-donating portion of the ligand, whereas LUMO conveys the tendency to accept the electron from the interacting amino acids of the protein. The electronic excitation energy was computed from the band-gap energy (HOMO-LUMO energy gap) to assess the reactivity, stability of the drug candidates (Zhan et al., 2003;Zheng et al., 2013). Lastly, the location of HOMO-LUMO over the catalytic pocket residues was identified to explain the binding mechanism at the electronic level. The Cubegen utility of Gaussion09 was employed to generate the cube files for visualizing the HOMO-LUMO orbitals. The MOs were visualized in Gauss-View (v-6.0) GUI for further analysis.

ADME/Tox properties calculation
In the drug discovery process, the ADME properties of the hit compounds are believed to play a crucial role to ensure safe drug administration. The QikProp (Jorgensen et al., 2006) module of Schr€ odinger is employed to compute druglikeness of highly active 8-MG derivatives for competitive inhibition of DHPS. The octanol-water partition coefficient (logP o/w ), H-bond donor, H-bond acceptor, molecular weight (MW), dipole moment, blood-brain barrier and membrane permeability, etc properties of the hit compounds are calculated to evaluate their drug-like properties.

The MD simulation protocol
All the MD simulations were performed in GROMACS-2016.5 (Abraham et al., 2015) using AMBER99SB force field (Hornak et al., 2006;Nguyen et al., 2014). The details of the atomistic-MD simulations of protein-ligand systems are given in Table  S9. The topology of the selected inhibitors was prepared using ACPYPE (Sousa da Silva & Vranken, 2012) software package. All the four protein-ligand systems were solvated separately in cubic boxes with TIP3P water molecules (Mark & Nilsson, 2001). A 1 nm buffer space was maintained between the surface of the protein and the edge of boxes for sufficient solvation. The approximate volume of each box was set to be 421.875 nm 3 . The overall charge of the system was neutralized by adding six Na þ ions and 0.15 M NaCl concentration was maintained. In order to remove bad contacts during solvation, energy minimization was performed using the steepest descent algorithm using 50,000 steps until a gradient threshold of 10 kJ/mol was reached. The bonds connected to hydrogen atoms in the protein are constrained by LINCS (Hess et al., 1997) algorithm, and SETTLE algorithm (Miyamoto & Kollman, 1992) was used to constrain the geometry of the water molecules. We have used leap-frog integrator algorithm to integrate the equation of motion with a time step of 2fs. All the systems were initially heated up to 300 K for 1 ns and further equilibrated for 4 ns. The temperature of the systems was maintained by V-rescale (Bussi et al., 2007) coupling algorithm (s t ¼ 0.1 ps). After temperature equilibration, the system condition is changed into an isothermal-isobaric ensemble from the canonical ensemble for 1 ns and further equilibrated up to 9 ns. During NPT equilibration the box volume was allowed to fluctuate isotropically to attain a steady-state and the system pressure was maintained by Parrinello-Rahman (Parrinello & Rahman, 1981) coupling algorithm (s p ¼ 2 ps). After attaining appropriate density of the system, we removed the restrains to allow the free or unbiased movement of protein-ligand complex for 100 ns. The cut-off radius for calculating short range electrostatic and van der Waals interaction was set to 1.2 nm. The long-range electrostatic interactions were calculated by the Particle Mesh Ewald method (Darden et al., 1993;Essmann et al., 1995). Further, the trajectories were saved at every 10 ps for our analysis. The mutation in the DHPS catalytic pocket was introduced employing Pymol software and the simulation protocol was same with the wild type protein-ligand simulations. To check the reproducibility, we further repeated the MD simulation of each system twice.

Free energy calculation
All the molecular processes such as protein folding, proteinprotein/ligand association, are known to be driven by free energy. Therefore, accurate determination of binding free energy of the inhibitors is believed to be crucial in the drug discovery process (Christ et al., 2010;Kumari et al., 2014;Srinivasan et al., 1998). The endpoint free energy method MM/PBSA (Kumari et al., 2014) was employed to compute the binding affinity of the compounds selected for MD simulation. It is evident from the literature that the thermodynamic condition for enzyme inhibition is isothermal -isobaric i.e. NPT. Hence, the last 40 ns production MD trajectory of protein-ligand complexes were considered for MM/ PBSA calculation. In general, the binding free energy of the ligand (L) at the catalytic pocket of receptor molecule (R) to form a complex (RL) can be defined as Further, the Equation (1) can be decomposed as Where, DE MM is the changes in gas phase molecular mechanics energy, which includes internal energies (DE int ) arising from the bond, atom, dihedral terms, electrostatic interaction (DE ele ), and van der Waals (DE vdw ) interaction terms. The solvation free energy term (DG sol ) consists of polar and nonpolar solvation energies between the solute and continuum solvent. The polar solvation energy term is generally computed by the Poisson-Boltzmann equation for obtaining electrostatic potential of solute in the solution while the nonpolar solvation is obtained from the linear equation of the solvent-accessible surface area of the solute (Gilson & Honig, 1988;Wang et al., 2019). The change in conformational entropy i.e. -TDS can be calculated by normal mode analysis (Srinivasan et al., 1998;Wright et al., 2014). However, due to the high computational cost, the change in conformational entropy of similar ligands is often neglected.

Umbrella sampling
Umbrella Sampling (US) simulations have long been used to explore the mechanism of ligand-protein unbinding and the determination of binding free energy (Huang et al., 2017;Miao et al., 2018;Sun et al., 2015;Zhou et al., 2015). However, the unbinding pathway modeled by such enhanced sampling methods must mimic the natural molecular movement to generate unbiased Potential of mean force (PMF). Therefore, identification of correct pulling direction is important to generate initial conformations for the US (You et al., 2019). In the present case, we employed steered molecular dynamics, an enhanced sampling method to pull the highly active 8MG derivative, and compared the unbinding free energetics with one inactive compound. The pull force (biasing potential) can be defined as Where k i is the spring constant and x i is the reference position of the ligand. x(t) is the time-dependent position of the ligand. In order to define RC, we selected the vector connecting the C b of Leu 215 and the central sulfur atom of 8-MG compounds. The reason behind selecting Leu 215 as a reference point of dissociation is due to it's smallest root mean square fluctuation (RMSF) value. The rate of ligand movement was set to 10 nm/ns and the ligand was pulled up to 4 nm from the reference point. The elastic constant k i was set to 1000 kJ/mol.nm 2 . The starting conformation for steered MD (SMD) simulation is taken from the last frame of our 100 ns conventional MD simulation performed above. The SMD simulations were carried out at 300 K, 1 bar using Nose-Hoover temperature and Parrinello-Rahaman pressure coupling algorithm (Lemkul & Bevan, 2010). Further, conformations are taken at every 0.2 nm distance as a starting conformation for the US, and 10 ns US was performed at each window. The unbiased PMF was constructed from each US window by WHAM (S. Kumar et al., 1992) and MBAR protocol (Shirts & Chodera, 2008).

Ligand based identification of structural determinants for DHPS inhibition
We have used pharmacophore modeling and 3 D-QSAR studies to primarily identify the crucial structural counterparts of 8-MG based molecules required for DHPS inhibition. A Total of thirty pharmacophore hypotheses are found to be predicted by PHASE and documented in Table S2. The DDHRR_1 ( Figure 2A) is chosen as the best pharmacophore hypothesis due to its highest Phase Hypo score of 1.34 and survival score of 5.81 which indicate good reliability and predictability of the model. This pharmacophore hypothesis consists of five pharmacophoric features such as two hydrogen bond donors (D), one hydrophobic center (H), and two aromatic ring (R) that are found to be present in all the active compounds. The inactive compounds are found to lack the hydrophobic features that are important for the competitive inhibition of DHPS. The distance and angle between the pharmacophoric sights are believed to be a crucial attribute in competitive inhibitions (Table S4). The active and inactive molecules are aligned on the pharmacophore hypothesis and shown in Figure 2B,C. It is evident from Figure 2B,C that the reason behind the difference in the active and inactive compounds is the intrinsic distances and angles of pharmacophoric sites. Further to assess the ability of the DDHRR_1 model to discriminate the bioactive compounds for DHPS inhibition, we tested the model against 252 decoy dataset compounds obtained from the decoy finder server. The pharmacophoric model is found to retrieve 100% hit compounds from the decoy dataset. The contribution of the active molecule in enrichment, the robust initial enhancement (RIE) value was calculated for all the pharmacophore models. The RIE for the DDHRR_1 model is found to be 18.14 indicating the superiority of the model ranking over the random distribution.
Lastly, taking the DDHRR_1 model as a 3 D-structural query, we performed screening to retrieve bioactive compounds as DHPS inhibitors from the database (containing 81,549 hits) made by Phase. Database screening yields 4128 hits that can act as DHPS product analogs.
Next, the 3 D-QSAR model is analyzed to unveil the effect of functional group substitution on the biological activity of 8-MG derivatives for DHPS inhibition. It is believed that atom-based 3 D-QSAR modeling is more efficient in predicting structure-activity relationships compared to pharmacophore-based 3 D-QSAR modelsjo. To reduce the number of descriptors, the original variables correspond to sixty-two 8MG compounds are converted to latent variables or factor scores. A nonhierarchical method such as k-means clustering (Everitt et al., 2011) was employed on the factor scores of all the compounds and five clusters are formed. The deviation of the compounds to their corresponding cluster is depicted in Table S1. The test set compounds are selected based on proximity to the centroids in each cluster and the remaining compounds are considered as the training set. The biological activity of the 8-MG compounds predicted by our model is shown in Table S2. In addition, the scatter plot of both training set and test set compound illustrated in Figure 2D,E indicates an excellent linear correlation and moderate difference for the biological activity of 8-MG compounds predicted by our 3 D-QSAR model. The PLS regression summary is depicted in Table 1.
The lower values of RMSE (0.21) and standard deviation (0.1349) at the sixth PLS imply that the data used for the model is reliable for 3 D-QSAR analysis. The predictive power of our 3 D-QSAR model is assessed by the cross-validation coefficient (Q 2 ) which is found to be 0.8185, indicating excellent predictive power of the model. In addition, the relevance of the model is reflected by the regression coefficient of 0.9617 for the training set compounds. The large value of F (158) and the range of stability from 0.958 to 0.063 on the maximum scale of 1 supports the statistical significance of this regression model. Furthermore, the degree of confidence of the 3 D-QSAR model is supported by the lower value of P (2.38e À25 ) (Bhole et al., 2021). The predictive power of the 3 D-QSAR model is further validated by the external test set (Dennis et al., 2014) (Table S5). The scatter plot of experimental activity vs. predicted activity and predicted activity vs. residual activity are depicted in Figures 2F and S1 respectively. The regression coefficient (R 2 ) of 0.79 indicates the good predictability of the activity of the compounds that are not considered in model development. It is believed that a 3 D-QSAR model with an R 2 value above 0.5 is a good model to predict the activity of unknown compounds (Golbraikh et al., 2003;Golbraikh & Tropsha, 2000). It is evident from Figure S1 that there are no outlier compounds predicted by our 3 D-QSAR model which further indicates the stability of the model. Furthermore, contour plot analysis is performed for 8MG derivatives to inspect the relevance of spatial arrangement of structural determinants to explain their activity. In Figure 2 the favorable zone (blue cubes) and the unfavorable zone (red cubes) for the biological activity of the compounds considered here are shown by applying the 3 D-QSAR model at 6th PLS on most active compound 62 in the training set. The blue cubes located at the -CONH group and -COOH group linked to the H11 pharmacophoric site of compound 62 (pK D ¼ 6.409) and compound 61 (pK D ¼ 6. 377) respectively indicate the preference of hydrogen bond donor group in this place ( Figure 2G). This is further supported by the presence of -CONH group in blue cube region for highly active such as compound 50 (pK D ¼ 5.319), compound 53 (pK D ¼ 5.347), compound 55 (pK D ¼ 5.538) and compound 56 (pK D ¼ 5.688). Moreover, the absence of such hydrogen bond donor in compound 59 (pK D ¼ 4.523), 60 (pK D ¼ 4.721) and 32-34 (pK D values 4.795-4.824) decrease their activity group. The red cubes located near the R12 pharmacophoric region  indicate the unfavorable region for the hydrogen bond donating group. The fact is supported by the lesser activity of compound 55 (pK D ¼ 4.721) and compound 56 (pK D ¼ 4.721) compare to the highest active compound 61. For hydrophobic interaction attributes, the blue cubes located at the benzene ring of compound 62 and compound 57 (pK D ¼ 5.690) indicate a favorable region for hydrophobic interaction with DHPS ( Figure 2H). This is further supported by the highly active compound 53 (pK D ¼ 5.347) and 54 (pK D ¼ 5.409). The location of the hydrophobic group attached with the H11 pharmacophoric site is found to be important in enhancing the activity of the compound. The placement of hydrophobic moieties greater than 4.7 Å from the sulfur atom is found to increase the activity of 8-MG compounds ( Figure S2). The red cubes located within the 4.7 Å of the sulfur atom indicate the unfavorable region of hydrophobic groups. The activity of compound 5 (pK D ¼ 3.796) and compound 9 (pK D ¼ 4.149) is decreased due to the presence of hydrophobic groups at this position. Furthermore, the red cubes located near R12 pharmacophoric site indicate the presence of bulky hydrophobic group which decreases the activity of 8-MG based compounds. This assumption is supported by the low activity of compound 1 (pK D ¼ 3.337), compound 2 (pK D ¼3.456) and compound 6 (pK D ¼ 3.796). Further, we analyzed the effect of the spatial arrangement of electron-withdrawing groups on the activity of 8-MG compounds ( Figure 2I). The blue cubes located at the benzene ring indicate the preference of electron-withdrawing groups in this region. The presence of Cl at 3rd position, Br in 4th position and -SO 2 -group in 4th position of benzene ring for compound 40 (pK D ¼ 5.004), compound 51 (pK D ¼ 5.337), compound 57 (pK D ¼ 5.609) respectively increases the activity. The red cubes located near the H11 pharmacophoric site indicate the unfavorable region for the electron-withdrawing group. It is evident from the 3 D-QSAR contour map that the presence of -CONH groups or ketone groups as a linker between the benzene ring and H11 pharmacophoric site enhances the activity of 8MG-based compounds. In addition, the presence of the electron-withdrawing group increases the activity of the compounds for DHPS inhibition.

Binding mode of 8MG-compounds to the DHPS catalytic pocket
Molecular docking studies of the 8MG-compounds retrieved from both literature and database screening were carried out at the DHPS catalytic pocket for visual inspection of crucial amino acid interactions of the inhibitors responsible for DHPS inhibition. We validated our docking results by measuring the RMSD between the crystal and redocked structure of the ligand (PDB ID: 5U10, 2.04 Å) which is found to be 0.04 nm ( Figure S3). RMSD value lesser than 0.2 nm is considered to be best for predicting the accurate binding pose of unknown inhibitors (Kramer et al., 1999). The docking score and Glide predicted binding free energy (Glide Emodel) of all the compounds are depicted in Table S6a. The best binding conformations of highly active and inactive 8MG derivatives are shown in Figure 3.
It is evident from our docking results that hydrogen bonding, p-cation, and p-stacking interactions are key interactions that anchor the ligand inside the catalytic pocket. All the interactions are found to be located between the amino acid stretch of Pro64 to Met139 and Asp185 to Arg 225 indicating the pronouncing binding pocket of DHPS. The highly active compounds found to be involved in hydrogen bond interaction with Asn96, Asn115, Gly189, Asn115, Asp185, Lyn 221 that reside inside the catalytic pocket of DHPS. The R12 and R13 pharmacophoric features of 8-MG rings are found to form p-p stacking interaction with Thr62 and Arg225 residues which play an important role in stabilizing the inhibitors in DHPS active site ( Figure S4). We also generated the map of hydrophobic and hydrophilic fields for highly active inhibitors and depicted in Figure S5. It is evident from Figure S5 that the DHPS pocket mostly consists of hydrophilic amino acids. However, small hydrophobic portions are also present in the active site. The benzene rings of highly active inhibitors are found to reside near the hydrophobic part of the surface which is found to be correlated with the 3 D-QSAR contour plot of the hydrophobic feature. The guanine ring of the 8MG derivatives is found to be buried inside the hydrophilic part located at the pterin binding site of the pocket. The docked pose of highly active compounds inside the DHPS pocket is shown in Figures 3A-F and S4A-F. The imidazole ring of compound 62 is found to be involved with p-p stacking interaction with Phe 190. There is also an equal probability of forming p-cation interaction with the guanidinium of Arg255 (Gallivan & Dougherty, 1999). Moreover, the nitrogen atoms substituted in the 8MG ring are found to accepted hydrogen bonds from Asn115 (8MG-N-H-NH-Asn115, 2.51 Å) and donate hydrogen bond to Asn115(Asn115-C ¼ O-H-8MG, 1.68 Å), Asp185(C-O À1 --H-8MG, 1.79 Å), Lyn221(8MG-NHCO-HN-Lyn221, 2.48 Å). One water-mediated hydrogen bond is found to anchor the C ¼ O of 8-MG ring to Gly217, Gly187, and Asp185. Another highly active compound 51 is found to be involved in hydrogen bonding interaction with Ser222 (NH-O ¼ S-8MG, 2.63 Å), Gly189(C ¼ O-HN-8MG, 2.00 Å). The imidazole ring of compound 51 is found to involved hydrogen bond interaction with Asp96 (C ¼ O-HN-IME-8MG, 2.78 Å) ( Figure 3C). We also perform XP docking of 4128 compounds (obtained from database screening) with respect to DDHRR_1 hypothesis and found five active compounds which can inhibit DHPS. The structures of the compounds are provided in Figure  S6. Compound S1 (2-((2-amino-6-oxo-6,9-dihydro-1H-purin-8yl)thio)-N-(2,3-dimethylphenyl)acetamide) and compound S3 (2-((2-amino-6-oxo-6,9-dihydro-1H-purin-8-yl)thio)-N-(4-((2-methylpiperidin-1-yl)sulfonyl)phenyl)acetamide) is found to obtain highest docking score ( Table S6b). The 8-MG rings of these compounds are found to have similar interactions as found in the case of highly docked compounds discussed above ( Figure  3E,F). In contrast, the inactive compound 2 and 11 are found to bind at the surface of the DHPS catalytic pocket and their guanine ring is found to interact with Gly 189 and Lys 221 ( Figure 3G,H). Such nonspecific interactions decrease the inhibitory effect of the inactive compounds and destabilize them inside the DHPS catalytic pocket. Therefore, it is evident from docking studies that the marcaptoguanine ring of the inhibitors is primarily stabilizing the base of the ligand, situated deep inside the pterin-binding pocket and the conjugation from the sulfur atom restricts the sallow region of pyrophosphate and pABA binding pocket. Such binding pose of the inhibitors is found to cover all the binding pockets of the DHPS catalytic cavity and circumvent the sulfa-resistance mutations often located at the pABA binding pocket.
The HOMO, LUMO energy is found to be varied significantly in vacuum and aqueous environment which are calculated in two different basis-set such as 6-31 G ÃÃ (d,p), 6-311 G ÃÃ (d,p). It is observed from Table 2 that HOMO, LUMO energies of 8MG compounds are lowest in the aqueous solution of 6-311 G ÃÃ (d,p), thus considered for further analysis. The HOMO energy ranges from À5.971 eV to À6.672 whereas the LUMO energy ranges from À0.903 eV to À3.125. It can be seen that there is not much difference in the HOMO energy but LUMO energy differs for each molecule. Hence, LUMO energy can be considered as the principal significator for the biological activity of the 8-MG inhibitors. Further, the band gap (energy difference between HOMO and LUMO) is calculated to assess the kinetic stability of the inhibitors (Bharathi et al., 2016). It is clear from Table 2 that inactive compounds have lower kinetic stability compare to active compounds due to a higher band-gap (E g > 5). Among the highly active inhibitors, compound 51 is found to have the lowest band-gap indicating the highest stability among the inhibitors. It can be further noted that the band-gap for pteroic acid is much lower compared to all the active inhibitors. Thus, 8-MG derivatives are chemically more inert compare to the product (pteroic acid) at the DHPS catalytic pocket which is helpful for better competitive inhibition of DHPS. The HOMO and LUMO regions are known to be directly proportional to the electron-rich and deficient regions of the molecule. Therefore, these regions indicate the reactive regions of inhibitor molecules for stable interaction with catalytic pocket residues. The occupancy of the HOMO, LUMO orbitals on the 8MG derivatives are depicted in Figures S7-S9. It can be found that the location of LUMO significantly differs in the case of active compounds. The electron reaches HOMO located at the guanine ring (R12 and R13 pharmacophoric site), which can provide strong hydrogen bonding and cation-p interaction. While the LUMO electron density majorly delocalized at the substitutions associated with the central sulfur atom ( Figure S7). The nonoverlapping position of HOMO-LUMO in 8-MG derivatives may provide better stability and interaction with the DHPS catalytic pocket. Similar observation is obtained from compound S1 and S3 ( Figure S8A,B). Unlike the active compounds, the location of HOMO-LUMO is found to be completely opposite in the case of pteroic acid ( Figure S8C). This can be the reason behind the higher chemical reactivity of this compound compare to 8MG analogs which help in maintaining kinetic equilibrium between substrate and product at DHPS active site. Furthermore, the HOMO-LUMO location is found to be placed in the guanine ring of compound 2, 6 and 11 which may hinder the stability and nonbonded interaction at the protein's active site, thereby decreasing the activity of these compounds ( Figure S9).
The global reactivity descriptors such as dipole moment, chemical potential (m), electronegativity (v), hardness (ղ), softness (s) and electrophilicity index (x) were calculated to explain the stabilization and interaction of 8-MG derivatives at DHPS catalytic pocket. It can be seen from Table S7 that compound 51, S3 is found to have a higher dipole moment due to higher number of polar atoms, thereby increasing their activity. Among highly active inhibitors, the chemical hardness of compound 51 is much lower, indicating higher biological activity of the compound . The hardness of pteroic acid is the lowest compare to all the inhibitors considered for DFT studies which further validates that the chemical reactivity of the inhibitors is comparatively lesser than the actual product of DHPS. It can be noted from Table S6 that the electronegativity and electrophilicity index of the active inhibitor are higher compare to inactive inhibitors. Therefore, the active inhibitors have well-defined electron-rich and deficient regions in their molecular geometry to stabilized itself at the active site of DHPS.

ADME/Tox profiling
In order to avoid poor pharmacokinetics, oral bioavailability, and side effects, the determination of ADME/Tox properties of the drug candidates is believed to be crucial in the drug discovery process (Llorach-Pares et al., 2018;Nolte et al., 1998). The highly active 8MG compounds (pK D > 5.0) with high docking scores are considered to be evaluated by Lipinski's rule of five, the percentage of human oral absorption. The ADME/Tox data of the 8-MG compounds are shown in Table S8. The molecular weights of all the 8MG derivatives are found to be ranging from 250-350 Da and other factors of Lipinski's rule (logP, MW, HBA, HBD) are also found to be in the permissible range. Moreover, the percent oral absorption values of the compounds are above 25%. This signifies good oral absorption of the hits in the human body. Further, reduced values of QPlogS of the 8-MG compound validate the higher solubility and drug-likeness of the compounds. It can be observed that some of the active drug candidates have the probability to block HGRT K þ channel and poor Caco-cell line permeability. However, the most active compounds (Compound 51, Compound 62, Compound 61 and Compound S1, Compound S3) have considerable values for the above two pharmacokinetic factors. Further, the bloodbrain barrier permeability values of most of the active compounds exhibited positive results which signifies that 8-MG compounds have the least possibility to affect the central nervous system. The QPlogKhsa values of active compounds indicate that the drug candidates have lesser interaction with human serum albumin. Lastly, the polar surface area (PSA) values of the active compounds indicate good permeability across the cell monolayer.

Dynamic stability of protein-ligand complex
In order to obtain further insights into the dynamic stability and conformational change of the protein-inhibitor complex, 100 ns long MD simulation of the highly active, inactive inhibitor bound DHPS and apo-DHPS was performed. We also carried out MD simulation of DHPS-pterolic acid (product) complex for the same time scale to compare it's dynamic behavior with inhibitors bound state. The docked conformations of the compounds are taken as the initial state for MD simulations. The root-mean-square deviation (RMSD) of the protein backbone and inhibitor compounds are shown in Figure 4 to determine the systematic deviation of protein-inhibitor complexes.
The upper range of backbone RMSD for DHPS in all simulations is found to be < 0.3 nm which indicates structural stability of the protein during motion (Carugo & Pongor, 2001). The RMSD of highly active compound 62 is stabilized at an average value of 0.13 nm throughout the course of simulation ( Figure 4A) while compound 61 maintained the average RMSD value of 0.1 nm up to $74 ns and then converged at a higher value of 0.17 nm ( Figure 4B) for the rest of the simulation. The RMSD of compound 51 is found to have fluctuated till $55 ns for exploring stable conformation and then reached steady-state with a lower RMSD value of 0.08 nm ( Figure 4C). Highly active compounds obtained from our pharmacophore-based virtual screening are found to reach a maximum RMSD value of $ 0.2 nm then stabilized at the end of the simulation ( Figure 4D,E). It can be seen from Figure 4 that compound 62 and 61 have the most stable conformation compares to other active compounds. In contrast, the RMSD of the highly inactive compound 6 showed several fluctuations throughout the simulation which indicates conformational instability of the inhibitor at the catalytic pocket of DHPS ( Figure 4F). Next, we calculated the RMSD of the pteroic acid-DHPS complex to compare it's stability with highly active inhibitors. The RMSD of pterolic acid is found to be stabilized at 0.1 nm which indicates that the conformational dynamics of the product is similar to the active 8-MG derivatives. It can be noted from Figure 4H that the backbone RMSD profile of the apo-DHPS is higher compare to inhibitor bound or product bound state. This indicates higher conformational stability of DHPS in their ligand-bound state. To assess the convergence in MD simulation, we calculated RMS average correlation (Galindo-Murillo et al., 2015) and given in Figure S10. It is evident from Figure S10 that the correlation profile for compound 51, compound 61, and compound S3 bound DHPS backbone approach quickly to the final structure compared to other inhibitor bound states. This further indicates the convergence of the RMS deviation for the above-mentioned 8-MG bound DHPS backbone. Further, the RMSF profiles of DHPS backbone residues in both ligand-bound and apo-state are displayed in Figure  5A,B to depict the highly mobile or flexible segments of the protein molecule. The average RMSF profile of the DHPS backbone is found to be almost similar in all highly active inhibitor-bound states. However, the residue stretches from 62 to 77 is found to have a comparatively higher RMSF value (0.4 nm) with respect to the rest of the DHPS backbone, indicating the higher flexibility of this region during motion ( Figure 5A). Apart from this most of the backbone regions of the DHPS showed an RMSF value less than 0.1 nm due to the compact tertiary structural arrangement. The degree of fluctuation at the catalytic pocket residues (Asp 96 -Met 139 & Asp 185 -His257) is found to be < 0.2 which helps the formation of the stable protein-inhibitor complex. In the case of inactive compounds, the residue starch 62-77 of DHPS is found to fluctuate more ($0.58 nm) compare to the active bound state (below 0.42 nm). Therefore, the flexibility of regions 62-77 is essentially restricted by active inhibitor binding for a stable protein-ligand complex. Further, it can be noted that the RMSF value of the residue stretches 222-239 is also stabilized by the active 8MG compounds. The C-terminal end of DHPS is found to fluctuate significantly in both active-bound and inactive bound states due to the solventexposed unstructured loops.
The compactness of ligand-bound and apo-state of DHPS-8MG complex was evaluated by calculating the gyration radius (rGyr) throughout the simulation time ( Figure 5C). The rGyr value of the DHPS backbone in all the simulations ranges from 1.72 nm to 1.79 nm with an average rGyr value of 1.75 nm. This confirms the fact that inhibitor binding does not cause the unfolding of the DHPS enzyme over the simulation period. The dynamic stability of the protein-ligand complexes correspond to simulation_batch 2 and simula-tion_batch 3 are given in Figures S11-S13. The trend of the results is found to be almost similar compare to simulation batch 1.

Non-bonded interactions
The stability of the docking predicted interaction profile of the inhibitors is assessed by MD simulation studies. The average number of hydrogen-bonding interactions involved in stabilizing the inhibitors inside the DHPS catalytic pocket was determined and depicted in Figure 5D. We also determined the percentage of hydrogen bond occupancy (HBO) for the donor and acceptor pairs (Table 3). It is found that compound 51 formed the highest hydrogen bonds with DHPS catalytic pocket compare to other 8-MG derivatives indicating the substitution of sulfonamide group enhances the electrostatic interaction to stabilize the outer part of the compound at the flexible loop of DHPS catalytic pocket. Similarly, compound S3 is found to form higher hydrogen bond compare to the S1 which further indicates the importance of the sulfonamide group. Compound 61 is found to have less hydrogen bond interaction compare to other active compounds. The residues such as Asn 115, Asp 185, Lys 221, Ser 222, and Arg 255 are found to be common which contribute hydrogen bonds to stabilize the active compounds. Additionally, Leu 184 (HBO: 98.0%), Leu 215 (HBO: 64.8%), Gly 215 (HBO: 47.8%) for compound 62, Asp 96 (HBO: 33.1%), Gly 189 (HBO: 25.7%) for compound 51 and Thr 62 (HBO: 381%) for compound S3 stabilizes them at the catalytic pocket of DHPS. Therefore, it is evident from our MD study that most of the docking predicted interactions were present throughout the trajectory. The water-mediated hydrogen bonds between the active site residues and ligand are also crucial for anchoring the ligand in the catalytic pocket. The visual inspection of our MD trajectory suggested that Glu 60, Asp 116, Asp 165, Asp 185, Lys 201, and Leu 215 are involved in water-mediated hydrogen bonds with the active inhibitor molecules.
Apart from hydrogen bonding, the ligand or drug molecules are believed to be stabilized at the enzyme catalytic pocket by a delicate balance of week interactions such as p-p stacking (Chen et al., 2018;Thakuria et al., 2019) and cation-p (Ma & Dougherty, 1997;Scrutton & Raine, 1996) interaction. In the present scenario, the contribution of such non-covalent interactions to stabilize the highly active inhibitors is assessed from our MD trajectories.
The distance criteria for p-p stacking interaction are reported to be 3-8 Å, while the angle between the normal vector of two aromatic rings should be 0 -120 (Thakuria et al., 2019). The Phe 190 is found to be involved in p-stacking interaction with the guanine ring of the highly active inhibitors buried deep inside the DHPS catalytic pocket. This is assessed by the time evolution of the distance and the angle between the aromatic ring of Phe 190 and the marcaptoguanine ring of active inhibitors (Figures S14 and S15). In all the cases the average distance is found to be ranging from 0.58 nm to 0.75 nm and the angle between the normal vectors is 90 -120 . In the case of compound 62, p stacking interaction is absent between the marcaptoguanine ring and Phe 190. However, Phe 188 is found to make such interaction with the substitute CONH-C 6 H 6 ring and stabilizes the solvent-exposed part of the ligand. Further, the marcaptoguanine ring of the active inhibitors is displayed cation-p interaction with Lys 221 and Arg 255 (K. Kumar et al., 2018). The time evolution of the distance between the marcaptoguanine ring and the positively charged atom of both the residues is found to be below 0.7 nm (Figures S16 and S17) indicating pronounced cation-p interaction to add stability for protein-inhibitor complex (Xiu et al., 2009).

Principal component analysis (PCA)
It is often challenging to identify and understand the important motions of a biomolecular system from the large dataset obtained during the course of MD simulation. Therefore, to reduce the directionality of the dataset and filter the largescale motion within the explored timescale, PCA is used (Amadei et al., 1993;Berendsen & Hayward, 2000;David & Jacobs, 2014). To calculate the collective motion of the protein, the mass-weighted matrix of principal components can be defined as (Sittel et al., 2014) < > is indicating the average of instantaneous structures sampled within the course of simulation and x i is a massweighted atomic coordinate. In the present scenario, we use the last 40 ns of the trajectory for analyzing the principal motions of the protein in an inhibitor-bound state. The eigenvectors are calculated and given in Figure 6A to comprehensively analyze the effect of inhibitor binding on protein motion. It is found that the total motion of inhibitor bound and apo DHPS monomer structure are dispersed over 789 eigenvectors. The rapid decrease of the eigenvectors towards the Eigen index indicates that the initial 5 eigenvectors are mostly contributing to the collective motions of DHPS-8MG complexes during MD. It is found that first eigenvector corresponds to 42.65%, 43.79%, 45.78%, 61.05%, 64.94%, 49.2%, 44.51% and 48.3% motions for DHPS bound with compound 62, 61, 51, S1, S3, 6, Pteroic acid and Apo-DHPS respectively ( Figure 6A). Thus, it can be inferred that the inhibitor-bound DHPS systems are stable compare to the DHPS bound with highly active compounds obtained from database screening. Additionally, the lower collective motions of the Apo-DHPS backbone indicate that the amalgamation of product analogs increases the principal motion crucial for minute conformational transition.
The first principal movement (PC1) of the ligand-bound and apo-DHPS are shown in Figure 6C-J. The length of the cone represents motion magnitude and the direction is shown by the pointing of the arrow. It is found that essential motions of apo and product-bound DHPS are only associated with flexible loop regions ( Figure 6B) that are crucial for ligand recognition. However, the entire DHPS structure is stable in both cases. In contrast, the flexible loops of inhibitor bound DHPS (except compound 62) showed closing motion which induces enhanced movements in the rigid domains of the protein backbone. Especially, the amplitude of the motions in the rigid domain is significant when the DHPS is bound with compound S1 and S3. This observation suggested that the binding of highly active 8MG compounds induces the shrinking of catalytic pockets to restrict natural substrate binding and product release. The essential space of the overall dynamics of inhibitor bound DHPS and apo-DHPS can be visualized by the 2 D-projection of the first two principal components such as PC1 and PC2. It is evident from the projection plot that the extent of sampling is remarkably different in each case ( Figure S18). DHPS bound with compound S1 and S3 is found to span larger conformational space compare to other inhibitor bound states, indicating highly systematic movement of flexible loops located at the surface of the DHPS catalytic pocket. This can be related to the porcupine plot shown in Figure 6F,G. In contrast, DHPS bound with highly active Compound 62 and 61 have confined conformational space which implies that intrinsic flexibility of DHPS reduces due to the presence of such compounds. A similar conformational landscape is found for the product-bound state of DHPS. Further, the larger essential space of apo-DHPS provides a hint regarding the recovery of the flexibility of protein backbone, gained due to the removal of the ligand from the catalytic pocket. Similar phenomena were found for the inactive ligand-bound state of DHPS as a result of nonspecific inhibitor binding. Lastly, the cosine content of PC1 is measured to confirm the absolute convergence of our MD trajectory (Papaleo et al., 2009). The cosine content corresponds to PC1 from all MD trajectories is determined to be under 0.5 indicating adequate sampling of MD trajectories.

MM/PBSA binding free energy
The binding energies of highly active 8MG derivatives such as Compound 62, 61, 51, and top-scored compounds S1 & S3 obtained from pharmacophore-based virtual screening are computed by MM/PBSA (Kumari et al., 2014) algorithm and the results are given in Table 4. It is found that compound 51 has comparatively lower binding energy (-67.26 kJ/mol) among the highly active 8MG analogs. The free energy decomposition shown in Table 4, indicates that the electrostatic energy (-154.25 kJ/mol) is the principal contributor to the better stabilization of compound 51 at the DHPS catalytic pocket. The compound 62 is found to have a stronger van der Waals component (-179.56 kJ/mol) due to the presence of hydrophobic benzene ring. It can be observed that the Compound-S3 obtained from database screening has the lowest binding energy of À71.29 kJ/mol due to the most stable electrostatic contribution of À260.65 kJ/mol. Further, the weak binding affinity of inactive Compound 6 is reflected by it's comparatively higher binding energy of À28.19 kJ/mol. However, the binding affinity of the product i.e. 7,8-dihydropteroic acid is lower compare to the 8MG compounds which indicates better competitive inhibition of DHPS by such product analogs. Notably, in all the cases the polar solvation term is found to be unfavorable in the context of stabilizing the protein-inhibitor complex. The binding free energies of the compounds correspond to simulation_batch 2 and simulation_batch 3 are given in Tables S10 and S11  respectively. The trend is found to be similar to simulation_ batch 1. Furthermore, the total binding energies in each protein-ligand complex are decomposed into the contribution made by each DHPS residue and shown in Figure S19. This enables a comparison between the relative contribution of the interacting residues towards the overall binding affinity. The residue-inhibitor interaction spectrum is found to be almost similar for the compounds which showed lower binding energy towards DHPS. It is found that interaction of Glu 105, Asn 115, Phe 188, Phe 190, Lys 221, and Arg 255 is common to the highly active 8-MG derivatives considered for free energy calculation. The amino acid region Thr53-Thr62 also stabilizes Compound 61, 51, and S3 at the catalytic pocket of DHPS. In the case of compound 6, the nonspecific residues such as Leu 201, Gly 217, Arg 220, and Ser 239 have contributed towards the binding free energy. It can be observed that amino acid residues with hydrophilic side chains have more contribution to stabilizing the inhibitors. However, compound 62 is mainly stabilized by hydrophobic residues.

Effect of mutation
In order to verify the mutation-resistant activity of the 8-MG derivatives, single point mutations were introduced at residue number 64 (Pro64Ser) and 221 (Lys221Gln) of the DHPS backbone. These two mutations are reported to destabilize both the pterin, sulfa-based inhibitors against DHPS, and the strongest drug-resistant mutation ever found for DHPS . The effect of the mutation was tested on highly active Compound 51 and S3 and compared with the wildtype Protein-ligand complex. The RMSD of compound 51 bound DHPS backbone is found to be stabilized with an average value of 0.25 nm while the RMSD of compound 51 is stabilized with an average value of 0.15 ( Figure S20A). It can be seen that the single point mutation did not cause a significant change in ligand conformation. There are five hydrogen bonds are found between mutated DHPS and compound 51 throughout the simulation. Similarly, in the case of compound S3, the backbone RMSD is found to be stabilized at 0.3 nm for the last 30 ns of the simulation. The compound S3 is found to be stabilized with a higher RMSD value of 0.23 nm for the last 10 ns of the simulation ( Figure  S20B). The total number of hydrogen bonds between compound S3 and mutated DHPS is found to be 4. Further, the RMSF values for the amino acid residues are slightly increased in the case of mutated DHPS. However, the change in RMSF is not much significant to destabilize the highly active 8-MG compounds from DHPS active site ( Figure S20C, D). It is evident from the gyration radius profile that the point mutations are not responsible for changing the overall shape of DHPS ( Figure S20E,F). Lastly, we calculated MM/ PBSA calculation to assess the binding affinity of the compounds at the mutated catalytic cavity of DHPS. The binding energy of compound 51 and S3 is found to be À52. 38 kJ/ mol and À69.49 kJ/mol (Table S12) which indicates that 8-MG compounds have a promising ability to circumvent DHPS resistant mutations.

Drug unbinding by umbrella sampling
In order to study the dissociation pathway of the 8MG-derivatives from the DHPS catalytic pocket and to determine their free energy of inhibitor dissociation, umbrella simulation studies have been employed. Due to the high computational cost, we have selected highly active compound 51, compound S1, and inactive compound 6 to study the dissociation process. The reaction coordinate is extended from 0 to $5 nm to ensure complete dissociation of the ligand from the binding pocket. It can be observed from Figure 1A that substrate entry and product release can occur through narrow, solvent-exposed, and metastable unstructured loops which acts as the gateway to enter a comparatively stable pterin binding pocket of DHPS. Consequently, there is a major dissociation energy barrier to be crossed by productanalogs to escape from narrow catalytic pocket. Therefore, an external biasing potential has been applied to steer the ligand towards the bulk solvent. The structure of the DHPS binding pocket is restrained by applying harmonic potential on the C a atoms of binding pocket residues. Application of pulling force allows the perturbation of dynamic equilibrium of the system, thus hindering the calculation of thermodynamic quantities directly from the steered molecular dynamics (SMD) trajectories with minimum errors. In such a scenario, the weighted histogram analysis method (WHAM) is commonly used to compute equilibrium data from nonequilibrium SMD trajectories. The force acts on the ligand are found to be increasing sharply due to constant velocity pulling until a breaking point is reached ( Figure S21). At this point, the critical interactions are disrupted which allows the inhibitors to dissociate from the catalytic pocket of DHPS. The events leading up to the uncoupling of each ligand considered here are found to be differed significantly, although the direction of pulling is the same in each case. Therefore, the force vs time curves, computed from COM pulling is not essentially comparable due to the path dependency of the dissociation process. In the case of compound 51, the hydrogen bonds with Asp 96, Asn 115, Asp 185, and water-mediated hydrogen bonds with Cys 137, Phe 188, Gln 226 are found to be broken at the point of maximum force. Similarly, the hydrogen bond between Asn 115, Asp 185, and compound S1 are disrupted at the breaking point. It is evident from the above discussion that Asp 115 and Asp 185 are critical to the stability of the inhibitors inside the DHPS catalytic pocket. However, compound 6 is found to be formed only water-mediated hydrogen bonds with nonspecific Arg 233 and Arg 235 during the dissociation process. Further, the free energy profile along the reaction coordinate (RC) is analyzed and the PMF profile for each studied system considered for umbrella sampling is shown in Figure 7.
It can be seen from the graph that the PMF profile for the inhibitor dissociation follows a similar trend. The PMFs are found to drop to a minimum value and then increase to a stable state when RC reaches to $3 nm-4 nm. The flat region of the PMF indicates that the inhibitor is in isotropic bulk or the interaction between the protein and ligand is completely disrupted. Free energy of dissociation is determined from the difference between the highest and smallest values of the PMF curve. It is evident from Figure 7 that compound S3 has the lowest binding energy of À93.29 kJ/mol indicating the highest energy barrier to dissociate from the DHPS catalytic pocket. Further, compound 51 is found to show lower binding energy of À73.42 kJ/mol, while compound S1 has a comparatively lesser binding affinity of À49.36 kJ/mol. This result indicates that the stabilization of 8-MG compounds containing sulfonamide or sulfonyl fragments, at the flexible pABA binding pocket enhances the barrier for inhibitor dissociation from DHPS catalytic pocket. In the case of compound 6, the dissociation energy barrier is found to be less due to the highest binding energy of À31.79 kJ/mol. It can be noted here that the binding free energy (i.e. dissociation energy) computed by umbrella sampling is significantly smaller compare to binding energy estimated by MM/PBSA or MM/GBSA method (Sun et al., 2015). This is due to the ignorance of conformational entropies and dependencies of polar solvation equation schemes (Genheden & Ryde, 2015) in binding free energy calculation. However, the trend of the binding free energy of the 8MG derivatives is found to be correlated in both US and MM/PBSA methods.

Conclusion
It is extensively recognized that drug discovery with respect to DHPS, in general, is a time taking, expensive, and highly complex process. In-silico studies have become an indispensable tool for probing unknown chemical space to identify novel bioactive molecules with improve effectivity. In this article, in-silico techniques such as pharmacophore mapping, atom-based 3 D-QSAR modeling, high throughput virtual screening, molecular docking, MD simulation with a pair of free energy techniques, namely MM/PBSA & umbrella sampling were employed to identify crucial counterparts of DHPS competitive inhibition with 8-MG derivatives. The contributions of various atomic features, like hydrogen bond donor, hydrophobic interaction, and electron-withdrawing groups on the 8MG analogs, towards DHPS inhibition were highlighted. 8 MG-compounds containing phenylacetamide (Compound S1) and phenyl sulfonyl (Compound S3) were identified as the highest-ranked hits among the organic molecules obtained from database screening. Docking study revealed that hydrogen bonding, p-p stacking, and p-cation interactions of product analogs with active site residues provide the stability of the inhibitors. It is found that most of the non-bonded interactions are present between the guanine ring of the inhibitors and the active site residues located deep inside the active site. The active site residues such as Asp 96, Asn 115, Asp 185, Phe 188, and Ser 222 are mainly involved in hydrogen bonding while Phe 190, Lys 221, and Arg 255 are associated with p-p stacking and p-cation interactions. Further, the DFT study reveals that the overlapping of HOMO-LUMO orbitals on common guanine ring decreases the activity of 8MG analogs. Moreover, the highly active 8-MG compounds were found to have higher water solubility as predicted by ADMET calculations. The dynamic stability of the inhibitors at the DHPS catalytic pocket was confirmed by atomistic MD simulation studies. The active compounds are found to be most stable compared to the inactive one due to their highly specific binding in DHPS catalytic pocket similar to pteroic acid. The binding of highly active compounds such as Compound 51, 62, and S3 was found to reduce the flexibility of residues 62-77 which is crucial for ligand stability. The non-bonded interaction predicted by docking calculations was found to have similarities with the MD simulation studies. PCA studies indicate that binding of the top-scored compounds induces closing motion in the flexible loops thereby causes prolonged inhibitions of DHPS. Further, MM/ PBSA free energy studies indicate that electrostatic energy is the principal contributor to the stability of the inhibitors. Lastly, the free energy of dissociation of the product analogs reveals that water-mediated hydrogen bonding is also important for the stability of the inhibitors. The overall study summarizes that the substitution of the sulfonamide class of fragments at the central sulfur atom of 8-MG compounds is crucial to enhance the activity and stabilize them at the flexible pABA binding loops. The information obtained from this study offers a basis for understanding the conformational dynamics of DHPS with highly active 8-MG analogs and inspire medicinal chemist to synthesize newer generation of antimicrobials.