ACTs: Artemisinin-Based Combined Therapies; ADMET: Absorption, Distribution, Metabolism, And Excretion– Toxicity; DHF: Dihydrofolate; MOE: Molecular Operating Environment; MTX: Methotrexate; P. Falciparum: Plasmodium Falciparum; Pfdhfr- TS: Plasmodium Falciparum Dihydrofolate Reductase-Thymidylate Synthase; PTX: Pemetrexed; RTX: Raltitrexed; THF: Tetrahydrofolate; TMX: Trimetrexate; TS: Thymidylate Synthase

Malaria is still prevalent in tropical and sub-tropical countries with a continued rate of morbidity and mortality [1]. Efforts to generate a vaccine against malaria have always been in grey. The antimalarial drugs provide suitable treatment but often resistance against these drugs develops. Different stages of malarial parasite viz., sexual, merozoite and sporozoite provide different stages of prevention against malaria. One of the stages involves folate proteins, such as Plasmodium falciparum dihydrofolate reductase (PfDHFR) which offers another pathway for intervention into the malaria [2]. The folate interference results into malaria arrest on folate synthesis pathway.

DHFR is an important factor in folate biosynthesis both in human and P. falciparum. In P. falciparum DHFR is associated with thymidylate synthase (TS) and works as a bifunctional enzyme while in human DHFR and TS exist as two independent proteins [3,4]. This distinguishing feature enhances the selectivity of this chemotherapeutic target [5]. In P. falciparum, a single polypeptide contains both DHFR and TS domains at N- and C-terminal, respectively, encoded by dhfr-ts genes [6]. The length of polypeptide for DHFR and TS varies among different protozoans. In P. falciparum, the DHFR domain comprises 231 amino acids (1-231 residues) and TS domain is of 288 amino acids (321-608 residues), which are separated by a junctional peptide (linker region) of 89 amino acids (232-320 residues) [4,7]. The native protein exists in homodimeric form with 140 kDa molecular weight (each monomer is of 70 kDa) [7].

Inhibition of PfDHFR effects methionine synthesis and reduced the level of thymidylate synthase that blocks DNA replication and kills the parasite [2,3]. Several antimalarials such as pyrimethamine and cycloguanil [8] have been utilized as potential inhibitors of PfDHFR. Antifolates mainly target the DHFR domain of PfDHFR-TS and remained effective until the development of resistance in Plasmodium strains due to point mutations in the target protein [8-11]. In resistant clinical isolates, different combinations of four point mutations N51I, C59R, I164L and S108N have been reported in the DHFR domain of PfDHFR-TS [2,12]. The S108N was the first mutation reported to cause a decrease in binding affinity due to steric clash between the heavy side chain of the mutated residue and the p-chlorophenyl substituent of pyrimethamine. Additional reduction in binding affinities of antifolates to the DHFR region is due to three other mutations (N51I, C59R and I164L). The effects of the single, double, triple and quadruple mutations on the binding affinity of antifolates have been extensively studied and described the mutation-induced resistance to antifolates [4,8,12-14].

Recently artemisinin-based combined therapies (ACTs) have been recommended as antimalarials. Although ACTs resulted in effective malaria control but lately, ACT resistant isolates have been reported in various endemic regions [15-17]. An effective antimalarial drug against resistant strains of Plasmodium is required. The development of potent drugs against resistant target proteins is cumbersome. Computational methods are helpful to predict the potency of drug candidates and to establish the preliminary information about interaction between active site amino acids and drug candidates. Current therapeutic strategies have been remarkably facilitated by in silico studies and resulted in effective drugs which were then confirmed to be excellent therapeutic agents [18-20].

Structurally rigid inhibitors like pyrimethamine failed to bind to quadruple mutant PfDHFR-TS, due to structural alteration in active site. However, structurally flexible inhibitors have been reported to bind in resistant and wild type strains of P. falciparum [21-23]. The methotrexate (MTX), TMX (trimetrexate), PTX (pemetrexed) and raltitrexed (RTX) are well known antifolates [23,24] with relatively flexible chemical structures. The 2-amino-4-oxo-pyrimidine, 2,4-diaminopteridine, 2,4-substituted quinazoline and glutamic acid are the important pharmacophores in chemical structures of MTX, TMX, PTX and RTX. The MTX, TMX, PTX and RTX contain a 2,4-substituted moiety and three of them share another chemical structure in the form of glutamic acid linked to the central part through carbonyl group. The 2,4-diamine is an important pharmacophore among antifolates and an essential part of the antifolate activity to DHFR [3,25,26].

In this study, potential of structural analogues of MTX, TMX, PTX and RTX to inhibit DHFR and TS domains of quadruple mutant PfDHFR-TS have been studied and evaluated. The analogues of MTX and TMX were designed to fit in the quadruple mutant DHFR domain of PfDHFR-TS, while analogues of RTX and PTX were constructed to target the TS domain. Analogues of four selected antifolates were designed, characterized and docked with 3D structure of PfDHFR-TS to search for appropriate drug candidates. Our work suggests that charge mediated interactions between receptor residues and the glutamic acid scaffold effectively contribute to determine the antifolate potential along with the flexibility of the linker region holding 2,4-diamino heterocyclic part of the ligand.

Selection and pre-processing of receptor protein

The X-ray crystallographic structure of the quadruple mutant Pf DHFR-TS was retrieved from RCSB Protein Data Bank (PDB). The PDB ID of selected crystal structure of PfDHFR-TS is 4dp3 and it is the most recently reported quadruple mutant (N51I, C59R, S108N, and I164L) structure in the databank [4]. It exists in homodimeric (chain A and B) form of PfDHFR-TS, but for docking analyses chain B was removed from the original coordinate file with the help of edit tools of Discovery Studio 3.5 Visualizer [27]. Additional pre-docking steps such as, removal of water molecules, co-factor and co-crystallized ligands were also performed by utilizing Discovery Studio 3.5 Visualizer.

Generation of inhibitors dataset

In next step, a dataset of 33 compounds containing 2,4-diamino heterocycles, substituted pyrimidine or glutamic acid scaffold (analogues of MTX, TMX, RTX and PTX) were designed to investigate the antifolate potential against mutated PfDHFR-TS. The 2D structures for all designed ligands (by addition or substitution in parent compound structures) are shown in Figures 1-4. The ChemDraw [28] was utilized for drawing the chemical structures of analogues. Validation of the compounds was performed by screening the drug databases PubChem [29], eMolecules [30], and ChemSpider [31] (Supplementary Table 1). All compounds were subjected to Lipinski’s Rule of Five and ADMET for their drug-like characteristics (Supplementary Table 2).

Table 1: Docking results of selected ligands.

Table 2: Comparative Analysis of top scoring ligands.

Figure 1: Methotrexate (MTX) analogues. X and Rⁿ (n=1, 2, 3) are the points of modification. The table describes the modifications at each position.

Figure 2: Trimetrexate (TMX) analogues. X and Rⁿ (n=1, 2, 3, 4, 5, 6) are the points of modification. The table describes the modifications at each position.

Figure 3: Raltitrexed (RTX) analogues. X, Y and R are the points of modification. The table describes the modifications at each position.

Figure 4: Pemetrexed (PTX) analogues. X and Rⁿ (n=1, 2) are the points of modification. The table describes the modifications at each position.

Molecular docking simulation

The Molecular Operating Environment (MOE) program [32] was used for ligand-receptor virtual docking and binding energy calculations. The 3D structures of ligands were protonated with default parameters by Protonate 3D tool of MOE. The energy minimization for each protonated ligand was performed by utilizing MOE. Finally, the four new databases for each group of compounds that served as ligand dataset were also generated to run docking simulation by MOE.

The primary preparation (as mentioned earlier) of receptor protein was performed in Discovery Studio Visualizer. The addition of hydrogen atoms and energy minimization of receptor protein was also performed by using MOE. This step was performed for the relaxation of newly added hydrogen atoms to adjust all non-hydrogen atoms. After performing all these steps, the receptor protein was prepared for docking simulations. The default parameters of MOE-Dock program were used for the molecular docking of the ligands. The Site Finder tool was utilized to specify the binding pocket of DHFR and TS domains.

The ligands were allowed to be flexible in search for the correct conformations of the ligands and to obtain minimum energy structures. After the generation of the docking results, the best conformations against both DHFR and TS domains were analyzed for their binding interactions by using the ligand interaction tool of MOE. Ligplot [33] was also generated for selected top scoring complexes to validate the interaction pattern purposed by MOE.

Post docking comparative analysis

The selected drug candidates for DHFR and TS domains were also docked with relevant human proteins. The 3D structure of Human DHFR (hDHFR) (PDB ID: 1drf), and Human TS (hTS) (PDB ID: 1hvy) proteins were retrieved from RCSB PDB. Both proteins were prepared for docking as it has been described for PfDHFR-TS and docked with lead compounds, MTX04 and PTX03 against DHFR and TS domains respectively. Ligand interaction diagrams were generated to analyze the binding pattern of proposed drug candidates with human analogue proteins.

PfDHFR-TS have various wild types and (single, double, triple and quadruple) mutant clinical isolates but the most resistant strain of P. falciparum is quadruple mutant. In the current study, the quadruple mutant (N51I, C59R, I164L and S108N) crystal structure (4dp3) with resolution of 2.40 Å, was selected as a receptor protein. The selected crystal structure covers 89% of the amino acid sequence (UniProtKB ID: D9N170) of both DHFR and TS domains.

The PfDHFR-TS is a well-known target of antifolates. Several point mutations in active site of DHFR domain developed the resistance against these antifolates. Previously effective antifolates like pyrimethamine, cycloguanil, WR99210, chlorcycloguanil and trimethoprim share some common functional groups like 2,4-diamino scaffold [2]. In the current study, the choice of 2,4-substituded heterocycles and glutamic acid is based on the observation that the conformation of 2,4-substituted heterocycles and the polar nature of glutamic acid plays an essential role to antagonize the activity of folic acid or N5,N10–methylene-5,6,7,8-tetrahydrofolate (mTHF) [4,34]. The similarity of these chemical structures with respective original enzyme substrates is the main reason of their binding to the active site. MTX and DHF only differ by a single substituent methyl group while RTX and PTX closely resemble mTHF.

The quadruple mutant PfDHFR domain inhibitors were designed by taking MTX and TMX as backbone structure. The analogues of MTX and TMX were considered to attain the required 3D conformation against mutant PfDHFR. The analogues of RTX and PTX were designed to target the TS domain of PfDHFR-TS. The drug-like properties of chemical compounds can be predicted by calculation of ADMET and Lipinski’s rule of five. Lipinski’s rule of five was established to estimate the absorption or permeation of drug like compounds [35,36]. According to this rule, if a compound violates two or more than two rules, it may consider as a drug candidate with poor bioavailability. In the present dataset of ligands, only five compounds (MTX02, MTX05, MTX06, MTX09 and MTX12) were found to violate three, while, thirteen compounds were found to violate two Lipinski’s rule of five. In this study, the Lipinski’s rule of five is used to estimate the oral bioavailability of the drug candidates, rather than screening of the compounds dataset.

The docking analyses of PfDHFR-TS with 33 designed ligands were performed by MOE. Out of 33 compounds, 22 (MTX01-12 and TMX01-10) were docked against DHFR while other 11 (RTX01- 05 and PTX01-06) were docked against TS domain based on their chemical structures. Docking simulation produced 10 possible binding conformations for each ligand and ranked them according to their predicted values of binding energies. The top scoring conformations for each ligand with minimum value of free binding energy (S-value) has shown in (Supplementary Figure 1). Screening of 330 complexes based on free binding energy (S-value) selected the best complex for each compound. The 33 top scoring complexes further screened by analysis of their interaction pattern with critical amino acids and the S-value of the complexes. It resulted in the selection of two compounds from each group based on maximum number of potential hydrogen bonds with critical amino acids and as well as the minimum the S-value. For example, among the selected complexes of MTX analogues, the single compound (MTX04) has shown hydrogen bond potential with four critical residues and it also has comparatively low S-value (-7.6275Kcal/mol) (Supplementary Figure 1). The MTX09 has shown hydrogen bond with three critical amino acids but it has slightly high S-value (-6.8873Kcal/mol) as compared to the other analogues of the MTX. Furthermore, each of the four compounds (MTX02, MTX03, MTX05 and MTX07) have shown interaction with two critical residues but MTX03 has the minimum value of free binding energy (-8.2940Kcal/mol) (Supplementary Figure 1). Therefore, the finally selected compounds from the MTX group were MTX04 and MTX03. The similar process was performed on the other three groups. The eight best selected complexes, two from each group include MTX03, MTX04, TMX01, TMX09, RTX01, RTX05, PTX03 and PTX04 (Table 1).

The 2D ligand interaction diagrams of selected (eight) ligands were generated by using the ligand interaction tool of MOE (Supplementary Figure 1). Further selection (among eight) was done by comparative analysis of interaction pattern with known drugs like pyrimethamine and WR992012. Among MTX03, MTX04, TMX01 and TMX09, the MTX03 shown the minimum value of binding energy (-8.2940Kcal/ mol) but it forms potential hydrogen bond with only two critical amino acids (Ile14 and Arg122). The MTX04 has shown more number of hydrogen bonds with critical residues and also has comparable value of binding energy (-7.6275Kcal/mol). The RTX01, RTX05, PTX03 and PTX04 were screened for TS domain. The PTX04 shown the minimum binding energy (-7.2820Kcal/mol) but PTX03 shown additional cation-π interaction with Arg345 with considerably similar value of binding energy (-6.9362Kcal/mol). Therefore, the MTX04 and PTX03 have selected as appropriate lead compounds among designed data set against DHFR and TS domain respectively. The Ligplot analysis validated the interaction pattern for MTX04 and PTX03. The Ligplot graphical representation reproduced all potential hydrogen bonds for PTX03 and three out of five hydrogen bonds for MTX04 (Figure 5). The ligand binding interaction diagram indicate that MTX04 have potential to form hydrogen bonds with Asp54, Arg59, Ser120, Arg122 and Leu164 (Figure 6) while other residues exposed towards ligand with the probability to interact are Cys15, Ala16, Ile14, Trp48, Ile51, Met55, Phe58, Phe116, Leu119, Ile112 (Figure 6). MTX04 has shown the potential to interact with critical residues (Asp54, Arg122) of PfDHFR active site as reported in complexes with other known drugs like pyrimethamine, WR99210 and with substrate. MTX04 has also shown potential for hydrogen bonding with two mutated residues (Arg59, Leu164) to provide appropriate steric complementarity to fit in active site of PfDHFR (Figure 6).

Figure 5: Ligplot analysis of MTX04-PfDHFR and PTX03-PfTS complex. (A) The 2D representation of interaction pattern for MTX04 in quadruple mutant PfDHFR domain. (B) The 2D representation of interaction pattern for PTX03 in PfTS domain.

Figure 6: Binding pose of the MTX04 in DHFR domain. The 2D (left side with white background) and 3D (right side with black background) views of MTX04 docked in DHFR domain. In 2D, the dotted arrows are showing potential hydrogen bonds with Asp54, Arg59, Ser120, Arg122 and Leu164. In 3D, the MTX04 is in green color while the active site residues are in element style.

The MTX04 ligand has a polar head group in the form of glutamic acid, substituted benzyl linker and 2,4-diaminopteridinyl group (Figure 1). The substituted, amino group interacts with carbonyl oxygen of Asp54 through charge mediated hydrogen bond while carbonyl oxygen of Leu164 act as donor to form hydrogen bond with other substituted amino group hydrogen of pteridine ring (Figure 6). The side chain amino groups of Arg59 and Arg122 act as donor to form hydrogen bond with glutamic acid carbonyl oxygen and hydroxyl group oxygen respectively (Figure 6). The other hydroxyl group oxygen of glutamic acid part is suspected to form hydrogen bond with electropositive nitrogen of Ser120 (Figure 6).

As opposed to DHFR, the TS domain has not been considered as primary target of anti-malarial antifolates. In the current study, the active site of the TS domain was also targeted by analogues of such inhibitors that have already been tested against the TS domain of other protozoa. Out of 220 complexes generated from 22 ligands, four top ranked ligands were further screened and only PTX03 was selected as a drug candidate to inhibit the TS domain. PTX03 has potential to form hydrogen bonds with Arg510, His551, Ser524 and Arg345. The Leu487, Phe520, Asn521, Ser511, Tyr553, Glu382, Asn407, His491, Ile403 and Trp404 are exposed towards the ligand and involved in ligand recognition (Figure 7). PTX03 interact with key residues in the active site of the TS domain as reported in crystal structure investigations and the flexibility of the PTX03 structure made it more suitable to fold and mimic the orientation of the substrate. Arg345 has potential to form hydrogen bonds at two different oxygen atoms (one from carbonyl group and another from hydroxyl group) of glutamic acid moiety of PTX03 (Figure 7). Arg345 also has potential to form a cation-π interaction with aromatic ring of PTX03 (Figure 7). Hydroxyl group oxygen of Ser524 acts as hydrogen bond donor to the amino substituent of dihydroquinazoline (Figure 7). The imidazole side chain of His551 has potential to make a hydrogen bond with carbonyl oxygen of PTX03 while the side chain amino group of Arg510 may form a hydrogen bond with the carbonyl oxygen of PTX03 (Figure 7). All interacting residues lie within the reported catalytic site of TS domain of P. falciparum [9] and play vital role to bind with cofactor

Figure 7: Binding pose of the PTX03 in TS domain. The 2D (left side with white background) and 3D (right side with black background) view of PTX03 docked in TS domain. In 2D, dotted arrows are showing potential hydrogen bonds with Arg345, His510, Ser524 and His551. The Arg345 is also showing cation-π interaction with PTX03. In 3D view the PTX03 is in green color while the TS domain is in element style.

Comparative analysis of binding complexes revealed that MTX04 and PTX03 interact with Human DHFR and TS in completely different ways as compared to the respective Plasmodium analogue protein (Figures 8 and 9). The minimum binding energy value for the hDHFR and MTX04 complex is -5.9421kcal/mol while minimum binding energy value for the hTS and PTX03 complex is -7.5608kcal/mol (Table 2). The binding energy values for the complex of PTX03 with TS domain of PfDHFR-TS and hDHFR do not show significant difference (Table 2). To explore the underlying reason, multiple sequence alignment has been done for TS domain of PfDHFR-TS and hTS (Supplementary Figure 1b). It revealed that conservation of sequence that could be the reason of similar binding potential shown by the TS inhibitor. The considerable difference between binding energies (Table 2), primary sequence (Supplementary Figure 1a) and 3D conformation of MTX04 (Figure 6) strongly suggests that this drug candidate has more potential to bind with DHFR domain of PfDHFR-TS then hDHFR.

Figure 8: Docking poses of the MTX04 in PfDHFR and hDHFR. The 3D view of the conformation of MTX04 in (a) DHFR domain of PfDHFR-TS and (b) hDHFR protein. Surface potential around ligand was calculated by Discovery Studio Tools. The residues with hydrogen bond donating potential and hydrogen bond accepting potential are in purple and green color respectively.

Figure 9: Docking poses of the PTX03 in TS domain of PfDHFR-TS and hTS. The 3D view of the conformation of PTX03 in (a) TS domain of PfDHFR-TS and (b) hTS protein. The residues with hydrogen bond donating potential and hydrogen bond accepting potential are in purple and green color respectively.

The docking simulation of designed dataset reveals that MTX04 and PTX03 have greater tendency to bind with the quadruple mutant DHFR domain and TS domain of PfDHFR-TS respectively. The current study could be helpful to develop antifolates against mutated targets by considering MTX04 and PTX03 as lead compounds. A successful synthesis of these compounds will lead to design experimental protocols for complex formation of PfDHFR-TS with these lead compounds followed by functional analyses of the complexes.

This research was funded by Pakistan Academy of Sciences and EMROWHO.