Study of novel pyrazolo[3,4-d]pyrimidine derivatives as selective TgCDPK1 inhibitors: molecular docking, structure-based 3D-QSAR and molecular dynamics simulation

Abstract

Toxoplasma gondii calcium-dependent protein kinase 1 (TgCDPK1) is a promising drug target to treat toxoplasmosis, and the selective TgCDPK1 inhibition over human kinases is significant for the development of potent antiparasitic TgCDPK1 inhibitors. In the present study, a total set of 46 pyrazolopyrimidine-based compounds possessing TgCDPK1 and Src inhibitory activity were studied by using a molecular modeling approach combining molecular docking, three dimensional quantitative structure–activity relationship (3D-QSAR) and molecular dynamics (MD) simulations. The best comparative molecular field analysis (CoMFA) models were established with satisfactory robustness and predictability, with R² = 0.968, q² = 0.666 and R_pred² = 0.745 for TgCDPK1 and R² = 0.970, q² = 0.581 and R_pred² = 0.635 for Src, respectively. Other tests on additional validations further confirmed the satisfactory predictive power of the models. The key residues impacting the interactions and the probable binding modes between inhibitor and enzymes (TgCDPK1 and Src) were identified by docking and further verified by MD simulations. Computational results demonstrated that bulky or electronegative substituents on R₂, and certain bulky groups attached to the terminal of R₁ may increase the potency of TgCDPK1 and TgCDPK1/Src selectivity. Finally, six new compounds showing high TgCDPK1 potency and TgCDPK1/Src selectivity were designed. We hope this study can be helpful for further development of novel potent TgCDPK1 selective inhibitors.

---

Introduction

Toxoplasma gondii is a pathogenic microorganism which can infect humans and other warm-blooded animals.¹ In most instances, the parasite invades the human body through the intake of infected water, food or undercooked meat with tissue cysts. Generally, the invasion process can be divided into two steps, the rapid proliferation of tachyzoites and the spread to other organisms of bradyzoites.² As a systemic disease, toxoplasmosis is distributed all over the world and most people are recessive carriers. However, the ​pregnancy infection with T. gondill can be vertically transmitted, which can cause birth defects or miscarriage.^3,4 What’s more, in recent years, toxoplasmosis has been recognized as a fatal infection for HIV patients. Currently, a combination regimen of pyrimethamine associated with sulfadiazine has achieved great success for toxoplasmosis therapy.⁵ Due to the teratogenicity of pyrimethamine for pregnant women and the sensitivity to sulfa drugs for patients, some second-line therapies are usually used in clinical treatment, but the effectiveness is not satisfactory.^6,7 So it is necessary to find new Toxoplasma gondii antiparasitic drugs with efficiency and low toxicity.

Calcium dependent protein kinase (CDPK) is an enzyme target that is uniquely found in plants and protozoans, controlling calcium-regulated biological processes, such as host cell invasion, gliding motility and exocytosis.^8,9 Many reports have shown that the over-activation of CDPK1 may cause parasites to invade host cells and then lead to the proliferation of T. gondii, demonstrating that controlling the TgCDPK1 activity is very important to treat toxoplasmosis. Thus, TgCDPK1 becomes a promising target for the development of an antiparasitic agent and several ATP-competitive inhibitors of TgCDPK1 have been developed. Recently, a pyrazolopyrimidine (PP) scaffold has been proved to be efficacious to the TgCDPK1 enzyme, which contains a glycine gatekeeper residue permitting large hydrophobic substituents for shapes complementary to an adjacent hydrophobic pocket.^10,11 At the same time, these pyrazolopyrimidine-based inhibitors also exhibit potential inhibition against related Src, which contains the smallest threonine gatekeeper residue in human kinases. Thus the Src kinase can act as a filter for probing minimally with potential off-target human kinases in the search for selective TgCDPK1 inhibitors.¹²

Nowadays, the in vitro assessment of selective TgCDPK1 inhibitors is still a time-consuming and labor-intensive task. In order to predigest the drug discovery process, it is necessary to carry out detailed theoretical studies on the mechanisms of ligand–receptor interactions.¹³ Computer-aided drug design technology, including molecular docking, molecular dynamics (MD) simulations and 3D-QSAR (quantitative structure activity relationship), etc., is the most effective approach to analyze the detailed binding modes between inhibitors and receptors.^14–16

In the current study, a total of 46 pyrazolopyrimidine-based derivatives which had dual TgCDPK1 and Src inhibitory activity were used to apply the molecular modeling study. The 3D-QSAR models using the popular comparative molecular field analysis (CoMFA) were built to explore the useful structural information influencing inhibitors’ activities and the ability to predict the optimum models was checked with internal and external validation. Moreover, to elucidate the rational conformations and interaction mechanisms for these compounds and two receptors, molecular docking and molecular dynamics simulations were also carried out. Finally, the binding free energy calculation by the MM-PBSA method was employed to verify the contributions of key residues to ligand potency and TgCDPK1/Src selectivity. Based on the above studies, some new compounds with higher TgCDPK1 inhibitory activity and TgCDPK1/Src selectivity were designed and they are waiting for experimental verification.

Materials and methods

Date set

A total of 75 pyrazolopyrimidine-based compounds with their inhibitory activities against TgCDPK1 and/or Src were collected from literature.¹⁷ The molecules that did not have extract IC₅₀ values in both receptors were ignored and the remaining 46 compounds were chosen for the present study. The general structural formulae and template compounds are shown in Fig. 1. In an approximate ratio of 3 : 1, the collected derivatives were divided into a training set (34 compounds) and a test set (12 compounds) for 3D-QSAR studies (ESI Table S1†).¹⁸ The principle of selecting test set compounds was to guarantee that not only are their IC₅₀ values randomly but uniformly distributed in the range of those for the whole set, but also their structures are as diverse as possible for the dataset. In this study, the IC₅₀ values of these compounds were converted to pIC₅₀ values (−log IC₅₀), which were used as dependent variables in the 3D-QSAR analyses.

Fig. 1 General structural formula and numbering of pyrazolo[3,4-d]pyrimidine derivatives (a) and template molecule (b) (compound 23).

3D-molecular structures of all compounds were built in Sybyl 6.9 software and then the energy minimization was carried out using the Powell gradient algorithm and the Tripos force field with a convergence criterion of 0.001 kcal mol⁻¹ Å⁻¹ and a maximum of 1000 iterations. MMFF94 charges were assigned to each compound. The energy minimized structure was used for molecular docking.

Molecular docking

The Surflex module in the Sybyl 6.9 software package was applied to dock these inhibitors into the binding pockets of receptors.^19,20 The crystal structures of TgCDPK1 and Src were retrieved from the Protein Data Bank with ID 3SXF and 3UQF, respectively. Prior to docking, it is necessary to do the preparatory work. We analyzed the proteins and filled the missing atoms and residues using the Biopolymer module implemented in Sybyl. Then, all water molecules and ligands were removed, the polar H-atoms were added and Kollman charges were assigned to the protein.²¹ In the present study, the ligand-based mode was used to generate the Protomol. The two Protomol extent parameters, threshold and bloat, were set to 0.5 and 0, respectively. By default, 20 conformations per ligand were received from the docking method. The best conformation score for each ligand was selected for further studies.²²

3D-QSAR studies

Active conformation selection and molecular alignment are the key procedures in 3D-QSAR analysis.²³ In the present study, the docked conformation of the most active compound 23 was used as the template alignment. The common fragment (the atoms numbered from 1 to 9) of TgCDPK1 and Src was selected for alignment. The two aligned superimposition results are depicted in ESI Fig. S1.†

The CoMFA steric and electrostatic fields were calculated based on Lennard-Jones and coulombic potentials, with standard 30 kcal mol⁻¹ cutoffs. All the aligned compounds were calculated by using a sp³ carbon atom as the probe and a +1.00 unit charge at regular grid point spaces.^24,25

The partial least squares (PLS) method was adopted for the regression analysis to build the 3D-QSAR equations. The optimal number of components and the cross-validation correlation coefficient (q²) were determined by leave-one-out (LOO) cross-validation procedures.²⁶ Then, the non-cross-validated analysis was performed to calculate the non-cross-validation correlation coefficient (R²), standard error of estimates (SEE), and F value. The bootstrapping analysis for 100 runs was applied to further assess the robustness and statistical confidence of the derived models.²⁷

The predictive abilities of two CoMFA models were identified by a test set of 12 compounds, and their pIC₅₀ values were predicted. The predictive correlation coefficient (R_pred²) was obtained using the following eqn (1):

where SD is the sum of the squared deviation between the biological activities of the compounds in the test set and the average activity in the training set, and PRESS is the sum of squared deviation between actual and predicted activities of compounds in the test sets.

Model validation. In order to examine the true predictive abilities and robustness of models, further statistical external validations were applied.^28,29

The regression of y_i against ỹ_i (or ỹ_i against y_i) through the origin, y^ro_i = kỹ_i (or ỹ^ro_i = k′y_i), slope k and k′ was calculated as follows:

where ntest is the total number of compounds in the test set, ȳ_tr is the average activity in the training set, y_i and ỹ_i are the observed and predicted activities in the test set.

While y^ro_i (or ỹ^ro_i) can be determined by r_o² (or r′_o²), which was obtained by the following formula:

r_o² could be used to calculate essential parameter r_m² as follows:

where, r² was the non-cross-validated correlation coefficient obtained from the PLS process.

The 3D-QSAR models are acceptable when the following parameters are all satisfied: q² > 0.5, R_pred² > 0.6, r_o² or r′_o² close to r², i.e. [(r² − r_o²)/r²] < 0.1 or [(r² − r′_o²)/r²] < 0.1, with the corresponding 0.85 ≤ k ≤ 1.15 or 0.85 ≤ k′ ≤1.15, r_m² > 0.5.³⁰

In addition, the Y-randomization test was used to test the robustness and statistical significance for the 3D-QSAR models. In this test, the Y-variable was shuffled while the X-matrix was kept unaltered to guarantee that a good calibration result was not due to chance correlation.^31,32 The whole Y-randomization test was repeated over several trials and expected to have low r_{Y-randomization}² and q_{Y-randomization}² values, which confirm that the developed QSAR models are robust.³³

Domain of applicability. It is necessary to limit the domain of application before the generated theoretical model is used for screening new compounds and it can be considered reliable for only the compounds that fall into this domain.^31,34 The applicability domain (APD) of the 3D-QSAR model was defined by similarity measurements, which were determined on the basis of Euclidean distances between all pairs of training and test set compounds. The APD was calculated by using the following expression:

APD = 〈d〉 + Zσ (5)

where 〈d〉 and σ are the average and standard deviation of all distances, respectively. Z is an empirical cutoff with a default value of 0.5.³⁵ The calculation of 〈d〉 and σ was performed as follows: first, the initial average of all Euclidean distances between the training set compounds was calculated. Then, the set of distances that were lower than the average was formulated.

Molecular dynamics simulations

Based on the docking results, the MD simulations were performed with the AMBER 9.0 software package.³⁶ The docked complexes of 3SXF and 3UQF with the best active compound 23 and the lowly active compound 4 were picked as initial conformations. The B3LYP density functional method and 6-31G(d) basis set was carried out for electrostatic potential calculations, and atomic charges for ligands were fitted by calculating the RESP procedure.³⁷ The potential parameters of the protein and ligand were assigned using the standard AMBER 03 force field and the general amber force field (GAFF).^38,39 Then, all the complexes were neutralized by adding counter ions and immersed in a TIP3P water model box with a margin distance of 12 Å.⁴⁰

In order to eliminate unfavorable steric clashes, energy minimizations were performed in two steps. Firstly, all complexes were restrained with a force of 2.0 kcal mol⁻¹ Å⁻² and minimized by 2000 steps of steepest descent (SD) followed by the conjugated gradient method for 3000 steps. Secondly, all the complex systems were performed by the steepest descent minimization for 5000 steps and conjugated gradient for 5000 steps without any constraint. Then, the complex system was gradually heated from 0 to 300 K over 200 ps at a constant volume with a force constant 1.0 kcal mol⁻¹ Å⁻², and equilibrated at 300 K for another 500 ps with a constant pressure of 1 atm. Finally, the MD simulations were run over 10 ns in the NPT ensemble with the target temperature at 300 K and at 1.0 atm. During the MD simulations, the time step was defined as 2 fs, and the particle mesh Ewald (PME) method was performed to deal with the long-range coulombic interactions with a cutoff distance of 8.0 Å.⁴¹ The SHAKE algorithm was applied for all hydrogen-heavy atoms. The coordinate trajectories were written every 1 ps and the last stable MD trajectory was extracted to apply a further binding free energy calculation analysis.

Binding free energy calculations

A MM-PBSA methodology was used to calculate the binding free energies.^42,43 For each system, a total of 200 snapshots of each complex were extracted from the last 2 ns of stable MD trajectory and used for calculations. The binding free energy was computed by following eqn (6):

ΔG_bind = G_complex − (G_protein + G_ligand) = ΔG_MM + ΔG_sol − TΔS (6)

where ΔG_MM is the molecular mechanics free energy, ΔG_sol is the solvation free energy, and TΔS is the entropy contribution. The molecular mechanics free energy (ΔG_MM) was calculated by the van der Waals and electrostatic interaction:

ΔG_MM = ΔG_vdw + ΔG_ele (7)

The solvation free energy ΔG_sol is the sum of two terms:

ΔG_sol = ΔG_PB + ΔG_SA (8)

where the solvation free energy (ΔG_PB) is calculated by solving the Poisson–Boltzmann (PB) equation or by the generalized Born (GB) model for MM-PBSA. ΔG_SA was determined according to eqn (9):

ΔG_SA = γ × SASA + β (9)

where γ is the surface tension proportionality constant which was set to 0.0072 kcal mol⁻¹ Å⁻², and β represents the offset value, which was 0 here. SASA is the solvent accessible surface area. As the entropy calculation was time-consuming and had low predictive accuracy, the entropy term was not calculated.^44,45 In order to find the difference of binding modes and the key residues of these complexes, the MM/GBSA model was also applied to decompose the binding free energy into each residue. There are four terms we can get from decomposition results: van der Waals contribution (ΔG_vdw), electrostatic contribution (ΔG_ele), polar solvation contribution (ΔG_ele,sol) and non-polar solvation contribution (ΔG_nonpol,sol).

ΔG_{inhibitor-residue} = ΔG_vdw + ΔG_ele + ΔG_ele,sol + ΔG_nonpol,sol (10)

Here the van der Waals (ΔG_vdw) and electrostatic interactions (ΔG_ele) were calculated with the Sander program. The polar contribution ΔG_ele,sol was determined by the generalized Born (GB) model (GB^OBC, igb = 2). The non-polar contribution ΔG_nonpol,sol was computed using the solvent accessible surface area (SASA).^46,47

Results and discussion

Validation of docking reliability

It is indispensable to verify the reliability of the docking project before docking analysis. The X-ray structures of TgCDPK1 (PDB code: 3SXF) and Src (PDB code: 3UQF) in complexes with their native ligands were obtained from the Protein Data Bank. The two native ligands were redocked into their corresponding binding pocket, whose RMSD values were 0.597 Å and 0.697 Å, respectively. These results indicated the docking procedure was highly trustworthy, which could be extended to search the binding modes of other inhibitors.

Docking results

All studied inhibitors were docked into the binding sites of TgCDPK1 and Src kinases, and the optimal conformations of these compounds were identified. To further elucidate the interaction mechanism, we selected the most potent TgCDPK1 inhibitor 23 to perform deeper docking studies and for discussion.

The binding models of compound 23 with respect to TgCDPK1 and Src are displayed in Fig. 2 and 3. It is noted that compound 23 is situated at the ATP binding site and bound to TgCDPK1 and Src in similar conformations, with the naphthalene ring being almost perpendicular to the pyrazolo[3,4-d]pyrimidine core. The substituent R₁ on the naphthalene ring was oriented toward the gatekeeper residue Gly128 in TgCDPK1 (or Thr338 in Src), and positioned in a large hydrophobic pocket formed by Ala78/293, Val79/323, Val100/Ile336, Met112/314, Leu114/Ala311, Leu126/Ala337, Phe196/307 and Leu198/Ala403. At the same time, the R₂ substituent, isopropyl, surrounded by Leu57/273, Gly58/274, Gly60/344 and Val65/281, was extended to the solvent area. In addition, in the static structures, the amino group at the C₆-position of pyrazolo[3,4-d]pyrimidine and N₁ atom can form two hydrogen bonds with the carbonyl oxygen of Glu129/339 (C=O⋯NH, 3.1 Å/2.9 Å) and NH backbone of Tyr131/Met341 (N–H⋯N, 2.9 Å/3.1 Å) for TgCDPK1 or Src, respectively. The amino group at the C₆-position can make another weak hydrogen bond with Thr338 (C=O⋯NH, 3.4 Å) for Src, as the gatekeeper Gly128 in TgCDPK1 is smaller than Thr338, such that H–B cannot be found at this site. Finally, the naphthalene ring can form a π–cation interaction with the NH₃⁺ group of Lys80 and form a p–π interaction with the O atoms of Asp195 for TgCDPK1, but due to the far distances of Lys295 and Asp404, such interactions in Src–23 were greatly reduced (Fig. 3). Comparing TgCDPK1 with Src, the hydrogen bond and multiple other interactions have changed; these differences may be responsible for the distinct activities, which will be discussed in detail in the following section.

Fig. 2 Docking structure of the highly potent compound 23 and corresponding surface of TgCDPK1 (a) at the active binding site, in which the red and blue regions represent oxygen and nitrogen atoms, respectively, whereas white regions indicate carbon or hydrogen atoms. (b) Interactions between the TgCDPK1 active binding site and compound 23. Hydrogen bonds are depicted as red dotted lines.

Fig. 3 Docking structure of the highly potent compound 23 and corresponding surface of Src (a) at the active binding site, in which the red and blue regions represent oxygen and nitrogen atoms, respectively, whereas white regions indicate carbon or hydrogen atoms. (b) Interactions between the Src active binding site and compound 23. Hydrogen bonds are depicted as red dotted lines.

3D-QSAR analyses

The 3D-QSAR models for both TgCDPK1 and Src kinases were generated from CoMFA analyses and their statistical parameters are listed in Table 1. For a reliable predictive model, the cross-validated coefficient q² should be greater than 0.5.

Table 1 Statistical results of the CoMFA model

Statistical parameters	CoMFA
	TgCDPK1	Src
R^2	0.968	0.970
N	5	6
q^2	0.666	0.581
SEE	0.107	0.110
F	135.026	143.844
Rbs^2	0.983	0.987
SDbs	0.012	0.007
Rpred^2	0.745	0.635
k	1.003	0.998
r^2	0.894	0.868
ro^2	0.901	0.878
(r^2 − ro^2)/r^2	0.070	0.085
rm^2	0.712	0.676
rY-randomization^2	0.074–0.310	0.087–0.369
qY-randomization^2	−0.211 to 0.193	−0.146 to 0.076

---

In our study, two satisfactory CoMFA models with high R² = 0.983, q² = 0.666, F = 135.026, SEE = 0.107 for TgCDPK1 and R² = 0.970, q² = 0.581, F = 143.844, SEE = 0.110 for Src analyses were obtained, showing the models had optimal internal predictive capability. The R_bs² of 0.983 and 0.970, SD_bs of 0.012 and 0.007 for these two models from the bootstrapping analysis (100 runs) further revealed the robustness and effectiveness of the constructed models.

To evaluate the true predictive power of the established models, both CoMFA models were subjected to systemic external validation using 12 compounds in the test set, giving the predictive correlation coefficient R_pred² of 0.745 for TgCDPK1 and 0.635 for Src and other parameters r² of 0.894, r_o² of 0.901, [(r² − r_o²)/r²] of 0.070, k of 1.003 and r_m² of 0.712 for TgCDPK1, and r² of 0.868, r_o² of 0.878, [(r² − r_o²)/r²] of 0.085, k of 0.998 and r_m² of 0.676 for Src. These parameters agree well with the above statistical criteria, further certifying that the established models are satisfactory.

The APD values of the test set molecules (1.210–1.855) were also observed to be within the limit range of the calculated APD (APD = 2.455) of the training set molecules, which indicates that the established models were reliable for the prediction of new compounds (ESI Table S1†).

In addition, the 10 random shuffles for the Y-randomization test gave r_{Y-randomization}² and q_{Y-randomization}² values in ranges of 0.074 to 0.310 and −0.211 to 0.193 for the TgCDPK1 model, as well as 0.087 to 0.369 and −0.146 to 0.076 for the Src model; the low r_{Y-randomization}² and q_{Y-randomization}² values mean the established models were not correlated.

The higher contributions of the steric field (70.8% and 51.2%) compared to the electrostatic field (29.2% and 48.8%) demonstrate that the steric field is more dominant than electrostatic field to the inhibitory activity for two models.

The calculated pIC₅₀ and the residual values of compounds for two CoMFA models are listed in Table 2. The correlations between the experimental pIC₅₀ and the calculated ones for two models are displayed in Fig. 4. We can judge that the models have good quality because most points are distributed well along the line Y = X.

Table 2 The actual pIC₅₀, predicted pIC₅₀ and the residuals of the studied compounds for CoMFA analyses

No.	TgCDPK1			Src
	Act. pIC50	Pred.	Res.	Act. pIC50	Pred.	Res.
Training set
1	8.347	8.324	0.023	6.886	6.951	−0.065
4	7.398	7.397	0.001	5.292	5.456	−0.164
5	8.301	8.275	0.026	6.187	6.368	−0.181
6	8.387	8.353	0.034	6.745	6.607	0.138
7	7.921	7.847	0.074	5.824	5.834	−0.010
8	8.041	8.022	0.019	5.921	5.984	−0.063
9	8.444	8.360	0.084	6.538	6.511	0.027
11	8.602	8.767	−0.165	5.678	5.563	0.115
13	7.509	7.502	0.007	5.658	5.642	0.016
15	8.102	8.222	−0.120	6.921	6.696	0.225
16	7.620	7.634	−0.014	6.699	6.774	−0.075
18	8.222	8.299	−0.077	6.161	6.186	−0.025
19	7.699	7.686	0.013	5.699	5.767	−0.068
20	8.301	8.357	−0.056	6.420	6.568	−0.148
21	7.482	7.298	0.184	5.658	5.526	0.132
22	8.000	8.102	−0.102	6.260	6.176	0.084
23	9.222	8.915	0.307	6.699	6.541	0.158
24	8.495	8.776	−0.281	6.523	6.458	0.065
25	8.268	8.240	0.028	6.585	6.648	−0.063
26	6.699	6.798	−0.099	5.959	6.045	−0.086
29	8.638	8.650	−0.012	5.509	5.530	−0.021
30	7.959	7.903	0.056	6.553	6.494	0.059
31	8.886	8.847	0.039	5.745	5.560	0.185
32	8.260	8.248	0.012	6.553	6.594	−0.041
33	9.155	9.159	−0.004	6.420	6.471	−0.051
35	9.155	9.148	0.007	5.796	5.760	0.036
36	8.585	8.600	−0.015	5.301	5.324	−0.023
37	8.620	8.616	0.004	5.252	5.348	−0.096
38	8.469	8.491	−0.022	5.061	5.017	0.043
40	8.215	8.305	−0.090	5.796	5.863	−0.067
42	8.721	8.703	0.018	5.482	5.556	−0.075
43	8.658	8.562	0.096	5.284	5.199	0.085
45	8.523	8.460	0.063	5.174	5.155	0.019
46	8.367	8.388	−0.021	5.086	5.148	−0.062
Test set
2	7.721	7.872	−0.151	5.886	6.114	−0.228
3	8.137	8.300	−0.163	6.237	6.155	0.082
10	7.886	7.988	−0.102	6.284	5.923	0.361
12	6.886	7.227	−0.341	5.056	5.641	−0.586
14	8.301	7.655	0.646	7.187	6.767	0.420
17	8.509	8.227	0.282	5.387	5.914	−0.527
27	9.046	8.925	0.121	6.108	6.516	−0.408
28	9.097	8.962	0.135	6.387	6.211	0.176
34	8.398	8.824	−0.426	6.569	6.486	0.083
39	8.367	8.353	0.014	6.137	5.859	0.278
41	8.553	8.902	−0.349	6.092	5.495	0.597
44	8.569	8.496	0.073	5.201	5.137	0.064

---

Fig. 4 The plot of the predicted versus actual pIC₅₀ values for the CoMFA model: (a) for TgCDPK1, (b) for Src, in which the black triangles means the training set compounds and the red dots expressed the test set compounds.

CoMFA contour maps analysis

To visualize the results of the CoMFA model, 3D color-coded contour maps were generated. Fig. 5a and b show the steric and electrostatic contour maps for the TgCDPK1 model with compound 23 as a reference structure.

Fig. 5 CoMFA contour maps based on the most active compound 23: (a) steric contour map of TgCDPK1 kinase inhibition, (b) electrostatic contour map of TgCDPK1 kinase inhibition, (c) steric contour map of Src kinase inhibition and (d) electrostatic contour map of Src kinase inhibition.

In Fig. 5a, one big green and three small yellow contours surrounding the terminal propenyl of the R₁ substituent, suggest that a moderately-sized group in this site is favorable, because a group that is too large may make collisions with Val100, Leu198 and Leu126. This is illustrated by experimental fact that compounds 23 and 24 with –CH₂CHCH₂ and –CH₂CH₂CH₃ as the terminal group of R₁, respectively, exhibit higher potency than compounds 20 and 25 with –CH₂CH₃ and –CH₂C(CH₃)₃ at the same position. In addition, compound 28, which showed better activity than 30 and 34, can also be explained by this finding. The other two yellow regions located under and above the naphthalene ring show that no substituents in the naphthalene ring except for the 6′-position are favorable, which may be blocked by the gatekeeper Gly128. This may be the reason why compound 12, with 1′-methylnaphthalene as R₁, is almost the most inactive compound. A yellow contour is shown in the vicinity of amino, which would be blocked by the near residues Glu129 and Val130. Similarly, a yellow region is found near the –CH₃ group of substituent R₂, which is located at the entrance of the solvent area, and in docking (Fig. 2b) with one methyl group may be blocked by the near residue Leu57, while another –CH₃ has no steric hindrance, thus certain bulky R₂ would be beneficial for the activity. The compound activity order 14 (8.3010) > 15 (8.1024), 18 (8.2218) > 19 (7.6990), 20 (8.3010) > 21 (7.4815) and 25 (8.2676) > 26 (6.6990) were all corresponding examples.

On the TgCDPK1 electrostatic contour map (Fig. 5b), several red contours near the O atom, the terminal of R₁ and the amino group at C₆-position, illustrate that electronegative groups at these regions may improve the binding activity, which is in agreement with the experimental conclusion that compounds 31, 33, 35 with an electronegative –Cl group at the terminal of R₁ reveal better biological activity than 29 with a H atom at this site. Two blue maps above naphthalene ring mean a positively charged group is favorable, which may result in electrostatic interactions with the electron-rich O atoms of Glu129 when docking, hence, most of the excellent compounds 23, 27–29, 31, 33, 35, 42–43 all possess a naphthalene ring. Furthermore, there are a large blue contour inserted between two big red maps near the substituent R₂, revealing that certain long alkyl groups in this blue region would enhance the activity, which can form many hydrophobic contacts with the residues Leu57, Gly58 and Gly60 at the entrance of the solvent region in docking. In addition, the terminal of R₂ bearing a relatively negative charged group would be advantageous to increase the activity. The docking study also shows that the substituent R₂ is extended outside to the solvent area, so the R₂ with a negative atom should have some good physicochemical properties, such as solubility, permeability, etc., to increase the binding affinity to the plasma protein. That may be the reason why the piperidine-containing compounds all confer potent inhibition of TgCDPK1 activity.

The CoMFA maps for Src are rather similar, while only slight differences were found compared to TgCDPK1 (Fig. 5c and d). The steric field shows two unfavorable contours near substituents R₂, while only one yellow contour near this position in TgCDPK1. It suggests the requirement of a small R₂ group to increase the activity. Most of the excellent inhibitors 1, 6, 14–16, 23, 25, 28, 30, 32 and 34 all bear a relatively small i-Pr or t-Bu group as R₂, and compounds 36–38, 43–46 possessing bulky substituent R₂ are the most inactive compounds. Moreover, a steric favorable green place is situated near the 3′- and 4′-positions of the naphthalene ring in Src, replacing the yellow contour at this site in TgCDPK1, which is in accordance with the steric results by comparing compound activity; 5 (-m-CH₃) > 4 (–H). In the docking experiment, we can find that gatekeeper Thr338 is close to the second phenyl of the naphthalene in Src, allowing certain group linking to the 3′-, and 4′-positions, while Gly128 is near the first phenyl in TgCDPK1. Meanwhile, the size of the green contour near the propenyl of R₁ is smaller than that of TgCDPK1, implying a certain bulky group in this site is more favorable for TgCDPK1 than Src. As for the electrostatic field, two blue plots above and below the C₆-position call for a positively charged group or atom to improve activity in Src, instead of the small red contour near the amino group in TgCDPK1. In docking, the pyrimidine ring may form π–π stacking with the near residue Tyr340, and the positive charge of C₆ can produce an electrostatic interaction with the electron-rich O atom of Tyr340. Fig. 5d shows that the blue regions near the naphthalene ring have disappeared; the red plots are still there though their size is enlarged. This indicates the importance of a negative substituent for Src over TgCDPK1, which further explains why compounds 31, 33, 35 are more active than compound 29. Meanwhile, the blue contour at the terminal of R₂ in TgCDPK1 is also missing, only a small red plot is seen here. Combining the above steric analysis, a small group with electronegative atoms at the terminal of R₂ should be good for Src.

MD simulations

MD simulations were carried out for 4 complexes of two series (23-3SXF, 4-3SXF, 23-3UQF and 4-3UQF) to validate the dynamic interactions between ligands and receptors. The RMSD values of the backbone atoms of the proteins with respect to their initial conformations are analyzed and displayed in Fig. 6.

Fig. 6 RMSD of backbone atoms of the studied complexes during MD simulations.

As shown in Fig. 6, the four complexes reached equilibrium at about 0.5 ns and the mean RMSD values of 23-3SXF and 4-3SXF complexes are significantly lower than those of 23-3UQF and 4-3UQF complexes. It indicates that the complexes are stable during the MD progress and the compounds bind more tightly to TgCDPK1 than to Src, which coincides with the much higher activities of 23 and 4 against TgCDPK1 than Src. Meanwhile, the superimposition of the carmine stick average structure from the last 2 ns trajectory and the yellow stick initial docked structure is displayed in ESI Fig. S2.† The obvious appearance of this picture is that the docked structures and the average structures are well overlapped at the active sites with only little declination, which verifies the reliable docking results.

Moreover, the analyses of root-mean-square fluctuation (RMSF) vs. the residue number for complexes are plotted in Fig. 7. The protein structures of the two complexes in each profile display similar distributions and trends of dynamics properties (Fig. 7a). It is noted that the 4-bound systems exhibit higher flexibility than 23-bound systems, especially the regions around Leu57, Lys80, Met112, Leu114, Leu126 and Gly128 for TgCDPK1 and Leu273, Val281, Ile294 and Met341 for Src. On this basis, we can deduce that 23 might have a more stable interaction with receptors than 4. The important hydrogen bond interactions (Table 3) also show that both compounds 4 and 23 can form two stable hydrogen bonds with Tyr131/Met 341 and Glu129/339 for TgCDPK1 or Src, respectively, with occupancies of over 85% during the MD simulations. Whereas, the supernumerary hydrogen bond interaction with Thr338 in the 23–Src docking result disappeared.

Fig. 7 RMSFs of each residue of the protein for the four systems. (a) Complexes 23-3SXF and 4-3SXF and (b) complexes 23-3UQF and 4-3UQF.

Table 3 H-bond analysis from MD

System	Donor	Acceptor	Occupancy (%)	Distance (Å)	Angle (°)
4-3SXF	Ligand@N1	Tyr131@N–H	94.44	3.179	18.64
	Glu129@O	Ligand@N–H	90.79	2.886	28.98
	Ligand@N3	WAT@O–H	18.40	3.103	23.71
23-3SXF	Glu129@O	Ligand@N–H	99.69	2.924	25.59
	Ligand@N1	Tyr131@N–H	95.05	3.044	29.70
	Ligand@N3	WAT@O–H	38.49	3.031	25.10
	Ligand@N8	WAT@O–H	23.30	3.201	34.15
4-3UQF	Glu339@O	Ligand@N–H	99.40	2.947	24.05
	Ligand@N1	Met341@N–H	97.60	3.103	19.40
23-3UQF	Ligand@N1	Met341@N–H	96.35	3.123	19.28
	Glu339@O	Ligand@N–H	85.26	3.128	34.18
	Ligand@N8	WAT@O–H	35.20	3.047	21.90

---

To explore the influence of water on the binding of inhibitors in the active sites of receptors, the stereoview of the proteins with 23 and water molecules is displayed in Fig. S3 (ESI†). Fig. S3a† shows that two water molecules can form weak hydrogen bonds with N₃ and N₈ atoms of compound 23 in 23-3SXF, with occupancies of 38.49% and 23.30% during the MD simulations, respectively, implying the water molecules may be favorable for the binding of 23 to TgCDPK1. From Table 3, we can also find that, in 4-3SXF, the N₃ atom of compound 4 forms a hydrogen bond with one water molecule, which is weaker than those in 23-3SXF, meaning water molecules can provide a greater effect for the binding force of 23 over 4 in TgCDPK1, although they are all weaker interactions. Meanwhile, Fig. S3b† shows that one water molecule can form a hydrogen bond with a N₈ atom of 23 in the 23-3UQF complex, whereas no hydrogen bond can be found for 4 with water molecules in the 4-3UQF complex. Thus, the water molecules would help to stabilizing the inhibitor binding to two receptors, which will be in good correlation with the activity.

Binding free energy analysis results

The MM/PBSA program was used to calculate the binding free energies of the four systems and the results are shown in Table 4. As can be seen, the calculated ΔG_bind values of the TgCDPK1 bound systems (23-3SXF and 4-3SXF) remain higher than those of the Src bound systems (23-3UQF and 4-3UQF), meaning that the binding affinities of the inhibitors to TgCDPK1 are stronger than those to Src. This is in agreement with the experimental activities. We can also ascertain from Table 4, that ΔG_vdw, ΔG_ele and ΔG_nonpol,sol terms are favorable for binding and that the positive ΔG_ele,sol displays adverse effects for all systems. Comparing 23- with 4-bound systems, it is noted that the sum of the non-polar term (ΔG_vdw + ΔG_nonpol,sol) of 23 (−61.15 kcal mol⁻¹ for TgCDPK1 and −54.35 kcal mol⁻¹ for Src) is obviously stronger than that of 4 (−52.27 kcal mol⁻¹ for TgCDPK1 and −43.58 kcal mol⁻¹ for Src). The polar contributions (ΔG_ele + ΔG_ele,sol) in two TgCDPK1 systems are similar (20.02 kcal mol⁻¹ for 23 and 20.79 kcal mol⁻¹ for 4), while they are different in two Src systems (23.23 kcal mol⁻¹ for 23 and 16.40 kcal mol⁻¹ for 4). Therefore, we can conclude that the non-polar interactions play a crucial role in the stronger binding affinities of 23 over 4 in both series and the polar interaction, especially ΔG_ele, also influences the binding to Src.

Table 4 The binding free energy of eight systems^a

System	Polar contributions		Non-polar contributions		ΔGbind	ΔGexp	Ki (μM)
	ΔGele	ΔGele,sol	ΔGvdw	ΔGnonpol,sol
a All energies are in kcal mol^−1. ΔGexp is the experimental binding free energy estimated from Ki values by ΔG ≈ −RT ln Ki.
23-3SXF	−20.09	40.11	−54.65	−6.50	−41.13	−14.04	0.00006
4-3SXF	−19.45	40.24	−46.98	−5.29	−31.48	−10.16	0.04
23-3UQF	−14.21	37.44	−47.68	−6.67	−31.12	−9.20	0.2
4-3UQF	−16.45	32.85	−38.14	−5.44	−27.18	−7.27	5.1

---

Furthermore, the binding free energy was decomposed to individual residues to gain more detailed interactions between residues and the inhibitor. Fig. 8 displays the results of the decomposition energy for each key residue of the four systems. As we can see in Fig. 8a, TgCDPK1 series residues Leu57, Val65, Ala78, Lys80, Met112, Glu129, Val130, Tyr131, Leu181 and Ile194 make the most favorable contributions to the binding of 23 (or 4) in TgCDPK1. From Fig. 8b, it can be also seen that the interactions between Src and 23 (or 4), mainly depend on residues Leu273, Val281, Lys295, Ile336, Thr338, Glu339, Tyr340, Met341 and Leu393. It is found that most of the important residues are hydrophobic, implying that van der Waals interactions play a key role for binding to two receptors. Meanwhile, residues Glu129/339 and Tyr131/Met341 can make hydrogen bonds with two compounds, suggesting that the hydrogen bond interactions are also indispensable. On the other hand, it is noted the residues with more favorable contributions to the better activities of 23 over 4 in TgCDPK1 are Lys80, Met112, Leu114 and Leu126. Comparing compound 23 with 4, substituent R₁ is the difference; the 4-propenyl-naphthalene group in 23 is replaced by 4-methoxybenzene in 4, whose binding conformations were almost the same in TgCDPK1 (Fig. 9a). Fig. 9a shows that the naphthalene ring of 23 can make a stronger π–cation interaction with Lys80 than the phenyl of 4. Meanwhile, the terminal propenyl group of R₁ in 23 is sandwiched between Leu114 and Leu126 and makes a number of hydrophobic contacts, but such interactions do not exist in compound 4. Through the above fact, the interactions of 23 with TgCDPK1 are stronger than those of 4, which is consistent with the higher activities of 23 over 4. The structure alignment of 23-3SXF and 23-3UQF systems also shows that the binding modes of 23 in two receptors are very similar, except for a slight offset, resulting in the spatial orientation of the terminal 3-propenyl (Fig. 9b). From Fig. 9b, we find that Glu129 binds closer to the amino group and makes a stronger hydrogen bond interaction with the amino group than Glu339, which is in line with the associated N–H⋯N distance increase from 2.924 Å in the 23-3SXF complex to 3.128 Å in 23-3UQF. Meanwhile, the positive part of the naphthalene ring can take electrostatic interaction with the electron-rich S atom of Met112 in TgCDPK1, but such interactions are not present in Src. Finally, the hydrophobic contacts of terminal 3-propenyl with near Leu126 and Leu114 in TgCDPK1 are also stronger than those with Ile336 in Src. The strengthening of the hydrogen bond, electrostatic and hydrophobic interactions are responsible for the higher activity of 23 against TgCDPK1 than Src.

Fig. 8 Free energy decomposition plots for the four systems. (a) Complexes 23-3SXF and 4-3SXF, (b) complexes 23-3UQF and 4-3UQF.

Fig. 9 (a) Superimposition of compound 23 (magenta) and 4 (green) in TgCDPK1. (b) Superimposition of 3SXF-23 (magentas) and 3UQF-23 (green).

The results confirm that the hydrophobic interactions play an important role in distinguishing the binding affinity of these inhibitors and also reveal the reliability of docking models.

Selectivities for TgCDPK1 over Src and design of new selective TgCDPK1 inhibitors

From the above analyses of 3D-QSAR and docking studies, some significant structural information for the selectivity is pointed out. Firstly, the size of substituent R₂ provides a primary determinant for the selectivity of TgCDPK1 over Src. It indicates that a bulky R₂ substituent, such as 4-piperidinemethyl or other bigger 4-piperidinemethyl-containing groups, may result in the greater selectivity for TgCDPK1 over Src. It can be easily revealed by the result that compounds 29, 36–46 exhibit much higher activities for TgCDPK1 inhibition than Src. Secondly, certain bulky groups introduced to the terminal of R₁ can improve TgCDPK1 selectivity, which is well demonstrated by the fact that compounds 31 and 35 with chlorobenzene groups attached to the terminal of R₁, had higher potency and selectivity for TgCDPK1 over Src. Furthermore, enhancing van der Waals interactions and avoiding the steric collision with gatekeeper Gly128 also benefit the activity.

Based on the above conclusions, six new compounds were designed using the most active TgCDPK1 compound 23 and the most selective compound 38 of TgCDPK1 over Src as a reference. Table S2† shows their structures and corresponding predicted activity, finding that the designed six new compounds exhibit higher TgCDPK1 activities and TgCDPK1/Src selectivity than the reference compounds. The calculated APD values of the new designed compounds were also within the limit range of the calculated APD (2.455) (ESI Table S2†). These results further validate the correctness of our models.

Conclusions

In the present work, a multiplex computational approach by molecular docking, 3D-QSAR and MD simulations was used to characterize the interactions and selectivity of a series of pyrazolopyrimidine-based inhibitors to TgCDPK1 and Src. Molecular docking produced possible poses of inhibitors for the TgCDPK1 and Src receptor, respectively. The graphical analyses of CoMFA models afford the critical structural features influencing the inhibitory activities and the TgCDPK1/Src selectivity. The docking study showed that the hydrogen bond interaction of residue Glu129 of TgCDPK1 was the crucial feature for TgCDPK1/Src selectivity. In addition, the van der Waals interactions, especially the important contacts with Leu114 and Leu126, play a key role on the distinct activities of inhibitors. The MD simulation and MM-PBSA calculations confirmed the reasonable binding modes of these complexes and the key interaction features. Based on the computational results, six compounds with higher TgCDPK1 potency and TgCDPK1/Src selectivity were designed. These results may provide some insights into the mechanisms between inhibitors and receptors and valuable clues to the further design of selective TgCDPK1 inhibitors.

Acknowledgements

The authors gratefully acknowledge support of this research by the Science and Technology planning Project of Guangzhou (No. 2013J4100071). We also heartily thank the computation environment support by the Department of Biochemistry, College of Life Sciences, SunYat-Sen University for SYBYL 6.9.

Notes and references

1. L. M. Weiss and J. P. Dubey, Int. J. Parasitol., 2009, 39, 895–901.
2. S. Lourido, C. Zhang, M. S. Lopez, K. Tang, J. Barks, Q. Wang, S. A. Wildman, K. M. Shokat and L. D. Sibley, J. Med. Chem., 2013, 56, 3068–3077.
3. G. Pappas, N. Roussos and M. E. Falagas, Int. J. Parasitol., 2009, 39, 1385–1394.
4. J. Dubey, E. Lago, S. Gennari, C. Su and J. Jones, Parasitology, 2012, 139, 1375–1424.
5. L. H. Bosch-Driessen, F. D. Verbraak, M. S. Suttorp-Schulten, R. L. van Ruyven, A. M. Klok, C. B. Hoyng and A. Rothova, Am. J. Ophthalmol., 2002, 134, 34–40.
6. E. J. Goldstein, J. G. Montoya and J. S. Remington, Clin. Infect. Dis., 2008, 47, 554–566.
7. M. Antczak, K. Dzitko and H. Długońska, Biomed. Pharmacother., 2016, 82, 677–684.
8. O. Billker, S. Lourido and L. D. Sibley, Cell Host Microbe, 2009, 5, 612–622.
9. S. Lourido, J. Shuman, C. Zhang, K. M. Shokat, R. Hui and L. D. Sibley, Cah. Rev. The., 2010, 465, 359–362.
10. K. K. Ojo, E. T. Larson, K. R. Keyloun, L. J. Castaneda, A. E. DeRocher, K. K. Inampudi, J. E. Kim, T. L. Arakaki, R. C. Murphy and L. Zhang, Nat. Struct. Mol. Biol., 2010, 17, 602–607.
11. S. Lourido, G. R. Jeschke, B. E. Turk and L. D. Sibley, ACS Chem. Biol., 2013, 8, 1155–1162.
12. E. T. Larson, K. K. Ojo, R. C. Murphy, S. M. Johnson, Z. Zhang, J. E. Kim, D. J. Leibly, A. M. Fox, M. C. Reid and E. J. Dale, J. Med. Chem., 2012, 55, 2803–2810.
13. W. L. Jorgensen, Science, 2004, 303, 1813–1818.
14. S. Ma, S. Tan, D. Fang, R. Zhang, S. Zhou, W. Wu and K. Zheng, RSC Adv., 2015, 5, 81523–81532.
15. N. Singh, S. Tiwari, K. K. Srivastava and M. I. Siddiqi, J. Chem. Inf. Model., 2015, 55, 1120–1129.
16. D. Zhou, J. Chen and Y. Xu, RSC Adv., 2016, 6, 51716–51724.
17. S. M. Johnson, R. C. Murphy, J. A. Geiger, A. E. DeRocher, Z. Zhang, K. K. Ojo, E. T. Larson, B. G. K. Perera, E. J. Dale and P. He, J. Med. Chem., 2012, 55, 2416–2426.
18. S. K. Verma and S. Thareja, RSC Adv., 2016, 6, 33857–33867.
19. M. Hao, Y. Li, Y. Wang, Y. Yan and S. Zhang, J. Chem. Inf. Model., 2011, 51, 2560–2572.
20. E. Cichero, C. Brullo, O. Bruno and P. Fossa, RSC Adv., 2016, 6, 61088–61108.
21. S. Ma, G. Zeng, D. Fang, J. Wang, W. Wu, W. Xie, S. Tan and K. Zheng, Mol. BioSyst., 2015, 11, 394–406.
22. F. Cao, X. Li, L. Ye, Y. Xie, X. Wang, W. Shi, X. Qian, Y. Zhu and H. Yu, Environ. Toxicol. Pharmacol., 2013, 36, 626–635.
23. J. Verma, V. M. Khedkar and E. C. Coutinho, Curr. Top. Med. Chem., 2010, 10, 95–115.
24. R. Safavi-Sohi and J. B. Ghasemi, Med. Chem. Res., 2013, 22, 1587–1596.
25. D. Han, M. Su, J. Tan, C. Li, X. Zhang and C. Wang, RSC Adv., 2016, 6, 27594–27606.
26. P. Kamsri, A. Punkvang, S. Hannongbua, P. Saparpakorn and P. Pungpo, RSC Adv., 2015, 5, 52926–52937.
27. V. Srivastava, S. Gupta, M. Siddiqi and B. Mishra, Eur. J. Med. Chem., 2010, 45, 1560–1571.
28. G. Alexander and T. Alexander, J. Mol. Graphics Modell., 2002, 20, 269–276.
29. E. Vrontaki, G. Melagraki, T. Mavromoustakos and A. Afantitis, Methods, 2015, 71, 4–13.
30. V. D. Mouchlis, G. Melagraki, T. Mavromoustakos, G. Kollias and A. Afantitis, J. Chem. Inf. Model., 2012, 52, 711–723.
31. S. Zhang, A. Golbraikh, S. Oloff, H. Kohn and A. Tropsha, J. Chem. Inf. Model., 2006, 46, 1984–1995.
32. G. Melagraki and A. Afantitis, RSC Adv., 2014, 4, 50713–50725.
33. P. S. Ambure, R. P. Gangwal and A. T. Sangamwar, Mol. Diversity, 2012, 16, 377–388.
34. G. Melagraki and A. Afantitis, Curr. Top. Med. Chem., 2014, 15, 1827–1836.
35. H. Jasuja, N. Chadha, M. Kaur and O. Silakari, Mol. Diversity, 2014, 18, 253–267.
36. R. Salomon-Ferrer, D. A. Case and R. C. Walker, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2013, 3, 198–210.
37. Y. L. Tian, M. Lv, J. J. Li, T. Xu, H. L. Zhai and X. Y. Zhang, Eur. J. Pharm. Sci., 2015, 76, 138–148.
38. K. Wichapong, A. Rohe, C. Platzer, I. Slynko, F. Erdmann, M. Schmidt and W. Sippl, J. Chem. Inf. Model., 2014, 54, 881–893.
39. H. Sun, Y. Li, M. Shen, S. Tian, L. Xu, P. Pan, Y. Guan and T. Hou, Phys. Chem. Chem. Phys., 2014, 16, 22035–22045.
40. J. P. do Nascimento, J. R. A. Silva, J. Lameira and C. N. Alves, Chem. Phys. Lett., 2013, 583, 175–179.
41. J. Shraberg, S. W. Rick, N. Rannulu and R. B. Cole, Phys. Chem. Chem. Phys., 2015, 17, 12247–12258.
42. T. Hou, J. Wang, Y. Li and W. Wang, J. Chem. Inf. Model., 2010, 51, 69–82.
43. G. Rastelli, A. D. Rio, G. Degliesposti and M. Sgobba, J. Comput. Chem., 2010, 31, 797–810.
44. U. Bren, V. Martínek and J. Florián, J. Phys. Chem. B, 2006, 110, 12782–12788.
45. X.-L. Shen, M. Takimoto-Kamimura, J. Wei and Q.-Z. Gao, J. Mol. Model., 2012, 18, 203–212.
46. C. Gao, M. Grøtli and L. A. Eriksson, J. Mol. Model., 2015, 21, 1–12.
47. M. Bello, R. Campos-Rodriguez, S. Rojas-Hernandez, A. C.-M. de Oca and J. Correa-Basurto, Immunol. Res., 2015, 62, 3–15.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra20277b

---