Pharmacophore and 3D-QSAR characterization of 6-arylquinazolin-4-amines as Cdc2-like kinase 4 (Clk4) and dual specificity tyrosine-phosphorylation-regulated kinase 1A (Dyrk1A) inhibitors.

Cdc2-like kinase 4 (Clk4) and dual specificity tyrosine-phosphorylation-regulated kinase 1A (Dyrk1A) are protein kinases that are promising targets for treatment of diseases caused by abnormal gene splicing. 6-Arylquinazolin-4-amines have been recently identified as potent Clk4 and Dyrk1A inhibitors. In order to understand the structure-activity correlation of these analogs, we have applied ligand-based pharmacophore and 3D-QSAR modeling combined with structure-based homology modeling and docking. The high R(2) and Q(2) (0.88 and 0.79 for Clk4, 0.85 and 0.82 for Dyrk1A, respectively) based on validation with training and test set compounds suggested that the generated 3D-QSAR models are reliable in predicting novel ligand activities against Clk4 and Dyrk1A. The binding mode identified through docking ligands to the ATP binding domain of Clk4 was consistent with the structural properties and energy field contour maps characterized by pharmacophore and 3D-QSAR models and gave valuable insights into the structure-activity profile of 6-arylquinazolin-4-amine analogs. The obtained 3D-QSAR and pharmacophore models in combination with the binding mode between inhibitor and residues of Clk4 will be helpful for future lead compound identification and optimization to design potent and selective Clk4 and Dyrk1A inhibitors.

Introduction

Cdc2-like kinases (Clk) and dual specificity tyrosine-phosphorylation-regulated kinases (Dyrk) both are CMGC family of protein kinases.[1,2] They are responsible for phosphorylation of serine-arginine-rich (SR) proteins and are important for regulation of basic cellular processes.[1,3,4] Specifically, the cdc2-like kinases promote phosphorylation within spliceosome, therefore regulating alternative splicing of mRNA isoforms.[5] Because abnormal gene splicing is the cause of many pathological conditions including cancers,[6,7] modulation of Clk may represent a promising approach for treatment of such diseases. Dyrk1A is the most ubiquitously expressed isoform of Dyrk family.[1] Located on the Down Syndrown critical (DS) region of chromosome 21, it has increased expression in DS patients,[8−10] and has shown involvement in growth and mental retardation and neurodegeneration.[1,9,11] Therefore, inhibition of Dyrk1A may be a strategy for development of drug candidates for these disorders. Some compounds have been identified as both Clk and Dyrk inhibitors, such as, 6-arylquinazolin-4-amine analogs,[5,12,13] leucettines,[14] bauerine C derivatives,[15] a benzothiazole analog,[16] and natural product extracts.[2] However, development of potent and selective Clk and Dyrk inhibitors is still yet to be explored.[12,13]

Pharmacophore and QSAR are ligand-based molecular modeling techniques based on the notion that compounds interacting with the same target could share similar structural or physicochemical properties. Structural properties such as hydrophobic, aromatic, and hydrogen-bond donor and acceptor could be featured by a pharmacophore model, which is used for characterization of structurally diversed compounds targeting the same protein.[17−22] In combination of virtual screening, pharmacophore modeling has been proved as an effective strategy for lead compound identification.[20,23] Compared to pharmacophore modeling, 3D-QSAR is also based on 3D-conformers but considers the overall force field around a molecule, instead of focusing on group features in a single region.[24−28] Typical programs that generate 3D-QSAR models include comparative molecular field analysis (CoMFA),[26] comparative molecular similarity indices analysis (CoMSIA),[29] and phase.[30,31] The force fields calculated by 3D-QSAR may be steric, electrostatic, hydrophobic, and hydrogen-bond donor and acceptor.[24] Because 3D-QSAR is best used when ligands share the same structural scaffold, it can be applied in lead optimization for rational drug design.[32,33]

The ligand-based pharmacophore and 3D-QSAR models may shed light on the design of novel Clk and Dyrk inhibitors and may help with issue of selectivity among Clk and Dyrk members. Previous publications have not identified pharmacophore or 3D-QSAR models for Clk and Dyrk ligands. Recently a series of 6-arylquinazolin-4-amines were reported as Clk and Dyrk inhibitors.[5,12,13] In the present study, we developed pharmacophore and 3D-QSAR models based on their activities against Clk4 and Dyrk1A by using the phase package of Schrodinger.[34] The obtained 3D-QSAR models have shown good predictive capabilities, according to the statistical validation based on training and test set compounds. Further, the binding mode between active ligands and the target Clk4 and Dyrk1A have been proposed based on docking program Glide.[35] The obtained ligand–protein interactions agree with the force field contours obtained via QSAR analysis and help to understand the protein–ligand interactions that are responsible for the biological activities on the molecular level when targeting Clk4 and Dyrk1A. The developed models give useful information of lead optimization for future rational design of Clk4 and Dyrk1A inhibitors and will be helpful in development of selective inhibitors between these two targets.

Methods

Pharmacophore Modeling

The pharmacophore and 3D-QSAR models were generated based on 56 recently published 6-arylquinazolin-4-amine analogs that were tested for their inhibition effects against Clk4 and Dyrk1A.[5,12,13] Training and test set compounds were selected in such way that they covered a similar range of biological activities. Their structures are shown in Table 1. The molecular structures were sketched and built with Maestro.[36] The pharmacophore models were generated with the “Develop Pharmacophore Model” module of phase. Observed activities (IC50) were converted to form of negative logarithm (pIC50) before pharmacophore generation. Multiple conformers were generated for each molecule followed by energy minimization based on OPLS-2005 force field.[37] The conformational space was explored by ConfGen,[38] with 100 conformers per rotatable bond and 1000 maximum of conformers per structure. A distance-dependent dielectric was applied for solvation treatment. The pharmacophore models were developed with the most active training set compounds, which are defined as “active ligands” for pharmacophore generation. Features of hydrogen bond acceptor and donor, hydrophobic, negative, positive, and aromatic rings were located in the pharmacophore models. Pharmacophores with five features that match to all active ligands were generated by using a tree-based partitioning technique[34] with maximum tree depth of five. The generated pharmacophore hypotheses were scored with default parameters, except that the weight of reference ligand activity is set to 0.3. The top two hypotheses were selected for further generation of 3D-QSAR models. All molecules were aligned in according to selected pharmacophore models.

Table 1

Molecular Structures of Training and Test Set Compounds and Their Clk4 and Dyrk1A Inhibitory Activity

Compounds are from ref (12) other compounds without footnote a are from ref (13).

The most active drugs used for generation of Clk4 pharmacophore model.

The most active drugs used for generation of Dyrk1A pharmacophore model.

The most inactive drugs used for evaluation of Clk4 pharmacophore model.

The most inactive drugs used for evaluation of Dyrk1A pharmacophore model.

Belong to test set for Clk4 QSAR model.

Belong to test set for Dyrk1A QSAR models.

pIC50 values calculated from IC50 data.

pIC50 values predicted based on 3D-QSAR model.

These compounds have IC50 > 10,000nM. pIC50 not available due to lack of exact IC values.

3D-QSAR Modeling

Atom-based 3D-QSAR is advantageous over pharmacophore-based 3D-QSAR in that the former considers the entire molecular space while the latter does not involve area beyond the pharmacophore model.[34,39] In this study, atom-based 3D-QSAR models were generated with training set compounds based on the molecular alignment obtained by pharmacophore generation. In the atom-based model, each atom is represented by a sphere with the van der Waals radius, in accordance to the atom type assigned to each atom. Training set molecules are covered with a regular grid of cubes, with each cube represented with up to six “bits”, representing six different classes of atoms. The atom types are hydrogen-bond donor (D), hydrophobic or nonpolar (H), negative ionic (N), positive ionic (P), electron-withdrawing (includes hydrogen-bond acceptors, W), and miscellaneous (X).[34] The 3D-QSAR partial least-squares (PLS) models were built with three maximum PLS factors in regression model and 1 Å length of the sides of cubic volume elements. The 3D-QSAR models were validated with test set compounds.

Homology Modeling

The crystal structure of Clk4 has not been published yet. A homology model of Clk4 was generated with template of Clk1 by using Prime, Schrodinger.[40] The sequence of human Clk4 was retrieved from the Protein Database at NCBI (http://www.ncbi.nlm.nih.gov/protein). Search of homologous proteins in the NCBI Protein Database (PDB) and sequence alignment were performed through remote access to the BLAST service at NCBI, a function imbedded in Prime. The initial alignment by BLAST was rectified by the second structure prediction (SSP) program SSpro (bundled with Prime), followed by refined alignment obtained via Prime. The homologous model was generated by including template ligand into the model. The initial model was refined with the refinement procedure of Prime. The quality of the final model was accessed by procheck.

Preparation of Receptor and Ligand Molecules for Docking

Low-energy conformations of ligands that were used for docking program Glide were generated via Ligprep[41] of Schrodinger. New structures were produced based on force field OPLS_2005, with protonation states generated at target PH 7.0 ± 2.0. Thirty-two stereoisomers computed by retaining specified chiralities were allowed for each ligand. Protein structures for use by Glide were prepared with the Protein Preparation Wizard[42] of Schrodinger. The structures were first preprocessed with bond order assignment, hydrogen addition, metal treatment, and deletion of all waters in the crystal structures. Hydrogen bonding network and orientation of Asn, Gln, and His residues were optimized based on hydrogen bond assignment. The states of histidine (HIS, HIE, or HIP) were assigned after optimization. Finally, the proteins were minimized to RMSD 0.3 Å based on force field OPLS2005.

Receptor Grid Generation and Docking

Docking is based on a grid represented by physical properties in the receptor volume that is searched for ligand–receptor interaction during docking process. Grid files were prepared with the “Receptor Grid Generation” panel of Glide.[43−45] Grid points were calculated within a region or an enclosing box defined with the centroid of the bound ligand and the size of a docked ligand with length ≤20 Å. To study possible hydrogen bonding interactions with docked ligands, constraints were applied on some Clk4 atoms, i.e., the backbone hydrogen of Leu242, according to the participation of its corresponding residues in hydrogen bonding in crystal structures of Clk1 (PDB ID: 1Z57) and Dyrk1A (PDB IDs: 3ANQ, 3ANR, 2WO6, and 2VX3). Docking was performed by Glide[43−45] of Schrodinger. The score function of Glide, or Glidescore,[43] a modified and expanded version of ChemScore,[46] was used for binding affinity prediction and ligand ranking. The docking can be on the level of either standard (SP) or extra precision (XP). The improvement of XP over SP includes the addition of large desolvation penalties to both ligand and protein, assignment of specific structural motifs that contribute significantly to binding affinity, and expanded sampling algorithms required by scoring function improvement.[44] The XP scoring function comprises four components: Ecoul (Coulomb energy), Evdw (Van de Waals’s energy), Ebind (items favoring binding), and Epenalty (items hindering binding).[44] In this study, extra precision docking was applied, with ligand conformations being generated during docking process. Although the protein keeps rigid, the surface of a ligand is “softened” by scaling the van der Waals radii of nonpolar atoms in order to decrease penalty caused by close contacts. The scaling factor was 0.8, while the partial charge cutoff was 0.15.

Results and Discussion

Pharmacophore Models

Pharmacophore models of Clk4 and Dyrk1A inhibitors were generated with five of the most active compounds. Table 1 showed that these two targets share overlapping but not exactly same active ligands. For example, compound 1 (denoted as comp-1; other compounds represented in the same way) has the highest activity against both enzymes, and comp-2 and comp-7 are active ligands that were used for both model generations. In addition, comp-4 and comp-5 were for generation of Clk4 models, while comp-3 and comp-19 were for Dyrk1A models. A total of 30 and 37 five-point hypotheses were generated for Clk4 and Dyrk1A inhibitors, respectively, by requiring all active ligands matched to the generated hypotheses. The initial hypotheses were evaluated by scoring both active and inactive ligands. Although inactive ligands were not involved in model generation, they were used to eliminate hypotheses that do not distinguish between active and inactive compounds, which is especially useful when all active ligands share common structural skeleton. Figure 1 shows the pharmacophores with the highest adjusted scores mapped to the most active compound 1. Both models are represented with AAARR, indicating they have three hydrogen-bond acceptors and two hydrophobic groups. It is not surprising that the models associated with Clk4 and Dyrk1A have features located at almost the same positions, considering both active sets have common scaffolds. For both models, two acceptors and one hydrophobe are matched to the quinazoline ring, which is shared among all tested compounds. The other two features, or one acceptor and one hydrophobe, are mapped to the R3 substituent 1,3-benzodioxol, which is shared among all active ligands.

Figure 1

Pharmacophore models generated for Clk4 (A) and Dyrk1A (B). Orange: aromatic rings. Rose: hydrogen-bond acceptor.

Atom-Based 3D-QSAR Models

The Clk4 and Dyrk1A inhibitors used for atom-based 3D-QSAR generation were aligned based on the above-mentioned pharmacophore models (Figure 2). One third of tested compounds were assigned to the test set, with training and test set compounds covering the same range of inhibition activities. There were 35 and 17 drugs in the training and test sets associated with Clk4, respectively, while 31 and 15 drugs were associated with Dyrk1A, respectively. Compounds with ICs above 10,000 nM[12,13] (denoted as NA in Table 1) were only used for evaluation of pharmacophore models but were excluded from 3D-QSAR modeling due to the lack of exact IC values. Statistically significant 3D-QSAR models were generated via partial least-squares (PLS) (three factors) based on the training set followed by validation with test set compounds. Table 2 shows the results of atom-based 3D-QSAR models. The correlation coefficients based on training (R2 in Table 2) and test set compounds (Q2) were 0.88 and 0.79 (factor 3) for Clk4, respectively, while those with Dyrk1A were 0.85 and 0.82 (factor 2), respectively. Figure 3 represents the observed versus predicted activities of training and test set compounds regarding Clk4 (Figure 3A) and Dyrk1A (Figure 3B) 3D-QSAR models.

Figure 2

Superimposing of training and test set compounds associated with Clk4 (A) (n = 52) and Dyrk1A (B) (n = 46).

Figure 3

Observed and predicted activities of training (red) and test (blue) set compounds associated with Clk4 (A) and Dyrk1A (B).

The 3D-QSAR models with combined effects of hydrogen bond donor, hydrophobic/nonpolar, electron-withdrawing, and other features were visualized in Figure 4. Figure 4A and B indicate the cubes generated with the most (comp-1) and least active (comp-52) compounds in training set regarding Clk4. The blue regions indicated favorable features contributing to the ligand interactions with target enzyme, while the red ones indicated unfavorable features. The effect of the hydrogen bond donor is revealed by the red region on the substitute R1 of comp-52, indicating a hydrogen on the amine group located on 4-position of the quinazoline ring is unfavorable for activity. The results are supported by the evidence that when R1 is changed from an alkyl group to a hydrogen atom, the activity decreases, as can be seen when comparing comp-1 with comp-14, comp-2/3 with comp-13, comp-10 with comp-20, comp-5 with comp-18, and comp-6 with comp-9.

Figure 4

Atom-based Clk4 (A and B) and Dyrk1A (C and D) 3D-QSAR models visualized with the most and least active compounds.

The comparison between the hydrophobic effects of the most and least active compounds can be seen at position of substituent R3. There is a big blue region at the five-member ring of the benzodioxol group of comp-1, indicating oxygen or hydrophilic atoms could be favorable at this region. On the contrary, for comp-52, there is a red area around the methyl group of the 3-methylphenyl group on the R3 substitution, indicating a hydrophobic group attached to the phenyl ring is unfavorable for activity. This observation is consistent with trend of activity that when the benzodioxol ring in comp-13 is replaced with less hydrophilic groups, such as methylphenyl (comp-52), methoxybenzene(comp-45), and chlorophenyl group (comp-54 and comp-51), the activities decreased dramatically. These results were also consistent with the identified pharmacophore feature, characterized by a hydrogen bond acceptor located at the benzodioxol group. The other difference between the hydrophobic contours regarding comp-1 and comp-52 is that there is a blue area at the 2-position of thiazole ring on R2 substituent of comp-1, indicating a hydrophobic substitute on a meta-position might be favorable for activity. The observation is supported by the fact that comp-4, with a methyl group on the furan ring, is more active than comp-20, which does not have a substituent on the furan ring. The oxygen atoms of the benzodioxol ring also contributed to the electron withdrawing features, indicated with a blue area around 1,3-dioxol group of comp-1.

The combinational effects of hydrogen bond donor, hydrophobic/nonpolar, electron-withdrawing, and other features regarding Dyrk1A are visualized in Figure 4C and D. Similar to the Clk4 activity data, comp-1 is the most potent inhibitor against Dyrk1A. Although comp-52 is among the compounds with lowest activity against Dyrk1A, it was not used in the training or test sets due to lack of exact activity value. Instead, Figure 4D represented the energy fields around comp-51, the compound with lowest activity among training set molecules. The roughly similar activity trends between Clk4 and Dyrk1A account for similar patterns of energy fields occupied by comp-1 and comp-51, compared with their Clk4 counterparts. Similar to the volume regarding Clk4 model, the one occupied by comp-1 (Figure 4C) had three big blue regions: those around the oxygen atoms of benzodioxol ring as R3 substituent, around the 2-methyl-thiazole ring of R2 substituent, and around the methyl group as R1 substituent, indicating a hydrophilic and electron-withdrawing group attached to phenyl ring of R3 substituent, a hydrophobic group attached to the substituting ring at R2 substituent, and a tertiary amine with bulky hydrophobic R1 substitute, could benefit the inhibitory activity. In contrast, the red regions in Figure 4D around the hydrogen of R1 substituent, the thiophen ring of R2 substituent, and the chlorophenyl ring of R3 substituent represent that a hydrogen on R1 and a hydrophobic group beside R2 and R3 could be harmful for the inhibitory effects.

A homology model of Clk4 was generated in previous publications.[5,13] Here, a homology model of Clk4 complexed with ligand 10Z-hymenialdisine was generated with template of Clk1 complexed with the same ligand (PDB ID: 1Z57) by using Prime. Initial sequence alignment was obtained by remote access to the NCBI BLAST service, leading to identification of PDB structures with high sequence identities as the query sequence. The structure with highest sequence identity 86% (Clk1, PDB ID: 1Z57) was chosen as the template to build the structure model for Clk4. Because residues 1–145 and 481 do not have corresponding residues in the template, only a homologous model with residues 146–480 was generated. The high sequence identity between the template and query structures account for a high level of alignment without leaving a gap among matched residues. The initial alignment was adjusted by Prime in terms of comparison between matched residues and secondary structure prediction (Figure S1, Supporting Information). Because all residues in the generated model found their corresponding residues in the template, loop refinement was omitted in the structure refinement procedure. Atoms with homology status of 1 indicate that their side chain coordinates are not taken from the template. For such atoms, coordinates were refined with the “predict side chain” tool of Prime. The refined model was compared with the template to ensure that side chains belonging to the binding site have same orientation as those of the template residues. The quality of the homologous model was assessed with Procheck (Figure S2, Supporting Information).

Binding Mode Identified by Docking

After the characterization of ligand–protein interaction by ligand-based pharmacophore and 3D-QSAR models, it was of interest to explore the interaction in a structure-based approach. The docking of inhibitors 1, 29, and 52 into the Clk4 ATP binding domain was performed with Glide.[35] Figure 5 demonstrated the binding modes obtained from docking without any hydrogen bond constraints imposed on protein atoms. Superimposing of the ligands in Figure 5A showed that they adopted similar poses in the binding pocket, with the R3 substituent at the hydrophilic entrance of the binding cleft sided by residues Asp248, Ser245, Glu290, and backbone of Leu165, Gly166, and Glu167, the quinazoline core overlapping at the bottom of the binding pocket, and R2 substituent fitting into a hydrophobic pocket surrounded by Leu165, Val173, Ala 187, Leu 241–242, and Leu293. In the crystal structures of Clk1[15,47] and Dyrk1A,[16,48] a hydrogen bond between a ligand and Lys191 in Clk1 (Lys188 in Dyrk1A; Figure S3, Supporting Information for sequence alignment among Clk1, Clk4, and Dyrk1A), a residue from a β sheet on one side of the ATP binding cleft is important for ligand–protein interaction. Similarly, binding model obtained by docking ligands to the ATP binding domain of Clk4 (Figure 5B, C, D) indicated that the corresponding residue Lys 189 in Clk4 formed hydrogen bonds with all three ligands.

Figure 5

Obtained binding modes between Clk4 and ligands. (A) Superimposing of compounds 1, 29, and 52 docked to homology model of Clk4. (B) Binding mode between compound 1 and Clk4. (C) Binding mode between compound 29 and Clk4. (D) Superimposing of Dyrk1A (blue, PDBID: 3ANQ) and Clk4 backbones with compound 29. All compounds are represented with ball and sticks, all protein residues in lines, except D248 of Clk4 and D247 of Dyrk1A in (D), which are represented in sticks. All residues are shown with only backbone atoms except D248 of Clk4 and D247 of Dyrk1A in (D).

Compound 1 has the highest inhibition activity among all tested compounds. Above mentioned 3D-QSAR model indicated that a hydrophobic R1 substitute on the position-4 amine is favorable. Figure 5B represented that the methyl group on compound 1 is oriented into a hydrophobic pocket surrounded by the side chains of residues Val173, Ala187, and Phe239, which could increase the van der Walls interaction between compound 1 and Clk4. Compound 29 was selected as a chemical probe for Clk4 that has selectivity of Clk4 against other Clk and Dyrk.[12,13] Figure 5C showed that there is a hydrogen bond between the hydroxyl group on the R3 substituent of compound 29 and the side chain of Asp248, which could contribute to the selective inhibitory effects of this compound against Clk4. The superimposing between structures of Dyrk1A (PDB ID: 3ANQ) and Clk4 is shown on Figure 5D. Compared to the side chain of Asp248 of Clk4, the corresponding atoms of residue Asp247 in Dyrk1A are moving away from the binding pocket by about 2 Å, which could account for the high selectivity of this compound between Clk4 and Dyrk1A.The interaction between Clk4 and ligands identified by docking agreed with the results from ligand-based pharmacophore and 3D-QSAR models. The hydrogen bond between side chain of Lys189 and the nitrogen of quinazoline ring of compounds 1, 29, and 52 was consistent with the hydrogen bond donor feature located on the position-1 nitrogen of quinazoline core identified by the pharmacophore model featured in this study (Figure 4). The orientation of the hydrophilic R3 substituent of compound 29 and 1 to the hydrophilic pocket of Clk4 was supported by the contour maps obtained via 3D-QSAR model indicating that hydrophilic and electron-withdrawing groups were favored in this area. By contrast, the unfavorable interaction between the hydrophobic methylphenyl group of compound 52 and the hydrophilic pocket could account for its much lower inhibitory activity than compound 1 and 29.

It is noticed that there were two hydrogen bond donors featured on the nitrogen atoms of quinazoline ring by the pharmacophore model but only one of them participated in the hydrogen bonding interaction with Clk4. Pharmacophore characteristics are indications of structural properties of ligands interacting with a receptor but do not necessarily identify key features that are responsible for ligand–protein interaction. Because all 6-arylquinazolin-4-amine analogs involved in this study, active or inactive, have the same quinazoline core, the identified two hydrogen bond acceptors associated with this core may not be essential for ligand recognition. For further study, a training set with structurally diverse compounds might be needed to explore the structural space of interacting ligands.

Correlation between Binding Energies and Protein–Ligand Interaction

The docking scores associated with the obtained Clk4 interaction with compounds 1, 29, and 52 were −8.63, −8.61, and −7.72 kcal/mol, respectively. The comparison among binding energies is consistent with the activities of these compounds, with compounds 1 and 29 being much stronger Clk4 inhibitors than compound 52, in terms of their IC50 values (observed pIC50 were 4.96, 3.87, and 2.10, respectively). Compound 1 had a slightly lower binding energy than the compound 29. Although the latter has one more hydrogen bonding interaction than the former, the lower binding energy of compound 1 could be attributed to its favorable hydrophobic interaction of R1 substitute (methyl group in compound 1 vs hydrogen in compound 29) and the more favorable electrostatic interaction of its R3 substitute, which fitted to the hydrophilic pocket sided by residues Asp248, Ser245, and Glu290. The very close binding energies between compounds 1 and 29 could be due to the overestimation of the hydrogen bond effect between Asp248 and the hydroxyl group on R3 of compound 29. The predicted pIC50 values of compounds 1, 29, and 52 were 5.13, 3.75, and 2.25, respectively. Compared with docking scores, the QSAR analysis seemed more effective in distinguishing the inhibitory activities of compounds 1 and 29 against Clk4.

Comparison with Previous Binding Modes between Clk4 and Its Inhibitors

The binding mode between Clk4/Dyrk1A and compounds 1 and 29 was discussed in previous publications.[5,12,13] Homology modeling of Clk4 and docking of 1 and 29 to the ATP binding domain of Clk4 were performed with different programs.[5,12,13] Similar to the current docking results, the previous binding mode between Clk4 and compound 29 indicated a hydrogen bond between the side chain of Asp 248 of Clk4 and the hydroxyl group of compound 29.[13] Previous superimposing of the homology model of Clk4 and crystal structure of Dyrk1A suggested that unfavorable backbone shift of residue Asp247 in Dyrk1A (counterpart residue Asp248 in Clk4) could be responsible for the decreased activity of compound 29 against Dyrk1A than Clk4, which is also confirmed in the present study. However, the difference between the current and previous binding mode is significant. Observed in the current ligand–enzyme interaction, the orientation of the quinazoline core and the R2 substituent attached to the 4-amine group flipped almost a 180 degree from previous position. Therefore, the current mode represented a hydrogen bond between a quinazoline nitrogen and the side chain of Lys 189, instead of between the nitrogen and the backbone of Leu242, a residue located on the hinge region of the ATP binding pocket, in the previous model. Both Lys189 and Leu242 are critical for ligand characterization with Clk4, according to the importance of their corresponding residues in the crystal structures of Clk1 and Dyrk1A. A hydrogen bond between ligand and the counterpart Lys (Lys 191 in Clk1 and Lys 188 in Dyrk1A) is present in all identified crystal structures of Clk1 (PDB ID: 1Z57 and 2VAG) and part of the crystal structures of Dyrk1A (PDB IDs: 3ANQ and 3ANR). By contrast, the involvement of the residue at the same position as Leu242 in the hydrogen bond interaction is only available at one of the Clk1 structures (PDB ID: 1Z57) but is available at all identified Dyrk1A structures (PDB IDs: 3ANQ, 3ANR, 2WO6, and 2VX3).

To further study the interaction between Clk4 and the ligands, alternative binding modes with hydrogen bonding interaction between Leu242 and Compounds 1, 29, and 52 were obtained by imposing H-bond constraint on backbone hydrogen of Leu242, requiring at least one hydrogen bond involving the constrained atom in the protein–ligand complex obtained from docking. The binding modes with a hydrogen bond involving Leu242 of Clk4 and compounds 1, 29, and 52 are shown in Figure S4 of the Supporting Information. The docking scores associated with the above protein–ligand interactions were −7.62, −7.67, and −7.55 kcal/mol, respectively. Compared with the binding modes obtained without any H-bond constraint (−8.63, −8.61, and −7.72 kcal/mol, respectively, for compounds 1, 29, and 52), the docking scores regarding those with Leu242 H-bond interaction are higher, indicating the binding modes with Lys189 participation in hydrogen bond interaction might be more favorable than those with Leu242 participation.

As can be seen from the Clk4-compound 1 complex (Figure S4, Supporting Information), the N3 of the quinazoline ring participated in the hydrogen bonding with Leu242 located at the hinge region. By contrast, a previous publication proposed two hydrogen bonding interactions between both quinazoline nitrogen atoms and the hinge region of Clk4.[13] Except for one hydrogen bond involving the amide NH of Leu242, the other one was marked between the backbone carbonyl oxygen of Glu240 and the N3 of the quinazoline core.[13] Because there is no other hydrogen bond donor close to Leu242 that is pointing to the active site of Clk4, it seems hard that both nitrogen atoms on the quinazoline ring can be involved in hydrogen bonding interaction with the hinge region. Although the current binding mode seems more favored than the previously published one in terms of their docking scores, further study is yet to be explored in order to identify the drug-target interaction associated with arylquinazolines and Clk4/Dyrk1A.

Insights into Design of New Clk4 and Dyrk1A Inhibitors with Higher Affinity and Specificity

The major goal of QSAR analysis and docking is to design new ligands with higher potency and selectivity. Both binding modes (Figure 5) and QSAR analysis demonstrated that a hydrophobic R1 group could be favorable for the inhibition of Clk4. Binding modes indicated that R1 group and the α-carbon of substitute R2 attached to the 4-amino of quinazoline ring were surrounded by a hydrophobic pocket formed with residues Phe239, Val223, Leu242, Val173, and Leu293. Therefore, modification on these two locations with hydrophobic groups might be a means of improving inhibitory activities against Clk4. QSAR prediction based on Clk4 pharmacophore model (Figure 1) indicated that an addition of methyl group to the α-carbon of group R2 of compound 1 could lead to an Clk4 inhibitor with pIC50 of 5.61, higher than the predicted 5.13 of compound 1. QSAR prediction also indicated that substitution of the hydrogen atom with methyl group on the R1 of compound 29 might increase pIC50 value by 0.49, compared with the predicted pIC50 of compound 29, or 3.75. Because compound 29 is a selective inhibitor and a chemical probe of Clk4 over other Clk and Dyrk,[12] the compound with a methyl modulation as R1 could represent a better probe that explores the phenotype particularly down-regulated by Clk4.

Conclusion

6-Arylquinazolin-4-amines have been recently identified as potent Clk and Dyrk1 inhibitors.[5,12,13] Characterization of ligand–protein interaction through ligand-based 3D-QSAR and pharmacophore models combined with structure-based docking will be of great help in future lead compound identification and optimization of novel Clk and Dyrk1 inhibitors. The comparison between the interaction features associated with Clk4 and Dyrk1A might shed light on the design of selective Clk4 and Dyrk1A inhibitors. In the present study, we have developed pharmacophore and atom-based 3D-QSAR models for the Clk4 and Dyrk1A inhibitory effects of a series of 6-arylquinazolin-4-amines. The high R2 and Q2 (0.88 and 0.79 for Clk4, 0.85 and 0.82 for Dyrk1A, respectively) based on validation with training and test set compounds suggested that the generated 3D-QSAR models are reliable in predicting novel ligand activities against Clk4 and Dyrk1A. Integrating molecular docking with ligand-based SAR models allows us to use structural information to further investigate ligand–protein interaction. The interactions identified through docking ligands to the ATP binding domain of Clk4 were consistent with the structural properties and energy field contour maps characterized by the pharmacophore and 3D-QSAR models and gave valuable hints regarding the structure–activity profile of 6-arylquinazolin-4-amine analogs, suggesting that the obtained protein-inhibitor binding mode is reasonable. The 3D contour maps obtained through atom-based 3D-QSAR modeling in combination with the binding mode between inhibitor and residues of Clk4 obtained with docking provide valuable insights into the rational design of novel Clk4 and Dyrk1A inhibitors, especially 6-arylquinazolin-4-amine analogs.