Rational Drug Design of Axl Tyrosine Kinase Type I Inhibitors as Promising Candidates Against Cancer

The high level of Axl tyrosine kinase expression in various cancer cell lines makes it an attractive target for the development of anti-cancer drugs. In this study, we carried out several sets of in silico screening for the ATP-competitive Axl kinase inhibitors based on different molecular docking protocols. The best drug-like candidates were identified, after parental structure modifications, by their highest affinity to the target protein. We found that our newly designed compound R5, a derivative of the R428 patented analog, is the most promising inhibitor of the Axl kinase according to the three molecular docking algorithms applied in the study. The molecular docking results are in agreement with the molecular dynamics simulations using the MM-PBSA/GBSA implicit solvation models, which confirm the high affinity of R5 toward the protein receptor. Additionally, the selectivity test against other kinases also reveals a high affinity of R5 toward ABL1 and Tyro3 kinases, emphasizing its promising potential for the treatment of malignant tumors.


INTRODUCTION
Receptor tyrosine kinases are transmembrane proteins, which consist of several domains that are activated upon ligand binding to their extracellular regions, triggering downstream signaling cascades (Robinson et al., 2000;Myers et al., 2016). They are involved in various regulatory processes, such as cell survival, growth, differentiation, adhesion, proliferation, and motility (Robinson et al., 2000;Ségaliny et al., 2015;Myers et al., 2016). Impaired gene functions by mutations or deletions may cause the abnormal expression of protein kinases, which, in turn, entails tumor formation and progression (Blume-Jensen and Hunter, 2001;Zhang et al., 2008).
One of the frequently identified kinases involved in the formation of various types of tumors is Axl receptor tyrosine kinase (Craven et al., 1995;Sun et al., 2003). Axl belongs to the TAM family receptors, which also includes Tyro3 and Mer (O'Bryan et al., 1991;Li et al., 2009). The kinase structure comprises an extracellular part with two immunoglobulin (Ig)-like domains responsible for ligand binding, a transmembrane region, and an intracellular domain (O'Bryan et al., 1991;Lemke and Rothlin, 2008). The growth arrest-specific 6 (Gas6) protein precursor and protein S are primarily responsible for kinase activation as their ligands (Stitt et al., 1995;Varnum et al., 1995;Li et al., 2009). Both ligands share a similar domain composition. In particular, they include two sex-hormone-binding globulin domains at the C-terminus, both with the laminin G1 and G2 proteins necessary for the subsequent binding to the Ig-like domain of the receptor, causing their dimerization and activation (Lemke and Rothlin, 2008). Close to the N-terminal, there are epidermal-growth-factor-like repeats and, the so-called, Gladomain that consists of gamma-carboxyglutamic acid, which is necessary for binding to phosphatidylserine of the apoptotic cell membrane in a vitamin-K-dependent reaction (Hasanbasic et al., 2005;Sasaki et al., 2006;Li et al., 2009).
Axl overexpression has been detected in a majority of human cancers, including acute myeloid leukemia (Rochlitz et al., 1999;Hong et al., 2008), breast cancer (Berclaz et al., 2001;Zhang et al., 2008;Gjerdrum et al., 2010), gastric (Wu et al., 2002) and lung cancer (Shieh et al., 2005), melanoma (Quong et al., 1994), osteosarcoma (Han et al., 2013), renal cell carcinoma (Gustafsson et al., 2009), etc. Therefore, targeting the Axl to inhibit its function might be a promising strategy for the treatment of various malignant tumors. Different strategies of targeting the Axl have already been considered. For instance, Rankin and Giaccia (2016), in their review, highlight the three classes of Axl inhibitors directed on cancer therapy. The first class includes small-molecule tyrosine kinase inhibitors that block Axl kinase activity (Rankin and Giaccia, 2016). The second class consists of anti-Axl antibodies (Rankin and Giaccia, 2016) that block Axl activation, which is triggered by the Axl-Gas6 interaction, and the third class comprises soluble Axl decoy receptors (Rankin and Giaccia, 2016) that serve as a trap for Gas6, hence, preventing the Axl-Gas6 binding.
Different experimental and computational techniques have been developed and applied in the last decades for rational drug design and discovery (Baldi, 2010;Ou-Yang et al., 2012;March-Vila et al., 2017). For instance, computational and experimental approaches focused on in silico design and organic synthesis of the Axl kinase inhibitors have already been performed by Mollard et al. (2011). In their research, the authors constructed a homology model for the active site of the Axl kinase and performed docking experiments for the designed compounds. Recently, the three-dimensional (3D) structure of the Axl kinase in a complex with its inhibitor (macrocyclic compound 1) has been successfully solved by Gajiwala et al. (2017) using differential scanning fluorimetry and hydrogen-deuterium exchange mass spectrometry. This 3D structure, as a tetrameric configuration, consists of two active (B and D chains) and two inactive (A and C) motifs in a complex with a small ATPcompetitive inhibitor. The active and the inactive states are characterized by the DFG (Asp-Phe-Gly) loop-in and loopout conformations.
According to the mode of binding, all tyrosine kinase inhibitors have been divided into different types. In their review, Zhang et al. (2009) distinguishes four basic types of inhibitors. According to this classification, the type I and the type II inhibitors bind to the DFG-in and DFG-out motifs, respectively. Additionally, the type III inhibitors interact with the protein outside the highly conserved ATP-binding pocket, representing allosteric binding and, therefore, named as allosteric inhibitors (Wu et al., 2015). Finally, the type IV inhibitors bind the active site irreversibly, forming covalent bonds within the binding pocket. A slightly expanded classification of the protein kinase inhibitors is suggested by Roskoski (2016).
In the current study, we introduce systematic computational analysis for the DFG-in conformation of the Axl kinase to inhibit its activity using different sets of molecular docking algorithms and molecular dynamics simulation techniques to handle cancerrelated diseases.

3D Structure Derivation
The 3D coordinates for the Axl kinase in a complex with the macrocyclic compound 1 have been retrieved from protein data bank (PDB)

GOLD
Molecular docking using the GOLD software was performed using version 5.5 (Jones et al., 1997). The binding site residues were defined by specifying the crystal structure ligand coordinates bound to the protein and using the default cutoff radius of 6 Å, with the "detect cavity" option enabled. The GOLD docking experiments were performed using the ChemPLP scoring function. For each compound, 50 complexes were generated. The highest-scoring compounds were selected as the most appropriate ones.

MOE
Docking has been performed by selecting the default "Rigid Receptor" protocol. As a binding site, the coordinates of co-crystalized ligand atoms have been selected. The ligand placement was performed using the Triangle Matcher protocol. The top 30 poses were ranked by London dG scoring function, FIGURE 1 | The superposition of the crystal structure of macrocyclic compound 1 (shown in forest green) with its docked pose performed by the GOLD (A, orange), MOE (B, deep pink), and AutoDock (C, yellow) software. The protein is represented as a secondary structure and highlighted in cyan. Some of the pocket residues are shown in licorice and colored in purple. and the resulting five poses were identified using the generalized-Born-volume-integral/weighted-surface-area function. The conformations with the more negative final score were considered as favorable.

AutoDock
Molecular docking using AutoDock software (Morris et al., 2009) has been performed using version 4.2.6 (available at: http://autodock.scripps.edu). The AutoDock tools were used to generate the input parameter files for docking. For the current study, the receptor was considered as a rigid molecule, while the ligands contained rotatable bonds. Pure protein was applied for the docking, while all non-protein moieties were removed. Additional hydrogen atoms were added to the receptor, and the new PDB coordinates were saved. The ligand PDB file was modified by the addition of groups representing the Gasteiger charges. The volume of the grid box was set as 50 × 50 × 50 Å, with 0.375-Å spacing. The center of the grid box was placed so that it coincided with the center of the co-crystalized structure of the compound (macrocyclic compound 1). A genetic algorithm was selected to set the search parameters. The number of docking runs was fixed to 50. The conformations with the lowest binding energies have been selected for further analysis.
The figures were prepared using the UCSF Chimera (Pettersen et al., 2004;

Molecular Dynamics
All molecular dynamics (MD) simulations were performed using the AMBER 16 package (Case et al., 2005) with the FF99SB and GAFF force fields for the Axl protein and its ligands. The systems were solvated with the TIP3P water models and neutralized by adding the Na + ions using the tLEap input script available from the AmberTools package. Long-range electrostatic interactions were modeled via the particle-mesh Ewald method (Essmann et al., 1995). The SHAKE algorithm (Miyamoto and Kollman, 1992) was applied to constrain the length of covalent bonds, including the hydrogen atoms. Langevin thermostat was implemented to equilibrate the temperature of the system at 300 K. A 2.0-fs time step was used in all of the MD setups. For the minimization and equilibration (NVT ensemble) phases, 10,000 steps and 1-ns time period were used, respectively. Finally, 50-ns classical MD simulations, with no constrains as NPT ensemble, were performed for each of the protein-ligand complexes using the molecular mechanics combined with the Poisson-Boltzmann (MM-PBSA) or generalized Born (MM-GBSA) augmented with the hydrophobic solvent-accessible surface area term (Kollman et al., 2000;Shityakov et al., 2017). The MM-PBSA/GBSA solvation models were applied as a post-processing end-state method to calculate the free energies of molecules in the solution by means of the optimized python script (MM-PBSA.py).

Validation of the Binding Poses
To validate the poses of the ligands, we performed docking of the co-crystalized ligand-macrocyclic compound 1 (Gajiwala et al., 2017) to the ATP-binding pocket of the Axl kinase using the GOLD, MOE and AutoDock simulation software. Figure 1 demonstrates the superposition of the crystal structure of the ligand (always shown in forest green) and the docked pose of the same ligand performed by the GOLD (orange, Figure 1A), MOE (deep pink, Figure 1B), and AutoDock (yellow, Figure 1C) software. The calculated RMSD values for this ligand are below the commonly accepted threshold of 2.0 Å, indicating the validity of the above-mentioned docking engines for the prediction of the ligand binding pose. The most accurate results were reproduced by the GOLD and AutoDock software, with the calculated RMSD values of 0.2 and 0.5 Å, respectively (Figures 1A,C). A less accurate result was shown by the MOE software, where the RMSD between the crystal structure and the docked pose was 1 Å (Figure 1B).

Docking of ATP-Competitive Type I Inhibitors
The first round of in silico testing was performed on the commercially available ATP-competitive Axl kinase type I inhibitors using the molecular docking protocols. From Table 1, it is clear that the scoring functions from the molecular docking results indicate a better affinity for the analyzed drugs in comparison to the ATP, which was used here as a reference. The top scores belong to macrocyclic compound 1, crizotinib, and R428 according to the three independent molecular docking protocols. The last two substances were chosen as the core components for the structure similarity search and in silico chemical modifications based on their ability to block specific kinases that are highly expressed in malignant cells (Solomon et al., 2014;Chen et al., 2018).

Molecular Docking of R428-Based Analogs
Next, the drug R428 was subjected to the structure similarity search using the PubChem search algorithm (Kim et al., 2016). The results came up with 136 patented analogs, which are suitable for the GOLD/MOE molecular docking to determine the top 10 hit molecules ( Table 2). Within this group the best four compounds (  Figure 3A) due to the bigger number of rings in The CIDs of the compounds corresponding to the best scores are shown in italics.
the molecule, thus reducing its flexibility. The RMSD value in this case is 5.7 Å.

Molecular Docking of Crizotinib-Based Analogs
A similar structural search was applied to crizotinib to ensure the identification of the crizotinib-based analogs. The search results provided 26 analogs for further GOLD/MOE screening. Finally, the top five compounds were determined ( Table 4), but none of them was subjected to in silico chemical modifications due to their lesser binding affinity than crizotinib itself. Therefore, the parental structure was modified to improve its binding properties toward the Axl kinase.

R428_1 Modifications
To enhance the binding affinity of the identified top compounds (Table 3), we devised the structural modification scheme for R428_1. Its invariant part was estimated from the best proteinligand binding mode, forming the "scaffold" to adjust the ligand conformation inside the binding pocket (Figure 4). Therefore, we decided to design the new compounds by implementing this part as a template and adding the molecular extensions (X, Y, and Z) into the "scaffold" at the locations indicated by the arrows. So, the first six compounds (R1-R6) were designed: the compound R1 has been derived by adding piperidine bound to a tri-cyclic moiety (see Figure 4), the compound R2 was formed by adding a triazole-like ring also connected with a tri-cyclic moiety, the compound R3 was derived by adding a histidinebound triazole ring, the compound R4 was formed by adding two repeats of triazole-like rings, the compound R5 was achieved by extending the template structure with a tyrosine-bound triazolelike ring, and the compound R6 has been derived by adding to the template an NH 2 group. The next four compounds (R7-R10) have been derived from the above-mentioned designed R6 and R3 compounds with a slight alteration of the template structure. Thus, the compound R7 was obtained by replacing a hydrogen  Frontiers in Chemistry | www.frontiersin.org 5 February 2020 | Volume 7 | Article 920 The superposition of the top four out of best 10 R428 analogs according to the MOE and GOLD docking results. The compound R428_1 is shown in red, the compound R428_2 is shown in orange, the compound R428_3 is shown in yellow, and the compound R428_4 is shown in green. The backbone of the kinase is represented as a cartoon and shown in cyan; the key residues inside the pocket are highlighted in purple. The hydrogen bond between the R428_2, R428_3, and R428_4 compounds and D690 is marked as a blue line. (B) Superposition of the R428_2, R428_3, and R428_4 compounds for a better representation of their conformational coincidence. atom with a hydroxyl group in the compound R6. The compound R8, in turn, was obtained by replacing an amino group with a nitroso group in the compound R6. The compound R9 represents a pure template (see the upper-left part of Figure 4). Finally, the compound R10 has been derived by replacing an aromatic amino group with a nitroso group in the compound R3.
The newly designed compounds were tested by docking to the ATP-binding site of the Axl kinase. Improved binding results according to the GOLD, MOE, and AutoDock software were obtained for the designed R3, R5, and R10 compounds (see Table 5). In particular, they score higher than the top R428 derivatives (see Table 2). The superposition of R3, R5, and R10 inside the ATP-binding pocket of the Axl kinase shows a similar level of spatial occupancy (see Figure 5). However, shifts as well as differences in orientations of the tri-cyclic rings are observable. The compound R3 (Figure 5, deep pink) is oriented so that it obtains a "horseshoe-like" conformation inside the binding pocket and points its tri-cyclic moiety toward N677, the compound R5 (Figure 5, orange) obtains a slightly extended conformation and exposes its tri-cyclic ring to D690, while the compound R10 (Figure 5, olive green) is oriented opposite from the R3 and the R5 direction, pointing an amino group toward D690.
All of the newly designed compounds, except R4 and R8, form hydrogen bonds (H-bonds) with D690 (see Table 6). The compounds R1-R3 and R10 form an H-bond with N677, while the compounds R2-R4, R6, and R8 form an H-bond with M623. Additionally, the compound R2 forms an H-bond with G626, the compound R3 forms an H-bond with K624, the compound R4 forms an H-bond with A562, the compound R5 forms an H-bond with H625, and the compound R8 forms an H-bond with P621. Van der Waals interactions have been observed as well between the pocket residues and the designed compounds. So, the compounds R1-R5, R7, R9, and R10, besides the above-mentioned residues, additionally make interactions with L542. The compounds R6 and R8, in turn, make interactions with A565.

R428_2, R428_3, and R428_4 Modifications
Next, we tried to obtain new compounds with improved binding based on the R428_2, R428_3, and R428_4 structural characteristics. Since these compounds share an almost identical structure and acquire a quite similar conformation inside the pocket (see Figures 2, 3B), we selected a common structural feature that is shared among these compounds (see the upperleft part in Figure 6) and used it as a template for further chemical modifications. Figure 6 describes a stepwise process of the design of the new (R1 ′ -R11 ′ ) compounds. So, the compound R1 ′ has been derived by the addition of but-2-enylazanium to the part indicated with an arrow (upper-left part of Figure 6). The compound R2 ′ has been derived by adding an amino group, the compound R3 ′ has been formed by adding an aromatic ring, the compound R4 ′ has been derived by replacing a single aromatic ring (indicated with a dashed green circle) in a template by double aromatic rings, and the compound R5 ′ has been formed by the replacement of two aromatic carbons with double-bonded nitrogen atoms in the structure of R4 ′ . The compound R6 ′ has been obtained by adding an extended aliphatic chain to the template, and the compound R7 ′ has been derived by adding FIGURE 4 | The design strategy of the new compounds starting from the R428_1 compound (see Table 3). The upper-left part is an invariant part of the compound used as a template. The upper/lower right and lower left parts are extensions/modifications of the template.
a chain with an arginine-like termination. The compound R8 ′ has been derived similarly to compound R7 ′ , however with a carbonyl termination. The compound R9 ′ has been formed by adding to the aliphatic extension an aromatic ring containing a hydroxyl group, the compound R10 ′ has been obtained similarly to the compound R9 ′ , however with the absence of a hydroxyl group, and, finally, the compound R11 ′ contains an extension with a tryptophan-like termination. Table S1 reports the docking scores for the above-mentioned modifications. The compounds did not show significant improvements of the binding scores compared to the original compounds-R428_2, R428_3, and R428_4 (please refer to Table 2 and Table S1). However, improvements were noticed for the designed compounds R9 ′ and R10 ′ . For both compounds, the scores were ∼85 (by GOLD) and −8 (by MOE) (see Table S1). We, therefore, did not perform a further analysis on these compounds as they did not demonstrate binding affinities that are stronger than the ones of R428_1 modifications (please refer to Table 5). The strategy of the design is based on the following scheme: each consecutive compound is derived from a previous one via the modification of a single group present in the structure. The parts that underwent modifications are indicated with dashed circles. Thus, the compound C1 is derived by the replacement of a fluorine atom, indicated with a violet dashed circle, with a hydroxyl group; the compound C2 is derived from compound C1 by the replacement of an amino group, indicated with a red dashed circle, with a hydroxyl group; the compound C3 is derived from compound C2, where the pyrazole ring, indicated with a dashed green circle, is replaced by a pyran ring; the compound C4 is derived from the compound C3 by the replacement of a chlorine atom, indicated by a purple ring, with a hydroxyl group; the compound C5 is derived from the compound C4 by the replacement of an aminocyclohexane group, indicated with a dashed magenta circle, with a cyclopropane group; the compound C6 is derived from the compound C5 by the replacement of a cyclopropane group with cyclooctane; the compound C7 is derived from the compound C6 by the replacement of cyclooctane with 3-hydroxy-5-amino-pentane, the compound C8 is derived from the compound C7 by the replacement of two cyclohexanes, indicated with a dashed cyan circle, with a triple aromatic network; the compound C9 is derived from the compound C8 by adding to a triple aromatic network one more hydroxyl group, and, finally, the compound C10 is derived from the compound C9 by double modifications: first, by the replacement of a hydroxyl group in a triple aromatic network with fluorine atoms and, second, by the replacement of 3-hydroxy-5-amino-pentane with 1-cyclopenta-1,3-dien-1ylbutan-2-ol. Table S2 shows the docking results for the crizotinib modifications-compounds C1-C10. The possible alterations in the structure of crizotinib did not result in a significant binding improvement. Therefore, we did not consider these new derivatives for further analysis.

Validation of Docking Results
To estimate the strength of our best-designed compound R5, we have compared the calculated inhibition constants (K i ) as well as The top four scores are highlighted in italics.
Frontiers in Chemistry | www.frontiersin.org 8 February 2020 | Volume 7 | Article 920 The compound R3 is shown in deep pink, the compound R5 is shown in orange, and the compound R10 is shown in olive green. The backbone of the kinase is represented as a secondary structure and highlighted in cyan; some of the key residues in the binding pocket-D690, N677, and L542-are shown as licorice and indicated in purple. the ligand efficiencies (LE) for ATP, R428, R428_1, and R5. The data are reported in Table 7. As the results show, the highest value of inhibition constant belong to ATP (566 µM). R428 shows the value of K i to be 12 times lower than that of ATP, indicating a highly competitive binding of the R428 over ATP. The compound R428_1 exhibits an even lower K i than the R428 itself, and, finally, the newly designed compound R5 shows the lowest value of K i (2.3 nM), which is 19 times lower than that of the R428. The calculated LE for R428, R428_1, and R5 compounds are showing almost the same values; the LE for the ATP is slightly lower than for the above-mentioned compounds. Taken together, these data indicate the high efficiency of the newly designed compound R5.
To validate further the molecular docking results, the free energies of binding ( G bind ) based on the implicit solvation models were calculated for R428, its best derivative-compound R428_1, the best-designed compound R5 and ATP as a reference molecule, establishing the best ligand affinity to the Axl kinase. The MM-PBSA/GBSA calculations (Tables 8, 9), using the 50 ns MD trajectories, confirm the previous data: the compound R5 has a much higher affinity to the Axl protein in MM-PBSA (Table 8). However, the MM-GBSA test provides the same affinity to Axl for this structure as the R248 derivative form ( Table 9). This outcome might be explained in such a way that the MM-PBSA model is the more "optimized, " providing the data, which are more in agreement with molecular docking experiments. This "hysteresis" phenomenon has been previously observed in the MD experiments for cyclodextrin-based complexes to assess their drug delivery potential at the blood-brain barrier (Shityakov et al., 2016a,b).
On the other hand, the entropy-enthalpy compensation analysis revealed significant differences only in the Axl-ATP thermodynamics, being either an exothermic (MM-GBSA) or an endothermic (MM-PBSA) binding reaction (Table 10). Nevertheless, the experimental binding enthalpy from the previous study was found to be negative for nucleotide binding to creatine kinase (Forstner et al., 1999), indicating the possibility of an exothermic binding event for Axl-ATP. Similarly, the R428-based compounds exhibited more negative H values, suggesting the enthalpy-driven binding process.

Selectivity Test
To establish how selective the best-designed compound R5 toward Axl kinase is, we further performed a selectivity test against a set of other kinases such as ALK5 (Gellibert et al., 2009), ABL1 (Pemovska et al., 2015), FYN (Kinoshita et al., 2006), JAK1 (Siu et al., 2017), MET (Porter et al., 2009), Tyro3 (Powell et al., 2012), and Mer (Gajiwala et al., 2017). Table 11 reports the docking results for the compound R5 to the kinases according to the three different software-GOLD, MOE, and AutoDock. The best binding scores are highlighted in italics (see Table 11). According to the docking results, the compound R5 has shown quite good binding scores for the studied kinases, which means that the compound R5 can be further applied to target the above-mentioned kinases. Among the studied kinases, improved binding has been observed toward ABL1 and Tyro3. Interestingly, the docking of the compound R5 to the Tyro3 binding site resulted in the highest score among all kinases considered, including the primary target-the Axl kinase. This result indicates a high selectivity of the compound R5 not only toward Axl (please refer to Table 5) but also to another TAM family member-the Tyro3 kinase.

DISCUSSION
The expression of Axl kinase was found in various types of tumor cells (see Table 12), making this protein a promising target for novel anti-cancer pharmaceuticals. Despite the fact that there are some Axl inhibitors in the market (Myers et al., 2016), the selectivity of drugs might be a serious obstacle to overcome in  order to improve their efficacy. Therefore, the novel drug-like molecules with an improved selectivity to a certain type of kinase might potentiate the therapeutic effect against malignant tumors.
To follow this path, molecular docking and molecular dynamics simulation approaches were devised and implemented to compare the binding affinities of known Axl kinase  inhibitors (Myers et al., 2016) and to design some novel drug candidates. We purposely did not implement here a completely new simulation strategy, but rather followed a wellvalidated modeling path, which has been successfully applied by experimental scientists in the design of new active compounds Mollard et al., 2011;Dakshanamurthy et al., 2012;Heifetz et al., 2012). It is worth to note that docking, despite its limitations, is a well-established and an experimentally approved technique for the prediction of potential active substances. For The best binding scores are indicated in italics.

Cancer type References
Acute myeloid leukemia Rochlitz et al., 1999;Hong et al., 2008Astrocytoma Vajkoczy et al., 2006Keating et al., 2010 Breast cancer Berclaz et al., 2001;Meric et al., 2002;Gjerdrum et al., 2010Colorectal Craven et al., 1995Dunne et al., 2014 Esophageal adenocarcinoma Hector et al., 2010;Hsieh et al., 2016 Glioblastoma multiforme Hutterer et al., 2008 Gastrointestinal stromal tumor Mahadevan et al., 2007 Gastric cancer Wu et al., 2002 Head and neck squamous cell carcinoma Ovarian adenocarcinoma Chen et al., 2013;Rea et al., 2015 Pancreatic ductal adenocarcinoma Koorstra et al., 2009 Renal cell carcinoma Gustafsson et al., 2009 Pleural mesothelioma Pinato et al., 2013 Urothelial carcinoma Hattori et al., 2016 Prostate cancer Sainaghi et al., 2005 Thyroid cancer Ito et al., 1999 Uterine endometrial cancer Sun et al., 2003 instance, Li et al. (2011) used a docking method to test 4,621 approved drugs from DrugBank against the crystal structure of MAPK14 to identify a potential anti-inflammatory drug for the treatment of chronic myeloid leukemia. The study revealed a potent inhibitor-the drug nilotinib, with an in vitro IC 50 of 40 nM . Dakshanamurthy et al. (2012) showed another successful application of computational modeling in drug discovery. In their methodology, the drug-target interaction among 3,671 FDA-approved drugs and 2,335 human protein crystal structures was predicted with 91% accuracy. In addition, Dakshanamurthy et al. discovered that the anti-parasitic drug mebendazole also revealed anti-cancer properties, with a strong inhibition activity against vascular endothelial growth factor receptor 2. The results were confirmed experimentally to the extent that effective treatment of different medulloblastoma models could be shown by applying mebendazole, including a clear impact on tumor angiogenesis (Bai et al., 2015).
Purely computational studies targeted on the understanding of agonists/antagonist interaction with a certain type of receptor for the design of new medications were performed recently. In particular, Zhu et al. (2019) applied a combination of in silico tools such as the 3D-QSAR, molecular docking, molecular dynamics, and free energy calculation to clarify a selectivity mechanism of glycogen synthase kinase 3β toward ATP-competitive inhibitors. They identified some key selective residues that might play an important role in the design of the novel ATP-competitive inhibitors (Zhu et al., 2019). Martínez-Campos et al. (2019) in their work employed molecular docking: firstly the QSAR method to analyze all known gamma aminobutyric acid (GABA B ) receptor agonists and later a structure-based drug design strategy for the design of the novel compounds. They came up with six potent baclofen analogs targeted on the activation of the GABA B receptor (Martínez-Campos et al., 2019). Liu et al. (2019) used virtual screening, molecular docking, and molecular dynamics tools to identify the potential methylguanine-DNA methyltransferase (MGMT) inhibitors. Their research resulted in two potent leads, the ZINC000008220033 and ZINC000001529323 compounds, which can be further optimized against the MGMT protein .
There are existing computational studies targeted on the development of the Axl kinase inhibitors. For instance, Fatima et al. (2017) performed a docking simulation to understand the compatibility of curcumin derivatives against a homology model of the Axl kinase active site. Similarly, Mollard et al. (2011) used an in silico approach to design the Axl kinase inhibitors. They also performed a homology model of the binding site for such a purpose. Moreover, these results were subsequently confirmed by experimental data (Mollard et al., 2011). In both of the above-mentioned cases, a homology model of the protein was employed. In our study, we rely instead on the experimentally solved 3D structure of the Axl kinase binding pocket (Gajiwala et al., 2017), and for more accuracy, besides the extensive molecular docking calculations, we additionally performed the MMPBSA/GBSA simulations (see Wang et al., 2019 for evaluation) to identify the best possible drug candidates. As the highest affinity toward Axl kinase was defined for the crizotinib and R428 molecules (see Table 1), they were further considered as the parental "scaffolds" for in silico structural modifications.
Crizotinib was already applied as a first-line therapy for the treatment of lung cancer (Awad and Shaw, 2014). This medication is also active against the ALK and hepatocyte growth factor receptor as proto-oncogene c-Met, especially in patients with the ALK-rearranged non-small cell lung cancer after oral administration (Awad and Shaw, 2014). However, most of the crizotinib-based modifications represented only a slight binding improvement to the Axl protein (Table S2; compounds C1-C10). Crizotinib and the related compounds, in line with this, work only for some time, and then there is a resistance observed in the clinic (Awad and Shaw, 2014).
On the other hand, R428 had proven to be a selective Axl kinase inhibitor (Holland et al., 2010), with the inhibition constant in the nanomolar (nM) range (Myers et al., 2016), as well as the blocker of other Axl-associated events such as autophosphorylation, metastasis development in breast cancer and proinflammatory cytokine production (Holland et al., 2010). In addition, R428 was found to induce apoptosis in many types of cancer cells (Chen et al., 2018). Our findings concerning R428-analog-based modifications demonstrate more promising results than that of crizotinib. In particular, the designed drug candidates-compounds R3, R5, and R10 (see Table 5) -possess an improved binding property toward Axl compared to the known type I inhibitors (see Table 1). Besides this, they are also involved in the interaction with the D690, N677, M623, and H625 residues (see Table 6) that are conserved for all TAM family receptors according to the sequence alignment as shown by Gajiwala et al. (2017). To evaluate the compatibility of the compound R5 to Axl kinase, in Figure 8, for comparison, we demonstrate ATP, R428, R428_1, and R5 inside the binding pocket of Axl kinase. As the figure shows, ATP does not fit the pocket well due to its small size (see Figure 8A) and, for this reason, does not result in strong binding. R428 fits much better than ATP (see Figure 8B), and, therefore, the binding score is much higher than that of ATP. Compound R428_1, in turn, fits even better than the compound R428 (see Figure 8C), resulting in the highest score among the R428 analogs considered. Finally, our best-designed drug candidate-compound R5 (see Figure 8D)-shows the best fit and the strongest binding affinity toward kinase.
The MD simulations are in strong accordance with the molecular docking results, indicating a consistency between the results obtained by two different methodologies. In particular, the best binding molecular docking results were obtained for the R5 compound compared to the parental R428 drug and the best R428 patented analog-compound R428_1. These findings were in agreement with the MM-PBSA calculations, which were more optimized for these systems analyzed than with the MM-GBSA approach ( Table 8). The latter method determined the second best affinity value for R5 to the Axl kinase pocket after Axl-R428_1 ( Table 9). The discrepancies in the MM-PBSA/GBSA performances were previously emphasized in the literature, showing that MM/PBSA performed better in calculating absolute, but not always binding free energies than MM-GBSA (Hou et al., 2011). Furthermore, the calculated K i for R5 has shown the lowest value (Table 7), indicating its high inhibition potency.
Additionally, the molecular-docking-based selectivity test was performed by screening a set of kinases (ALK5, ABL1, FYN, JAK1, Met, Tyro3, and Mer) to evaluate the binding affinity of the compound R5 to them. Interestingly, this compound has shown the strongest binding properties toward Tyro3, which is another receptor of the TAM family. We tried to explain this phenomenon by performing a comparative analysis for the binding pockets of the high-scoring (Axl, Tyro3, and ABL1) and the relatively low-scoring (MET) kinases (see Table 11). It FIGURE 8 | The ATP binding cavity on the surface of the Axl kinase and the drug accommodation inside the pocket. The ATP is shown in cyan (A), the compound R428 is shown in green (B), the compound R428_1 is shown in yellow (C), and our best designed compound R5 is shown in red (D). The pocket is represented as a surface, the protein chain is represented as a cartoon, the hydrophobic patches of the pocket are shown in light green, the hydrophilic patches are shown in purple, and the hydrogen bond acceptors are shown in red. The compounds in the pocket are represented as a licorice.
is clear from the molecular docking poses (see Figure 9) that R5 exhibits a "stretched" conformation inside the binding sites of Axl, Tyro3, and ABL1 and a "shrunk" form in the MET. There are four residues (D690, E546, D627, and G626) of the Axl pocket involved in the interaction with R5 ( Figure S1A). Compared to the Axl pocket, five residues in the Tyro3 (K597, M596, R512, D663, and K540) and the ABL1 (D325, E329, L248, T315, and E316) pockets are involved in the interaction with R5 (see Figures S1B,C). In contrast, only three residues of the MET pocket (K1161, Y1230, and D1164) are involved in the interaction with R5 ( Figure S1D). Therefore, we think that the 3D shape of the binding sites as well as the number of interacting residues influence the change in the binding scores. Furthermore, the Axl and Tyro3 kinases were already investigated as important drug targets for various types of cancers (Duan et al., 2016;Dantas-Barbosa et al., 2017), expanding the possibility of R5 application against different malignant tumors.
TAM activation can, thus, be efficiently inhibited by the smallmolecule inhibitors for this interaction. For the TAM activation blocking compounds RU-301 and RU-302, Kimani et al. (2017) suggested and experimentally observed that these may act as pan-TAM inhibitors as they suppress the H1299 lung cancer tumor growth in the mouse xenograft NOND-SCIDγ model tested. Individual studies are, of course, limited to only either modeling or experiments in vitro or in cell culture. Similarly, inhibiting only one cancer pathway may not be the strongest strategy, particularly in view of clinical application. However, there is now consensus that, due to their well-established and clear inhibitory effect on several cancer pathways such as cell survival, invasion, migration, chemo-resistance, and metastasis, such TAM inhibitors hold great promise for cancer treatment (Davra et al., 2016). Moreover, further studies have shown that the TAM family expression often correlates with poor clinical outcomes, and there are several classes of promising TAM inhibitors according to the accumulated experimental data (Davra et al., 2016). We reveal here a particular promising TAM inhibitor derived from the experimentally validated compound R428; this is compound R5. We have, hence, good reasons, including the solid experimental validation of closely related compounds and using only well-established and experimentally validated methods and compound scoring schemes, to suggest R5 as a novel anti-cancer compound. However, this is, of course, only the starting point for further research. To be sure, this needs even more experimental data and subsequent compound-specific tests and clinical trials. As a bioinformatics group, we hope to stimulate more experimental research toward this end. Moreover, the molecular approach of the TAM-receptor inhibition implemented here is also a valuable strategy not only in cancer but also in hemostasis and thrombosis and platelet function (Law et al., 2018), such that the compound R5 should be explored and tested further and also analyzed for such indications.

General Considerations and Caveats
It is, of course, possible to get a very high false-positive rate in the virtual screening, but this can mainly happen if you do not do any scoring or accept all targets equally well. In particular, as the reader can see in our manuscript, the chosen strategy strongly reduces the amount of false-positives. We always start from the best possible leads, and that is not only according to the scoring but also, more importantly, according to all available experimental data.
Nevertheless, we think that it is necessary to validate the proposed compounds with actual experiments, and that is the reason why we want to publish our results here, considering other purely theory articles such as those of Fu et al. (2019), Duan et al. (2019), and Zheng et al. (2019). We acknowledge that we highly value the hard work and the final proof of a direct experiment, and the whole point of our publication is to exactly simulate this. We also have a non-trivial result, i.e., several strong pharmacological leads, including experimental validation and their actual use in clinical settings, are incorporated and now further optimized by us, yielding a top compound, among the best candidates for an important and well-known cancer target, that looks very promising by three different scoring schemes.

CONCLUSION
In summary, systematic computational studies were performed to establish the best drug candidates for targeting Axl kinase to treat cancer. The first compound selection has shown that the most potent drugs for such a purpose are crizotinib, macrocyclic compound 1, and R428. Further searches were performed based on the analogs of pharmaceutically available drugs, namely, crizotinib and R428. The docking results for crizotinib analogs did not demonstrate scores higher than the original drug itself. Therefore, we focused on crizotinib by modifying it further. The R428 analog docking resulted in four high-scoring compounds, which were further utilized for structural modifications. The structural modifications of crizotinib did not show any significant Axl binding improvements, while modifications of the R428 analogs resulted in a new potential drug candidate-compound R5. Apart from the docking tests, the potency of this newly designed compound as an Axl kinase inhibitor has been confirmed by a molecular dynamics simulation in terms of protein-ligand complex stability. Furthermore, an in silico selectivity test against other kinases has also shown a high affinity of R5 toward ABL1 and Tyro3. Our in silico results suggest the application of the newly designed compound R5 particularly as an anti-cancer agent. Direct proof by experiment is now the next important step. Further synthesis and in vitro tests of this compound against various cancer cell lines, where Axl, ABL1, and Tyro3 play a significant role in proliferation, division and metastases are further steps to be taken, and further indications for TAM inhibition may follow.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
ES planned the study, performed the in silico screening and docking of the compounds, designed the new compounds, performed a selectivity test, analyzed the data, made corresponding figures and tables, and drafted the manuscript. SS performed molecular dynamics simulation, analyzed and reported the corresponding data, made tables, and participated in the manuscript drafting and editing. TD supervised the study, analyzed its data, and participated in the manuscript drafting and editing. All authors finalized the manuscript together.

FUNDING
This work was supported by the Land Bavaria and the Federal Ministry of Education and Research (BMBF), Era-Net grant 01KT1801 (project CLEARLY). The Land provides our Core Funding. Neither the Land nor the BMBF had any role in the study, they are unbiased governmental funding agencies that definitely did not influence any of the outcomes or data reported or had any role in the preparation of the study.