Identification of Novel Potential Inhibitors of Pteridine Reductase 1 in Trypanosoma brucei via Computational Structure-Based Approaches and in Vitro Inhibition Assays

Pteridine reductase 1 (PTR1) is a trypanosomatid multifunctional enzyme that provides a mechanism for escape of dihydrofolate reductase (DHFR) inhibition. This is because PTR1 can reduce pterins and folates. Trypanosomes require folates and pterins for survival and are unable to synthesize them de novo. Currently there are no anti-folate based Human African Trypanosomiasis (HAT) chemotherapeutics in use. Thus, successful dual inhibition of Trypanosoma brucei dihydrofolate reductase (TbDHFR) and Trypanosoma brucei pteridine reductase 1 (TbPTR1) has implications in the exploitation of anti-folates. We carried out molecular docking of a ligand library of 5742 compounds against TbPTR1 and identified 18 compounds showing promising binding modes. The protein-ligand complexes were subjected to molecular dynamics to characterize their molecular interactions and energetics, followed by in vitro testing. In this study, we identified five compounds which showed low micromolar Trypanosome growth inhibition in in vitro experiments that might be acting by inhibition of TbPTR1. Compounds RUBi004, RUBi007, RUBi014, and RUBi018 displayed moderate to strong antagonism (mutual reduction in potency) when used in combination with the known TbDHFR inhibitor, WR99210. This gave an indication that the compounds might inhibit both TbPTR1 and TbDHFR. RUBi016 showed an additive effect in the isobologram assay. Overall, our results provide a basis for scaffold optimization for further studies in the development of HAT anti-folates.


Introduction
African trypanosomes are flagellated hemo-parasites, transmitted by Tsetse flies, and cause zoonotic infection in mammalian hosts [1]. In animals the disease is known as Nagana while in humans it is known as Human African Trypanosomiasis (HAT) [2,3]. Acute HAT is caused by Trypanosoma brucei rhodesiense (Tbr) while chronic HAT is caused by Trypanosoma brucei gambiense (Tbg).
humans it is known as Human African Trypanosomiasis (HAT) [2,3]. Acute HAT is caused by Trypanosoma brucei rhodesiense (Tbr) while chronic HAT is caused by Trypanosoma brucei gambiense (Tbg). This neglected tropical disease (NTD) remains of considerable public health and animal production concern [4,5].
Trypanosomes are unable to synthesize folates and pterins de novo [6]. Reduced folate and pterin cofactors are essential for parasite survival where they are critical in pathways such as protein and nucleic acid biosynthesis [7]. In order to survive, trypanosomes scavenge extracellular folate and pterin precursors from their hosts [8,9]. Hence, the pathway is an interesting drug target. Drugs targeting the folate pathway have been used in the treatment of several infections, most notably in the treatment of bacterial and malarial infections [10]. However, their use in the treatment and management of HAT has not been successful to date.
The key enzymes involved in trypanosome folate metabolism are dihydrofolate reductase (DHFR) and pteridine reductase 1 (PTR1) (Figure 1) [11][12][13]. DHFR (EC 1.5.1.3) is an NADPHdependent enzyme that catalyzes reduction of folate to dihydrofolate (H2F), and H2F to tetrahydrofolate (H4F) (Figure 1) [11,12]. Folate is essentially a pteridine that has been conjugated to p-aminobenzoic acid (pABA) that is glutamylated (Figure 1) [14]. DHFR is a validated and primary target of most anti-folate drugs [12]. However, the use of traditional anti-folates against DHFR in trypanosomatids such as Trypanosoma and Leishmania has been largely unsuccessful [12,13,15,16]. PTR1 (EC 1.5.1.33), which is a short-chain dehydrogenase reductase family member and an NADPH-dependent enzyme, is unique to trypanosomatids [8]. It is important in the reduction of biopterin to dihydrobiopterin (H2B), and of H2B to tetrahydrobiopterin (H4B) (Figure 1). PTR1 also reduces folate to H2F, and H2F to H4F ( Figure 1) [8]. In trypanosomatids, PTR1, which is less susceptible to traditional anti-folate inhibition, contributes about 10% to total folate metabolism [13]. It is important to note that studies have shown that under DHFR inhibition PTR1 is over-expressed, thus promoting anti-folate resistance in Leishmania major and Trypanosoma cruzi [8,13,15,16]. This has been proposed as the key mechanism by which trypanosomatids are able to resist anti-folates targeting DHFR [8,13,15,16]. Gene knock down and knock out studies in T. brucei have shown that PTR1 is essential for parasite survival. As such, its inhibition alone might be sufficient to negatively impact parasite survival [17,18]. PTR1 (EC 1.5.1.33), which is a short-chain dehydrogenase reductase family member and an NADPH-dependent enzyme, is unique to trypanosomatids [8]. It is important in the reduction of biopterin to dihydrobiopterin (H 2 B), and of H 2 B to tetrahydrobiopterin (H 4 B) ( Figure 1). PTR1 also reduces folate to H 2 F, and H 2 F to H 4 F ( Figure 1) [8]. In trypanosomatids, PTR1, which is less susceptible to traditional anti-folate inhibition, contributes about 10% to total folate metabolism [13]. It is important to note that studies have shown that under DHFR inhibition PTR1 is over-expressed, thus promoting anti-folate resistance in Leishmania major and Trypanosoma cruzi [8,13,15,16]. This has been proposed as the key mechanism by which trypanosomatids are able to resist anti-folates targeting DHFR [8,13,15,16]. Gene knock down and knock out studies in T. brucei have shown that PTR1 is essential for parasite survival. As such, its inhibition alone might be sufficient to negatively impact parasite survival [17,18].
There are several studies that have reported successful combination of PTR1 and DHFR inhibitors in order to achieve synergistic inhibition of the trypanosomatid folate pathway in T. cruzi, L. major and T. brucei [18][19][20][21][22]. However, the identification of a single inhibitor motif that can target both enzymes has remained largely elusive. It has been hampered by poor selectivity against human DHFR as has been the case with PTR1 inhibitors that contain functional groups derived from DHFR inhibitors, such as 2,4 diaminoquinazoline, 2,4 diaminopteridine, or 2,4 diaminopyrimidine moieties [18,20,21]. Further, the current drugs used to treat HAT are old, toxic, and reducing in efficacy due to resistance [23,24]. A recent development in African HAT chemotherapy is the promising oral drug fexinidazole that is currently in clinical testing for the treatment of late stage chronic HAT (Tb gambiense) [25].
In this study, we sought to identify novel T. brucei PTR1 (TbPTR1) inhibitors that can be used in conjunction with known DHFR inhibitors or single inhibitors that target both enzymes with minimal human toxicity. Here, we performed structure-based virtual screening of 5742 small ligand molecules against TbPTR1 and its orthologs from T. cruzi (TcPTR1), L. major (LmPTR1) and H. sapiens (HsDHRS4). In silico docking experiments identified 18 compounds preferentially bound to the trypanosomatid PTR1s and not the human DHRS4 ortholog. These promising 18 potential hits which complexed with TbPTR1 were then subjected to molecular dynamics (MD) simulations, molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) binding free energy calculations, and dynamic residue network (DRN) analysis. Based on their computational binding modes, selectivity, dynamic stability during MD simulations, DRN analysis, free energy of binding, and commercial availability, 13 compounds were subsequently subjected to a blood stream form (BSF) T. brucei in vitro inhibition assay and an H. sapiens in vitro cytotoxicity assay. Five out of the 13 compounds, named RUBi004, RUBi007, RUBi014, RUBi016 and RUBi018, exhibited anti-trypanosomal activities against trypanosomes in culture with IC 50 values of 9.6 ± 3.2 µM, 34.9 ± 17.1 µM, 14.6 ± 9.9 µM, 25.4 ± 4.7 µM, and 12.7 ± 3.7 µM, respectively. Compounds RUBi007, RUBi016, and RUBi018 showed no significant human cell cytotoxicity at 100 µM while RUBi004 and RUBi014 had cytotoxicity IC 50 s of 23.6 ± 5.8 µM and 32.9 ± 2.2 µM, respectively.
Compounds RUBi004, RUBi007, RUBi014, and RUBi018 showed reduced inhibition when used in combination with a known DHFR inhibitor (WR99210), which may be suggestive of competitive inhibition of TbDHFR confounded by the upregulation of TbPTR1 expression as a result of simultaneous inhibition. Compound RUBi016 showed an additive effect when used in combination with WR99210, suggesting that it preferentially inhibits TbPTR1 and the addition of WR99210 further contributes to the reduction of available reduced folates, resulting in reduced parasite viability. In summary, the current study reports five compounds with inhibitory activities at low µM levels, and their scaffolds may be further optimized to design safe and effective HAT chemotherapeutics targeting the folate pathway.

Results
Please note that unless otherwise indicated, TbPTR1 numbering is used throughout the article.

Overview of PTR1 Structure and Conservation
PTR1 is a short-chain dehydrogenase/reductase (SDR) family member that has the ability to catalyze the NADPH-dependent two-stage reduction of biopterins to their 7,8-dihydro and 5,6,7,8-tetrahydro forms as well as folates to their H 2 F and H 4 F forms [8] (Figure 1). PTR1 is a homotetramer with four equivalent active sites and four NADPH binding sites (Figure 2A). Each single α/β-domain subunit is constructed around an NADPH binding Rossman-fold repeat that is composed of seven parallel β-sheets that are between three α-helices on either side (Figure 2A) [16].
Multiple sequence alignment (MSA) of the TbPTR1, TcPTR1, LmPTR1, and HsDHRS4 ortholog protein sequences showed conservation of several key residues within the SDR family signature (residues ASP161-ALA193) as well as in the substrate binding loop (residues SER207-GLU215) present in trypanosomatids ( Figure 2B). Pterin and folate substrates, along with inhibitors, interact with PTR1 complexes quite similarly, often via binding in a π-sandwich between the NADPH nicotinamide ring and residue PHE97 [19,26]. The NADPH cofactor is known to be essential in creating both the substrate binding site as well as the catalytic center [19,27]. ARG14, SER95, PHE97, ASP161, and TYR174 are important residues that interact with the folate and pterin substrates, and are well conserved among the trypanosomatids [19,26] (Figure 2B). In T. cruzi, TcPTR1 and TcPTR2 are isoforms that show very high sequence homology but also display varied enzymatic activity [28]. TcPTR1 in comparison to TcPTR2 shows higher activity with biopterin and folate than with H 2 F or H 2 B [28]. TcPTR1 has no crystal structure, so for this study a structure was calculated using homology modeling with TcPTR2 as the template [29].

Figure 2.
A cartoon representation of the Trypanosoma brucei pteridine reductase 1 (TbPTR1) protein structure (PDB: 2X9N) and a multiple sequence alignment of T. brucei PTR1, T. cruzi PTR1, T. cruzi PTR2, L. major PTR1, and H. sapiens dehydrogenase/reductase (SDR family) member four (DHRS4). (A) The protein is colored by chain, with the NADPH cofactor colored blue and the co-crystallized ligand cyromazine colored orange. TbPTR1 is a tetramer and the α/β single domain subunit (chain A) is shown in green. The substrate binding loop is colored red and was composed of residues SER207-GLU215, while the SDR family signature, which is colored brown, was composed of residues ASP161-ALA193. (B) The multiple sequence alignment (MSA) showed notable conservation in the SDR family signature, as shown by the sequence logo of the extracted motif. The MSA also showed that within the substrate binding loop there was a four residue deletion that was present only among the trypanosomatids.

Eighteen Potential Hits Out of 5742 Compounds are Identified via Virtual Screening
TbPTR1 (PDB: 2X9N) [26], LmPTR1 (PDB: 1E92) [16], HsDHRS4 (PDB: 304R), and a homology model of TcPTR1 were used in in silico virtual screening investigations. All structures included the NADPH cofactor, which is essential for arrangement of the substrate binding site and catalytic center. 5742 compounds were docked against four proteins using Autodock Vina, as is described in the Materials and Methods section, in order to identify potential hits. Top compounds were selected Figure 2. A cartoon representation of the Trypanosoma brucei pteridine reductase 1 (TbPTR1) protein structure (PDB: 2X9N) and a multiple sequence alignment of T. brucei PTR1, T. cruzi PTR1, T. cruzi PTR2, L. major PTR1, and H. sapiens dehydrogenase/reductase (SDR family) member four (DHRS4). (A) The protein is colored by chain, with the NADPH cofactor colored blue and the co-crystallized ligand cyromazine colored orange. TbPTR1 is a tetramer and the α/β single domain subunit (chain A) is shown in green. The substrate binding loop is colored red and was composed of residues SER207-GLU215, while the SDR family signature, which is colored brown, was composed of residues ASP161-ALA193. (B) The multiple sequence alignment (MSA) showed notable conservation in the SDR family signature, as shown by the sequence logo of the extracted motif. The MSA also showed that within the substrate binding loop there was a four residue deletion that was present only among the trypanosomatids.

Eighteen Potential Hits Out of 5742 Compounds are Identified via Virtual Screening
TbPTR1 (PDB: 2X9N) [26], LmPTR1 (PDB: 1E92) [16], HsDHRS4 (PDB: 304R), and a homology model of TcPTR1 were used in in silico virtual screening investigations. All structures included the NADPH cofactor, which is essential for arrangement of the substrate binding site and catalytic center. 5742 compounds were docked against four proteins using Autodock Vina, as is described in the Materials and Methods section, in order to identify potential hits. Top compounds were selected based on their Autodock Vina energy score of <−8.0 kcal/mol and their hydrogen bonding profiles. Eighteen compounds showed good selectivity for trypanosomatid PTR1 as presented in Table 1. Out of 18 compounds, only RUBi006 bound to the HsDHRS4 active site but with a weaker binding energy than the trypanosomatids. The docked complexes were analyzed using PyMOL [30] and Discovery Studio [31]. The docking energy scores were also evaluated using Xscore [32]. A summary of the compounds with the top binding modes and their corresponding energies are shown in Table S1. The top binding modes involved ligand interactions with residues that are known to be of catalytic importance, i.e., ARG14, SER95, PHE97, ASP161, and TYR174 (Table S1) [16,19]. Residues ARG14, SER95, PHE97, and ASP161 were conserved among all the trypanosomatids, while TYR174 was conserved in all the PTR1 orthologs ( Figure 2). Furthermore, residues ASP161 and TYR174 are located within the SDR family signature, which is important in the NADPH cofactor and substrate binding ( Figure 2). The IUPAC names of the ligands are shown in Table S2. based on their Autodock Vina energy score of <−8.0 kcal/mol and their hydrogen bonding profiles. Eighteen compounds showed good selectivity for trypanosomatid PTR1 as presented in Table 1. Out of 18 compounds, only RUBi006 bound to the HsDHRS4 active site but with a weaker binding energy than the trypanosomatids. The docked complexes were analyzed using PyMOL [30] and Discovery Studio [31]. The docking energy scores were also evaluated using Xscore [32]. A summary of the compounds with the top binding modes and their corresponding energies are shown in Table S1. The top binding modes involved ligand interactions with residues that are known to be of catalytic importance, i.e., ARG14, SER95, PHE97, ASP161, and TYR174 (Table S1) [16,19]. Residues ARG14, SER95, PHE97, and ASP161 were conserved among all the trypanosomatids, while TYR174 was conserved in all the PTR1 orthologs ( Figure 2). Furthermore, residues ASP161 and TYR174 are located within the SDR family signature, which is important in the NADPH cofactor and substrate binding ( Figure 2). The IUPAC names of the ligands are shown in Table S2. based on their Autodock Vina energy score of <−8.0 kcal/mol and their hydrogen bonding profiles. Eighteen compounds showed good selectivity for trypanosomatid PTR1 as presented in Table 1. Out of 18 compounds, only RUBi006 bound to the HsDHRS4 active site but with a weaker binding energy than the trypanosomatids. The docked complexes were analyzed using PyMOL [30] and Discovery Studio [31]. The docking energy scores were also evaluated using Xscore [32]. A summary of the compounds with the top binding modes and their corresponding energies are shown in Table S1. The top binding modes involved ligand interactions with residues that are known to be of catalytic importance, i.e., ARG14, SER95, PHE97, ASP161, and TYR174 (Table S1) [16,19]. Residues ARG14, SER95, PHE97, and ASP161 were conserved among all the trypanosomatids, while TYR174 was conserved in all the PTR1 orthologs ( Figure 2). Furthermore, residues ASP161 and TYR174 are located within the SDR family signature, which is important in the NADPH cofactor and substrate binding ( Figure 2). The IUPAC names of the ligands are shown in Table S2. based on their Autodock Vina energy score of <−8.0 kcal/mol and their hydrogen bonding profiles. Eighteen compounds showed good selectivity for trypanosomatid PTR1 as presented in Table 1. Out of 18 compounds, only RUBi006 bound to the HsDHRS4 active site but with a weaker binding energy than the trypanosomatids. The docked complexes were analyzed using PyMOL [30] and Discovery Studio [31]. The docking energy scores were also evaluated using Xscore [32]. A summary of the compounds with the top binding modes and their corresponding energies are shown in Table S1. The top binding modes involved ligand interactions with residues that are known to be of catalytic importance, i.e., ARG14, SER95, PHE97, ASP161, and TYR174 (Table S1) [16,19]. Residues ARG14, SER95, PHE97, and ASP161 were conserved among all the trypanosomatids, while TYR174 was conserved in all the PTR1 orthologs ( Figure 2). Furthermore, residues ASP161 and TYR174 are located within the SDR family signature, which is important in the NADPH cofactor and substrate binding ( Figure 2). The IUPAC names of the ligands are shown in Table S2. based on their Autodock Vina energy score of <−8.0 kcal/mol and their hydrogen bonding profiles. Eighteen compounds showed good selectivity for trypanosomatid PTR1 as presented in Table 1. Out of 18 compounds, only RUBi006 bound to the HsDHRS4 active site but with a weaker binding energy than the trypanosomatids. The docked complexes were analyzed using PyMOL [30] and Discovery Studio [31]. The docking energy scores were also evaluated using Xscore [32]. A summary of the compounds with the top binding modes and their corresponding energies are shown in Table S1. The top binding modes involved ligand interactions with residues that are known to be of catalytic importance, i.e., ARG14, SER95, PHE97, ASP161, and TYR174 (Table S1) [16,19]. Residues ARG14, SER95, PHE97, and ASP161 were conserved among all the trypanosomatids, while TYR174 was conserved in all the PTR1 orthologs ( Figure 2). Furthermore, residues ASP161 and TYR174 are located within the SDR family signature, which is important in the NADPH cofactor and substrate binding ( Figure 2). The IUPAC names of the ligands are shown in Table S2. based on their Autodock Vina energy score of <−8.0 kcal/mol and their hydrogen bonding profiles. Eighteen compounds showed good selectivity for trypanosomatid PTR1 as presented in Table 1. Out of 18 compounds, only RUBi006 bound to the HsDHRS4 active site but with a weaker binding energy than the trypanosomatids. The docked complexes were analyzed using PyMOL [30] and Discovery Studio [31]. The docking energy scores were also evaluated using Xscore [32]. A summary of the compounds with the top binding modes and their corresponding energies are shown in Table S1. The top binding modes involved ligand interactions with residues that are known to be of catalytic importance, i.e., ARG14, SER95, PHE97, ASP161, and TYR174 (Table S1) [16,19]. Residues ARG14, SER95, PHE97, and ASP161 were conserved among all the trypanosomatids, while TYR174 was conserved in all the PTR1 orthologs ( Figure 2). Furthermore, residues ASP161 and TYR174 are located within the SDR family signature, which is important in the NADPH cofactor and substrate binding ( Figure 2). The IUPAC names of the ligands are shown in Table S2. based on their Autodock Vina energy score of <−8.0 kcal/mol and their hydrogen bonding profiles. Eighteen compounds showed good selectivity for trypanosomatid PTR1 as presented in Table 1. Out of 18 compounds, only RUBi006 bound to the HsDHRS4 active site but with a weaker binding energy than the trypanosomatids. The docked complexes were analyzed using PyMOL [30] and Discovery Studio [31]. The docking energy scores were also evaluated using Xscore [32]. A summary of the compounds with the top binding modes and their corresponding energies are shown in Table S1. The top binding modes involved ligand interactions with residues that are known to be of catalytic importance, i.e., ARG14, SER95, PHE97, ASP161, and TYR174 (Table S1) [16,19]. Residues ARG14, SER95, PHE97, and ASP161 were conserved among all the trypanosomatids, while TYR174 was conserved in all the PTR1 orthologs ( Figure 2). Furthermore, residues ASP161 and TYR174 are located within the SDR family signature, which is important in the NADPH cofactor and substrate binding ( Figure 2). The IUPAC names of the ligands are shown in Table S2. based on their Autodock Vina energy score of <−8.0 kcal/mol and their hydrogen bonding profiles. Eighteen compounds showed good selectivity for trypanosomatid PTR1 as presented in Table 1. Out of 18 compounds, only RUBi006 bound to the HsDHRS4 active site but with a weaker binding energy than the trypanosomatids. The docked complexes were analyzed using PyMOL [30] and Discovery Studio [31]. The docking energy scores were also evaluated using Xscore [32]. A summary of the compounds with the top binding modes and their corresponding energies are shown in Table S1. The top binding modes involved ligand interactions with residues that are known to be of catalytic importance, i.e., ARG14, SER95, PHE97, ASP161, and TYR174 (Table S1) [16,19]. Residues ARG14, SER95, PHE97, and ASP161 were conserved among all the trypanosomatids, while TYR174 was conserved in all the PTR1 orthologs ( Figure 2). Furthermore, residues ASP161 and TYR174 are located within the SDR family signature, which is important in the NADPH cofactor and substrate binding ( Figure 2). The IUPAC names of the ligands are shown in Table S2. based on their Autodock Vina energy score of <−8.0 kcal/mol and their hydrogen bonding profiles. Eighteen compounds showed good selectivity for trypanosomatid PTR1 as presented in Table 1. Out of 18 compounds, only RUBi006 bound to the HsDHRS4 active site but with a weaker binding energy than the trypanosomatids. The docked complexes were analyzed using PyMOL [30] and Discovery Studio [31]. The docking energy scores were also evaluated using Xscore [32]. A summary of the compounds with the top binding modes and their corresponding energies are shown in Table S1. The top binding modes involved ligand interactions with residues that are known to be of catalytic importance, i.e., ARG14, SER95, PHE97, ASP161, and TYR174 (Table S1) [16,19]. Residues ARG14, SER95, PHE97, and ASP161 were conserved among all the trypanosomatids, while TYR174 was conserved in all the PTR1 orthologs ( Figure 2). Furthermore, residues ASP161 and TYR174 are located within the SDR family signature, which is important in the NADPH cofactor and substrate binding ( Figure 2). The IUPAC names of the ligands are shown in Table S2.
Molecules 2018, 23, x FOR PEER REVIEW 5 of 26 based on their Autodock Vina energy score of <−8.0 kcal/mol and their hydrogen bonding profiles. Eighteen compounds showed good selectivity for trypanosomatid PTR1 as presented in Table 1. Out of 18 compounds, only RUBi006 bound to the HsDHRS4 active site but with a weaker binding energy than the trypanosomatids. The docked complexes were analyzed using PyMOL [30] and Discovery Studio [31]. The docking energy scores were also evaluated using Xscore [32]. A summary of the compounds with the top binding modes and their corresponding energies are shown in Table S1. The top binding modes involved ligand interactions with residues that are known to be of catalytic importance, i.e., ARG14, SER95, PHE97, ASP161, and TYR174 (Table S1) [16,19]. Residues ARG14, SER95, PHE97, and ASP161 were conserved among all the trypanosomatids, while TYR174 was conserved in all the PTR1 orthologs ( Figure 2). Furthermore, residues ASP161 and TYR174 are located within the SDR family signature, which is important in the NADPH cofactor and substrate binding ( Figure 2). The IUPAC names of the ligands are shown in Table S2.

Drug Likeness
A principal component analysis (PCA) of the compounds based on their molecular descriptors, as is listed in the Materials and Methods section, showed that the top docking compounds clustered well with known Food and Drug Administration (FDA)-approved central nervous system (CNS)permeable drugs (Figure 3) occupying the same chemical space as FDA-approved drugs as well as FDA-approved CNS-permeable drugs. The compounds are of 'drug-like' desirability and are likely to cross the blood-brain barrier (BBB), making them good candidates for HAT chemotherapeutics [33].
RUBi006 was the only compound that bound to the human ortholog.

Drug Likeness
A principal component analysis (PCA) of the compounds based on their molecular descriptors, as is listed in the Materials and Methods section, showed that the top docking compounds clustered well with known Food and Drug Administration (FDA)-approved central nervous system (CNS)permeable drugs (Figure 3) occupying the same chemical space as FDA-approved drugs as well as FDA-approved CNS-permeable drugs. The compounds are of 'drug-like' desirability and are likely to cross the blood-brain barrier (BBB), making them good candidates for HAT chemotherapeutics [33].
RUBi006 was the only compound that bound to the human ortholog.

Drug Likeness
A principal component analysis (PCA) of the compounds based on their molecular descriptors, as is listed in the Materials and Methods section, showed that the top docking compounds clustered well with known Food and Drug Administration (FDA)-approved central nervous system (CNS)permeable drugs (Figure 3) occupying the same chemical space as FDA-approved drugs as well as FDA-approved CNS-permeable drugs. The compounds are of 'drug-like' desirability and are likely to cross the blood-brain barrier (BBB), making them good candidates for HAT chemotherapeutics [33].
RUBi006 was the only compound that bound to the human ortholog.

Drug Likeness
A principal component analysis (PCA) of the compounds based on their molecular descriptors, as is listed in the Materials and Methods section, showed that the top docking compounds clustered well with known Food and Drug Administration (FDA)-approved central nervous system (CNS)permeable drugs (Figure 3) occupying the same chemical space as FDA-approved drugs as well as FDA-approved CNS-permeable drugs. The compounds are of 'drug-like' desirability and are likely to cross the blood-brain barrier (BBB), making them good candidates for HAT chemotherapeutics [33].
RUBi006 was the only compound that bound to the human ortholog.

Drug Likeness
A principal component analysis (PCA) of the compounds based on their molecular descriptors, as is listed in the Materials and Methods section, showed that the top docking compounds clustered well with known Food and Drug Administration (FDA)-approved central nervous system (CNS)permeable drugs (Figure 3) occupying the same chemical space as FDA-approved drugs as well as FDA-approved CNS-permeable drugs. The compounds are of 'drug-like' desirability and are likely to cross the blood-brain barrier (BBB), making them good candidates for HAT chemotherapeutics [33].
RUBi006 was the only compound that bound to the human ortholog.

Drug Likeness
A principal component analysis (PCA) of the compounds based on their molecular descriptors, as is listed in the Materials and Methods section, showed that the top docking compounds clustered well with known Food and Drug Administration (FDA)-approved central nervous system (CNS)-permeable drugs (Figure 3) occupying the same chemical space as FDA-approved drugs as well as FDA-approved CNS-permeable drugs. The compounds are of 'drug-like' desirability and are likely to cross the blood-brain barrier (BBB), making them good candidates for HAT chemotherapeutics [33].

Five Hit Compounds Show Anti-Trypanosomal Activity in Vitro
As a next step, a total of 18 TbPTR1-ligand complexes were subjected to 200 ns all-atom MD simulations followed by MM-PBSA free energy calculations. All the compounds showed linear stable MD trajectories as observed in RMSD (Figures S1-S7), radius of gyration (Rg) calculations (Figures

Five Hit Compounds Show Anti-Trypanosomal Activity In Vitro
As a next step, a total of 18 TbPTR1-ligand complexes were subjected to 200 ns all-atom MD simulations followed by MM-PBSA free energy calculations. All the compounds showed linear stable MD trajectories as observed in RMSD ( Figures S1-S7), radius of gyration (Rg) calculations ( Figures S1-S6 and S8), and promising hydrogen bonding features ( Figures S1-S6 and S9). Root mean square fluctuation (RMSF) values (Figures S1-S6 and S10) and binding free energies ( Figure S11) were also calculated. PCA was also carried out to investigate the overall dynamics of the protein systems [34]. Overall, in silico analysis indicated that all 18 compounds could be further studied for in vitro analysis. However, only 13 of these were commercially available and tested for anti-trypanosomal activity against T. brucei BSF in culture. Compounds RUBi003, RUBi006, RUBi009, RUBi013, and RUBi017 could not be purchased, and as such were not used in the in vitro inhibition assays, even though they showed similar binding modes to folate, pterins, and known TbPTR1 inhibitors. In spite of good binding modes and stable MD trajectories, compounds RUBi001, RUBi002, RUBi005, RUBi008, RUBi010, RUBi011, and RUBi012 did not show anti-trypanosomal activity in in vitro experiments when used at 20 µM. The five compounds shown in Figure 4 that showed significant inhibition of T. brucei viability at this concentration were subjected to dose-response assays to derive their IC 50 values against T. brucei ( Figure 5), and are discussed further below.  RUBi016 inhibited blood stream form trypanosome growth in culture with an IC50 value of 25.4 ± 4.7 μM ( Figure 5) and it was the only compound which showed an additive effect when used in combination with WR99210 ( Figure 6A). RUBi004, RUBi007, RUBi014, and RUBi018 inhibited BSF trypanosome growth in culture with IC50 values of 9.6 ± 3.2 μM, 34.9 ± 17.1 μM, 14.6 ± 9.9 μM, and 12.7 ± 3.7 μM, respectively ( Figure 5). When used in combination with the known DHFR inhibitor WR99210, these compounds showed antagonism, with RUBi018 to the least and RUBi014 to the highest extent ( Figure 6).
TbPTR1 and TbDHFR both catalyze the reduction of folate to H2F and H4F; the substrate binding sites, however, are differently ordered [35]. Additionally, TbDHFR is capable of undergoing significant conformational changes when in complex with thymidylate synthase (TS), while TbPTR1 is more rigid [36][37][38]. All the compounds appear to bind in similar patterns to pterin, folate, and known TbPTR1 inhibitors interacting with the NADPH cofactor and the TbPTR1 protein [19] (Figure  7). As a next step, all the hit compounds were docked against TbDHFR active sites to see potential binding activity. RUBi004 ( Figure 8A), RUBi007 ( Figure 8B), RUBi014 ( Figure 8C), and RUBi018  Figure 5) and it was the only compound which showed an additive effect when used in combination with WR99210 ( Figure 6A). RUBi004, RUBi007, RUBi014, and RUBi018 inhibited BSF trypanosome growth in culture with IC 50 values of 9.6 ± 3.2 µM, 34.9 ± 17.1 µM, 14.6 ± 9.9 µM, and 12.7 ± 3.7 µM, respectively ( Figure 5). When used in combination with the known DHFR inhibitor WR99210, these compounds showed antagonism, with RUBi018 to the least and RUBi014 to the highest extent ( Figure 6).
The use of TbPTR1 inhibitors in combination with TbDHFR inhibitors has long been proposed as a viable avenue for the generation of new anti-trypanosomal anti-folate drugs [19]. It is important to note, however, that studies have shown that under DHFR inhibition, PTR1 is over-expressed, thus promoting anti-folate resistance in Leishmania major and Trypanosoma cruzi [8,13,15,16]. This likely contributes to the antagonism of the RUBi compounds' activity by WR99210, an effect that might be exacerbated by the additional possible inhibition of TbDHFR by RUBi004, RUBi007, RUBi014, and RUBi018. Figure 5. Dose-response plots to derive IC50 values for RUBi004 (red), RUBi014 (orange), RUBi007 (green), RUBi016 (blue), RUBi018 (purple), WR99210 (maroon), and pentamidine (black). Pentamidine was used as the positive control anti-trypanosomal compound. Parasite % viability was determined using the resazurin method. Error bars indicate the standard deviation of % viability obtained in three replicate wells. The above graph represents one of three independent experiments that were performed on separate occasions. The combined IC50 values (± standard deviation) of RUBi004, RUBi007, RUBi014, RUBi016, and RUBi018 were determined to be 9.6 ± 3.2 μM, 34.9 ± 17.1 μM, 14.6 ± 9.9 μM, 25.4 ± 4.7 μM, and 12.7 ± 3.7 μM, respectively. The IC50 values of WR99210 and pentamidine were determined to be 0.5 ± 0.4 μM and 0.014 μM, respectively. (green), RUBi016 (blue), RUBi018 (purple), WR99210 (maroon), and pentamidine (black). Pentamidine was used as the positive control anti-trypanosomal compound. Parasite % viability was determined using the resazurin method. Error bars indicate the standard deviation of % viability obtained in three replicate wells. The above graph represents one of three independent experiments that were performed on separate occasions. The combined IC 50 values (± standard deviation) of RUBi004, RUBi007, RUBi014, RUBi016, and RUBi018 were determined to be 9.6 ± 3.2 µM, 34.9 ± 17.1 µM, 14.6 ± 9.9 µM, 25.4 ± 4.7 µM, and 12.7 ± 3.7 µM, respectively. The IC 50 values of WR99210 and pentamidine were determined to be 0.5 ± 0.4 µM and 0.014 µM, respectively. as a viable avenue for the generation of new anti-trypanosomal anti-folate drugs [19]. It is important to note, however, that studies have shown that under DHFR inhibition, PTR1 is over-expressed, thus promoting anti-folate resistance in Leishmania major and Trypanosoma cruzi [8,13,15,16]. This likely contributes to the antagonism of the RUBi compounds' activity by WR99210, an effect that might be exacerbated by the additional possible inhibition of TbDHFR by RUBi004, RUBi007, RUBi014, and RUBi018. Figure 5. Dose-response plots to derive IC50 values for RUBi004 (red), RUBi014 (orange), RUBi007 (green), RUBi016 (blue), RUBi018 (purple), WR99210 (maroon), and pentamidine (black). Pentamidine was used as the positive control anti-trypanosomal compound. Parasite % viability was determined using the resazurin method. Error bars indicate the standard deviation of % viability obtained in three replicate wells. The above graph represents one of three independent experiments that were performed on separate occasions. The combined IC50 values (± standard deviation) of RUBi004, RUBi007, RUBi014, RUBi016, and RUBi018 were determined to be 9.6 ± 3.2 μM, 34.9 ± 17.1 μM, 14.6 ± 9.9 μM, 25.4 ± 4.7 μM, and 12.7 ± 3.7 μM, respectively. The IC50 values of WR99210 and pentamidine were determined to be 0.5 ± 0.4 μM and 0.014 μM, respectively.    TbPTR1 and TbDHFR both catalyze the reduction of folate to H 2 F and H 4 F; the substrate binding sites, however, are differently ordered [35]. Additionally, TbDHFR is capable of undergoing significant conformational changes when in complex with thymidylate synthase (TS), while TbPTR1 is more rigid [36][37][38]. All the compounds appear to bind in similar patterns to pterin, folate, and known TbPTR1 inhibitors interacting with the NADPH cofactor and the TbPTR1 protein [19] (Figure 7). As a next step, all the hit compounds were docked against TbDHFR active sites to see potential binding activity. RUBi004 ( Figure 8A), RUBi007 ( Figure 8B), RUBi014 ( Figure 8C), and RUBi018 ( Figure 8E) molecular TbDHFR dockings showed high binding affinity towards the protein, with −9.4 kcal/mol, −9.5 kcal/mol, −8.7 kcal/mol, and −7.9 kcal/mol, respectively. RUBi016 bound to TbDHFR with a binding affinity of −7.6 kcal/mol ( Figure 8D). Interestingly, TbDHFR compound binding energy scores were in agreement with isobologram results. Three compounds (RUBi004, RUBi007, and RUBi014) with the greatest antagonist activities also had the highest binding energies towards TbDHFR. The possibility that RUBi004, RUBi007, and RUBi014 also inhibit DHFR may explain why they show more marked antagonism when used in combination with DHFR inhibitor; WR99210 and each of these three compounds both inhibit TbDHFR, thus possibly enhancing the upregulation of TbPTR1.
Compound RUBi018, with the least antagonist activity, had a similar TbDHFR binding energy to RUBi016 that showed an additive effect with WR99210. RUBi018 also displayed the highest binding energy for TbPTR1 among all five compounds (−8.4 kcal/mol; Table 1). Interestingly, with the exception of RUBi018 which was more active than predicted, the anti-trypanosomal IC 50 values of the remaining four compounds correlated with their TbPTR1 binding energies. RUBi004, with the highest affinity (−10.4 kcal/mol), displayed the most potent activity (with an IC 50 of 9.6 ± 3.2 µM).
The use of TbPTR1 inhibitors in combination with TbDHFR inhibitors has long been proposed as a viable avenue for the generation of new anti-trypanosomal anti-folate drugs [19]. It is important to note, however, that studies have shown that under DHFR inhibition, PTR1 is over-expressed, thus promoting anti-folate resistance in Leishmania major and Trypanosoma cruzi [8,13,15,16]. This likely contributes to the antagonism of the RUBi compounds' activity by WR99210, an effect that might be exacerbated by the additional possible inhibition of TbDHFR by RUBi004, RUBi007, RUBi014, and RUBi018.

Detailed in Silico Compound Analysis
Overall, our in silico rational-based drug discovery approach was able to identify five compounds that showed anti-trypanosomal in vitro activity. Compounds RUBi007, RUBi016, and RUBi018 showed no significant human cell cytotoxicity at 100 μM, while RUBi004 and RUBi014 had cytotoxicity IC50s of 23.6 ± 5.8 μM and 32.9 ± 2.2 μM, respectively ( Figure S12). RUBi014, also known as eriodictyol, is a flavonoid that has previously been reported to be selectively anti-protozoal with activity against T. brucei in culture [39,40]. RUBi018 is related to 4-aryl-2-(1-substituted ethylidene)

Detailed in Silico Compound Analysis
Overall, our in silico rational-based drug discovery approach was able to identify five compounds that showed anti-trypanosomal in vitro activity. Compounds RUBi007, RUBi016, and RUBi018 showed no significant human cell cytotoxicity at 100 µM, while RUBi004 and RUBi014 had cytotoxicity IC 50 s of 23.6 ± 5.8 µM and 32.9 ± 2.2 µM, respectively ( Figure S12). RUBi014, also known as eriodictyol, is a flavonoid that has previously been reported to be selectively anti-protozoal with activity against T. brucei in culture [39,40]. RUBi018 is related to 4-aryl-2-(1-substituted ethylidene) thiazoles that have been previously reported to show anti-bacterial activity [41]. To our knowledge, there has been no specific activity reported in the literature for the rest of the compounds.

Molecular dynamics:
The RMSD calculations revealed that the protein backbone and the NADPH cofactors in the TbPTR1-ligand complexes showed conformational changes relative to their initial structures, deviating by 0.3 nm-0.5 nm (Figures S1a, S2a, S3a, S4a, S6a and S7a). Rg analysis showed that the binding of the ligand led to increased compactness in TbPTR1-RUBi004 and TbPTR1-RUBi007 (Figures S2b and S3b). In the TbPTR1-RUBi004 complex we observed slight increases in the flexibility of residues LYS13, PHE97, TYR98, and TYR174 via RMSF calculations ( Figure S2d). Loop residues MET169 and ALA170 showed increased flexibility, while helix residues ALA188 and ALA189 showed reduced flexibility ( Figure S2d). The substrate binding loop SER207-GLU217 was stable ( Figure S2d). In the TbPTR1-RUBi007 complex we observed decreased flexibility in ARG14, PHE97, GLY205, PRO210, ALA212, and TRP221 ( Figure S3d), and increased flexibility in GLY247, SER248, and ALA249 ( Figure S3d). In the TbPTR1-RUBi014 and TbPTR1-RUBi016 complexes we observed a significant increase in the flexibility of the substrate binding loop and the small α6 helix (Figures S4d and S5d). Lastly, in the TbPTR1-RUBi018 complex there was also an increase in the flexibility of the substrate binding loop (PRO210-MET213), but to a lesser extent ( Figure S6d).
There were no large motions observed in the substrate binding loop in the TbPTR1-RUBi007 complex and the largest motions were in the α6 helix, the C-terminal residues HIS267-ALA268, the modelled loop 1, and the modelled loop 2 ( Figure 9C).
The TbPTR1-RUBi014 and TbPTR1-RUBi016 complexes showed very strong alteration of the motion of the substrate binding loop, the α6 helix, and the C-terminal residues HIS267-ALA268, as shown by the RMSF analysis and PCA ( Figure 9D,E, Figures S4d and S5d). They appear to make the active site widen and extend ( Figure 9D,E) when compared to the apo protein ( Figure 9A).
Lastly, PCA showed that the TbPTR1-RUBi018 complex had the largest motions in the substrate binding loop, the α6 helix, the CYS160-TY174 loop region, the C-terminal residues ALA268, the modelled loop 1, and the modelled loop 2 (residues LYS143-SER151) ( Figure 9F).
Hydrogen bond analysis: RUBi004 formed two hydrogen bonds at a frequency of 0.02 during 200 ns of the MD simulation with TbPTR1 residue LYS13 and the NADPH cofactor ( Figure 7A and Figure S2c). RUBi007 formed two hydrogen bonds at a frequency of 0.27 with the protein residues ARG14 and GLY205 ( Figure 7B and Figure S3c). RUBi014 formed an average of four hydrogen bonds at a frequency of 0.002 with the protein residues ASP161, ASN175, PRO204, and GLY205 ( Figure 7C and Figure S4c). RUBi016 formed three hydrogen bonds at a frequency of 0.12 with the protein during the MD simulation, and these were with ALA96, LEU208, and the NADPH cofactor ( Figure 7D and Figure S5c) while RUBi018 formed hydrogen bonds with GLY205 and TRP221 at a frequency of 0.07 ( Figure 7E and Figure S6c).
Binding free energy: All five compounds bound stably to the protein throughout the MD simulation. A summary of the free binding energies of all compounds in the study is shown in Table 2. RUBi004, RUBi007, RUBi014, RUBi016, and RUBi018 bound to TbPTR1 with free binding energies of −63 ± 14 kJ/mol, −88 ± 10 kJ/mol, −56 ± 12 kJ/mol, −23 ± 10 kJ/mol, and −88 ± 12 kJ/mol, respectively ( Table 2). The per residue energy contribution offers an interesting insight into the study compounds ( Figure 10). Many of the protein-ligand complexes had residues contributing to binding that were of catalytic importance and involved in the binding of known TbPTR1 inhibitors ( Figure 10) [16,19,26,27]. Residues LYS13, ARG14, PHE97, MET163, TYR174, and substrate binding residues LEU209-TRP221 appear to enhance ligand interactions as they contribute the most favorably energetically to ligand binding ( Figure 10A). RUBi004 and RUBi007 showed NADPH cofactor binding residues LYS13 and ARG14 contributing favorably to binding (Figure 10(Bi,ii)). Active site residue ASP161 generally had poor interactions with the ligands where it gave unfavorable energy contributions to ligand binding ( Figure 10). ASP161 forms hydrogen bonds with MET163 and TYR174 and is important in proton transfer to the substrate during catalysis [35]. Notably in RUBi016 where there wasn't any interaction with residues MET163 or TYR174 as shown in the docking analysis, ASP161 did not give an unfavorable energy contribution ( Figures 7D and 10(Biv)). This work provides insight into important TbPTR1 protein-ligand interactions that can be used in rational-based drug design to characterize potential inhibitors with the end goal of designing and optimizing HAT anti-folate drugs. Dynamic residue network analysis: Average shortest path (Average L) and average betweenness centrality (Average BC) metrics over a MD trajectory [43] were calculated for comparative DRN analysis between ligand-bound and unbound TbPTR1. These metrics were also compared to RMSF data ( Figure S13a). As shown previously by Penkler et al., a general trend between Average BC, Average L -1 and RMSF -1 has been observed [44,45]. Pearson correlation coefficient values are presented Dynamic residue network analysis: Average shortest path (Average L) and average betweenness centrality (Average BC) metrics over a MD trajectory [43] were calculated for comparative DRN analysis between ligand-bound and unbound TbPTR1. These metrics were also compared to RMSF data ( Figure S13a). As shown previously by Penkler et al., a general trend between Average BC, Average L −1 and RMSF −1 has been observed [44,45]. Pearson correlation coefficient values are presented in Table S3. Overall, only very slight changes, especially in the substrate binding loop (SER207-GLU215) and the small α6 helix (GLY214-VAL225), were observed for Average L between the apo protein and the TbPTR1-ligand complexes (Figures S1e, S2e, S3e, S4e, S5e, and S6e).
Residues THR9, GLY205, and SER207 were conserved among all the PTR1 orthologs while SER95, CYS160, MET163, ASP165, LEU208, PRO210, and SER233 were conserved among the trypanosomatid PTR1 orthologs only (Figure 2). ILE15 and ALA238 were conserved in all the PTR1 orthologs except LmPTR1, where they were replaced by a leucine and serine residue, respectively ( Figure 2). Residue VAL206 was only present in TbPTR1 while the other PTR1 orthologs had a leucine instead (Figure 2).

BC-related residue interaction analysis:
A close examination of the residues THR9, SER95, and ALA238, which were shown to have high Average BC values in the TbPTR1 apo protein, indicated that they are functionally important ( Figure S1e). Analysis of their residue interaction networks in the apo protein showed that residue THR9 formed vdW interactions with ALA11, ASN92, and ALA94, and hydrogen bonds with VAL91 and ASN93 ( Figure S14(1a)). SER95 formed hydrogen bonds with ALA96, ASN127, and the NADPH cofactor ( Figure S14(1b)). Furthermore, SER95 formed vdW interactions with TYR174 and PHE97, both of which are involved in ligand binding (Figure 7). Lastly, ALA238 formed an alkyl interaction with ALA18 that is hydrogen bonded to ILE15 ( Figure S14(1c)). These interactions with functionally important residues indicate that Average BC is helpful in identifying residues crucial to communication flow within the protein-ligand dynamic network.
In addition, the binding of ligands appears to alter the flow of information across the protein network. In the TbPTR1-RUBi004 complex, residues VAL164, SER172, and SER207, all of which were close to residues involved in ligand and NADPH cofactor binding, showed changes in Average BC ( Figures S2f and S14(2a)). SER207 formed a hydrogen bond with the NADPH cofactor, while RUBi004 formed vdW interactions with residues MET163, TYR174 and VAL206 ( Figure 7A and Figure S14(2c)).
These interactions with functionally important residues indicate that Average BC is helpful in identifying residues crucial to communication flow within the protein-ligand dynamic network.

Discussion
In this study, structure-based molecular docking was used to screen 5742 selected compounds against trypanosomatid PTR1s and a human homolog (HsDHRS4) to identify potential hits for HAT. Eighteen compounds showed good selectivity for trypanosomatid PTR1s, and only compound RUBi006 bound to the HsDHRS4 active site but with a weaker binding energy than the trypanosomatids. MD simulations, DRN calculations, and MMPBSA free energy calculations indicated that all 18 compounds were potentially good hits. Of the 18, 13 commercially available compounds were tested for anti-trypanosomal activity using in vitro inhibition assays. Five compounds out of the 13 (RUBi004, RUBi007, RUBi014, RUBi16, and RUBi018) exhibited anti-trypanosomal activity against trypanosomes in culture, with IC 50 s of 9.6 ± 3.2 µM, 34.9 ± 17.1 µM, 14.6 ± 9.9 µM, 25.4 ± 4.7 µM, and 12.7 ± 3.7 µM, respectively. RUBi007, RUBi016, and RUBi018 showed no significant human (human cervix adenocarcinoma (HeLa)) cell cytotoxicity at 100 µM (HeLa cell viability was >90% at this concentration) while RUBi004 and RUBi014 had cytotoxicity IC 50 s of 23.6 ± 5.8 µM and 32.9 ± 2.2 µM, respectively. When used in combination with WR99210, RUBi004, RUBi007, RUBi014, and RUBi018 displayed antagonistic effects, while RUBi016 showed an additive effect in the isobologram assay.
When anti-folate drugs that target trypanosomatid dihydrofolate reductase-thymidylate synthase DHFR-TS are used, PTR1 is over-expressed, allowing for a by-pass mechanism to ensure parasite survival [8,13]. This is the escape mechanism that has hampered the use of traditional anti-folates against trypanosomatids. PTR1 is an important drug target, as demonstrated by gene knock out in Leishmania [8] and knock down in Trypanosoma brucei [17,35] studies that show that the enzyme is essential for parasite survival. However, TbPTR1 is less susceptible to inhibition than TbDHFR [21,35,38]. Given the nature of the interaction between TbPTR1 and TbDHFR, a combination therapy would offer several advantages, especially against resistance problems as has been shown by anti-malarial combination treatment strategies [35,46].
Our experimental data allow us to draw the following conclusions. Compounds RUBi004, RUB007, RUBi014, and RUBi018 inhibited parasite growth with IC 50 s of 9.6 ± 3.2 µM, 34.9 ± 17.1 µM, 14.6 ± 9.9 µM, and 12.7 ± 3.7 µM when assayed on their own. When used in combination with WR99210 (with an IC 50 of 0.55 µM), which is a known TbDHFR inhibitor, each compound showed an antagonistic effect. From our molecular docking studies, we demonstrate that it is reasonable that the compounds RUBi004, RUB007, RUBi014, and RUBi018 can bind both TbPTR1 and TbDHFR with good binding affinities and in binding modes similar to those of traditional folates, pterins, and known inhibitors. Our molecular dynamics simulations also show that the ligands bind to TbPTR1 stably and with acceptable binding energies. In line with these results, we theorize that the compounds could be competing for the TbDHFR active site with WR99210, which would result in the observed antagonism. Further, the resulting over-expression of TbPTR1 would result in a further reduction in compound efficacy.
Compound RUBi016 inhibited parasite growth with an IC 50 of 25.4 ± 4.7 µM when assayed on its own. From our molecular docking results, it appears that RUBi016 binds TbDHFR with a lower binding affinity than TbPTR1, given the binding affinity values of −7.6 kcal/mol and −8.9 kcal/mol, respectively. During MD simulations it bound to TbPTR1 stably and had a good binding energy. Unlike the four compounds, RUBi016 showed an additive effect when used in combination with WR99210. We theorize that RUBi016 may be more selective for TbPTR1 compared to the other compounds, which may have a tendency to bind to both PTR1 and DHFR. The other compounds thus display antagonistic effects with WR99210, while RUBi016 has an additive effect. In the case of RUBi016, the addition of WR99210 inhibits TbDHFR, further impeding folate reduction and resulting in more inhibition of parasite survival, producing the observed additive effect. Furthermore, the fact that RUBi016 does not have the lowest IC 50 value does not exclude the possibility that it might be the most selective of the compounds. Interestingly, the ratio of its binding energy for PTR1 and DHFR is the highest (except for RUBi018 for PTR1), which means that docking scores would predict it to be the most selective. For example, RUBi004 has the lowest IC 50 value and the lowest binding energy, but it also has a very low binding energy for DHFR which means it is not really selective.
We note that compounds RUBi004 and RUBi014 contain pan-assay interference compounds (PAINS) features [47]. While this is a cause for concern, there are over 60 FDA approved drugs that have PAINS features [47,48]. The binding patterns of all five compounds are consistent with those of folates, pterins, and known inhibitors. Further analyses such as ligand structure activity relationship (SAR) analysis and protein-ligand co-crystallizations are required to validate these compounds as potential HAT anti-folate chemotherapeutics [49].

Ligand Library Preparation
The small-molecule ligands were obtained from the South African Natural Compounds database (SANCDB) [50] and the ZINC database (ZINC12) [51]. The compounds in the ZINC dataset already conformed to the Lipinski rules, and are commercially available drug-like compounds [51,52]. In the development of potential HAT chemotherapeutics, it is important to select for BBB permeability because in the second stage of the disease the parasites invade the CNS [53,54]. As such, the ligand filtering strategy was based on drug physicochemical properties that have been shown to be useful in predicting BBB permeability [55,56]. These included lipophilicity as calculated by XlogP, rotatable bonds, hydrogen donors, net charge, and molecular weight [55,56]. The values were determined in accordance with features common in CNS FDA-approved drugs and literature [56][57][58].
The ligand library was prepared by filtering 10,639,555 compounds from the ZINC Drugs Now subset [51,52] for compounds with XlogP ≤ 3, fewer than four rotatable bonds, at least two hydrogen donors, a net charge of zero, and a molecular weight ≤ 490. The ZINC subset was reduced to 5107 compounds after filtering, while the SANCDB contained 635 compounds. The final ligand library comprised of 5742 compounds.

Preparation of Protein-Ligand Complexes
A crystal structure of TbPTR1 that has a resolution of 1.15 Å was retrieved from the RCSB Protein Data Bank (PDB: 2X9N) [26]. Multiple sequence analysis using MUSCLE [59] was carried out to analyze TbPTR1 (Uniprot: O76290) and the homologues sequences, including T. cruzi (Uniprot: O44029), L. major (Uniprot: Q01782, PDB: 1E92), and the H. sapiens dehydrogenase/reductase SDR family member four (DHRS4) (Uniprot: Q9BTZ2, PDB: 3O4R). The crystal structure of TbDHFR was also retrieved and had a resolution of 2.2 Å (PDB: 3QFX) [42]. Homology modelling was done using in-house Python scripts to fix missing residues in 2X9N (residues GLN104-GLY113 and LYS143-SER151) as well as to generate a homotetramer TcPTR1 structure from its PTR2 isoform (Uniprot: Q8I814, PDB: 1MXH). The modelling was done by MODELLER (version 9.19, Departments of Biopharmaceutical Sciences and Pharmaceutical Chemistry, and California Institute for Quantitative Biomedical Research, San Diego, USA) using the 'automodel' class and included the NADPH cofactor [60,61]. For both TbPTR1 and TcPTR1, of the 100 models generated, the top models were validated using the ProSA [62] online server ( Figure S15). A table gathering a summary of the TbPTR1, TcPTR1, LmPTR1, HsDHRS4, and TbDHFR protein structures is presented in Table S4. The TbPTR1 structure was co-crystallized with cyromazine, which was used to validate the docking procedure.
We carried out blind docking of the ligand library against TbPTR1, TcPTR1, LmPTR1, and HsDHRS4 homotetrameric protein structures that included their NADPH cofactors using Autodock Vina (version 7.4, Scripps Research Institute, San Francisco, USA) [63]. Later, the compounds were blind docked to the TbDHFR (PDB: 3QFX) [42] dimeric structure that included its NADPH cofactors using Autodock Vina. The docking parameters used for each of the proteins are summarized in Table S5. Protein-ligand complexes were then evaluated based on if the ligand was located in the active site, as well as based on binding mode, selectivity, docking energy scores, and hydrogen bonding. The docking energies were further evaluated by re-docking the compounds to their protein targets and then using Xscore to give an independent energy score [32].

Prediction of Blood-Brain Barrier Permeability
To prioritize compounds that can cross the BBB, a PCA was carried out to identify which compounds occupied the same chemical space as known CNS-permeable drugs [64][65][66]. The PCA was based on the molecular descriptors of the top binding compounds, FDA-approved drugs, and FDA-approved CNS-permeable drugs [51]. The molecular descriptors used included XlogP, number of H-bond donors (HBD), number of H-bond acceptors (HBA), net charge (NC), topological polar surface area (tPSA), molecular weight (MWT), number of rotatable bonds (NRB), and polar and apolar desolvation. The first and second principal components were used to create a scatter plot that explained the largest percentage of the variance.

Molecular Dynamics
Eighteen protein-cofactor-ligand complexes were parametrized using the AMBER03 force field utilizing ACPYPE [67] and GROMACS (5.1.4) [68]. Each protein-ligand complex was solvated using a Simple Point Charge (SPC) water model in a cubic box of 5.07 × 5.18 × 5.16 nm with a 2.37 × 2.86 × 3.30 nm center. A minimum distance of 1.5 nm was allowed between any protein or ligand atom with the wall. The systems were then neutralized using Na + and Cl − counter ions. The MD systems also included simulating the protein in complex with the NADPH cofactor without the ligand. The MD simulations were performed using GROMACS 5.1.4 [68]. To correct for any structural distortions, the systems were minimized using a steepest descent algorithm using a 100 kJ/mol/nm tolerance value. This was followed by an equilibration using constant number of particles, pressure, and temperature (NPT) and constant number of particles, volume, and temperature (NVT) ensembles. This was finally followed by a 200 ns production run at 300 K without any restraints. Trajectories were generated every 2 fs with protein bonds involving hydrogens being constrained using a Linear Constraint Solver (LINCS) algorithm and saved after every 10 ps. The MD trajectory analysis included RMSD, Rg, RMSF, and PCA, using the GROMACS toolbox, Visual Molecular Dynamics (VMD) [69], and ProDy [70]. All MD simulations were performed at the Center for High Performance Computing (CHPC) (Cape Town, South Africa) using 240 cores (CPU: Intel ® Xeon ® ).
Graphs and diagrams were generated using, JChemPaint [63], PyMOL [29] and GRACE software (http://plasma-gate.weizmann.ac.il/Grace/). Protein-ligand complex structures were generated from the equilibrated trajectories at the end of the simulation. These structures were then used to analyze the protein-ligand interactions as well as residue interactions using Discovery Studio [31].

MM-PBSA Free Energy Calculations
The last 50 ns of the equilibrated MD trajectories were used to perform binding free energy (BFE) calculations of the ligand-protein complexes using the g_mmpbsa package (version 1.6, Jawaharlal Nehru University, New Delhi, India) [71]. The BFE calculation was based on the MM-PBSA method [72,73]. The BFE of the protein-ligand complexes was calculated using the equations below (in general terms): (1) The binding free energy of the protein-ligand complex in solvent (∆G binding ) where G complex is described as the total energy of the protein-ligand complex, G protein as the isolated free energy of the protein and G ligand as the isolated free energy of the ligand.
(2) The free energy of either the ligand, protein, or protein-ligand complex (G x ), where the average mechanical potential in a vacuum is described as <E MM >, TS as the entropic contribution (T is temperature and S is entropy), and G solvation as the free energy of solvation.
(3) The vacuum molecular mechanics potential energy (E MM ), where E bonded represents bonded interactions such as bonds, dihedrals, angles and improper interactions. The non-bonded interactions (E nonbonded ) are modelled using the Coulomb and Lennard-Jones (LJ) potential functions. They include electrostatic interactions (E electrostatic ) and van der Waals interactions (E vdW ).
(4) The energy required to transfer the protein-ligand solute from a vacuum into a solvent is described as the free energy of solvation (G solvation ), where G polar and G nonpolar describe the electrostatic and non-electrostatic energy contributions, respectively [71].
Furthermore, to determine the energy contribution of each protein residue that binds with the ligand, a free energy decomposition was carried out using g_mmpbsa [71]. This allowed for a better understanding of the protein-ligand interactions and helped identify PTR1 binding residues of functional significance.

Average Shortest Path, and Average Betweenness Centrality
We carried out dynamic network analysis on the equilibrated (after 50 ns) apo protein and protein-ligand MD trajectories using MD-TASK [43] in order to identify changes in the topological properties of the proteins brought about by the ligand interactions. This was used to glean the impacts of ligand binding on protein dynamics, function, and conformation. A cut off of 6.7 Å was used in the creation of the dynamic residue networks in MD-TASK. The average shortest path gave the density of shortest paths (L) between all node pairs [43]. The average betweenness centrality was used to identify residues in the dynamic network that were important for communication flow. Additionally, by comparing the apo protein and the protein-ligand complexes, we were able to use Average BC to assess how communication flow across the dynamic network was altered by ligand binding during the MD simulations. We generated equilibrated structures at the end of the simulations in order to map the interaction networks of any important residues identified using Discovery Studio [31].

Trypanosoma In Vitro Inhibition Assay
Compounds RUBi001, RUBi005, and RUBi015 were purchased from MCULE while RUBi002, RUBi004, RUBi007, RUBi008, RUBi010, RUBi011, RUBi012, RUBi014, RUBi016, and RUBi018 were obtained from MolPort (not all the compounds were commercially available). These compounds were assayed for trypanocidal activity by adding 20 µM of each compound to cultures of T. b. brucei (strain Lister 427) in 96-well plates. The parasites were maintained at 37 • C and 5% CO 2 in an IMDM medium containing 25 mM HEPES, 10% fetal bovine serum, 1 mM hypoxanthine, 0.05 mM bathocuproine disulfonic acid, 1.5 mM cysteine, 1.25 mM pyruvic acid, 0.09 mM uracil, 0.09 mM cytosine, 0.16 mM thymidine, and 0.014 % 2-mercaptoethanol. Parasites were diluted to 2.4 × 10 4 cells in a volume of 200 µL per well and incubated with the test compounds for 24 h. Parasite percentage viability was determined using the resazurin method [74]. Resazurin (0.5 mM in phosphate-buffered saline; 20 µL/well) was added to each well and incubation continued for a further 24 h, after which fluorescence (Ex 560 /Em 590 ) was read in a Spectramax M3 microplate reader (Molecular Devices). Trypanocidal activity of the compounds was reported as the percentage of viable parasites in the compound-treated wells when compared to untreated controls (% viability). Pentamidine, an FDA-approved trypanocidal drug, was used as the control drug standard [75]. For compounds that produced <20% viability, IC 50 values were subsequently determined. The assays were conducted as described above, except that parasites were incubated with 3-fold serial dilutions of the test compounds and IC 50 values derived from % parasite viability vs. log[compound] dose-response plots by non-linear regression analysis using GraphPad Prism (version 5.02, GraphPad Software Inc., San Diego, USA). IC 50 evaluations were carried out on three independent occasions and the mean IC 50 ± standard deviation is reported in the text. To assess compound interactions, the compounds were assayed for trypanocidal activity when used in combination with WR99210, a known TbDHFR inhibitor [42]. For combination assays, IC 50 values were determined for RUBi004, RUBi007, RUBi014, RUBi016, and RUBi018 as well as WR99210 alone using a starting concentration of 100 µM and 20 µM for the RUBi compounds and WR99210, respectively, and in combinations at ratios of 75:25, 50:50 and 25:75, respectively (thus starting with concentrations of 75 µM/5 µM, 50 µM/10 µM, and 25 µM/15 µM for the RUBi compounds/WR99210). For isobologram analysis, the fractional inhibitory concentrations of the RUBi compounds and WR99210 were calculated by dividing the IC 50 s obtained for the compounds at the various combination ratios with the IC 50 obtained for the compounds in the absence of the partner drug, and the FIC values plotted against each other, i.e., RUBi compound FIC versus WR99210 FIC.

In Vitro Human Cytotoxicity Assay
The compounds assayed for trypanocidal activity were also tested to determine if they caused adverse effects against human cells in vitro. For this assay, HeLa cells were used. The cells were cultured in DMEM supplemented with 10% fetal calf serum and antibiotics (penicillin/streptomycin/amphotericin B) at 37 • C in a 5 % CO 2 incubator. Cells were plated at a density of 2 × 10 4 cells/well and after an overnight incubation the compounds were assayed for cytotoxic activity by adding three-fold serial dilutions of each compound (with a 100 µM starting concentration) to the 96-well plates, followed by incubation for 48 h. Cell viability was determined using the resazurin method [74]. Resazurin (0.5 mM in phosphate-buffered saline; 20 µL/well) was added to the cells and, after a 2-h incubation, fluorescence was read in a Spectramax M3 plate reader at excitation and emission wavelengths of 560 nm and 590 nm, respectively. Fluorescence readings were converted to percentage cell viability relative to control wells untreated with compounds. Plots of % cell viability versus log[compound] were used to determine IC 50 values by non-linear regression using GraphPad Prism (version 5.02, GraphPad Software Inc., San Diego, USA). Three repeats of the experiment were carried out. Emetine, a drug that induces cell apoptosis, was used as a control [76], and produced an IC 50 of 0.013 µM.

Pan-Assay Interference Compounds Assay
Compounds that showed trypanocidal activity were also subjected to a PAINS assay using the web server located at http://www.cbligand.org/PAINS/. This was done to identify and flag any compounds that contained PAINS features [47,77].

Conclusions
Dual inhibition of TbPTR1 and TbDHFR is a promising approach to successfully developing safe and effective anti-folate based anti-trypanosomal chemotherapeutics. As shown in this study, computation-based approaches are useful in fast and rapid rational drug design. Furthermore, in the discovery of novel TbPTR1 inhibitors, when the compounds are assayed in combination with known DHFR inhibitors, careful interpretation of isobologram assays is required to obtain an optimal outcome. When used in combination with WR99210, a known TbDHFR inhibitor, compounds RUBi004, RUBi007, RUBi014, and RUBi018 showed moderate to strong antagonism as demonstrated by isobologram results, which would indicate that they might be binding to both TbPTR1 and TbDHFR. RUBi016, as shown by its additive effect and molecular docking results, appears to selectively bind to TbPTR1. The five compounds assayed showed anti-trypanosomal activity with no significant human cell cytotoxicity in vitro. The merging of these scaffolds could yield to the development of even more potent and selective TbPTR1 inhibitors.