|Home | About | Journals | Submit | Contact Us | Français|
Drug resistance acquired by Plasmodium falciparum (Pf) is a major problem in the treatment and control of malaria. One of the major examples of drug resistance is that caused by mutations in the active site of dihydrofolate reductase (DHFR) of Pf (PfDHFR-TS). A double mutation, A16V+S108T, is specific for resistance to the marketed drug cycloguanil. In this study, we used 58 cycloguanil (2,4-diamino-1,6-dihydro-1,3,5-triazine) derivatives to explore the relationship between various physico-chemical properties and reported binding affinity data on wild type and A16V+S108T mutant type. Using the Hansch 2D-QSAR method, we obtained a parabolic relationship of hydrophobicity of substituents at the N1-phenyl ring with the wild type binding affinity data. Hydrophobicity being a key property for wild type binding affinity data, we found steric factors to be crucial for A16V+S108T mutant resistance. We investigated FlexX, GOLD, Glide and Molegro virtual docking programs and 13 different scoring functions on 10 of the cycloguanil derivatives to evaluate which program was best for reproducing the experimental binding mode and correlating the docking scores with the reported binding affinity data. We identified GOLD using its GoldScore fitness function as the most accurate docking program for predicting binding affinity data of cycloguanil derivatives to DHFR and Molegro virtual docker with its template docking algorithm and MolDock [GRID] scoring function as most accurate for reproducing the experimental binding mode of a reference ligand that is structurally similar to the cycloguanil derivatives studied. We also report an interaction index which best describes the structure-activity relationships exhibited by these analogs in terms of PfDHFR-TS active site interactions.
Classified as one of the leading causes of death in the world, malaria is a major threat for public health, especially in the tropical and sub-tropical countries. Malaria has been declared endemic particularly in the sub-Sahara African countries, Asia, Oceania and Latin America.1 Malaria is caused by Plasmodium parasites, of which the most lethal is Plasmodium falciparum (Pf). The biology of this parasite is well understood with the help of genomics2 and proteomics projects3 and has aided identification of new molecular targets that could be used to aid rational design of drugs or vaccines. Drug resistance acquired by the malarial parasites is a major problem in the treatment and control of malaria. Dihydrofolate reductase (DHFR) in Pf (PfDHFR-TS; TS refers to thymidylate synthase, bound to DHFR in Pf), which catalyzes the reduction of dihydrofolate to tetrahydrofolate, is one of the most widely studied enzymes in antimalarial drug design due to its potential role in DNA synthesis.4, 5 DHFR is also considered to be a good target for other protozoal diseases like leishmaniasis, trypanosomiasis and Chagas’ diseases.6 Resistance of the parasite to antifolates has been observed and is mainly caused by mutations around the enzyme active site.7–9 Some of the resistance-causing mutations include single S108N, double C59R+S108N, triple N51I+C59R+S108N, C59R+S108N+I164L and quadruple N51I+C59R+S108N+I164L. A double mutation, A16V+S108T, is specific for resistance to the marketed drug cycloguanil. Over the past few years, several computational techniques have been applied to better understand PfDHFR-TS-ligand interactions.10 Currently there are four crystal structures of protozoal DHFR available in the Protein Data Bank, including the wild type PfDHFR-TS complexed with triazine inhibitor WR99210 (1j3i),11 a double mutant PfDHFR-TS complexed with pyrimethamine (1j3j),11 a quadruple mutant PfDHFR-TS complexed with WR99210 (1j3k)11 and wild type Plasmodium vivax DHFR-TS complexed with pyrimethamine (2bl9).12 Before the release of PfDHFR-TS crystal structures, researchers developed homology models to study the resistance mechanism of commonly used antifolate drugs.13 Recently, Fogel et al. performed molecular modeling studies including docking studies with the GOLD program on pyrimethamine derivatives into the 1j3k DHFR structure, evaluated several scoring functions, including those from GOLD, Molegro virtual docker, Discovery Studio and MOE and found Molegro virtual docker’s Protein-Ligand interaction score to show good correlation with the reported binding affinity data.14 The same authors performed docking studies of pyrimethamine derivatives into both wild type and mutant enzymes.15 QSAR studies using neural networks on this class of compounds have also been performed by the same authors.16 4D-QSAR studies for some of the antifolates have been published.17 A 3D pharmacophore model has been prepared using Catalyst software on a diverse set of PfDHFR-TS inhibitors, including cycloguanil and pyrimethamine derivatives.18 The Autodock program has also been shown to predict the correct binding modes of folic acid-competitive DHFR antagonists identified by Plouffe et al. through their cell-based high-throughput screening.19
In summary, previous modeling studies focused on studying the interactions of various small molecule inhibitors of DHFR with the active site residues of wild and mutant types using different arbitrarily chosen docking programs and QSAR studies. However there is need for new research on the importance of physico-chemical properties such as the linear and/or bilinear relationship of hydrophobicity with binding affinity data and on systematic attempts to find a good docking program/scoring function to study molecular interactions accurately for this highly interesting and therapeutically important class of compounds. Hence our report has two main objectives: (1) to gain new insights into the structure activity relationship (SAR) among the cycloguanil derivatives binding to wild and mutant type PfDHFR-TS and (2) to find an optimal docking program/scoring function combination to study the binding interactions between these compounds and wild type PfDHFR-TS that are not revealed by the SAR studies. For objective 1, we made use of QSAR to study a series of 58 cycloguanil (2,4-diamino-1,6-dihydro-1,3,5-triazine) derivatives20–22 which all originated from the same laboratory, since well-defined binding constants against wild and mutant types of PfDHFRTS have been reported for them. We chose to use traditional Hansch 2D-QSAR analysis for this series of compounds for two main reasons. First, 48 of the 58 compounds had an undefined stereogenic center at C6 of the triazine ring of cycloguanil, which rendered the data unsuitable for 3D-QSAR studies. Secondly, Hansch 2D-QSAR is well known for its simplicity and ease of mechanistic interpretation. For the docking studies, we used the only 10 compounds with well-defined stereochemistry to study interactions with the wild type PfDHFR-TS that could not be delineated from our 2D-QSAR studies. For objective 2, we evaluated FlexX, GOLD, Glide and Molegro virtual docking programs. We then investigated in detail the interactions of the analogs with the active site of PfDHFR-TS, using what we found to be the most suitable docking program for this class of protein-ligand interactions. The two objectives of this work fit together well since the methods of 2D-QSAR and docking give complementary information on the structure-activity relationship.
2D-QSAR models were constructed using the data set of 58 compounds.20–22 The compounds used for the present analysis along with their observed affinities for wild and mutant type PfDHFR-TS are listed in Table 1. We used two types of standard descriptors: indicator variables and physico-chemical constants. A complete list of all the descriptors used in 2D model generation is given in Table 2. An indicator variable designated as I with a relevant subscript was set to 1 if a particular substituent is present and to 0 if absent. R1 and R2 are both equivalently attached at C6, so to make sure the indicator variables were well defined, we ordered the R1/R2 substituents so that for any analog the less bulky substituent is at R1. Physico-chemical constants—hydrophobic (π); electronic (σ), including resonance effect (R) or field/inductive effect (f); molar refractivity (MR); and hydrogen bond donor (HDR) or acceptor (HAR)—were taken from the literature.23, 24 The substituent hydrophobicity constant (π) which we used from the literature is in fact measured experimentally. The hydrophobicity constant (πX) for any substituent X is calculated according to the equation, πX = logPX − logPH, where PX and PH are the partition coefficients of the reference compound with and without substituent, respectively. The Hammett electronic parameter (σ) is calculated based on the influence of substituent X on the ionization of benzoic acid. The molar refractivity (MR), defined in the Lorenz-Lorentz equation as , where n is the refractive index, MW is the molecular weight and d is the density, is considered to be a measure of the volume occupied by an atom or group of atoms according to Hansch and co-workers, since usually there is only very slight variation in the refractive index between a set of analogs.25 A correlation matrix was used to correlate the biological activity with physico-chemical and structural predictor variables. Descriptors with inter-correlation |r| > 0.6 were not included in the same QSAR model. Predictor variables with p > 0.05 were eliminated whilst deriving the QSAR models, in order to assure statistical reliability. Linear regression analysis was performed using Systat Version 11. The Student’s t-distribution was used to assess the significance of individual regression terms. Forward stepping regression was used to build the QSAR models. This method initially generates a QSAR model containing only one variable, which is chosen to be the one with the highest t-statistic, and subsequent variables are added based on their relative importance, also as determined by t-statistics. The QSAR models were evaluated using standard statistical parameters including the correlation coefficient (r), coefficient of determination (r2), standard error of estimate (s) and Fisher F-value. The figures within the parentheses following the coefficient terms are the standard error of regression terms and constants. A data point was considered as an outlier if its residual value exceeded twice the standard error of estimate of the model. Self-consistency of the derived models was verified using the leave-one-out (loo) process, and the predictability of each model was assessed using cross-validated r2, called q2. The best measure of reliability of a QSAR model is a high q2, not just a high r2 which could be a result of over-fitting to data. In general, a value of q2 > 0.5 is considered acceptable.26, 27
For docking studies, we used 10 of the 58 compounds (indicated in Table 1) to study the interaction of these compounds with the active site residues of PfDHFR-TS, information not uncovered from our QSAR studies. The 10 were the only compounds with well defined stereochemistry at C6 of the triazine ring among the 58 cycloguanil derivatives used in the QSAR studies. For all the docking studies we used wild type PfDHFR-TS complexed with WR99210 (1j3i).11 The structure of this reference ligand is shown in Table 1. The whole tetrameric protein was used. All water molecules and the co-crystallized ligand, WR99210, were removed before docking. All the docking analyses were confined to chain A of 1j3i. Hydrogen atoms were added to the protein for each type of docking. The protein was optimized to remove any bad contacts only for Glide docking (it is a part of the standard protein preparation in the Glide program). The 10 ligands were each constructed starting from an overlay with the bioactive conformation of WR99210 and then minimized using MMFF94 in Sybyl 7.2. We used two types of measures to evaluate the accuracy of the docking programs. First, we considered the ability of the docking program to reproduce the experimental binding mode of WR99210. For this we used the minimized structure of WR99210 and docked it along with the 10 other ligands during each docking run and then compared the RMSD between the predicted binding mode and the experimental structure of WR99210. We used an all-atom-in-place superposition using Schrödinger 8.5 suite of programs to calculate the RMSD. The second measure we used was R2 for correlation between the docking scores and the experimental binding affinity data.
FlexX Release 2.0 implemented in Sybyl 7.228 was used as one of the four docking programs investigated. We prepared an active site comprising of all complete amino acids containing any atom within 6.5 Å of any ligand atom of WR99210. We considered four scoring functions, DScore,29 GScore,30, 31 CScore,32 PMF Score33 and the composite FlexX Total scores. 30 poses for each ligand were generated and all were considered for the RMSD calculation, whereas the scores from only the top-most pose were used to correlate with the reported binding affinity.
We used the Glide program from Schrödinger 7.5 software, which is well known for its speed and accuracy.34–36 We used the standard precision (SP) module of Glide. The Glide GScore and EModel which are based on the ChemScore function of Eldridge et al.32 were considered. The protein was prepared using the protein preparation and refinement tool in Glide. For the active site, a grid box centered at the ligand, WR99210, was used, set to accommodate a maximum ligand length of 15 Å. We docked the ligands flexibly, allowing for the flip of 5- and 6-membered rings, writing out a maximum of 30 poses per ligand and also enabling the post-docking minimization of the ligands. We also evaluated the performance of Ruvinsky’s correction to GScore and EModel.37, 38 This uses a number of docked poses within each energy well and calculates a refined docking score. For Ruvinsky’s correction treatment, we generated 100 poses per ligand using the SP mode of Glide. These poses were then clustered and the entropy was calculated for each cluster. This was accomplished using the pose_entropy script (revision 1.8) available within the Schrödinger suite of programs.
We used Molegro virtual docker, which has been recently introduced and gained attention among medicinal chemists.39, 40 Due to the availability of the bioactive conformation of a structurally-related compound, WR99210, we used the template docking available in Molegro virtual docker and evaluated MolDock, rerank and protein-ligand interaction scores from MolDock and MolDock [GRID] options. Template docking is based on extracting the chemical properties like the pharmacophore elements of a ligand bound in the active site and using that information for docking structurally similar analogs. We used W99210 from 1j3i as the template. We used the default settings, including a grid resolution of 0.30 Å for grid generation and a 15 Å radius from the template as the binding site. We used the MolDock optimizer as a search algorithm, and the number of runs was set to 10. A population size of 50, maximum iteration of 2000, scaling factor of 0.50, crossover rate of 0.90 and a variation-based termination scheme for parameter settings were used. The maximum number of poses to generate was increased to 10 from a default value of 5. Since the Molegro virtual docker works by an evolutionary algorithm, consecutive docking runs do not give exactly the same fitness scores. To address this issue of inherent randomness, we used five consecutive runs and then used the top ranked score from each run to calculate the average score for each ligand to correlate with the experimental binding
We used the GOLD 4.0 docking program.41 The active site was defined as any atom that lies within a 15 Å radius of the ligand, WR99210. Since GOLD, like Molegro, uses a genetic algorithm to dock ligands, we performed five consecutive docking runs to obtain an average docking score for each ligand. Default GA settings that ensure 100% search efficiency were used for the genetic algorithm (GA). An early termination of the number of GA runs was allowed when the RMSDs of the top three GA solutions were within 0.5 Å. We evaluated GoldScore, ChemScore and Astex Statistical Potential (ASP) scoring functions.42 ASP is an atom-atom potential that is based on a database of protein-ligand complexes. Its accuracy is comparable to GoldScore and ChemScore.42
Hansch 2D-QSAR models are listed in Table 3, and the corresponding statistical parameters of the regression equations are given in Table 4. Models 1–2a were developed by considering the observed binding affinity data on the wild type enzyme as a dependent variable (pKi(wt)) and the derived physico-chemical and indicator variables as independent variables. Model 1 explains 80.2% variance in observed binding affinity data on wild type. We found 22 to be a serious outlier when deriving Model 1. Model 1a, prepared without 22, had significantly improved statistical parameters such as much higher r2 compared to Model 1, supporting that 22 is an outlier. Model 1a explains 90.4% variance in the observed binding affinity. The coefficients for the descriptors are nearly the same for Model 1a as those for Model 1, so we analyze in detail Model 1a. None of the descriptors used in Model 1a are auto-correlated with each other. Model 1a gave the order of importance of the variables as I4 > I7 > I2 > I5 > Σ;πR > I6 > I8 > I3. It is quite obvious that other than hydrophobicity, which is a very common feature for compounds acting on DHFR, certain structural features of these analogs are crucial for PfDHFR-TS inhibition. Published results show that substituents at the C6 position interact with Ala16 in the wild type enzyme.9 For example, one of the C6-methyl groups of cycloguanil showed hydrophobic interactions with this residue.9 The negative regression coefficient of the indicator variable I2 in Model 1a indicates that a hydrogen atom at both C6 positions on the triazine ring eliminates the attractive hydrophobic interactions which other substituents can make with the adjacent alanine side chain. The negative contribution of I4, showing the untoward effect of the isopropyl group at the R2 position, is probably due to a steric clash of this bulky hydrophobic group with the alanine side chain. On the other hand, the positive values of the regression coefficients for other descriptors such as I3 (Me at R1) or I5, I6, I7 and I8 [for 3-phenoxyphenyl, 3-benzyloxyphenyl, 4-phenoxyphenyl and 3-(4-chlorophenoxy)phenyl at R2, respectively] indicate that some bulky substituents increase the binding affinity of the ligands towards the wild type. This in part may be because the compounds in the data set with bulky aryloxy R2 substituents have only a hydrogen atom at the R1 position and have either 3- or 4-chloro substitution at the N1-aryl ring (R position) (note that cycloguanil possesses a 4-chloro group at this position).
The positive coefficient of the ΣπR term of Model 1a shows that hydrophobic substituents at R increase the binding affinity of the ligands to the wild type. The importance of the hydrophobicity of these anti-folates is well documented in the literature and has previously been studied using Hansch 2D-QSAR.43–45 We then tested for any non-linear relationship between hydrophobicity and the binding affinity data on the wild type. Model 2, developed using the square of ΣπR, shows that adding in this new descriptor is not statistically acceptable, since the standard error of the regression term for this descriptor is more than the regression coefficient, making a low t (two-tail) value considerably lower than the tabulated value at a 95% confidence interval. An explanation for this is that Model 2 suffers from the presence of 2 outliers, namely compounds 5 and 22. Model 2a was developed after eliminating 5 and 22 as outliers. Model 2a shows that the addition of the (ΣπR)2 term is statistically acceptable at greater than 95% confidence interval. The occurrence of a non-linear relationship is not very surprising as several bilinear relationships between hydrophobicity (using π) and DHFR binding affinity data for various species have already been reported in the literature.44, 45 The optimum π values for meta substituents of the N1-aryl ring of cycloguanil analogs towards chicken, human, rat, murine, E. coli, Leishmania casei, L. major and P. carinii were 5.00, 2.10, 5.00, 1.76, 3.00, 4.31, 4.54 and 4.00, respectively. According to our literature search, there is no report about optimal hydrophobicity requirements for cycloguanil analogs against PfDHFR-TS binding affinity. Hence, we calculated the optimum value of the hydrophobicity of substituents at R for better binding affinity towards the wild type using simple differentiation of the parabolic equation (Model 2a) with respect to ΣπR, treating contributions from all other indicator variables as constant, and found d[pKi(wt)]/d[ΣπR] = −0.349*2(ΣπR) + 1.303 = 0.61 as the optimum πR value. Most of the potent compounds possess either a 4-chloro or 3-chloro substituent (π = 0.71 for Cl); however most of the compounds containing di-chloro substituents are inactive (π = 1.42 for di-Cl because of the additive effect of hydrophobicity). Neither a methyl substituent (π = 0.56, which is low) nor the ethyl substituent (π = 1.02, which is too high) are optimal for binding affinity. Substituents like CF3, CH2Br, SCH3, SeCH3, CH=CH2, SCOC2H5, NHC3H7 with π values 0.88, 0.79, 0.61, 0.74, 0.62, 0.82, 0.64, respectively, would be some better choices, depending upon the synthetic feasibility of these analogs. Some thioalkyl, selenium and CF3 containing-derivatives have already been reported to show promising activity on Pneumocystis carinii and human DHFR binding affinity.44, 45 The low coefficient of this π term shows that there exists a total desolvation of substituents in the PfDHFR binding site. Since this descriptor is not a whole molecular property, it is difficult to say whether the hydrophobic contributions of the C6-position substituents are important or not; but certainly the hydrophobic contribution from the N1–aryl ring substituents are additive and an optimal value is required for enhanced binding affinity to the wild type.
Model 3 was obtained by considering binding affinity data reported on the A16V+S108T mutant form of the enzyme, pKi(mut), as the dependent variable and the physico-chemical and indicator variables as independent variables. Model 3a was obtained by removing 22 as an outlier from Model 3. The positive contribution of the indicator variable I1 in Model 3a suggests that having the R1 position unsubstituted is favorable for enhancing binding to the mutant form of the enzyme. This is consistent with the SAR observations of Vilaivan et al.,21 who showed that removal of one of the dimethyl groups of cycloguanil at the C6 position results in better binding affinity to the mutant form. The reason, as pointed out by that research group, is that one of the methyl groups experiences a steric clash with the Val16 side chain in the A16V+S108T mutant (referred to as the “steric constraint hypothesis”9). Analogous to our observation for the wild type, a bulky isopropyl group (reflected by the negative contribution of the indicator variable, I4) is detrimental for binding to the mutant form. The positive contributions of the indicator variables I5, I6, and I7 show that 3-phenoxyphenyl, 3-benzyloxyphenyl and 4-phenoxyphenyl at R2 are favorable for enhanced binding affinity. However the strong binding of these analogs is in part due to the presence of H as the other substituent (at the R1 position) of C6, a fact that was already described above for Model 1a. B1Para(R) is a measure of the width of the first atom of the R substituents. Its negative coefficient indicates that bulkier para substituents around the N1-aryl ring (R) are detrimental for binding to this mutant type. This is again because of steric constraints imposed by Ser108. Most of the active compounds for this mutant type possess, instead, meta substituents.
Compound 22, an outlier in Models 1–3, is the least active among the series, and possesses a bulky t-butyl group at the R2 position. Interestingly, a compound with a t-butyl group at the meta position of the N1-aryl ring of cycloguanil was noted as an outlier in the QSAR study of cycloguanil analogs for Pneumocystis carinii DHFR inhibition by Marlowe et al.45 It was pointed out by those authors that the t-butyl group’s interaction with Pneumocystis carinii or Leishmania major DHFR is more highly twisted than in bacterial and mammalian DHFRs. The probable reason suggested for such behavior is that certain bulky amino acids like Trp, Ile or Phe restrict the entry of cycloguanil analogs to the hydrophobic binding pocket of DHFR in the former species; the inactivity of 22 suggests that this also occurs for Pf DHFR. Since compound 5 does not have any unusual substitution pattern and its binding affinity can be captured adequately by the descriptors used in Model 2, the outlying behavior of this compound may be due to experimental error. The reported binding affinity data (Ki) for these compounds were an average of triplicated assay results and there was a relatively high standard error reported for compound 5 (Ki = 1.5±0.8) compared to for other compounds 9 (Ki = 1.4±0.2), 16 (Ki = 1.6±0.2), 17 (Ki = 1.5±0.3) and 57 (Ki = 1.4±0.3) with Ki’s in a similar range to that for compound 5. The major results derived from the Hansch 2D-QSAR models are summarized in Figure 1.
A summary of all the average RMSDs and the R2 between docking scores and reported binding affinities is given in Table 5. Here we present the results for the 4 docking algorithms.
FlexX docking scores, including the Total score, GScore, ChemScore, PMF Score and DScore, all resulted in a poor correlation of R2 < 0.3 with the reported binding affinity data. FlexX predicted consistent binding modes for 9 of the 10 cycloguanil derivatives. Compound 39, which is unsubstituted at R1 and R2 and with a 4-F substitution at R, was predicted to adopt a different binding mode than the other 9 compounds. In contrast, the only other compound with such unsubstituted R1 and R2 and with a 4-Cl substitution at R was predicted to bind in a similar way to the other 9 compounds. This shows that the predictions of FlexX are inconsistent because such a big change in the binding mode is not expected to be caused by the small steric change from the F to the Cl atom at the R position, considering the space available in the active site. The average RMSD between the 30 docked poses from FlexX and the crystal structure of W99210 was 4.971 Å. The RMSD between the top-ranked of the 30 poses and the experimental structure was 3.3 Å. This high RMSD shows that FlexX was not capable of reproducing the experimental structure. Together with the performance of the docking scores, it can be concluded that FlexX as implemented in Sybyl 7.2 is not an accurate docking program to study this class of protein-ligand interactions. The binding modes of the 10 cycloguanil derivatives and the reference ligand, W99210, are given in Supplementary Figure S1. The individual docking scores for the 10 cycloguanil derivatives from FlexX and the RMSD values are given in Supplementary Tables S1 and S2, respectively.
The Glide SP module predicted consistent binding modes for all 10 compounds, matching the reference ligand binding mode. Ruvinsky’s treatment of the docking scores also resulted in similar binding modes of 9 compounds, but one compound, 17, for no apparent reason failed to refine through this procedure. The correlation R2 between the Glide GScore and the reported binding affinity data for the 10 compounds was 0.416. The correlation R2 between the Glide EModel and the reported binding affinity data was 0.426. As can be seen from the correlation coefficient values, neither the Glide GScore nor the EModel was satisfactorily able to predict the binding affinities of these analogs. The Glide GScore and EModel refined using Ruvinsky’s correction were still worse than the corresponding unprocessed GScore and EModel from Glide. R2 between the reported binding affinity data and the Ruvinsky refined Glide GScore and EModel for the nine compounds were 0.167 and 0.196, respectively. The average RMSD between the experimental and 30 predicted binding poses of W99210 was 3.399 Å. But the RMSD of the top-ranked pose among the 30 poses of the reference ligand and the experimental structure was only 1.734 Å, and hence, it can be concluded that the top-ranked pose from Glide can be confidently used for further calculations, since RMSD < 2 Å. The top-ranked pose of each ligand and the reference molecule from a Glide SP run before and after Ruvinsky’s treatment are given in Supplementary Figure S2. The individual docking scores for the 10 cycloguanil derivatives from Glide SP and the RMSD values are given in Supplementary Tables S3 and S4, respectively.
The chemical properties used in the bound ligand, W99210, which was used for template docking using Molegro virtual docker are given in Figure 2. Four hydrogen bond donor atoms (purple), three hydrogen bond acceptors (green) and 12 ring atoms (yellow) of the reference ligand, as shown in Figure 2, were used to dock other compounds. Molegro predicted consistent binding modes for all 10 compounds, matching the binding mode of the reference ligand, just as was found with Glide. The predicted binding modes are given in Figure 3. Three different scores can be obtained using two different scoring methods, MolDock and MolDock [GRID], so a total of six sets of scores were investigated using Molegro, and the results are given in Supplementary Table S5. Each score in the table is an average of five consecutive runs. The correlation R2 between the reported binding affinity data of the 10 cycloguanil derivatives and the MolDock, rerank and Protein-Ligand interaction scores from the MolDock scoring method were 0.741, 0.740 and 0.740, respectively. The correlation R2 between the reported binding affinity data of the 10 cycloguanil derivatives and the MolDock, rerank and the Protein-Ligand interaction scores from the MolDock [GRID] scoring method were 0.735, 0.703 and 0.739, respectively. This shows that Molegro docking scores were approximately equally successful in predicting the binding affinity, and that for docking scores Molegro performed better than FlexX or Glide. Also, since hydrogen bonding interactions do not influence the variation in the binding affinity data exerted by these analogs (discussed in detail below, under active site interactions), the correlation obtained with the two scoring functions used, MolDock and MolDock [GRID], which differ only by taking hydrogen bond directionality into account, was similar. The RMSDs obtained using the two different scoring functions are given in Supplementary Table S6. The average RMSDs were 2.14 and 2.18 Å, respectively. We found an impressively low RMSD = 1.02 Å between the top-ranked pose and the experimental structure of W99210 using the MolDock [GRID] scoring function. This clearly demonstrates that the Molegro virtual docker was very accurate in reproducing the experimental binding mode and the software also gave a reasonably accurate prediction of binding affinity of the cycloguanil derivatives.
Using GOLD software, among the three available docking scores investigated, only GoldScore predicted a consistent binding mode for all 10 cycloguanil analogs which was similar to that of the reference ligand. The correlation R2 between the reported binding affinity and the GoldScore, ChemScore and ASP were 0.911, 0.855 and 0.626, respectively. Hence, the order of accuracy of the scoring functions for predicting the binding affinity was GScore > ChemScore > ASP. The average RMSD between the predicted binding modes of W99210 and the experimental structure using GoldScore, ChemScore and ASP were 2.99, 3.88 and 1.96 Å, respectively. The order of accuracy for reproducing the experimental binding mode was ASP > GoldScore > ChemScore, and this order was maintained when the RMSDs of only the top-ranked poses were considered. This data demonstrates that the GoldScore scoring function predicted the binding affinity of the cycloguanil analogs very well, whereas the ASP was most accurate for reproducing the experimental binding mode. The ChemScore performance was not satisfactory for predicting either the binding mode or the binding affinity when compared to the other two scoring functions. Since the GoldScore yielded consistent binding modes and also an excellent correlation with the experimental binding affinity data, we selected the top-ranked binding pose of each ligand from GoldScore to study in detail the interaction of these derivatives with the active site residues. The predicted binding modes from GoldScore are given in Figure 4, whereas results from the other two scoring functions, ChemScore and the ASP, are given in Supplementary Figure S3.
We analyzed the top-ranked poses of the 10 cycloguanil ligands from GOLD (GoldScore) for important hydrogen bond and hydrophobic interactions (Figure 4). Two polar interactions are conserved in all the 10 cycloguanil derivatives: the 4-amino group forms a hydrogen bond with the side chain carboxylic acid of Asp54, and the 2-amino group forms a hydrogen bond with the backbone carbonyl of Ile14. Except for 31 and 39, which are unsubstituted at the C6 position, all the compounds showed an additional hydrogen bond with Asp54, between the carboxylic acid OH of Asp54 and N5 of the triazine ring of the ligand.
C6 of the triazine ring of cycloguanil derivatives lies within 4.5–5.0 Å of the methyl of Ala16 (Dist1 in Table 6). C4′ of the N1 phenyl of the triazine ring lies within 4.8–5.3 Å of the side chain methyl of Ser108 (Dist2 in Table 6), allowing favorable hydrophobic interactions. Both of these hydrophobic interactions are roughly equally important for these analogs as we can see from the correlations, R2 = 0.793 and 0.725 between the reported wild type binding affinity and Dist1 and Dist2, respectively. To make one interaction index which would account for both these interactions we used the average of the two distances = (Dist1+Dist2)/2, which showed an excellent correlation of R2 = 0.874 with the binding affinity data (Figure 5). We also calculated a key dihedral angle C2-N1-C1′-C2′ (atom labels given in Table 1) involved in the various substitution patterns in these analogs, using the docked poses, and found it to show a good correlation of R2 = 0.873 with binding affinity data (Figure 5). This clearly demonstrates that hydrophobic interactions of these analogs with the active site of PfDHFR-TS and the C2-N1-C1′-C2′ dihedral angle are crucial for their binding affinity in the wild type enzyme.
Since the X-ray structure is unavailable for the double mutant relevant to these cycloguanil compounds, A16V + S108T, we did not proceed with docking into the mutant protein. It would be possible to mutate the wild type protein or another mutant X-ray structure and then perform docking but this would be less reliable than working directly from a crystal structure, as we were able to do for the wild type.
In the modern scenario of the drug design program, in addition to designing potent and safer drug candidates against malaria, necessary care has to be taken in addressing the emergence of resistance to anti-malarial drugs. In order to understand the importance of various physico-chemical and structural properties of cycloguanil derivatives towards their wild and mutant type PfDHFR-TS binding and to identify and use an optimal docking program to study the molecular interactions, we reported herein our extensive 2D-QSAR investigations followed by comparative docking analysis. Our 2D-QSAR investigations suggested that several structural features of the chosen ligands and the hydrophobicity of N1-phenyl substituents of 2,4-diamino-1,6-dihydro-1,3,5-triazine scaffold are crucial for effective binding with the wild type of the enzyme. This observation was analogous to our previously published results on a different antimalarial enzyme target which have shown the importance of ligands’ hydrophobicity to be crucial for their in vitro antimalarial activities.46 For binding to the A16V+S108T mutant form of the same enzyme, we determined that several structural variables such as the presence of 3-phenoxyphenyl, 3-benzyloxyphenyl and 4-phenoxyphenyl at the C6 position of the triazine ring of cycloguanil derivatives influence binding, and also found a negative steric tolerance for para-substituents on the N1-phenyl ring.
We studied the binding mode of the cycloguanil derivatives with DHFR-TS by focusing on the 10 reported analogs which had well-defined stereochemistry. Among the four different docking programs and 13 different scoring functions analyzed, we obtained the best correlation R2 = 0.911 between reported binding affinity data and docking scores for the GOLD program using the GoldScore fitness function. Also, among the programs and scoring functions analyzed, we found the Molegro virtual docker program, with its template docking algorithm and MolDock [GRID] score, to be very accurate for reproducing the experimental binding mode of ligand WR99210, with an RMSD of 1.02 Å. The average of two key hydrophobic distances and a key dihedral angle were found to be well correlated to the wild type activity exhibited by the 10 cycloguanil derivatives. Hence a combination of 2D-QSAR and docking analysis provided fresh insights into the existing SAR of cycloguanil derivatives towards their PfDHFR-TS inhibition.
This research was supported in part by funding from the University of Mississippi, including from its Faculty Research Program and the Office of Research and Sponsored Programs; from the National Center for Zoonotic, Vector-borne, and Enteric Diseases (CK) of the Centers for Disease Control and Prevention (U01/CI000211); from the National Science Foundation (EPS-0556308); and from the National Institute of Health’s National Center for Research Resources (P20 RR021929). SP is a Natural Products Neuroscience Fellow. PNT was a Summer Research Institute for Undergraduates participant. This investigation was conducted in a facility constructed with support from research facilities improvement program C06 RR-14503-01 from the NIH National Center for Research Resources.
Supporting Information: A complete list of different docking scores from FlexX, Glide, Molegro and GOLD and predicted binding modes using FlexX, Glide and GOLD for all 10 compounds used in the docking analysis are given in the Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org.