Origins of the specificity of inhibitor P218 toward wild-type and mutant PfDHFR: a molecular dynamics analysis

Molecular dynamics simulations were performed to evaluate the origin of the antimalarial effect of the lead compound P218. The simulations of the ligand in the cavities of wild-type, mutant Plasmodium falciparum Dihydrofolate Reductase (PfDHFR) and the human DHFR revealed the differences in the atomic-level interactions and also provided explanation for the specificity of this ligand toward PfDHFR. The binding free energy estimation using Molecular Mechanics Poisson-Boltzmann Surface Area method revealed that P218 has higher binding affinity (~ −30 to −35 kcal/mol) toward PfDHFR (both in wild-type and mutant forms) than human DHFR (~ −22 kcal/mol), corroborating the experimental observations. Intermolecular hydrogen bonding analysis of the trajectories showed that P218 formed two stable hydrogen bonds with human DHFR (Ile7 and Glu30), wild-type and double-mutant PfDHFR’s (Asp54 and Arg122), while it formed three stable hydrogen bonds with quadruple-mutant PfDHFR (Asp54, Arg59, and Arg122). Additionally, P218 binding in PfDHFR is stabilized by hydrogen bonds with residues Ile14 and Ile164. It was found that mutant residues do not reduce the binding affinity of P218 to PfDHFR, in contrast, Cys59Arg mutation strongly favors inhibitor binding to quadruple-mutant PfDHFR. The atomistic-level details explored in this work will be highly useful for the design of non-resistant novel PfDHFR inhibitors as antimalarial agents.


Introduction
Malaria is an important cause of death and illness in children and adults in the tropical and sub-tropical countries, causing 300-500 million clinical cases and 1.2 million deaths annually (Murray et al., 2012). Many malaria-control strategies exist, but none are appropriate and affordable in all contexts. Resistance to antimalarials has been observed for all currently used antimalarials (amodiaquine, chloroquine, mefloquine, quinine, and sulfadoxine-pyrimethamine, cycloguanil) and more recently, in artemisinin derivatives (Dondorp et al., 2009;Gregson & Plowe, 2005). The development and spread of drugresistant strains of malaria parasites have been identified as a key factor in this resurgence and is one of the greatest challenges for the design of an effective antimalarial compound.
Plasmodium falciparum Dihydrofolate Reductase (PfDHFR) is a well-known target for combating malaria. Its function is to catalyze the biochemical reaction which involves the reduction of 7,8-dihydrofolate (DHF) to 5,6,7,8-tetrahydrofolate (THF), using the reduced form of nicotinamide adenine dinucleotide phosphate (NADPH) as a cofactor. Hence, it plays a key role in the folate biosynthesis pathway and is responsible for the generation of the DNA base, deoxythymidine monophosphate. A difference between host (human) and parasite (P. falciparum) DHFR active sites was reported and this difference offers opportunities for the design of selectively binding drugs (Nirmalan, Sims, & Hyde, 2004;Zhang & Rathod, 2002). Two enzymes of the parasite (DHFR and TS (Thymidylate Synthase)) reside on the same polypeptide chain as a DHFR-TS bifunctional protein (P. falciparum DHFR-TS). As demonstrated by X-ray crystallographic structure of P. falciparum DHFR-TS, it contains 608 amino acids, the first 231 comprising the DHFR domain, the next 89 residues forming the junction region, which joins the remaining 288 residues of the TS domain (Ivanetich & Santi, 1990).
Antifolate drugs known to inhibit the normal function of PfDHFR are Pyrimethamine (Pyr) and Cycloguanil (Cyc) (Figure 1). The triazine proguanil was discovered in the mid 1940s as an effective antimalarial compound (Curd, Davey, & Rose, 1945). Proguanil is converted in vivo to its active metabolite, cycloguanil, which selectively binds to and inhibits PfDHFR (Carrington, Crowther, Davey, Levi, & Rose, 1951). Recently, quantum chemical method has been employed to study the possible mechanisms for cytochrome-mediated Proguanil to cycloguanil conversion (Arfeen, Patel, Abbat, Taxak, & Bharatam, 2014). Pyrimethamine is also a potent and selective inhibitor of PfDHFR, and has been used as an antimalarial agent since 1952. Genetic mutations in the active site of PfDHFR at positions, Ala16Val, Ser108-Asn, Cys59Arg, Asn51Ile, and Ile164Leu have considerably decreased the therapeutic activities of these antifolates (Bras & Durand, 2003;Gregson & Plowe, 2005;Peterson, Milhous, & Wellems, 1990;Peterson, Walliker, & Wellems, 1988;Rastelli et al., 2000;Vanichtanankul et al., 2012;Yuvaniyama et al., 2003). Later, a dihydrotriazine derivative WR99210 (Figure 1) was developed which showed activity in wild as well as double and quadruple-mutant PfDHFR. The flexible side chain of WR99210 enables it to adopt a conformation that can fit into the active site modified by the Ser108Asn mutation (Nzila, 2006;Rastelli et al., 2000). Unfortunately, low bioavailability and inadequate pharmacokinetic profile of WR99210, limits its use in malaria chemotherapy (Rieckmann, Yeo, & Edstein, 1996). All these known antifolates have common structural features which makes essential molecular recognition interactions with the PfDHFR active site residues, to inhibit its activity. The important molecular recognition centers of PfDHFR for these known antifolates include residues Ile14, Asp54, and Ile/Leu164. The carboxylate oxygen atoms of Asp54 form two hydrogen bonds with WR99210 ( Figure S1, Supplementary material), Cyc, and Pyr. The backbone oxygen atoms of Ile14 and Ile164 also form H-bond with the amino group of antifolates in the active site of PfDHFR.
Quantitative structure-activity relationship (QSAR) studies have also been extensively carried out on these known antifolates and their derivatives to identify structural features responsible for the differences in anti-plasmodial activities with respect to their electrostatic, steric, and hydrophobic nature (Adane & Bharatam, 2009;Basak & Mills, 2010;Basak, Mills, & Hawkins, 2011;Hecht, Cheung, & Fogel, 2008;Maitarad, Kamchonwongpaisan, et al., 2009;Maitarad, Saparpakorn, et al., 2009;Ojha & Roy, 2011;Santos-Filho, Mishra, & Hopfinger, 2001;Sivaprakasam, Tosso, & Doerksen, 2009). The studies reported slightly different binding alignment of the ligands to the mutant form of an enzyme as compared to the wild-type form. Asn108 of the mutant enzyme was found to cause steric clash with Pyr and Cyc and is considered as the cause of resistance to these drugs.
Recently, structure-based drug design approach was adopted to discover an antimalarial compound P218, currently in its preclinical trial ( Figure 1); to be a selective and efficient inhibitor of both wild-type and quadruplemutant form of PfDHFR . P218 is designed to have basic functional groups required for an inhibitor to bind in PfDHFR cavity as well as it includes pyrimidine side chain flexibility and a carboxylate group which can make hydrogen bonds with conserved Arg122 (similar to natural substrate DHF). Sequence and X-ray crystal structure analysis showed that Met55, Cys/Arg59, and Phe116 in PfDHFR are replaced at structurally equivalent positions in human DHFR by Phe31, Gln35, and Asn64 which causes significant difference in the active sites of human DHFR (hDHFR) and PfDHFR . These differences in the active sites of PfDHFR and hDHFR lead to the loss of P218 interaction with Arg70 (a conserved residue in hDHFR), as observed from the X-ray crystal structure of hDHFR. The selectivity of P218 for PfDHFR is clear from the significant difference in experimentally determined binding constant values of P218 for PfDHFR (wild-type .51 ± .03 nM and quadruple-mutant .53 ± .13 nM) and hDHFR (2841 ± 319 nM) . P218 is reported to have in vitro and in vivo antimalarial activity against wild-type and pyrimethamine-resistant PfDHFR, as well as it is reported to possess suitable pharmacological, metabolic, and safety profiles for its consideration for clinical development .
Structural details of PfDHFR became available for the past decade and provided clues for the interaction between the catalytic site of an enzyme and also the clues regarding the inhibition as well as the reasons for resistance by the mutated enzyme. As mentioned above, many studies using computer-aided methods have been carried out to understand the enzyme-inhibitor interaction in the past one decade. All these efforts (X-ray, molecular docking, and QSAR) are under static structural conditions. Only two reports which involved the dynamics of enzyme-inhibitor interaction appeared. In 2000, Rastelli et al. built the homology models of wild-type and mutant PfDHFR and carried out molecular docking of Pyr, Cyc, and WR99210 and performed molecular dynamics (MD) to study the enzyme-inhibitor interactions (Rastelli et al., 2000). More recently, MD studies were carried out to evaluate the flexible environment of PfDHFR (Mokmak et al., 2014). The authors observed that the flexibility given to the enzyme during dynamical conditions helped in elucidating the binding patterns of drug (Pyr) and the lead (WR99210) in the active sites of wild-type and quadruple-mutant PfDHFR. Similarly, the dynamical environment of enzyme needs to be taken into account to understand the behavior of various inhibitors and the origin of specificity of the inhibitors and the origin of resistance to these inhibitors. Particularly, the recently developed lead compound P218 and the origin of specificity of this toward wild-type and mutant PfDHFR as against hDHFR is required to be explored. Such study becomes important in the wake of the fact that no new antifolate was introduced in the recent past despite extensive studies. This information gain might help researchers to design potent and specific inhibitors of PfDHFR. Hence, present computational study explores the binding mode and interaction profile of P218 to PfDHFR (wild and mutant) and human DHFR using MD simulations approach, which could be helpful for the design and optimization of novel antimalarials.

Protein structure preparation
The X-ray crystallographic structures of wild-type PfDHFR (wtPfDHFR, protein data bank (PDB) ID: 4DPD, co-crystallized ligand: DHF), quadruple-mutant PfDHFR (qmPfDHFR, PDB ID: 4DP3, co-crystallized ligand: P218), double-mutant PfDHFR (dm1PfDHFR, PDB ID: 1J3J, co-crystallized ligand: Pyr), and doublemutant PfDHFR (dm2PfDHFR, PDB ID: 3UM5, cocrystallized ligand: Pyr) enzymes for this molecular modeling exercise were obtained from RCSB protein databank (www.rcsb.org). Only DHFR domain (1-231 amino acids) of chain-A (having bound ligand in its active site) was considered, for both wild and mutant enzymes and water coordinates were removed. The missing residues in the chain-A for wtPfDHFR (residue numbers 22-29 and 84-96); qmPfDHFR (residue numbers 86-95); dm1PfDHFR (residue numbers 86-95); and dm2PfDHFR (residue numbers 86-95) were identified. For a successful dynamics run, it is necessary to add these missing residues. Modeller9v7 (Fiser & Šali, 2003;Šali & Blundell, 1993) software was used to add these missing residues. In order to add missing residues, sequence of PfDHFR enzyme in FASTA format was obtained from Universal Protein Resource. The FASTA file was used as input in the alignment script of model-ler9v7 and an alignment file was created. Then the script for adding missing residues was run and the 3D model with the added residues was generated from alignment information. The model was then validated by comparing the root mean square deviation (RMSD) of models with the corresponding crystal structures. The all atom RMSD values for wtPfDHFR (.09 Å), qmPfDHFR (.05 Å), dm1PfDHFR (.03 Å), and dm2PfDHFR (.01 Å) as compared to their corresponding crystal structures were found to be in acceptable range which signify that the model did not vary much as compared to crystal structure. Further, protein preparation of the hDHFR (X-ray crystal structure) and of PfDHFR (missing residues added structures) was done using Maestro (Maestro, 2011). Hydrogen atoms were added during protein preparation wizard. Receptor Grid Generation Panel within Glide suite (Glide, 2011) was used to set up receptor grid for the prepared structures. The grid was defined by 10 Å inner-box and 30 Å outer-box on each side of the centroid of bound ligand. The centroid of the grid box was same as the centroid of the bound ligand. The OPLS_2005 force field was employed for the grid generation.

Molecular docking
The 3D structures of P218 (bound ligand in X-ray crystal structures of hDHFR and qmPfDHFR), DHF (bound ligand in X-ray crystal structure of wtPfDHFR) and pyrimethamine (bound ligand in X-ray crystal structure of dm1PfDHFR and dm2PfDHFR) were built using sketch module of SYBYL7.1 and subjected to 1000 cycle minimization with standard Tripos force field and .005 kcal/mol energy gradient convergence criterion using Powell's method (Powell, 1964). The 3D structure of P218 built as above adopts a linear structure whereas its bioactive conformation (PDB ID: 4DP3) adopts a bent shape bringing the C 2 H 5 and C 6 H 4 groups into close proximity ( Figure S2, Supplementary material). Optimization of the bioactive conformation leads to a macrocyclic structure involving two intramolecular hydrogen bonds between carboxylic group and diaminopyrimidine ring, in a zwitterionic state. Further, conformational analysis using molecular mechanics-based systematic search followed by quantum chemical analysis showed that P218 adopts two macrocyclic arrangements: chair and boat. The macrocyclic boat structure is more stable than the macrocyclic chair structure by 3.15 kcal/mol (see Supplementary material Figure S3 for further details). Three-dimensional structures of these compounds were then prepared using LigPrep module of maestro (LigPrep, 2011) implementing OPLS_2005 force field and ionic states (carboxylates) for the ligands at pH values of 7.0 ± 2.0 were generated. Docking was performed using Glide software (Grid-based Ligand Docking with Energetics), (Glide, 2011) with the standard precision (SP) mode to estimate protein-ligand binding affinities and static intermolecular interactions (Friesner et al., 2004(Friesner et al., , 2006. Initially, to check the reliability of docking protocol, the bound ligands (as given in X-ray crystal structures) were redocked in the active sites of respective enzymes and RMSD values of the docked poses were compared to the X-ray crystallographic conformation of the ligand. After the validation of docking protocol, P218 was docked flexibly into the active site of the wtPfDHFR, dm1PfDHFR, and dm2PfDHFR enzymes. A maximum of 10 docking poses per ligand were generated in each case and analyzed further for the binding mode and intermolecular interactions.

MD simulations
MD simulations were performed to get detailed insight into the stability of hDHFR-P218 and PfDHFR-P218 complexes (both wild-type and mutants), intermolecular interactions of ligands in dynamical conditions, and to estimate the binding free energy of P218 in all the complexes. Full atomistic MD simulations were performed at the molecular mechanics level using the Amber 11 program package with the AMBERff99SB force field (Case et al., 2005(Case et al., , 2010. For the purpose of MD simulations, X-ray co-crystal structures of hDHFR (PDB ID: 4DDR) and qmPfDHFR (PDB ID: 4DP3) and the best docking poses of P218 in wtPfDHFR (PDB ID: 4DPD), dm1PfDHFR (PDB ID: 1J3J), and dm2PfDHFR (PDB ID: 3UM5) were selected. Enzyme-ligand complex structures were solvated in a truncated octahedral box of TIP3P water model that extended up to 12 Å in each direction of the solute. To calculate the long-range electrostatic interactions, the cut-off distance was kept at 12 Å and they were treated using the Particle-Mesh-Ewald method. The appropriate ionization states of all the amino acids were taken into consideration during the simulation. Before performing MD simulations, the systems were relaxed by a series of steepest descent (SD) and conjugate gradient (CG) minimizations. The minimized system was then gradually heated from 0 to 300 K for about 50 ps under constant volume-constant temperature (NVT) conditions, while keeping the solute atoms fixed by a weak harmonic restraint of 2 kcal/mol/ Å 2 . Then, 200 ps density equilibration with the same restraints on the solute was carried out. Subsequently, constant pressure equilibrium (NPT) for 2 ns at 300 K was carried out (relaxing the restraints), followed by the production run for 35 ns at 300 K under periodic boundary conditions. Isotropic position scaling was used to maintain the pressure (NTP = 1) and a constant pressure periodic boundary with an average pressure of 1 atm (PRES0) and a relaxation time of 2 ps were used (TAUp = 2.0). The five complexes in the study were large, especially when the water molecules were added to solvate the complex. For reasons of computational time and resource use, the simulation for each of the complexes was limited to 35 ns to explore the binding difference of P218 in the studied enzyme systems. Verlet algorithm (Verlet, 1967) was used to solve the equations of motion, with an integration time step of 2 fs and the trajectory was recorded at 2 ps interval. SHAKE algorithm (Ryckaert, Ciccotti, & Berendsen, 1977) was used to constrain all the covalent bonds containing hydrogen atoms. Langevin thermostat (NTT = 3) was used to maintain the temperature of the system at 300 K.

MM-PBSA calculation and normal mode analysis
The binding free energy between the enzyme and P218 was estimated using the Molecular Mechanics-Poisson Boltzmann Surface Area (MM-PBSA) approach incorporated in AMBER software. MM-PBSA calculations were performed on the last 5 ns trajectories of the total simulation run (35 ns). By using snapshot frequency of 2, total 1250 representative snapshots were extracted and the binding energies were calculated by running the python script in parallel (eight processors) (Miller et al., 2012). MM-PBSA analysis was performed to evaluate the binding affinity of the ligand toward hDHFR and PfDHFR (wild-type and mutants), which corresponds only to the enthalpy contributions of binding free energy. Hence, the binding energy value obtained from MM-PBSA analysis is not the real binding free energy since entropy contribution (unfavorable) to binding is not estimated. Normal mode analysis was carried out to calculate the binding entropies of the ligands. The average of the results from normal modes for the complex, receptor and ligand, was carried out to obtain an estimate of the binding entropy using SANDER module of AMBER. Since these calculations are computationally expensive, only 100 snapshots from last 5 ns simulation run were selected for the binding entropy (TΔS) calculation of the ligands. The binding free energy was estimated by taking into account enthalpy (from MM-PBSA calculations) and entropy (from normal mode analysis) contributions to binding free energy.

Multiple sequence alignment
Multiple sequence alignment of human DHFR (hDHFR), wild-type PfDHFR (wtPfDHFR) quadruple-mutant PfDHFR (qmPfDHFR; Asn51Ile, Cys59Arg, Ser108Asn, and Ile164Leu) double-mutant PfDHFR (dm1PfDHFR; Cys59Arg and Ser108Asn), and double-mutant PfDHFR (dm2PfDHFR; Ala16Val and Ser108Thr) was carried out using ClustalX (Jeanmougin, Thompson, Gouy, Higgins, & Gibson, 1998;Thompson, Gibson, Plewniak, Jeanmougin, & Higgins, 1997). The sequence alignment was carried out for the DHFR domain (231 residues) of the bifunctional DHFR-TS enzyme. This analysis was carried out to study the degrees of similarity and identity among hDHFR and PfDHFR and also to identify structurally homologous residues in hDHFR corresponding to PfDHFR. The alignment results show that hDHFR and wtPfDHFR share 30% identity and the regions (highlighted in box) of non-conservative amino acid differences in hDHFR as compared to PfDHFR which are shown in Figure 2. It has also been reported that the active site residues Met55, Cys/Arg59, and Phe116 in PfDHFR are substituted by non-conservative amino acids Phe31, Gln35, and Asn64 in the active site of hDHFR (highlighted with #, Figure 2). These amino acids changeover is reported to lead significant change in the active site among the host and parasite DHFR, such that P218 binds differently to hDHFR and PfDHFR .

Molecular docking analysis
To study the binding mode and interactions of P218 with human and PfDHFR (wild-type and mutants), the protein-ligand complex having P218 bound in their active site was required. X-ray crystal structure for P218-hDHFR and P218-qmPfDHFR complexes is reported in PDB. Initially, to validate the docking exercise, P218 was redocked in the active site of hDHFR and qmPfDHFR, to calculate its RMSD with respect to pose observed in X-ray crystallography ( Figure S4(A)-(B), Supplementary material). The calculated RMSD for hDHFR (.31 Å) and qmPfDHFR (.29 Å) was found to be within acceptable range. P218 binding mode in hDHFR and qmPfDHFR was further analyzed to identify the key molecular recognition interactions required for P218 to bind in the active site of these enzymes. It was found that diaminopyrimidine ring of P218 interacts with Ile7 and Glu30 in hDHFR (Figure 3(A) and (B)), whereas in qmPfDHFR, diaminopyrimidine ring forms hydrogen bonds with Ile14, Asp54, and Leu164 (Figure 3(C)). The terminal α-carboxylic group in qmPfDHFR forms hydrogen bonds with Arg122 and Arg59, similar interaction of α-carboxylic group with hDHFR residue Arg70 (structural homolog of Arg122 in PfDHFR) does not occur in X-ray crystal structure (Figure 3(A)) as well as in the docked pose of P218 in the active site of hDHFR (Figure3(B)). P218 docked in the active site of hDHFR was found to show hydrogen bonding interaction (3.2 Å) between α-carboxylic group and Asn64.
After validation of docking protocol, P218 was docked in the active site of wild-type PfDHFR and double-mutant PfDHFR structures (dm1PfDHFR and dm2PfDHFR) using GLIDE, molecular docking software. Docking pose analysis was carried out to examine the binding mode and key protein-ligand interactions  and the pose with high docking score and which showed key molecular recognition interactions was considered for further studies. The diaminopyrimidine ring of P218 docked in the active site of wtPfDHFR, dm1PfDHFR and dm2PfDHFR showed hydrogen bond interaction with Ile14, Asp54, Ile164, and also α-carboxylic group of P218 formed hydrogen bonds with Arg122 (Figure 3(D)-(F)).
To study the structural features that lead to the selectivity of P218 for PfDHFR and also to study the affect of PfDHFR mutations on P218 binding, all the above-mentioned five complexes were superimposed. From the superposition of X-ray crystal structures of P218 in the active sites of hDHFR and qmPfDHFR and docked poses of P218 in the active-sites of wtPfDHFR, dm1PfDHFR, and dm2PfDHFR (Figure 4), it can be realized that P218 binding mode in hDHFR and PfDHFR is significantly different. The terminal carboxylic group of P218 makes hydrogen bond with Arg122 in the case of PfDHFR but it does not form hydrogen bond with Arg70 (structural homolog of Arg122 in PfDHFR) in hDHFR. The docking results also suggested that P218 binding to various forms of PfDHFR is not affected by existing mutations, as P218 binds in a similar fashion to wild-type and mutant PfDHFR (double and quadruple-mutant PfDHFR). To further analyze whether the interactions observed in all the five complexes were stable over a period of time, MD analysis was carried out.

MD analysis
Five different enzyme-ligand structures (i) X-ray crystal structure of human DHFR (PDB ID: 4DDR), (ii) X-ray crystal structure of quadruple-mutant PfDHFR (PDB ID: 4DP3), (iii) docked complex of wild-type PfDHFR, (iv) docked complex of double-mutant PfDHFR (dm1PfDHFR), and (v) docked complex of doublemutant PfDHFR (dm2PfDHFR), all carrying NADPH as the cofactor and P218 as the bound ligand in their respective sites were used for the MD analysis. MD simulations were carried out for 35 ns to ensure that the equilibration state of the complex lasted long enough for further hydrogen bonding and binding energy analysis. The convergence of all the simulations was analyzed in terms of the energy components, B-factor analysis, and RMSD from the initial X-ray crystal structure.
Initially, the stability of the complexes was evaluated by analyzing the kinetic energy, potential energy, and total energy for all the five complexes during the entire simulation run of 35 ns ( Figure S5, Supplementary material). The energy plots reveal that the changes in total energy, potential energy, and kinetic energy were smooth. The stability of all the five complexes was also analyzed by studying the complex RMSD values of the backbone atoms ( Figure 5(A)) relative to those in the starting structure. The trajectory analysis in Figure 5(A) shows that all the complexes attain equilibrium after 10 ns with acceptable RMSD fluctuations (<2 Å). From the RMSD plot of the ligand P218 ( Figure 5(B)), it can be inferred that P218 binds stably (has very low flexibility) in all the complexes.
P218 binding energies MM-PBSA and normal mode analysis was performed to evaluate the binding free energies (ΔG bind ) of P218 in all the above-mentioned complexes ( Table 1). The experimental binding constants K i for P218 in hDHFR, Figure 4. Superimposed X-ray crystallographic conformation of P218 in hDHFR (sky blue), qmPfDHFR (red), docked conformations of P218 in wtPfDHFR (green), dm1PfDHFR (orange), and dm2PfDHFR (yellow). Mentioned residues are the residues which forms hydrogen bonds with P218 in hDHFR and PfDHFR (wild-type and mutant). wtPfDHFR, and qmPfDHFR were reported in literature   (Table 1). The correlation of these experimental binding affinities (based on K i values) with the estimated binding free energies for three complexes (hDHFR, wtPfDHFR, and qmPfDHFR) showed that the experimental and calculated binding affinity values follow the same trend, i.e. binding affinity of wtPfDHFR ≈ qmPfDHFR > hDHFR. P218 binding energy for double-mutant PfDHFR (dm1PfDHFR and dm2PfDHFR) was also calculated (Table 1) and it was found that P218 binding affinity to dm2PfDHFR was comparable to wild-type and quadruple-mutant PfDHFR whereas a decrease in binding affinity of P218 for dm1PfDHFR was noticed. However, the binding affinity of P218 for all the PfDHFR complexes was estimated to be more than hDHFR. This suggests that P218 binding to PfDHFR is not affected by the existing mutations in the enzyme. Overall, from the ΔG bind values, it can be envisaged that there is a significant difference between P218 binding to human DHFR and PfDHFR. The binding free energies of P218 with wtPfDHFR, qmPfDHFR, and dm2PfDHFR are almost comparable (−34.90, −34.15 and −34.93 kcal/mol, respectively), while decrease in binding energy of P218 with dm1PfDHFR (−29.53 kcal/ mol) was noticed. In contrary, ΔG bind value of P218 in hDHFR was −22.14 kcal/mol. The difference between the overall binding energy for hDHFR as compared to wtPfDHFR (12.76 kcal/mol), qmPfDHFR (12.01 kcal/ mol), dm1PfDHFR (7.39 kcal/mol), and dm2PfDHFR (12.79 kcal/mol) is significantly high. As a result of this binding energy difference between hDHFR and PfDHFR, P218 can bind selectively and efficiently to PfDHFR and prevent it from executing its function of DNA synthesis in wild-type PfDHFR and mutant PfDHFR and subsequently kill the parasite, but does not affect the same process in human DHFR.  In the P218-PfDHFR complexes (wild-type and mutant), the contributions favoring P218 binding are VDW (van der Waals [VDW] contribution), ranging from −40 to −49 kcal/mol, and ELE (electrostatic energy), ranging from −179 to −243 kcal/mol. The nonpolar interaction with the solvent (PBSUR) yielded contributions of~−4 kcal/mol and the polar interaction with the solvent (EPB) yielded contributions in the range of 169-207 kcal/mol. Overall electrostatic contribution favors P218 binding in PfDHFR, with PBELE (sum of electrostatic solvation free energy and MM electrostatic energy) values ranging from −4 to −16 kcal/mol, dm1PfDHFR having the least value of~−4 kcal/mol. This reduction in PBELE value in the case of dm1PfDHFR may be the possible reason for the observed lower P218 binding affinity for dm1PfDHFR as compared to other PfDHFR structures. However, the binding free energy of P218 to hDHFR is comparatively less favorable (−22.14 kcal/mol). The main contribution leading to lower binding affinity in hDHFR originates from the unfavorable energy of desolvation of polar groups having PBELE value of 1.22 kcal/mol (as com-pared to −4 to −16 kcal/mol in PfDHFR). Additionally, the binding of P218 to hDHFR had the least favorable electrostatic interactions, as compared to P218-PfDHFR complexes; its ELE is −131.46 kcal/mol, whereas its non-polar contribution to solvation free energy (PBSUR) was found to be comparable to the PfDHFR complexes.
The VDW contribution is lower in qmPfDHFR than the other complexes but its electrostatic interaction (ELE) and PBELE contribute significantly more as compared to all the other four complexes, with least contribution in hDHFR (Table 1). This difference in electrostatic and VDW terms leads to difference in binding affinities among the studied complexes (Table 1).

Energy decomposition of each amino acid residue
To find the key residues that play important role in P218 binding in the active site of hDHFR and PfDHFR, and also to find out the residues responsible for decreased binding affinity of dm1PfDHFR as compared to other PfDHFR complexes, the interactions were quantified in terms of the pair interaction decomposition of free Figure 6. 'Per-residue' decomposition energy of residues in hDHFR and PfDHFR (wtPfDHFR, qmPfDHFR, dm1PfDHFR, and dm2PfDHFR) having binding energy with P218 larger than 1.0 kcal/mol or the residues with the difference between the complexes larger than .5 kcal/mol. Residues mentioned in human DHFR are structurally equivalent homolog residues corresponding to PfDHFR (as obtained from multiple sequence alignment, Figure 2). energy. The 'per-residue' binding energy of important amino acids that interact with P218, in all the five complexes is shown in Figure 6. From the figure, it can be seen that the patterns of the energy decomposition in wtPfDHFR and dm2PfDHFR are similar and their total binding free energies were also found to be almost equal ( Table 1). The patterns of the energy decomposition in wtPfDHFR and dm1PfDHFR differ by strong interactions of P218 with Asp54 and Cys15 in wtPfDHFR as compared to dm1PfDHFR. Residues Ile14, Arg122, and Ile164 slightly favor P218 binding to dm1PfDHFR and the mutated residue Arg59 significantly favors P218 binding to dm1PfDHFR as compared to wtPfDHFR, whereas, in all other cases, the 'per-residue' contribution is decreased for dm1PfDHFR. It can be noticed from the figure that in comparison to qmPfDHFR, the mutation at position 59 from Cys to Arg, does not strongly favor P218 binding in double-mutant PfDHFR (dm1PfDHFR). Visualization of binding mode of P218 in the active site of dm1PfDHFR revealed that Arg59 in dm1PfDHFR was oriented in such a way that it does not interact stably with P218 ( Figure S6, Supplementary material). This finding can also be related to comparatively low binding affinity of P218 in dm1PfDHFR and low electrostatic contribution to overall binding free energy (Table 1) as compared to other mutant and wild PfDHFR forms. Although the binding affinity of P218 in qmPfDHFR is similar to wtPfDHFR, differences are clear in the 'perresidue' contribution to total binding energy. The residue wise major difference in binding of P218 to qmPfDHFR and wtPfDHFR can be attributed to residues Ile14, Cys15, Cys/Arg59, Arg122, and Ile/Leu164, where on a comparative scale Ile14, Cys15, and Ile164 favor P218 binding in wtPfDHFR and residues Arg59 and Arg122 strongly favor P218 binding in qmPfDHFR. Unlike in dm1PfDHFR, the mutation at position 59 from Cys to Arg, strongly stabilize P218 binding to qmPfDHFR, whereas, mutation at position 164 from Ile-Leu in qmPfDHFR is not favorable and slightly destabilizes P218 binding to the enzyme. On the other hand, the patterns of the energy decomposition of P218-hDHFR are quite different from that of P218-PfDHFR complexes. Figure 6 shows the 'per-residue' decomposition energy of residues in hDHFR that are structurally equivalent homologs corresponding to PfDHFR, as studied from sequence alignment exercise (Figure 3). The residues Gln35, Pro60, and Leu67 in hDHFR corresponding to residues Cys59/Arg59, Pro113, and Leu119 of PfDHFR have significantly low 'per-residue' decomposition energy as compared to PfDHFR. The residues Phe31 and Figure 7. Hydrogen bond occupancies of important active site residues in hDHFR and PfDHFR (wtPfDHFR, qmPfDHFR, dm1PfDHFR, and dm2PfDHFR). Analysis was carried out on the last 5 ns (2500 frames) of simulation run. Angle and distance cut-off for hydrogen bonding was set to 120°(i.e. 120°-180°) and 3.2 Å.
Arg70 in hDHFR corresponding to residues Met55 and Arg122 in PfDHFR, disfavors and might destabilize P218 binding in hDHFR. This difference in 'per-residue' decomposition energy in hDHFR and PfDHFR (wildtype and mutant) can be correlated with poor overall binding affinity of P218 to hDHFR as compared to PfDHFR. These findings support the experimental results  and also provide additional insights into the unfavorable binding of P218 to hDHFR.

Hydrogen bond analysis
Binding affinity of a ligand in the active site of an enzyme can be correlated to hydrogen bonding interactions. Hydrogen bond (H-bond) occupancy between P218 and residues Ile14, Asp54, Cys/Arg59, Arg122, and Ile164 of PfDHFR and residues Ile7 and Glu30 of hDHFR, during the last 5 ns simulation run are shown in Figure 7. Hydrogen bond occupancies of >75%, 50-75% and <50% were used to define strong, medium, and weak H-bond interactions, respectively (Khuntawee, Rungrotmongkol, & Hannongbua, 2011;Mokmak et al., 2014). Any residue having hydrogen bond occupancy of >100% means that it can form more than one hydrogen bond at a time. In the case of PfDHFR (wild-type and mutant), P218 forms two strong hydrogen bonds with Asp54 (between the N7 amino group of P218 and O1 of Asp54 and between N1 of P218 and O2 of Asp54) having occupancy >150% and two hydrogen bonds with Arg122 (between the phenolic carboxylic O1 of P218 and N1H of Arg122 and between phenolic carboxylic O2 of P218 and N2H of Arg122) with occupancy >100% (Figure 7). The H-bond distances of the abovementioned hydrogen bonds during the last 5 ns simulation run ( Figure S7(A)-(D), Supplementary material) also depict that these bonds were stable during the simulation and thus strongly stabilize P218 binding in the active site of PfDHFR (wild-type and mutant). In case of wtPfDHFR, dm1PfDHFR and dm2PfDHFR, N8 amino group of ligand and backbone carboxylic group of Ile14 formed medium H-bond while it formed weak H-bond with the backbone carboxylic group of Ile164 (Figure 7) as also shown by hydrogen bond distances over last 5 ns simulation run ( Figure S7(E)-(F), Supplementary material). Quadruple-mutant PfDHFR formed weak hydrogen bond between N8 amino group and backbone carboxylic group of Ile14, while Leu164 does not form H-bond interaction with P218 ( Figure S7(E)-(F), Figure 8. 'Per-residue' VDW energy of residues in hDHFR and PfDHFR (wtPfDHFR, qmPfDHFR, dm1PfDHFR, and dm2PfDHFR) having absolute binding energy with P218 larger than 1.0 kcal/mol or the residues with the difference between the complexes larger than .5 kcal/mol. Residues mentioned in hDHFR are structurally equivalent homolog residues corresponding to PfDHFR (as obtained from multiple sequence alignment, Figure 2).
Supplementary material). Double-mutant dm1PfDHFR and qmPfDHFR have Cys59Arg mutation. Unlike dm1PfDHFR, qmPfDHFR forms additional strong and stable hydrogen bonds between phenolic carboxyl group and Arg59 ( Figure S7(G), Supplementary material). Cys15 is the other amino acid involved in weak hydrogen bond in case of PfDHFR. The pattern of strong H-bonds in all the P218-PfDHFR complexes is similar except for extra strong hydrogen bonds involving Arg59 in qmPfDHFR. For the P218-hDHFR complex, strong Hbonds are formed between the N7 amino group and Glu30 (structural homolog of Asp54 in PfDHFR), similar to the observation in PfDHFR (Figure 7 and S7(A)-(B), Supplementary material), whereas it forms medium hydrogen bond between N8 amino group and Ile7 (Figures 7 and S7(E), Supplementary material). Additionally, Asn64 (structural homolog of Cys59 in PfDHFR) was found to form weak hydrogen bond interaction with phenolic carboxyl group of P218 during the simulation ( Figure S7(H), Supplementary material). These results showed that the wild-type PfDHFR and dm2PfDHFR have similar H-bond profile which involves strong hydrogen bond with Asp54 and Arg122; medium hydrogen bond with Ile14 and weak hydrogen bond with Ile164, while dm1PfDHFR differs only in terms of extra weak hydrogen bond with Arg59. H-bonding profile of qmPfDHFR involves strong H-bond with Asp54, Arg59, and Arg122, weak hydrogen bond with Ile14, and no hydrogen bond with Leu164. Human DHFR forms strong hydrogen bond only with Glu30, medium H-bond with Ile7, and weak H-bond with Asn64 (this interaction was also demonstrated in the docked pose of P218 in the active site of hDHFR, Figure 3) during the simulation run; hence, the binding affinity of P218 in hDHFR is reduced as compared to PfDHFR. This result is in agreement with the contribution of individual energy components to binding energy (Table 1), where major difference to P218 binding to hDHFR and PfDHFR is due to electrostatic component. On a relative scale, it was noticed that hydrogen bond occupancy in dm1PfDHFR for residues Ile14, Asp54, Arg122, and Ile164 was lower than wtPfDHFR. The hydrogen bond occupancy in dm1PfDHFR for Arg59 was also significantly reduced as compared to qmPfDHFR. This may also account to decreased electrostatic contribution to binding free energy of P218 for dm1PfDHFR as compared to other PfDHFR complexes. Based on the hydrogen bond occupancies during the last 5 ns simulation, it was noticed that P218 formed three-four and five-six hydrogen bonds with hDHFR and PfDHFR, on an average, respectively (Figure S8, Supplementary material).

Hydrophobic contact
The VDW term of the MM-PBSA (Table 1) mostly represents the non-polar interaction: the hydrophobic contact between the molecules. This can also be found by the decomposed energy contribution. The non-polar contributions of selected amino acid residues in PfDHFR and their corresponding homologs in hDHFR are given in Figure 8. Residues Met55 and Phe58, in PfDHFR and their corresponding homologous residues Phe31 and Phe34 in hDHFR play an important role in hydrophobic contribution. While the residues in wtPfDHFR, qmPfDHFR, dm1PfDHFR, and dm2PfDHFR generally had comparable values but there exists small differences in the 'per-residue' hydrophobic contribution in hDHFR. Almost comparable 'per-residue' contribution to hydrophobic interaction justifies the above-mentioned results (Table 1), where hydrophobic contribution to binding energy in all the complexes was found to be similar. Overall, on a comparative scale, the hydrophobic interactions were slightly more prevalent in PfDHFR relative to hDHFR.

B-factor analysis
As indicated by the B-factors calculated over the last 5 ns trajectories (Figure 9), it was noticed that the flexible regions in all the five complexes correspond to loop regions. In the case of wtPfDHFR and qmPfDHFR, the most flexible regions correspond to the loop regions (residues 20-30 and residues 75-100). The flexibility is more pronounced in the case of qmPfDHFR as compared to wtPfDHFR. In the case of double-mutants, the flexibility again corresponds to the loop regions (residues 20-30, 75-100, and 125-140) with more pronounced flexibility in case of dm2PfDHFR.

Conclusions
The rise of resistance in PfDHFR is predominantly attributed to the steric constraints caused by specific mutations, the mutation Ser108Asn, being the prime cause of the loss in binding affinity of rigid antifolates like pyrimethamine, cycloguanil, and chlorcycloguanil. Recently, a new inhibitor P218, of resistant and wild-type PfDHFR strains, which does not inhibit human DHFR, has been reported and is in preclinical trials . To evaluate the detailed interaction profile of P218 with human DHFR and to find the affect of prevailing mutations in PfDHFR on P218 binding, MD simulations of P218 with human DHFR, wild-type PfDHFR, quadruple-mutant PfDHFR (Asn51Ile, Cys59Arg, Ser108Asn, and Ile164Leu), double-mutant PfDHFR, dm1PfDHFR (Cys59Arg and Ser108Asn), and double-mutant PfDHFR, dm2PfDHFR (Ala16Val and Ser108Thr) were carried out.
The binding mode analysis of P218 in the active sites of PfDHFR and hDHFR, as carried out by molecular docking studies revealed that P218 binds in a similar fashion to wild-type and mutant PfDHFR (quadruplemutant and double-mutants) but differently to hDHFR. It was noticed that diaminopyrimidine ring of P218 in the active site of PfDHFR forms hydrogen bonds with Ile14, Asp54, and Ile/Leu164, whereas it forms hydrogen bonds with Ile7 and Glu30 in the active site of hDHFR. In addition, the terminal carboxyl group of P218 showed hydrogen bond interactions with Arg122 in case of PfDHFR but no such hydrogen bond was observed between Arg70 of hDHFR (structural homolog of Arg122 in PfDHFR). This may lead to the lower binding affinity of P218 to hDHFR as compared to the affinity in PfDHFR. Quadruple-mutant PfDHFR and dm1PfDHFR have Cys59Arg mutation. In the case of qmPfDHFR, Arg59 showed strong H-bond interaction with terminal carboxyl group of P218 but such interaction was diminished in the case of dm1PfDHFR. As all these interactions were studied in static structures further to check whether these interactions were stable over a period of time, MD simulations for 35 ns were carried out.
The MD simulation studies confirmed the selectivity and efficiency of P218 against wild-type and quadruplemutant PfDHFR, in agreement with experimental data and also it was found that P218 binding affinity to the double-mutants is not much affected by the mutations. Calculated binding free energy of P218 with hDHFR (−22.14 kcal/mol) is significantly lower than wtPfDHFR (−34.90 kcal/mol); qmPfDHFR (−34.15 kcal/mol); dm1PfDHFR (−29.53 kcal/mol); and dm2PfDHFR (−34.93 kcal/mol). It was also noticed that majorly electrostatic contribution can be attributed to the differences in binding free energies among these complexes.
Detailed hydrogen bonding analysis during the last 5 ns simulation run revealed that P218 in the active site of PfDHFR binds stably, and in the entire simulation run, it was anchored at its both ends i.e. its diaminopyrimidine ring interacts with Ile14, Asp54, and Ile164 and its terminal carboxylic group forms strong hydrogen bonds with Arg122 (and Arg59 in case of qmPfDHFR). Whereas, in human DHFR, P218 binding is stabilized only by the interaction of its diaminopyrimidine ring Ile7, Glu30 at its one end, terminal carboxylic group at opposite end does not interact with Arg70 (structural homolog of Arg122 in PfDHFR). The mutations in PfDHFR cause no significant loss in interaction energy with P218 inhibitor, rather mutated residues additionally forms stable electrostatic and VDW interactions with P218. The atomic-level details studied in the present work, on binding pattern of P218 can be utilized to design novel and potent inhibitors of PfDHFR which binds specifically to PfDHFR and are not affected by mutations.

Supplementary material
The information on binding mode of WR99210 in wildtype and quadruple-mutant PfDHFR; difference between energy minimized and bioactive conformation of P218 and conformational analysis of P218; RMSD of the redocked ligands; kinetic energy, potential energy, and total energy graphs of all the five complexes during the simulation run; different orientation of Arg59 in qmPfDHFR and dm1PfDHFR; and detailed H-bond analysis in all the five complexes are given in supplementary material. The supplementary material for this paper is available online at http://dx.doi.10.1080/07391102.2014.979231.