Computer-aided design, synthesis and biological characterization of novel inhibitors for PKMYT1

Abstract
In the current work, we applied computational methods to analyze the membrane-associated inhibitory kinase PKMYT1 and small molecule inhibitors. PKMYT1 regulates the cell cycle at G2/M transition and phosphorylates Thr14 and Tyr15 in the Cdk1-cyclin B complex. A combination of in silico and in vitro screening was applied to identify novel PKMYT1 inhibitors. The computational approach combined structural analysis, molecular docking, binding free energy calculations, and quantitative structure–activity relationship (QSAR) models. In addition, a computational fragment growing approach was applied to a set of previously identified diaminopyrimidines. Based on the derived computational models, several derivatives were synthesized and tested in vitro on PKMYT1. Novel inhibitors active in the sub-micromolar range were identified which provide the basis for further characterization of PKMYT1 as putative target for cancer therapy.

Introduction
The membrane-associated inhibitory kinase PKMYT1 belongs to the Wee1-kinase family and regulates the cell cycle at G2/M transition [1]. PKMYT1 plays a major role in the inactivation of the second checkpoint G2 through its inhibitory phosphorylation of cyclin-dependent kinase (Cdk1) [1, 2]. Functional PKMYT1 kinase phosphorylates Thr14 in the Cdk1-cyclin B complex following a phosphorylation of Tyr15 by Wee1 kinase. This double phosphorylation of Cdk1 leads to a disassociation between cyclin B and Cdk1 and consequently deactivates its function [3]. As a result, the cell cycle is restricted until DNA damage is repaired [2]. Many cancer cells lack a functional p53 signaling, which abrogates the G1 checkpoint, resulting in increased DNA damage at the G2 checkpoint compared to normal cells [4]. Thus, cancer cells rely more on the G2 checkpoint than normal cells [5]. A new strategy for cancer treatments is to keep the cell in the cell cycle with unrepaired DNA damage in premature mitosis. The abrogation of the G2 checkpoint can be induced by pharmacological manipulation, resulting in mitotic catastrophe with a high level of unrepaired DNA damage and immediately causes apoptotic or non-apoptotic cell death [6, 7]. Several PKMYT1 inhibitors including the known tyrosine kinase inhibitors dasatinib and bosutinib [8], the pyridopyrimidine derivatives PD-0166285 [9], PD-173952 [10], PD-173955[8] and PD-180970 [9] have been identified by applying different approaches. Examples of the most active compounds from each chemotype are given in Figure 1.

In addition we recently screened the GSK kinase inhibitors sets (PKIS sets I and II) [11-13] and identified novel PKMYT1 inhibitors belonging to different chemical classes [13].Until recently, only one crystal structure of PKMYT1 in the apo form was available. In a previous study we applied ligand docking and QM/MM based binding free energy (BFE) calculations to virtually screen several compound databases and identified active inhibitors [10]. However, we recognized that the derived hits were either highly similar to already known PKMYT1 inhibitors or showed only weak inhibitory activity. To overcome the limitations of the previously developed models and inhibitors, we were interested in developing novel computational models with increased accuracy. For this purpose we included in the current work eight new crystal structures of PKMYT1 that were released very recently by Zhu et al. [14]. The new PKMYT1 structures were resolved in the active form and in complex with known kinase inhibitors and thus represent more suitable starting points for generating computational models for inhibitor design.In addition, we extended the ligand data set for establishing computational models for activity prediction. We considered all known PKMYT1 inhibitors including the ones that we identified from the GSK kinase inhibitors sets. (Supplementary data: Figure 1S, Table 1S). Based on docking and BFE derived models, we then applied a fragment growing approach to guide the structural optimization. The most promising compounds were subsequently synthesized and theirinhibitory activity on PKMYT1 was tested [15]. As a result, a set of novel diaminopyrimidines could be identified which should serve as valuable tools to further analyze the potential role of PKMYT1 inhibitors in cancer therapy.

Results and Discussions
Totally, nine crystal structures of PKMYT1 kinase are available in the protein Data Bank (PDB), eight of which were resolved in holo form (PDB IDs: 5VCV, 5VCW, 5VCX, 5VCY, 5VCZ, 5VD0, 5VD1, and 5VD3) and one structure in apo form (PDB ID: 3P1A) [14, 16]. Two structures (PDB IDs: 5VCX and 5VD3) were co-crystallized with saracatinib in different phosphorylation states of the cloned PKMYT1, where saracatinib was found to display an identical conformation.The ATP-binding pocket of PKMYT1 and other kinases is generally divided into a back pocket behind the gatekeeper (PKMYT1: Thr187), a front pocket within the solvent-accessible area, an adenine region close to the hinge region, a ribose-binding pocket and a phosphate binding region [16, 17]. The individual pockets are constructed by several amino acids which interact with ATP or other ligands occupying the pockets. Moreover, the surrounding amino acids and the accessibility of each pocket to the solvent impart distinct properties for the individual sub- pockets. The back-pocket in PKMYT1 is characterized as a hydrophobic pocket due to the surrounding hydrophobic amino acids: Val171, Leu185, Val124, and Lys139 (Figure 2). In addition, the back-pocket shows a limited accessibility to hydrophilic amino acids in the αC helix (His161 and Glu157). Meanwhile, the front-pocket is located in a solvent-accessible area and is surrounded by Leu116, Gly191 and Pro192 (Figure 2). The adenine region is located close to the hinge region to allow the embedded ligand/ATP to form conserved hydrogen bonds with the hinge region residues.

The general characteristics of the sugar pocket and phosphate region are solvent-accessible and hydrophilic. The phosphate region is surrounded at the bottom by Asp251 of the DFG motif and at the top Tyr121 in the P-loop.side chain of Phe240, targeting the back hydrophobic pocket by a substituted phenyl or a small chemical group, in addition to polar and van der Waals interactions with the front pocket residuese.g. Leu116, Pro191 and Gln196 [13]. The polar interactions with the phosphate and sugar regions residues e.g. Asp251 and Tyr121were also observed in several inhibitors [13].We found in our previous work that the docking program GOLD was performing well in redocking studies. GOLD gave the best results for all included inhibitors when the PKYMT1 structure bound to pelitinib (PDB ID: 5VCW) was taken for ligand docking [13]. Therefore, we docked all inhibitors under study into this protein conformation and used the derived poses to find a correlation between the computed and experimental affinities [13]. As training set we used 19 PKYMT1 inhibitors (whole data set shown in Supplementary data, Figure S1, Table S1) identified by us and the experimentally determined inhibitory activity from the activity assay (FPIA) and the binding assay (FPBA), respectively. Firstly, we investigated a correlation between the experimentally determined activities and several docking scores. Therefore, the obtained docking solutions were re-scored using implemented scoring functions in GOLD and Glide. Low correlation coefficients were observed using scoring function (R2 values ranged between 0.35 to 0.38, Supplementary data, Table S2) in agreement with our previous study [10].

Further studies were performed to predict the affinities of the ligands using different BFE calculation methods (MM-PB/SA and MM-GB/SA, as well as hybrid QM/MM methods). Rescoring was performed using either a single frame after short minimization or considering ten snapshots from shortmolecular dynamics (MD) simulations (100 ps). The solvation free energy was calculated using Poison Boltzmann (MM-PB/SA) and different models of the Generalized Born approach (MM(QM/MM)-GB/SA). These GB models were implemented in AMBERTools12 as GBHCT (igb= 1), GBOC1 (igb = 2), GBOC2 (igb = 5) and IGB8 (igb = 8) [18-23]. The QM region of the QM/MM hybrid mechanics was firstly applied for the ligand only and later for the ligand and selected residues from the hinge region, the P-loop, the hydrophobic pocket and DFG motif of PKMYT1, surrounding the ligand. Two models for QM were used: Parameterized Model number 3 (PM3) and Austin Model 1 (AM1) [24]. The performance of the binding free energy methods was evaluated by computing the correlation coefficients between the biological activity (Ki values) and the enthalpy score (∆H). The results are summarized in Table S3 (see Supplementary data).The best correlation coefficient for the training set was 0.68 (RMSE was 0.59 kcal.mol-1) by applying the Poison Boltzmann (PB) solvation model (MM-PB/SA) and using 10 snapshots from a short MD simulation (Figure 4, Supplementary data, Table S3). The leave-one-out cross- validated q2 was 0.60, and the RMSE 0.59 kcal.mol-1. Using only one frame after minimization the model showed slightly weaker correlation (R2 = 0.62, RMSE = 0.65 kcal.mol-1 and cross- validated q2 = 0.51. Using the GB solvation models resulted in poor correlation coefficients in all applied settings. Also the hybrid QM/MM-GB/SA models failed to correctly predict the experimental activity (Supplementary data, Table S3).Table 1 displays ∆H values using (MM-PB/SA method) for the inhibitors under study.

MM- PB/SA gave the lowest (favourable) value for PD-0166285 among the inhibitors -49.80 kcal.mol-1 with very low Z-score (Z score= 0.58). Meanwhile, for weakly active compounds e.g. tyrphostin and GSK3206866A unfavourable values were obtained (-24.12 and -27.17 kcal.mol-1, respectively). Thus, the applied MM-PB/SA model is able to accurately identify the high and weakly active compounds. The comparison between pyridopyrimidine derivatives (PD-series) emphasizes the latter observation. Furthermore, Z-scores for all inhibitors were less than 2 with the exception of saracatinib (Z-score of 2.11).Based on the MM-PB/SA-derived model, which showed a significant correlation with the experimentally determined data, further improvement of this correlation was analyzed by including two-dimensional (2D) molecular descriptors. The available 2D molecular descriptors implemented in Molecular Operation Environment System (MOE; Chemical Computing Group Inc., Montreal, QC, Canada) were tested. As a result, incorporating the 2D descriptor PEOE_VSA-4 improved the correlation and increased R2 of the MM-PB/SA model to 0.77 and cross-validated q2 to 0.70 (RMSE decreased to 0.50 and 0.58 kcal.mol-1, respectively) (Figure 4).PEOE_VSA-4 indicates electro-negative parts of an approximated accessible van der Waals surface area (in Å2) [25]. The following equation (Eq. 1) shows the generated QSAR model (named QSAR_1A) that includes two descriptors (MM-PB/SA enthalpy score ∆H and PEOE_VSA-4) to predict Ki values of PKMYT1 inhibitors.QSAR_1A: Pred_Ki = -0.11803* (MM-PB/SA) -0.03631* (PEOE_VSA-4) +2.81241 (1)In order to design novel inhibitors for PKMYT1, two different approaches were employed in the current work.

The hybrid approach, fragment growing and single point modification, was used to design diphenyl diaminopyrimidine derivatives. In the previous work, the conserved interactions and the relevant interactions of the PKMYT1 binding site have been discussed [13]. In the current work we used this information and the derived BFE/QSAR models to optimize selected hits.First, ligand efficiency (LE) and ligand lipophilicity efficiency (LLE) were calculated for the starting fragments (Figure 5). LLE does not take into account the size of the ligand, therefore, it can be used in the optimization process for the diaminopyrimidine series (LLE = 1.41) in order to improve the potency without solely increasing the lipophilicity [26-28]. Thediphenyldiaminopyrimidine scaffold was substituted to improve the potency from a predicted Ki value of 33.25 µM for the starting fragment to design selective compounds with predicted Ki in the nanomolar range. The visualization of the putative binding mode, structural analysis of the PKMYT1 binding site, and structure-activity relationship study of the diaminopyrimidine derivatives revealed that three vectors on the diaminopyrimidine scaffold could be expanded to target the adjacent pockets (Figure 5). Based on the derived QSAR_1A model, the biological activities of the designed compounds were predicted.The obtained results of the fragment-based growing approach on the diaminopyrimidine scaffold helped to select useful substituents to optimize the compounds. All synthesized compounds were first tested in vitro and IC50 values were determined (Table 5).

Since our QSAR model is based on Ki values, we used the Cheng-Prusoff equation [29] to calculate the corresponding IC50 values for the identified actives based on the predicted Ki values.Synthesis of the diaminopyrimidine derivatives was realized according to a protocol previously reported for the synthesis of potential aurora kinase-inhibitors [30] and an adapted protocol for microwave synthesis conditions [31]. Figure 9 shows the general synthetic scheme. The two step synthesis started with the reaction of the respective 5-substituted-2,4-dichloropyrimidine derivatives 1 with the respective aniline derivatives 2 (substituted R2-5) in ethanol in the presence of DIPEA to quench the released HCl. In the second synthesis step, the obtained intermediates, which were used without further purification, were reacted with morpholino-aniline 4 using HClas catalyst in 1-butanol. Purification of the final compounds was performed via preparative HPLC.First, a compound was synthesized to check whether the back-pocket is able to accommodate large groups. A benzyloxy moiety at position R1 was chosen to target the back-pocket, while using a 3-hydroxy-4-methylphenyl group and a morpholine ring to target the ribose/phosphate- and front-pocket, respectively. Therefore, compound 5a was synthesized and biologically tested, where it did not show any activity against PKMYT1 (displacement percentage less than 10%, Table 5). We hypothesized, that the back-pocket might not to be suitable to accommodate large groups. Hence, the analogous compound 5l, which contains a methoxy group instead of the benzyloxy moiety, was synthesized and tested. It noteworthy, that compound 5l showed a high displacement percentage 100%.

We thus decided to only use small hydrophobic groups at position R1 of the pyrimidine ring to target the back-pocket.The chosen substitutions at position R1 tried to mimic those found in active compounds and included chloro, bromo, methyl, methoxy and methylsulfanyl substituents. Methyl amide was proposed for position R2, whereas a hydroxyl group was introduced at position R3. Position R4 was substituted by either methoxy or methyl groups. Finally, a chloro substituent was used at position R5. So, novel hits presented in Table 5 were synthesized and biologically evaluated.The first set of compounds was substituted at position R1 by several small groups and position R3 by a hydroxyl group, whereas the remaining positions at the pyrimidine were unsubstituted. As a result, five compounds containing chloro, bromo, methyl, methylsulfanyl, methoxy groups (compounds 5f, 5h, 5b, 5d, and 5k) were biologically evaluated at 20 µM in order to determine the displacement percentage and assess their activities. Their displacement percentages ranged from 21.11% till 100% (Table 5). One exception is derivative 5b, containing a methyl group at the position R1, which could not displace ATP in PKMYT1 showing a displacement percentage less than 10% (Table 5).

The substituents at position R5 were evaluated: derivatives 5e and 5g contain a chloro substituent at position R5 and a hydroxyl or methyl group at position R3, whereas position R1 contains either chloro or methylsulfanyl moieties (Table 5). Compounds 5f, and 5d, containing chloro andmethylsulfanyl groups at R1 position, showed high displacement percentages 75% and 100%, respectively. It was observed that the chloro group at position R5 has a negative impact on the activity, where compound 5g showed a significant loss of activity (the displacement percentage is less than 10%). Compound 5e is also very weakly active (the displacement percentage is 13%) (Table 5). This can be due to the interfering of the substituent e.g. chloro at position R5 with the molecular surfaces of Gln196 and Leu116 (Figure 8). Thus, position R5 cannot be substituted due to the restricted space available nearby.The importance of position R4 was evaluated by using different small groups, which led in most cases to an enhancement of the inhibitory activity. Compound 5l, which contains an additional methyl group at R4 compared to compound 5k, showed a comparable activity (100% displacement compared to 5k). Meanwhile, using a methoxy at R4 had a negative impact e.g. 5m showing a displacement of 39.67% (Table 5).Further, compound 6 containing a methyl group as part of thiolane ring, was employed to assess the influence of a methylene group at position R0 (Figure 10). Compared to 5d the modified compound 6 showed loss of activity exhibiting a displacement percentage less than 10%.

This confirms the negative impact of introducing a second substituent at position R0 of the pyrimidine ring.The IC50 values of the compounds, which showed high displacements, were measured using both an activity (FPIA) and binding (FPBA) assay. The tested compounds showed an inhibitory activity in the low-micromolar and sub-micromolar range (Table 5). Then we tested whether the BFE and the derived QSAR model were able to predict the biological activities. Both models did not show a significant ability to predict the affinities of the novel compounds. However, some of the compounds predicted to be highly active (e.g. 5l and 5n) were experimentally confirmed to be sub-micromolar inhibitors of PKMYT1. Furthermore, both models predicted the negative impact of the chloro group at position R5, where the predicted pIC50 was decreased compared to their parent molecules (5d and 5e: predicted pIC50 = 5.99 and 5.16, 5f and 5g: 5.79 and 5.13,respectively, Table 5). Generally, the MM-PBSA model overestimated the affinity of the novel compounds compared to QSAR model 1A.In order to evaluate the predicted binding mode of the diaminopyrimidines (e.g. compound 5l), we searched for analogues that have been co-crystallized with other kinases. Thus, the database of ligands deposited in the PDB was screened using Swisssimilarity [32]. This resulted in 23 similar compounds showing a similarity score between 0.56 and 0.79. We found that diaminopyrimidine derivatives are found as inhibitors of the following kinases: Aurora A [33], Ephb4 [34], EGFR[35] and CDK2 [36]. A series of diaminopyrimidine analogues have been co-crystallized with Aurora A kinase (PDB IDs 3UNZ, 3UO5, 3UO6, 3UOD, 3UOH, 3UOJ, 3UOK, 3UOL,3UP2, 3UP7, 4DEB, 4DED, 4DEE) [37, 38]. In addition, other analogues were co-crystallized with CDK2 (PDB IDs 3UNJ, 3UNK [33, 37] and 1H08 [36]) and EGFR (PDB ID 5UWD [35]).

In Aurora A, the diaminopyrimidine compounds show the same binding mode as predicted for PKMYT1. This binding mode is characterized by a horseshoe shape conformation and two hydrogen bonds to the hinge region of the ATP pocket, whereas the phenyl ring occupies the ATP-ribose binding site (Figure 11A). Meanwhile, diaminopyrimidine compounds display two distinct conformations co-crystallized with CDK2; a horseshoe shape (PDB ID 1H08 [36]) and a conformation (PDB ID 3UNJ [33]) where the phenyl ring flips towards the back-pocket (Figure 11B). Also the analogues co-crystallized with EGFR (PDB ID 5UWD [35]) show a horseshoe shape binding mode at the ATP pocket.The diaminopyrimidine derivatives containing a substituent on the pyrimidine at position R1 show in the X-ray structures a horseshoe shape binding mode. The diaminopyrimidine analogue co- crystalized with CDK2 (PDB ID 1H08 [36]) contains a bromo group at position R1 and shows the horseshoe-shape binding mode (Figure 12). A second example is the crystal structure of EGFR (PDB ID 5UWD [35]), which contains an analogue with a trifluoromethyl attached at position R1 (Figure 12). Thus, substituting R1 with e.g. halogens, methoxy or sulfur are only found in the horseshoe shape binding mode, which was also predicted for compound 5l that contains a methoxy group at position R1.The binding mode of compound 5l displays the conserved hydrogen bonds with the hinge region residue Cys190 [13], which are critical for the inhibitory activity of the ATP-competitive inhibitors.

The interactions with Tyr121 and Asp251 in addition with the P-loop residue Gly117 (Figure 11C) might be the reason for the increase in potency as proposed by the computational fragment growing approach. The methyl group at R4 of compound 5l forms additional interactions with the P-loop residues, thus, it showed higher activity than compound 5k (Table 5). In contrast, a methoxy group at R4 (5m) decreases the activity, which might be due the steric hindrance with the P-loop residues. Substitution at position R5 by a chloro group decreases the activity of the diaminopyrimidine compounds 5e and 5g. This is due to the interference of the chloro group with Asn196 and Leu116, which is predicted to be unfavourable. Compound 5n containing a methyl amide at the position R2 showed decreased activity compared to 5l. This shows the importance of the interactions with the P-loop residues to boost the inhibitory activity. Substituting R1 using small groups (methoxy group in case of 5l) has a big impact on the inhibitory activity of the diaminopyrimidine derivatives. The substitution is favouring a horse shoe shape binding conformation that fits into the ATP-binding pocket of PKMYT1 (Figure 11C). On the other hand, using more bulky groups (benzyloxy group in 5a) results in loss of inhibitory activity due to the limited space in the back-pocket. Furthermore, a combination of the methyl group at R1 and the hydroxyl group at R3 might destroy the bioactive conformation (horse shoe shape), thus, compound 5b and 5c lose their activities (Table 5). The methylsulfanyl and methoxy groups of compounds 5d and 5l fit perfectly in the PKYMT1 back-pocket and interact with the gatekeeper Thr187. We found also that adding a second substituent on the pyrimidine ring e.g. compound 6 decreases the activity dramatically.

Conclusion
In the current study, we generated several computational models in order to predict the biological activity of novel PKMYT1 inhibitors. We extended our previously developed models by including novel X-ray structures of PKMYT1 and by increasing the ligand data set. We observed a good reproduction of the binding mode of recently reported PKMYT1 inhibitors. However, the docking and scoring failed to correctly predict the affinity of the ligands, where only a weak correlation was obtained between computed and experimental affinities. We were able to show that binding free energy calculations (especially the MM-PB/SA solvation model) were better to predict the activities of the compounds. Moreover, including 2D molecular descriptors based on the partial charges of the ligands slightly improved the quality of the prediction. The applied fragment growing allowed to design a first set of diaminopyrimidine derivatives, which are active in the sub-micromolar range (e.g. compound 5l). As a consequence, we were able to generate several active analogs and to assess first structure-activity relationships of the diaminopyrimidine derivatives. The derived results provide now the basis for further chemical modifications to design more potent PKMYT1 inhibitors.
Materials and Methods

1.Molecular docking
The inhibitors dataset was prepared using the LigPrep module in Schrödinger (LigPrep, Schrödinger, LLC, New York, NY, 2017-1). The preparation involved 3D protonation at pH=7.4 using Epik module [39, 40], minimization using the integrated Optimized Potentials for Liquid Simulations (OPLS_2005) force field [41] and generation of ten low energy conformations per ligand. Further preparation comprised conformer generation with OMEGA in OpenEye Scientific Software [42], where a maximum of twenty conformers for each compound were generated. Preparation of PKMYT1 X-ray structures Nine crystal structures of PKMYT1 are stored in the Protein Data Bank (PBD; www.rcsb.org) [43] (PDB IDs: 3P1A (no associated citation), 5VCV, 5VCW, 5VCX, 5VCY, 5VCZ, 5VD0, 5VD1 and 5VD3) [14]. In our previous work, we found that the structure of PKMYT1 co- crystalized with Pelitinib (PDB ID: 5VCW) was able to re-produce and predict the binding mode of native and non-native ligands [13]. Thus, it was used in the current work to dock the newly designed compounds. Preparation module in Molecular Operation Environment System (MOE) (Chemical Computing Group Inc., Montreal, QC, Canada) was used for preparing the PKMYT1 structures. The process encompassed 3D protonation, followed by energy minimization employing Amber10: Extended Huckel Theory (EHT) force field implemented in MOE [23, 44] with a tethering (0.5) Å and a gradient of 0.1 Kcal.mol-1. Å-2 for all atoms during the minimization. The inhibitors were docked into the prepared PKMYT1 structures using GOLD v5.2 [45]. The co-crystallized ligand (Pelitinib) within the PKMYT1 structures was defined as a center of the binding site. One hydrogen bond formed with the backbone-NH of Cys190 was used as hydrogen bond constraints. In addition, one conserved water molecule, which forms a water bridge with His161 and Glu157, was considered during the docking procedure. Generally, ten poses were generated for each docked ligand. Other docking options were kept as default. The obtained top- ranked docking poses were subjected to re-scoring by different docking scores implemented GOLD and Glide (Schrödinger 2017u1) [46-48].

2.Molecular Dynamics (MD) simulations
Molecular dynamics (MD) simulations were performed prior the calculation of the binding free energy. MD simulations were carried out using AMBER 12 [24] and employing AMBER03 force field [49-51]. Antechamber module was used to generate parameters for atom type and AM1- BCC [52] atomic charges for each ligand in general AMBER force field GAFF [53]. The preparation of the protein-ligand complexes was performed using LEaP module in AMBER. Water model TIP3BOX was used to solvate the system (protein-ligand complex and ions) in a box size of 10 Å [54] and counter ions were add to neutralize the system. Several steps were performed to prepare the system for MD simulations, starting with two consecutive steps of energy minimization. Firstly, minimization was carried out for water and ion molecules while keeping the ligand-protein complexes restrained to their initial coordinates with a force constant of 500 Kcal.mol-1.Å-2. The procedure helped to reduce unreal van der Waals interactions with the surrounding solvent molecules and rearrange of the inserted water and ions molecules. In this step 2000 iterations (beginning with 1000 steepest descent and followed by 1000 conjugate gradient) were performed. The second minimization step was applied for the whole system through 10000 iterations (first 5000 steepest descent and then 5000 conjugate gradients). The next step involved heating. The temperature was equilibrated from 0 to 300 K and occurred gradually over a timescale 2 fs of 100 ps to avoid problems with the hot solvent cold solute. Langevin dynamics was set for temperature control using a collision frequency of 1 ps-1 during the temperature equilibration [55]. The cutoff force constant of the restrained system was set at 10 Kcal.mol-1. Å-2. The output coordinates after the heating were inserted in a pressure equilibration routine. The third step was called density, which was used to set the constant pressure periodic boundary. The constant pressure periodic boundary was converted from 1 bar to 2 bar during 100 ps with a timescale of 2 fs at 300 K for the temperature and applying the same cutoff of force constant at 10 Kcal.mol-1. Å-2. During the MD simulation, a cutoff of 10 Å for non-bonded atoms was applied using Particle Mesh Ewald methods [56]. Moreover, all simulation were run using SHAKE to constrain hydrogen bonds [57]. Finally, 100 ps MD simulation was performed.

3.Binding free energy (BFE) calculations
The binding free energy calculations were performed for the putative binding modes described in the previous work [13]. AMBERTools12 [24] and MOE (Chemical Computing Group Inc., Montreal, QC, Canada) were used to calculate the binding free energies. The calculation was first carried out using one snapshot after a short energy minimization of the whole system in explicit water. Later, ten snapshots (frames) of ligand-protein complexes from MD simulation were considered for calculating the binding free energy. Different theories were applied using Generalized Born (GB) [58, 59] and Poison Boltzmann (PB) solvation models [60]. The molecular mechanics (MM) and quantum mechanics/molecular mechanics (QM/MM) hybrid were considered too [61, 62]. The QM region of the QM/MM hybrid mechanics was firstly applied for the ligand only and later for the ligand and selected residues from the PKMYT1 binding pocket, surrounding the ligand (In total 10 residues were considered: Leu116, Lys139, His161, Val171, Leu185, Thr187, Leu189, Cys190, Gln196, Asp251). Two models for QM were used: Parameterized Model number 3 (PM3) and Austin Model 1 (AM1) [24]. Finally, the binding free energy calculation of MM-PB(GB)/SA and MM/QM-GB/SA models were computed using MMPBSA.py [63] module in AMBERTools12[64-69]. In addition, the implemented MM-GB/SA model in MOE was used to estimate the binding free energy. The computed affinity of a ligand to a protein can be estimated by differences between the free energy of complex and the sum of the free energies of ligand and protein separately, eq. (2) [70]: ∆Gbinding = Gprotein/ligand – (Gprotein + Gligand) (2)
The absolute free energy for each part can be computed as a sum of the gas-phase free energy EMM plus the changing of the free energy because of solvation ∆Gsolvation and the entropic contributions -T∆S, eq. (3) [71]: Gmolecule = EMM + ∆Gsolvation – T∆S (3) The conformational entropy value -T∆S was not considered. Thus, the change of enthalpy (∆H) was considered to calculate the binding free energy and it will indicate to the binding free energy, eq. (4). ∆H = EMM + ∆Gsolvation (4)

4.Fragment-based drug design
Fragment-based drug design (FBDD) shows a high ability to sample the chemical space rather than other techniques [72]. Fragments are defined by Astex through ‘rule-of-3’ [73]. Fragment growing approach provides useful tools for the hit and lead optimization, which aims to build additional interactions with the adjacent pockets [74, 75]. Two metrics, namely ligand efficiency (LE) and ligand lipophilicity efficiency (LLE), are available to aid in the fragment hit selection and optimization. The ligand efficiency (LE) is used to rank fragments and to monitor the progress of optimization, whereas ligand lipophilic efficiency (LLE) is a metric used to monitor the lipophilicity with respect to in vitro potency of a molecule [26, 27, 76]. Moreover, the ligand lipophilicity efficiency does not consider the size of the ligand. During the optimization studies, the goal is to improve potency without increasing lipophilicity. However, the ideal values of LLE are located in the range of 5 to 7 or greater [28, 77, 78]. The following equations were used to estimate the ligand efficiency and ligand lipophilicity efficiency: LE = 1.4(-logKi)/HAC, HAC is heavy atoms count. LLE = pKi – cLogP, cLogP was calculated using MOE.

5.Biological screening of the compound collection
Aside from kinase inhibitors, all chemicals used were purchased from Sigma (Schnelldorf, Germany) unless stated otherwise. PD-0166285 was from Santa Cruz Biotechnology (Santa Cruz, CA, USA), PD-173952 and PD-180970 from Sigma (Schnelldorf, Germany) and MK-1775 from Selleckchem (Houston, TX, USA). All other tested kinase inhibitors were purchased from LC Laboratories (Woburn, MA, USA). The GSK published kinase inhibitor set (PKIS) I and II was kindly provided by Glaxo Smith Kline LCC (Research Triangle Park, NC USA). The monoclonal p-Tyr P-100 antibody was from Cell Signaling (#9411, Danvers, MA, USA). PKMYT1 full-length (FL) and PKMYT1 kinase domain (KD) were expressed and purified as described before [79]. (6-FAM)KI(pY)VV was from IKFZ (Leipzig, Germany). EFS247-259 and pY–EFS247-259 were synthesized via solid-phase peptide synthesis on Fmoc-Rink MBHA as basically described in the work of Rohe et al. [80] following a standard Fmoc/tBu based strategy. Processing and purification was performed as described in Rohe et al.
To test the activities of the compounds as PKMYT1 inhibitors, the synthesized compounds were screened with a fluorescence anisotropy-based binding assay, using the PKMYT1 domain (KD) and a fluorescently labeled ATP-competitive compound [15].

The initial screening was performed at 20 µM inhibitor concentration as triple measurement of each compound. Dasatinib (10 µM) was used as inhibitory control and 1% DMSO as vehicle control. Compounds with a tracer displacement >60% were retested at 20 µM and 5 µM concentrations in triple measurement. The IC50 and Ki values were determined for confirmed hits using the FP binding assay. To confirm the inhibitory activities of the hit compounds, we further used the newly developed activity assay (FPIA) [81], and determined the IC50 and Ki values for PKMYT1 (FL) inhibition. Similar to the FP-binding assay, 10 µM dasatinib and 1% DMSO were used as controls. The data processing was performed as described before. Because the FP does not behave additively, anisotropy (r) was calculated according to Rohe et al. [79, 80]. The ∆r was calculated as the difference between r of positive control (PKMYT1 reaction) and negative control (reaction without PKMYT1). Data is represented as means ± SEM. Kinetics, binding and inhibition curves as well as linear and nonlinear regression were generated using GraphPad Prism 5 (GraphPad Software, Inc., La Jolla, CA, USA).

6.Synthesis and characterization of novel compounds
General: All materials and reagents were purchased from Sigma–Aldrich Co. Ltd., abcr GmbH, ChemPUR Feinchemikalien und Forschungsbedarf GmbH, and Carbolution Chemicals. All solvents were analytically pure and dried before use. Synthesis was carried out under microwave conditions using Monowave 450 (Anton Paar, Graz, Austria). Thin-layer chromatography (TLC) was carried out on aluminum sheets coated with silica gel 60 F254 (Merck, Darmstadt, Germany). Purification of the final compounds was carried out with preparative HPLC (Shimadzu, Kyoto, Japan; LC-10AD, SIL-HAT auto sampler). Therefore a 7,8 x 300 mm XTerra RP-18 column (7 µM), Waters (Milford, MA, USA) with a UV-Vis-detector SPD-M10A VP PDA was used. Separation was achieved with a gradient of MeOH/H2O (starting at 95% H2O going to 5% H2O). Final compounds were confirmed to be of >95% purity based on HPLC. Purity was measured by UV absorbance at 254 nm. HPLC instrumentation consisted of an XTerra RP18 column (3.5 mm 3.9_100 mm; Waters, Milford, MA, USA) two LC-10AD pumps, an SPD-M10A VP PDA detector, and an SIL-HT auto sampler all from the manufacturer Shimadzu (Kyoto, Japan). The mobile phase was in all cases a gradient of MeOH/H2O (starting at 95% H2O going to 5% H2O). Mass spectrometry analyses were performed with a Finnigan MAT 710C (Thermo Separation Products, San Jose, CA, USA) for ESI-MS spectra, and with a LTQ (linear ion trap)-Orbitrap XL hybrid mass spectrometer (Thermo Fisher Scientific, Bremen, Germany) for HRMS-ESI (high- resolution mass spectrometry) spectra. For HRMS analyses, the signal for the RP-6306 isotopes with the highest prevalence was given and calculated for 35Cl and 79Br. 1H NMR spectra were taken on a Varian Gemini 2000 and a Varian Inova 500 using CDCl3 and [D6] DMSO as solvents. Chemical shifts (d, ppm) are referenced to the residual solvent signals.

General procedure for the conversion of the 2,4-dichloropyrimidine derivatives 1 with corresponding aniline derivatives 2: The 2,4-dichloropyrimidine derivative 1 was dissolved in EtOH and 1.1 equiv. of the aniline derivative 2 and 1.1 equiv. of DIPEA were added. The reaction was carried out for 50 min at 150 °C under microwave conditions. The reaction mixture was diluted with ethyl acetate and washed with water. The organic phase was evaporated and the crude intermediate was used for the following reaction step without purification. General procedure for conversion of the 2-chloro-4-anilinopyrimidine derivatives 3 with 4- morpholinoaniline 4: The crude 2-chloro-4-anilinopyrimidine derivative 3 was dissolved in 1-BuOH and 1.1 equiv. of 4-morpholinoaniline 4 and 1.1equiv. of 1 M HCl were added. The reaction was carried out for 20 min at 160 °C under microwave conditions. The reaction mixture was diluted with ethyl acetate and washed with 1 M aqueous NaOH. The mixture was evaporated and purified with preparative HPLC.