Abstract

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 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. 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 technique34 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.

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

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 Ligprep41 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 Wizard42 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 Glide43−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 ​Figure11 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.

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 ​(Figure2).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 nM12,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 ​Figure33 represents the observed versus predicted
activities of training and test set compounds regarding Clk4 (Figure ​(Figure3A)3A) and Dyrk1A (Figure ​(Figure3B)
3D-QSAR3B)
3D-QSAR models.

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 ​Figure4.4. Figure ​Figure4A4A 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.

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 ​Figure4C4C 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 ​Figure4D4D 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 ​(Figure4C)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 ​Figure4D4D 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.

Homology Modeling

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 ​Figure55 demonstrated
the binding
modes obtained from docking without any hydrogen bond constraints
imposed on protein atoms. Superimposing of the ligands in Figure ​Figure5A5A 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 Clk115,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 ​(Figure5B,5B, C, D) indicated that the corresponding residue Lys 189 in Clk4
formed hydrogen bonds with all three ligands.

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 ​Figure5B5B 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 ​Figure5C5C 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 ​Figure5D.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 ​(Figure4).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 ​(Figure5)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 ​(Figure1)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.

Acknowledgments

This research was supported by the Intramural
Research Program
of the National Institutes of Health (NIH), National Library of Medicine
(NLM). We thank Helix Systems, High-Performance Computing at the NIH
for making Maestro, Schrodinger available for this study.

Glossary

Abbreviations

Clk4

Cdc2-like kinase

Dyrk1A

dual specificity tyrosine-phosphorylation-regulated
kinase 1A

QSAR

quantitative
structure–activity
relationship

Funding Statement

National Institutes of Health, United States

Supporting Information Available

Supporting Information Available

Sequence alignment compared
with second-structure prediction between Clk4 and Clk1. The quality
of homologous model of Clk4 evaluated with Procheck. Sequence alignment
among Clk1, Clk4, and Dyrk1A obtained by using ClustaIW2. This material
is available free of charge via the Internet at http://pubs.acs.org.