Small Molecule Inhibitors of DYRK1A Identified by Computational and Experimental Approaches

Dual-specificity tyrosine phosphorylation-regulated kinase 1A (DYRK1A) is a protein kinase with diverse functions in cell regulation. Abnormal expression and activity of DYRK1A contribute to numerous human malignancies, Down syndrome, and Alzheimer’s disease. Notably, DYRK1A has been proposed as a potential therapeutic target for the treatment of diabetes because of its key role in pancreatic β-cell proliferation. Consequently, DYRK1A is an attractive drug target for a variety of diseases. Here, we report the identification of several DYRK1A inhibitors using our in-house topological water network-based approach. All inhibitors were further verified by in vitro assay.


Introduction
The phosphorylation of proteins catalyzed by protein kinases plays a major role in the regulation of cellular processes, such as cell proliferation, differentiation, apoptosis, and signal transduction. Abnormal protein phosphorylation has been implicated in several diseases. Consequently, protein kinases have emerged as major drug targets [1]. Dual-specificity tyrosine phosphorylation-regulated kinase 1A (DYRK1A) belongs to the DYRK family of kinases, which includes four other members (DYRK1B, DYRK2, DYRK3, and DYRK4). This kinase regulates critical cellular processes, including the proliferation and differentiation of neuronal progenitor cells [2]. DYRK1A is abnormally expressed in Down syndrome, Alzheimer s disease, Pick s disease [3], lung cancer, cervical cancer, gastrointestinal stromal tumors (GIST), glioblastoma, melanoma, acute megakaryoblastic leukemia, acute lymphoblastic leukemia, and acute myeloid leukemia [4][5][6][7]. Recently, DYRK1A was found to be involved in human pancreatic β-cell proliferation, making it a potential therapeutic target for the treatment of Type 1 and Type 2 diabetes [8][9][10][11][12]. Insufficient pancreatic β-cell mass or function leads to diabetes mellitus. Under high glucose conditions, β-cells increase intracellular calcium (Ca 2+ ) levels. This activates calcineurin, which in turn dephosphorylates the nuclear factor of activated T-cells cytoplasmic (NFATc) proteins. The dephosphorylation and activation of NFATc lead to nuclear import. Nuclear NFATc kinases, such as DYRK1A and GSK3β phosphorylate NFATc proteins, cause them to undergo nuclear export. The inhibition of DYRK1A and GSK3β blocks NFATc nuclear export and increases β-cell proliferation [9,13].

Results
We performed MD simulations on the DYRK1A and GSK3β structures in the apo state, and analyzed TWNs in their binding site. Staurosporine is a potent pan-kinase inhibitor [39]. The planar structure of staurosporine with few rotatable bonds allows it to occupy the adenosine triphosphate (ATP) binding sites of kinases. The co-crystal structure of human GSK3β, in complex with staurosporine (PDB code 1Q3D), is available [40]. However, the co-crystal structure of DYRK1A, in complex with staurosporine, had not yet been determined. We therefore obtained the crystal structure of human DYRK1A (PDB code 4YLL) [29] and docked staurosporine into its ATP binding site. For TWN analysis, we divided the ATP binding site into five regions (A-E), based on staurosporine's binding mode (Figure 1).

DYRK1A vs. GSK3β
DYRK1A is a dual-specificity kinase that possesses both serine/threonine and tyrosine kinase activities, while GSK3β is a serine-threonine kinase. Both DYRK1A and GSK3β have been implicated in diabetes [9,13]. As shown in Table 1, these kinases do not have a high total sequence similarity or binding site similarity. Total sequence similarity and binding site similarity were found to be 44.1% and 32.0%, respectively. Staurosporine exhibits comparable activity against these kinases. It showed IC 50 values of 19 and 15 nM against DYRK1A [20] and GSK3β [41], respectively. Sequence and binding site similarities were not able to account for the comparable activities of staurosporine against DYRK1A and GSK3β. Thus, staurosporine was extracted from the co-crystal structure of GSK3β (PDB code 1Q3D) and re-docked into the ATP binding site. DYRK1A and GSK3β docking results were compared. Although they exhibited similar IC 50 values, the binding energies differed by more than −10 kcal·mol −1 . Staurosporine displayed a binding energy of −68.9 and −79.0 kcal·mol −1 for DYRK1A and GSK3β, respectively. We therefore analyzed TWNs in the ATP binding sites of these kinases. This analysis revealed a high percentage of TWNs in the C regions (hinge region) of both kinases (Table 1 and Figure 2). DYRK1A exhibited 40.5% TWNs while GSK3β demonstrated 36.4% TWNs in the C region. The hinge regions of kinases are known to play a key role in ligand binding. Furthermore, we calculated the TWN-ligand shape similarity and observed comparable values of 54% and 61% for DYRK1A and GSK3β, respectively. Based on the TWN results, we anticipated that known GSK3β inhibitors with a high percentage of TWNs in the C region and reasonable TWN−ligand shape similarities could be repositioned as DYRK1A inhibitors.

TWN-Based Repositioning
Compounds AZD1080 [42] and SB-415286 [43] are known GSK3β inhibitors with IC 50 values of 31 and 78 nM, respectively. These compounds were docked into the binding site of DYRK1A, and TWNs were analyzed around the ligands. AZD1080 and SB-415286 showed binding energies of −58.3 and −83.9 kcal·mol −1 , respectively ( Table 2). Similar to staurosporine, both compounds displayed a high percentage of TWNs in the C region (hinge region) of DYRK1A. AZD1080 showed 61.8% TWNs, while SB-415286 showed 47.2% TWNs in the C region. Docking results revealed that both compounds formed hydrogen bonds with hinge residues Glu239 and Leu241, which occupy the C region. Additionally, they showed hydrogen bond interactions with Lys188 ( Figure 3). These residues are known to play important roles in DYRK1A kinase activity and in the binding of DYRK1A inhibitors [14,44]. We calculated TWN-ligand shape similarity for the compounds and obtained reasonable similarities of 36% and 44% for AZD1080 and SB-415286, respectively. Compounds AZD1080 and SB-415286 are available commercially. We purchased compounds AZD1080 and SB-415286 and verified their DYRK1A inhibitory activities by in vitro assay. In accordance with the TWN-ligand shape similarity values, SB-415286, with an IC 50 of 445 nM, was found to be more potent than AZD1080 (IC 50 = 2911 nM).

Docking and TWN-Based Screening
Our results suggest that compounds with a high percentage of TWNs in the C region and reasonable TWN-ligand shape similarity could inhibit DYRK1A. We docked our in-house compounds into the binding site of DYRK1A to identify potential ligands. Structures of the selected ligands (compounds 1−7) are displayed in Figure 3. Binding energy and TWN results are provided in Table 2.
The selected compounds showed a higher percentage of TWNs in the C region compared with the other regions. Except for compound 6, all compounds showed TWNs of more than 50% in the C region. We calculated TWN-ligand shape similarity for these ligands. Compounds 5 and 7 showed low similarity (<25%). Similar to AZD1080 (36%), compounds 1, 2 and 6 exhibited moderate similarity (25−40%). Similar to staurosporine (54%) and SB-415286 (44%), compounds 3 and 4 showed high similarity (>40%). We verified DYRK1A inhibitory activities of these compounds by in vitro assay. In accordance with the TWN-ligand shape similarity, compounds 5 and 7 displayed low DYRK1A inhibition (<40%) at 10 µM. Similar to AZD1080 (71%), compounds 1, 2 and 6 exhibited moderate inhibition (40-80%). However, in contrast to SB-415286 (91%), compounds 3 and 4 did not show high DYRK1A inhibition. Compound 3 demonstrated 82% inhibition, whereas compound 4 exhibited 71% inhibition. Structural analysis revealed that compounds 3 and 4 are smaller compared with the other compounds. They could not cover the A and D regions because of their small size. This might be the reason for their moderate inhibition values, despite the high TWN-ligand shape similarity. Based on the percentage inhibition data, IC 50 values against DYRK1A were calculated for all the selected compounds, except compounds 5 and 7. The activity values are provided in Table 2.

Discussion
DYRK1A phosphorylates NFAT, thereby inhibiting the effects of calcium signaling and maintaining NFAT in an inactive state [45]. DYRK1A is a negative regulatory factor of the cell cycle that promotes the G0 state or conversion to differentiation. In malignant cells, DYRK1A promotes cell survival by inhibiting pro-apoptotic signaling [46].
Recently, DYRK1A was found to be associated with human pancreatic β-cell proliferation [10,32]. Type 1 and type 2 diabetes are characterized by a decrease in pancreatic β-cell mass [9]. β-cells in high glucose conditions increase intracellular calcium levels. Activation of calcineurin leads to dephosphorylation of NFATs and subsequent nuclear translocation. As nuclear NFATc kinases, DYRK1A and GSK3β phosphorylate NFATc proteins to induce nuclear export. Both DYRK1A and GSK3β are negative regulators of the NFAT pathway, and this pathway is fundamentally important for β-cell proliferation [9,13].
Our research group has developed and is optimizing algorithms for determining TWNs. In our previous study, we performed MD simulations on 26 kinases in aqueous solution and analyzed TWNs in their ATP binding sites [35]. Other previous studies used TWN analysis to study protein hydration [38], kinase selectivity [33], kinase inhibitor design [35,36], and drug repositioning [35].
The purpose of this study was to identify DYRK1A inhibitors by screening our in-house library using TWN-ligand shape similarity based on known inhibitors of other kinases with TWN patterns similar to DYRK1A. Through multiple kinase TWN analyses, GSK3β with low binding site similarity, but high distribution of water molecules at the C region compared to DYRK1A, was selected. Notably, GSK3β participates in the same signaling pathway as DYRK1A, and their staurosporine inhibition values are similar. The GSK3β inhibitors AZD 1080 and SB-415286 were tested for their IC 50 values against DYRK1A. The binding energy and TWN-ligand shape similarity were also analyzed. The IC 50 values revealed better inhibition by SB-415286 than AZD1080 against DYRK1A. Compounds were selected from the in-house library based on the IC 50 values and TWN-ligand shape similarity of staurosporine, AZD 1080, and SB-415286 against DYRK1A. In the TWN analysis, compounds that account for more than 35% of the C region were selected through screening. The IC 50 values for compounds with TWN-ligand shape similarity values of less than 25% were not determined. Compounds with TWN-ligand shape similarities of 25-40% showed IC 50 values of 2.9-6.8 µM. Compounds with more than 40% TWN-ligand shape similarity showed IC 50 values of less than 0.5 µM.
However, compounds 3 and 4 deviated from these trends, with TWN-ligand shape similarity values over 50%, but IC 50 values of about 3.2 µM. We investigated the underlying cause of the deviation by analyzing the binding mode predicted through molecular docking studies. The analysis revealed that because the compound sizes were small, they could not extend in the direction of the hinge (Glu239, Leu241) or the key residue (Lys188), which is key to kinase competitive inhibition. These compounds nonetheless have inhibitory activity and are suitable for use as a scaffold. To develop more effective inhibitors from these compounds, optimization is required by adding functional groups in the direction of the hinge region as well as key residues. In a future study, we plan to optimize the TWN-ligand shape similarity range and compare the electro-shape similarities.
In this paper, we described how to identify a DYRK1A inhibitor from another known inhibitor using the TWN-ligand shape similarity method. As a computational drug discovery method, we propose the TWN-ligand shape similarity method through TWN analysis as a way to rapidly identify compounds amenable to drug repositioning. The TWN-ligand shape similarity method can be used to search for target compounds by acquiring scaffolds through high throughput screening (HTS) and prediction of biological activity.

Protein Preparation
The X-ray crystal structures of human DYRK1A (PDB code 4YLL) [29] and human GSK3β (PDB code 1Q3D) [40] were downloaded from the Protein Data Bank (PDB) and all non-protein molecules were discarded. These structures were further processed using the Prepare Protein module in Discovery Studio 2017 (BIOVIA, San Diego, CA, USA). This process included the identification of missing residues, addition of hydrogen atoms, assignment of bond orders, and formal charges. Protonation states were assigned under the assumption that the systems were at a pH of 7.4.

Molecular Dynamics (MD) Simulation
GROMACS is a freely available, versatile package to perform MD simulation. It is extremely fast and supports all the usual algorithms you would expect from a modern MD implementation. It provides extremely high performance compared to other programs and contains quite a few features that make it stand out from the competition. It is user-friendly, with a fully automated topology builder for proteins. There is plenty of consistency checking, and clear error messages are issued if something is incorrect. It can be run in parallel and can write MD trajectory data in a very compact way. There is no need to write any code to perform usual MD analysis, as it provides many flexible tools for the analysis [47]. Due to these reasons, we selected GROMACS for the MD simulation. MD simulations were performed with GROMACS 4.5.3 [47]. Protein topology was generated using CHARMM27 force field [48]. Protein was solvated in a cubic box using the TIP3P water model [49]. Counter ions were added to ensure the neutrality of the system. Then, the system was subjected to energy minimization using 500 steps of the steepest descent algorithm. This was followed by equilibration over two stages. Firstly, the system was equilibrated in the NVT ensemble (constant number of particles, volume, and temperature) for 0.1 ns, followed by equilibration in the NPT ensemble (constant number of particles, pressure, and temperature) for 0.2 ns. NVT equilibration was executed at a temperature of 300 K using V-rescale thermostat [50], whereas NPT equilibration was carried out at a pressure of 1 bar using the Parrinello-Rahman barostat [51]. Finally, a production run was performed for 10 ns with a time step of 1 fs. Periodic boundary conditions (PBC) were applied to the system. In MD simulation, PBC are usually applied to avoid problems with boundary effects. PBC make it possible to simulate a small system that is not terminated by a surface, as it is periodically repeated in all directions [47]. A linear constraint solver (LINCS) algorithm was used to constrain bonds [52]. The LINCS algorithm provides a high performance and it is faster and more stable than other constraint algorithms [47]. A cut-off distance of 1.2 nm was used for all short-range non-bonded interactions, while long-range electrostatics were calculated using the Particle mesh Ewald method [53]. After MD simulation, 100 trajectory files were extracted for TWN analysis.

Topological Water Network (TWN) Analysis
Water molecules form water-ring networks through hydrogen bonds, which the authors have termed topological water networks (TWNs) [33][34][35][36][37][38]. These networks include small rings, such as trimers (R3), tetramers (R4), pentamers (R5), and hexamers (R6). The potential functions considered in the TWNs involve a rigid TIP3P water model. The interactions between water molecules are conveniently modeled using Lennard-Jones and Coulomb potentials [49]. The interactive potential energy between two water molecules (a and b) is expressed by the Equation (1)  (1) where, v(a, b) = interaction potential energy r oo = distance between oxygen atoms q i = partial charge on the i site (−0.834e) q j = partial charge on the j site (0.417e) r ij = distance between q i and q j A = repulsive force between i and j (582,000 kcal·Å 12 mol −1 ) C = attractive force between i and j (595 kcal·Å 6 mol −1 ) Parameters were chosen in such a way that they produced reasonable structural and energetic results for liquid water. The energy criterion of −2.25 kcal·mol −1 was used to determine hydrogen bonding between water molecules. This value was selected as a criterion because it closely corresponds to the minimum value of the water-water pair potential energy distribution [49].

Binding Site Similarity
Binding site similarity was calculated using the geometric hashing method [54]. This method compares a set of binding sites quickly. The algorithm identifies equivalent heavy atoms between binding sites and matches them in the same relative spatial orientation. Binding site similarity is expressed by the following Equation (2): where R 3 represents the similarity score. It takes into account the total size of the two binding sites (n site1 and n site2 ). It is calculated analogously to the Tanimoto coefficient, and its value ranges from 0 to 1. A value of 1 indicates the self-comparison of a binding site. The n match denotes the number of atoms comprising the largest possible matching [55].

TWN-Ligand Shape Similarity
Shape similarity was calculated using the ultrafast shape recognition (USR) method [56]. This method is based on the assumption that the relative position of atoms defines the shape of a molecule. The molecular shape is described by a set of one-dimensional distributions with three-dimensional shape information. The USR method uses the distributions of all the atomic distances to four different reference locations: the molecular centroid (ctd), the farthest atom to ctd (fct), the closest atom to ctd (cst), and the farthest atom to fct (ftf ). The first three moments from each of the four one-dimensional distributions are considered to describe a molecule, as in Equation (3): Shape similarity is estimated by the Equation (4) below: where S qi is the similarity score function, are the vectors of shape descriptors for the query and the ith screened molecule, respectively.

Molecular Docking
Crystal structures of proteins were obtained and processed as described in the protein preparation section. Molecular docking studies were performed on the processed structures using the LigandFit module [57] of Discovery Studio 2017 (BIOVIA). The Prepare Ligand protocol was used to build and optimize ligands. Partial charges were assigned using the Momany-Rone partial charge method. Energy minimization was carried out with the CHARMM force field. The binding site was defined based on the co-crystallized ligand. For each ligand, 50 docked poses were generated and scored using scoring functions. Protein-ligand interactions were considered for selecting the binding modes of the ligands.

In Vitro Assay
Enzymatic assays were performed by Eurofins Scientific Inc. Korea (Brussels, Belgium). DYRK1A(h) was incubated with 8 mM MOPS pH 7.0, 0.2 mM EDTA, 50 µM RRRFRPASPLRGPPK, 10 mM MgAcetate, and (γ-33P-ATP (specific activity approx. 500 cpm/pmol, concentration as required). The reaction was initiated by the addition of the MgATP mix. After incubation for 40 min at room temperature, the reaction was stopped by the addition of 3% phosphoric acid solution. Then, 10 µL of the reaction was then spotted onto a P30 filtermat and washed three times for 5 min in 75 mM phosphoric acid and once in methanol prior to drying and scintillation counting. IC 50 was calculated for inhibitors, including staurosporine (from 10mM DMSO stock solution), depending on various final concentrations. All assays were performed in duplicate, and the average IC 50 value was reported.

Conclusions
In conclusion, we identified inhibitors of DYRK1A using a computational TWN-based approach, and we subsequently verified their inhibitory activity experimentally. More potent DYRK1A inhibitors can be developed through further optimization of these molecules.