**3.1 Homology modeling**

*L. major* trypanothione reductase was modeled with an appreciable degree of accuracy and validated via several bioinformatics tools (**Figure 1A**). The protein sequence of *L. major* with the accession number XP\_001687512.1, having 491 amino acids when searched amongst protein homologs in NCBI resulted in a list of similar structures with PDB codes 2JK6, 2X50, 6ER5, 1FEA from *Leishmania infantum* with their respective percentage identity of 95.72, 95.71, 95.70, and 78.78%, respectively. The similarity of structure and sequence between *Lm*TR protein sequence and the template 2JK6 was 95.72%. This favored its selection for the modeling process. The best model selection was based on the least discrete optimization potential energy (DOPE) score which informs about the energy of the protein.

## **3.2 Active site prediction**

In predicting the active site of the protein, results obtained showed 73 binding pockets. Several pockets were identified but the pocket with the largest volume and surface area of 595.278 Å3 and 924.887 Å<sup>2</sup> , respectively was selected. Larger pockets favor conformational rotation during virtual screening. A total of 80 amino acid active site residues were predicted from CASTp 3.0 (**Table A1**). The predicted binding site was visualized using PyMOL (**Figure 1B**).

#### **Figure 1.**

*(a) 3D structure of modeled protein in cartoon representation. A monomer of* Lm*TR with alpha helices represented in red, the beta sheets in yellow and the loops in green. (b) Surface representation of the* Lm*TR. Active site (in red) is the FAD binding region.*

The predicted active site was finally confirmed to be FAD binding site. TR has been studied to be a homodimer protein constituting FAD-binding, NADPH-binding, central, and interface domains [5]. The predicted site was concluded to be an FAD site in *Lm*TR by comparative studies between the FAD-2JK6 and FAD-*Lm*TR complexes (**Figure A1-A**). The study showed common hydrogen bonding residues including Ser14, Gly15, Arg287, and Thr335 interacting with FAD. These residues correlate with that which was predicted by CASTp 3.0 including other several common residues participating in hydrophobic interactions such as Gly13, Asp35, Ala46, Gly50, Thr51, and Cys52 (**Figure A1-A**).

#### **3.3 Structural validation and quality prediction**

PROSA was used to determine the quality of the structure by comparing it to protein structures that are experimentally solved in the PDB database. It validated the model based on the "quality score or *z*-score" with a value of 11.68. The *z*-score shows whether the predicted model is of X-ray or NMR quality with regards to the amino acid residues length. This *z*-score value showed that the modeled protein fell in the range of proteins solved experimentally by X-ray crystallography (**Figure A2-A**). Verify3D validated the 3D structure of the model. A good 3D structure is expected to have at least 80% of its amino acids to have scored greater than or equal to 0.2 in a 3D/1D profile. This model passed with an appreciable result of 91.65% of residues having a score greater or equal to 0.2 (**Figure A2-B**). Further validation with PROCHECK resulted in the generation of a Ramachandran plot (**Figure A3**). The plot described the rotations of the polypeptide backbone around the bonds between N-Cα (Phi, φ) and Cα-C (Psi, ψ). The plot allowed the viewing of the distribution of these torsional angles taking into consideration the allowed rotations and rotations that are unfavored which can result in a collision or steric hindrance. A protein with 90% of its residues in the most favorable region is considered a good model. ProQ predicted the quality of protein based on the MaxSub and LGscore scores [29]. The predicted LG score was 6.447 and MaxSub score was 0.520. LGscore >4 implied that the model was extremely good. Also, a MaxSub score > 0.5 implied that the model was very good [47].

#### **3.4 Pharmacophore modeling**

The active ligands used as training sets to develop a pharmacophore allowed features similar to the two compounds to be identified and combined into a single geometric function as the basis for the generation of the pharmacophore (**Figure 2**). Pharmacophore generated utilized features that contributed regions of hydrophobicity and hydrogen bond acceptors incorporated in the model for selective screening. The oxygen from nitrogen dioxide contributed to hydrogen bond acceptors with aromatic and alkene groups contributing to the hydrophobic region (**Figure A4**). A number of 10 hypotheses were developed for the model. The best hypothesis with a similarity of 58.12% had a score of 0.8537 and was selected based on the AUC score of 0.99 generated by LigandScout (**Figure A5-A**). The AUC score was used as a metric to validate how best the pharmacophore model created could distinguish rightly between active compounds and decoys. This intends to reduce false positives and negatives during the screening process.

The screening process was successfully completed with the model which identified 42 compounds that matched the pharmacophore model with a pharmacophore fit score ranging from 55.32 to 57.98 (**Table A2**).

#### **Figure 2.**

*2D structure of PNTPC (left) and CNQB (right) showing shared features. PNTPC and CNQB showed two hydrogen bond acceptors. Both compounds had an aromatic ring as a common feature. These served as active ligand for the training set.*

#### **3.5 Docking validation**

The docking system was validated using the five inhibitors and 250 decoys generated with DUD-E and further used their binding energies to plot an AUC with a value of 0.702 (**Figure A5-B**). AUC value of 1.0 verifies that the prediction of hits obtained from the hypothesis is perfect whereas values of 0.5 and less than 0.70 imply average and moderate random selection respectively [48].

Furthermore, the validation of the molecular docking was also undertaken by aligning the re-docked ligands with their respective co-crystallized complexes taken from PDB. The RMSD values of the alignment between the re-docked ligands and the co-crystalized ligands in complexes of 2YAU, 4APN, 5EBK, 6ER5, 6I7N were 1.483, 3.020, 1.920, 2.712, and 2.465 Å, respectively. Only two of the RMSD values (5EBK and 2YAU) of the alignments were below 2.0 Å, which is considered the threshold for good alignment.

Superimposition also validated the accuracy of docking at the predicted active site. The FAD molecule from the template selected for modeling, 2JK6 when extracted and docked in *Lm*TR, showed similar surrounding residues in the pocket of FAD-2JK6 and FAD-*Lm*TR complex. Ligand alignment of these two complexes gave an RMSD of 3.291 Å which is above the expected threshold. Nonetheless, the FAD-*Lm*TR docking simulated a pose that showed common residues such as Ser14, Gly15, Arg287, and Thr335 taking part in hydrogen bonding in these complexes (**Figure A1**). All these verified that the docking system performed very well in docking ligands to the active site.

#### **3.6 Virtual screening of pharmacophore hits**

When docking validation was verified, molecular docking was carried out. Molecular docking predicted various conformations of each ligand in the binding site of the *Lm*TR. Compounds were selected based on their binding energies. The binding energies gave a theoretical value that relates the affinity of the ligand to the protein model. The results of the respective pharmacophore fit scores and binding energies are well documented (**Table A2**). The 42 compounds obtained were

narrowed down to 11 by considering their binding energies and pharmacophore fit scores (**Table 1**). With respect to these studies to find inhibitors of the FAD binding site for which FAD is a prosthetic group that already binds tightly to the catalytic site, it was reasonable to select the compounds with binding energies below 9.0 kcal/mol or relatively closer to 9.0 kcal/mol which can provide a plausible binding advantage when acting as inhibitors at the site. The compound ZINC95486081 had the lowest binding energy of 9.8 kcal/mol and the highest was observed from evoxine as 6.2 kcal/mol (**Table A2**).



#### **Table 1.**

*The table shows parameters involved in the selection of lead compounds. This included pharmacophore fit score, the binding energy and the number of hydrogen and hydrophobic bond interacting residues.*

#### **3.7 Protein–ligand interaction**

Molecular interaction studies are important for understanding the mechanism of biological regulation at the molecular level and as such also provides a theoretical basis for drug design and discovery [49, 50]. Hydrogen and hydrophobic interactions are key players in stabilizing energetically favored ligands, in an open conformational environment of protein structures [29]. The intermolecular interaction and bond lengths of these 11 compounds were observed. The compound which showed the highest binding affinity, ZINC95486081 formed two hydrogen bonds with residues Lys60 and Ser178 with respective bond lengths of 2.92 Å and 2.87 Å. Five compounds including MTPA, Karatavicinol, Marmin, Betaxanthin and ZINC38658035 had four residues as the highest number of residues partaking in hydrogen bonding. MTPA formed hydrogen bonds with residues Thr51, Thr293, Asp327, and Ser14. Karatavicinol on the other hand formed hydrogen bonds with Thr51, Arg287, Thr293, and Asp327 (**Figure A6-A**). Marmin formed hydrogen bonds with Val34, Thr51, Thr160, and Thr335 (**Figure A6-B**). Betaxanthin also bonded with Ser14, Cys52, Gly127, and Thr335. Finally, ZINC38658035 formed hydrogen bonds with Tyr198, Val362, Tyr374, and Gly376. The compound 13 hydroxyfeselol was the only hit that formed three hydrogen bonds with Val362, Thr374, and Gly376. Taccalin and Pectachol formed only two hydrogen bonds with Ser14 and Ala365. On the other hand, Colladonin and Feselol formed the least hydrogen bond residues with Asn330. The shortest bond length of 2.83 Å was

exhibited by Marmin with Thr160. Betaxanthin showed the highest pharmacophore fit score of 57.51 followed by Pectachol (57.18). 13-Hydroxyfeselol showed the lowest fit score of 55.62.
