Identification of Novel AXL Kinase Inhibitors Using Ligand-Based Pharmacophore Screening and Molecular Dynamics Simulations

Abstract

AXL kinase is a promising target in novel drug discovery for cancer. A ligand-based pharmacophore model was generated with the Pharmit web server. Its inbuilt PubChem molecule database was screened and led to 408 candidate molecules. Docking of the AXL kinase active sites with the identified list of candidate molecules was carried out with Autodock Vina docking software. This resulted in four compounds selected for further investigation. Molecular dynamics simulation of two ligands (PubChem-122421875 and PubChem-78160848) showed considerable binding with AXL kinase. From the MM-PBSA binding free energies investigation, the PubChem-122421875 (G = −179.3 kJ/mol) and PubChem-78160848 (G = −208.3 kJ/mol) ligands had favorable protein-ligand complex stability and binding free energy. Hence, PubChem-122421875 and PubChem-78160848 molecules identified in this work could be a potent starting point for developing novel AXL kinase inhibitor molecules.

1. Introduction

Cancer is a prominent cause of death and an important barrier to rising life expectancy in every corner of the globe. Approximately half of all newly diagnosed cancers can be treated using the currently available treatment methods. Existing treatments, namely chemotherapy, radiotherapy, and surgeries or a combination of the above, are not the complete cure for cancer patients worldwide [1]. Along with low efficiency for some patients, these treatments show severe adverse effects. The development of novel treatment solutions is needed to mitigate the constant rise in cancer cases and death rates.

AXL is a receptor tyrosine kinase involved in the growth, variation, existence, and motility of numerous diverse cell types [2,3]. The Tyrosine kinase AXL receptor, belonging to the TAM family (Tyro3, Axl, Mer), is found in a variety of oncogenic processes. It is a promising target for anti-cancer therapy [4,5]. Initially, AXL gene expression was identified in chronic myelogenous leukemia patients [6,7], and it was later found overexpressed in breast cancer [3,8] and gastrointestinal stromal tumor cells [9]. Further, AXL overexpression has been identified in most human tumors such as prostate, kidney, pancreatic, and breast cancer [10,11,12,13,14]. Therefore, AXL kinase has become a target for potent small-molecule cancer inhibitors.

Previous studies identified various small molecules with the capacity to inhibit AXL kinase. The selective inhibitors Gilteritinib [15] and Bemcentinib (BGB324) [16] bind the active site of AXL kinase. Gilteritinib was approved by the US Food Drug Administration (FDA), and BGB324 is known as the first selective AXL inhibitor. Cabozantinib [17], another inhibitor chemical licensed by the FDA for different malignancies, has a 200-fold higher potency against the vascular endothelial growth factor receptor 2 (VEGFR2) than AXL [18]. Foretinib [19], merestinib [20], and CEP-40783 [21] are other AXL ligands undergoing clinical studies. Y. Wang et al. developed Gilteritinib derivatives using a popular side chain ring closure approach and found new molecules with good AXL inhibition activity [22]. To date, there are only a few selective inhibitors for AXl kinase, and non-selective inhibitors that also inhibit other kinases. The development of novel selective AXL inhibitors and validation of their safety and efficacy are required.

2. Methodology

2.1. Pharmacophore Modeling Based on Structure

Ligand-based, structure-based pharmacophore modeling and screening of databases were performed with the Pharmit server (http://pharmit.csb.pitt.edu/ accessed on 20 February 2022) [23,24,25]. Pharmacophore modeling involves structure-based modeling of the ligand-protein complex using a 3-dimensional structure of the residues of the active sites of the protein interacting with its crystal or the docked inhibitor molecule. Hence, pharmacophores found with this approach include information about the protein active site regions that help design or search for new ligands with specific features. The Pharmit web server is free, user-friendly, and interactive when searching for pharmacophores and screening inbuilt databases and user-created databases. The screening is carried out by energy minimization of conformations in the active sites [26].

The pharmacophore generated from Pharmit can pick the selective molecules that can show less strain and map with the active site residues’ complementary features. The main drawback of the program is that sometimes, mapping of the known inhibitors cannot map the pharmacophore’s features due to its selective conformational features. The PubChem database molecule library was selected for the screening of the three generated pharmacophore model results from the most similar molecules with pharmacophore mapping. The best-mapped molecules docking on the AXL kinase active site is performed using Vina software.

A decoy molecule database was prepared from the 26 different kinase decoy sets from the DUD-E database (http://dude.docking.org/ accessed on 1 August 2022). The molecules were uploaded to the Pharmit server as zip files, and the conformers were generated automatically. The pharmacophore generated from the protein complex can be validated with the decoy set of molecule conformations.

2.2. Virtual Screening

Virtual screening involves searching large chemical libraries. It supports finding potential leads for specific target proteins [27]. PDB coordinates (PDB_ID: 5U6B) [28] are the AXL kinase protein coordinates from the RCSB protein databank as well as identified molecules from the pharmacophore screening database. These were converted to a format suitable for autodock docking (pdbqt) using AutodockTools software [29]. Virtual screening of molecules from the database was performed against AXL kinase with Vina software [30]. Missing residues in AXL kinase were modeled with a 5U6B crystal structure template using the online SWISS-MODEL (https://swissmodel.expasy.org/ accessed on 20 February 2022). Docking input files and grid boxes around AXL kinase’s ATP binding site and the grid center (−20.518 × −4.68 × 10.651) with suitable dimensions (20 × 18 × 34 ų separated by 1.0 Å) were prepared in Autodock Tools. Each docking output with 20 conformers and their affinity scores were visualized using Pymol.

2.3. Molecular Dynamics (MD) Simulations and MMPBSA Free Energies

To understand the stability of the ligand in the active site, four complexes were chosen based on their molecular docking affinity and orientation. These were submitted to molecular dynamics (MD) simulations of 100 nanoseconds under NPT conditions [31,32,33,34]. Gromacs 2018 simulations package with Ubuntu-18.04 LTS OS was used to perform MD simulations of the complexes. The inbuilt GROMOS96 53a6 force field was assigned to the protein [35], and an online-based PRODRG server was used to add the GROMOS forcefield to the inhibitor molecules [36]. A 10 Å edge length cubic box was placed around the protein, and the box was loaded with SPC (simple point charge) water simulation models. These complex systems had negative charges and were neutralized by adding Na⁺ ions. MD simulations under NVT and NPT conditions were performed, and 100 ns of NPT (isothermal and isobaric ensemble) MD simulations were performed with a 0.002 ps time step. The temperature of the total system was set at 298 K by using the Thermostat with V-rescale [37], and a pressure of 1 bar was maintained by the Parrinello–Rahman method [38]. A LINear Constraint Solver (LINCS) algorithm was used to constrain the hydrogen bonds in the system [39]. The converged region and the last 20 ns trajectory of 100 ns simulations in total extracted coordinates were used in MM-PBSA free energy calculations. AXL kinase bound with ligand binding free energies were calculated using a free energy code mm_pbsa [40]. The equation below was applied when calculating the free energy change of ligands bound to the protein.

ΔG_binding = G_{(ligand + protein)} − (G_ligand + G_protein)

G = 〈E_MM〉 − TS + 〈G_solvation〉

Here, the average molecular mechanics potential energy = ⟨E_MM⟩; temperature in Kelvins = T; entropy = S; component of free energy of solvation = ⟨G_solvation⟩; component free energy of solvation G_solvation = G_non-polar + G_polar.

The addition of both polar and non-polar free energies gives free energy solvation. Solving the equation of Poisson Boltzmann yields G_polar, and the non-polar model of the solvent-accessible surface area yields G_non-polar.

3. Results and Discussion

3.1. AXL Kinase Active Site Pharmacophore Modeling

A ligand-based pharmacophore model was generated using a macrocyclic inhibitor with molecular formula-C₂₁H₂₂ClFN₆O₃ called (10R)-7-amino-11-chloro-12-fluoro-1-(2-hydroxyethyl)-3,10,16-trimethyl-16,17-dihydro-1H-8,4-(azeno)pyrazolo [4,3-h][2,5,11]benzoxadiazacyclotetradecin-15(10H)-one). The macrocyclic inhibitor present in the protein (PDB_ID: 5U6B) resulted in seven pharmacophore features, which include (purple, with a 1 Å radius) aromatic features; (green, with a 1 Å radius) two hydrophobic groups; (orange, with a 1 Å radius) two hydrogen bond acceptors; and (white, with a 1 Å radius) two hydrogen bond donors. Figure 1 depicts the pharmacophore features of macrocyclic inhibitors present in the AXL protein.

The generated pharmacophore model was used as input for structural-based pharmacophore screening of lead molecules in the PubChem database. PubChem is a popular chemical repository database of chemical compounds. The molecules matching with the pharmacophore were aligned with the crystal structure of the inhibitor, and during the database search, the compound’s conformations were modified to fit within the AXL kinase protein’s active site. The decoy set validation also resulted in no molecule mapped with the pharmacophore.

3.2. Virtual Screening Analysis

Virtual screening of the >400 molecule library on AXL kinase’s active site was performed using Autodock Vina software to obtain the best hit molecules based on binding affinities (kcal/mol). The validation of the docked complex comprises docked conformation and orientation of the ligand in the active site. From the screening of 408 molecules, the top 100 hits were filtered with the best binding affinity scores. Further evaluation of these hit molecules was carried out by manual structural binding interactions through visualization tools (Supplementary File, Table S1). Non-covalent interactions are also crucial in maintaining the stability of this protein system complex [41,42,43,44]. Based on the position of the docking of the molecule in the binding pocket and binding scores, the top four molecules were chosen from a total of 100 molecules. The top-scoring and the best ten binding-oriented molecules are shown in

Table 1. Two-dimensional structures of the best binding molecules to the active site of the AXL kinase with its binding energies (in kcal/mol).

S. No.	PubChem ID	Molecule Structure	Binding Score
1	PubChem-78160848		−8.3
2	PubChem-122421875		−7.5
3	PubChem-76973938		−8.2
4	PubChem-57928218		−7.9
5	PubChem-54673081		−8.3
6	PubChem-126726905		−7.9
7	PubChem-68651602		−7.9
8	PubChem-76167946		−7.8
9	PubChem-25224351		−7.8
10	PubChem-68650375		−7.7

. Molecular docking and the interactions between the molecules and the AXL kinase residues are shown in

Table 2. Molecular docking of four selected molecules in the active site of the AXL kinase and their interactions with the AXL kinase amino acid residues.

S.No	Molecule	Bonds between Protein and Ligand	Type of Bond	Distance in Å
1	PubChem-76973938	LIG:NH–Pro110 LIG:N–Met112 LIG–Met168 LIG–Leu31 LIG–Asp179 LIG–Met168 LIG–Arg165	Hydrogen Bond Hydrogen Bond Pi-Sigma Pi-Sigma Pi-Anion Pi-Sulfur Halogen	2.09 2.18 3.44 3.77 4.71 5.87 6.72
2	PubChem-57928218	LIG: NH–Met112 LIG:N–Met112 LIG:OH–Lys56 LIG:O–Lys56 LIG–Leu31 LIG–Met168 LIG–Ala54 LIG–Leu109	Hydrogen Bond Hydrogen Bond Hydrogen Bond Hydrogen Bond Pi-Sigma Pi-Sigma Pi-Alkyl Pi-Alkyl	2.69 2.41 2.49 2.87 4.00 3.48 4.02 4.93
3	PubChem-122421875	LIG:N–Lys56 LIG:NH–Pro110 LIG:NH–Asp179 LIG:NH–Leu31 LIG–Met168 LIG–Ala54 LIG–Val39 LIG–Leu109 LIG- Met112	Hydrogen Bond Hydrogen Bond Hydrogen Bond Hydrogen Bond Pi-Sigma Pi-Alkyl Pi-Alkyl Pi-Alkyl Pi-Alkyl	2.51 2.05 2.38 2.29 3.73 4.73 4.99 4.28 5.49
4	PubChem-78160848	LIG:NH–Pro110 LIG:NH–Lys113 LIG:N–Met112 LIG–Met168 LIG–Leu31 LIG–Asp179 LIG–Met168 LIG–Asn166 LIG–Asp179	Hydrogen Bond Hydrogen Bond Hydrogen Bond Pi-Sigma Pi-Sigma Pi-Anion Pi-Sulfur Fluorine Fluorine	2.08 2.60 2.24 3.45 5.45 4.55 5.47 4.23 3.89

. In the four molecules of the molecular dynamic simulations, two molecules presented stable binding in the active site.

The molecule identified as PubChem-76973938 becomes stabilized with both hydrophobic and hydrophilic nonbonding interactions (Figure 2A,B). The pyridinic NH2 group is involved in hydrogen bonding with the Pro110 residue. The pyridinic nitrogen is involved in a hydrogen bond with the Met112 residue. The NH of the five-member cyclic ring is part of a hydrogen molecule with the side chain of the Asp116 residue. The pyridine ring combines with Ala54 and Met168 residues via a CH-π interaction.

The docking of the molecule PubChem-57928218 has strong interactions with the AXL kinase active site residues. Its azaindole nitrogens are part of hydrogen bonds with the Met112 residue. The hydroxyl group forms a hydrogen bond with the side chain of the Lys56 residue (Figure 2C,D). With Ala54, Met168, and Leu313 residues, the azaindole ring forms a CH-π interaction. It was observed that the PubChem-122421875 molecule shows a very strong binding interaction in the docking studies with AXL kinase active site residues (Figure 3). The azaindole nitrogens of the PubChem-122421875 molecule are part of two hydrogen bonds with Pro110 and Met112 residues. The nitrogen of the pyridine ring interacts via hydrogen bonding with the Lys54 residue. The pyridinic NH forms a hydrogen bonding with the Asp179 residue. The oxygen substitution in the pyridine ring and the amide of the molecule form hydrogen bonds with the Leu31 residue. Molecular docking of the PubChem-78160848 molecule shows strong affinity toward the AXL active site (Figure 3). The amino group in the pyridine forms hydrogen bonding with the Met112 and Pro110 residues. The NH of the five-membered ring interacts via hydrogen bonding with the Lys113 residue’s main chain. The pyridine ring is stabilized with interaction with Ala54 CH-π, Leu313, and Met168 residues.

3.3. Molecular Dynamics Simulations

Molecular docking studies led to four molecules that were subjected to molecular dynamics simulations of 100 ns with AXL kinase to reveal the molecule binding and stability in the active sites. The two molecules PubChem-122421875 and PubChem-78160848 showed binding stability in the active site. Further, the results of PubChem-122421875 and PubChem-78160848 molecules are discussed. The Apo form of the AXL kinase simulations results in the protein RMSD with 0.5 nm at the end of the simulations indicating the protein stability, and the RMSF shows the fluctuations of the atoms at the regions Ser97-Phe102 (914–982 atoms), Tyr123-Gln129 (1178–1250 atoms), and Ala211-Tyr215 (2097–2150 atoms). The regions are the loop regions and are far from the active site (Figure 4A,B).

The RMSD of the AXL protein complex with PubChem-122421875 is converged around 0.45–0.5 nm. In the end, the restoration of the simulations was observed through a slightly high RMSD at 80 ns in the simulation (Figure 5A). The structure of the converged trajectory and RMSD plot shows the structural stability of the molecule. The RMSF plot indicates high fluctuations of the residue side chains (Figure 4B). The regions Met58-Glu67 (525–622 atoms), Ser97-Phe102 (914–983 atoms), and Arg125-Gln129 (1204–1256 atoms) at the N-terminal domain of the kinase showed high fluctuations, but these were in the loop regions, which are far from the binding site of the ligand. All the binding site residues are stable, and their side-chain fluctuations are limited. The stability of the ligand was observed from the average structure of the convergence for the last 10 nanoseconds. The molecules retained their hydrogen bond with the Pro110 residue. The NH from the N substituted pyridine ring created a new hydrogen bond with the His114 residue. Now, the amide carbonyl of the inhibitor interacts with the main chain of Ala54 residue via a hydrogen bond. Overall, from MD simulations, all residues in the binding site were stable, and their side-chain fluctuations were limited. The stability of the ligand was obtained from the average structure of the last 20 nanoseconds of convergence. These results reveal that the molecule PubChem-122421875 is strongly bound to the AXL kinase active site. The interactions of the ligand are depicted in Figure 6A,B.

The RMSD plot of the PubChem-78160848 and AXL protein complex convergence around 0.30 nm and the RMSD convergence observed from the last 15 ns of the simulation are shown in Figure 7A. The protein-ligand structure stability was also observed from the visualization of the trajectory structure. The RMSF plot indicates low fluctuations of the residue side chains except for a few regions (Figure 7B). The regions Ser97-Phe102 (914–983 atoms) with 0.7 nm fluctuations and Arg125-Gln129 (1204–1256 atoms) with 0.5 nm fluctuations are looped regions of the kinase showing high fluctuations and which are also far from the binding site of the ligand.

The molecules retain their hydrogen bonds with the Pro110 residue. The NH five-membered ring of the molecule forms a new hydrogen bond with the Glu41 residue side chain. The pyridine ring of the inhibitor shows a π–π bonding with the side chain of Phe111. On the whole, from MD simulations, all residues in the binding site are stable, and their side-chain fluctuations are limited. The stability of the ligand was noticed from the average structure of the last 20 ns of the convergence test. These results reveal that the molecule PubChem-78160848 is strongly bound to the AXL kinase active site, and the interactions of the ligand are depicted in Figure 8A,B.

Overall, the Apo and the other two molecules with protein complexes converged at 100 ns in the simulations. The RMSD plots of the protein show convergence at various points. The RMS fluctuations indicate the atomic fluctuations are at similar regions of the protein Apo form and the complex forms. The most common fluctuation regions Ser97-Phe102 and Tyr123-Gln129 are at the loop regions connecting the beta sheets of the protein.

3.4. MMPBSA (Molecular Mechanics Poisson–Boltzmann Surface Area) Free Energy Calculations

The AXL kinase and the ligand binding free energy were calculated using the g_mmpbsa theory [40]. The calculation by mmpbsa was conducted using a convergence of 10 selected structures using the last 10 nanoseconds of the convergence of a total simulation time of 100 ns. The binding free energy of PubChem-122421875 was estimated to have −179.3 kJ/mol, and its components were: van der Waals energy (estimated at −181.4 kJ/mol), electrostatic energy (−309.7 kJ/mol), polar solvation energy (332.7 kJ/mol) and SASA energy (−20.9 kJ/mol). The binding free energy of PubChem-78160848 was estimated at −208.3 kJ/mol, and its components were van der Waals energy (−192.4 kJ/mol), electrostatic energy (−170.4 kJ/mol), polar solvation energy (173.5 kJ/mol) and solvent accessible surface area (SASA) energy (−19.0 kJ/mol). Both molecules present prominent interactions with good binding free energies. The PubChem-122421875 molecule has an azaindole ring as core moiety, and PubChem-78160848 has a pyridine as a core moiety.

4. Conclusions

Ligand-based pharmacophore modeling and screening of the PubChem database were used to short-list molecules to be screened for their interaction with the AXL kinase active site. This yielded the four best ligands with high affinities for interacting with a pocket in the receptor. From molecular docking results, among the four AXL–molecule complexes, two showed particularly strong binding. These were considered for stability studies using molecular dynamics simulations. The molecule identified as PubChem-122421875 led to −179.3 kJ/mol, and the molecule identified as PubChem-78160848 presented −208.3 kJ/mol in MM-PBSA ligand binding free energy calculations. These molecules stayed stable and had low RMSD values during simulations. The binding energies were negative for key residues Ala54, Glu74, Pro110, Met112, Leu31, and Asp179 with both ligands PubChem-122421875 and PubChem-78160848 and played important roles in the stability of the complexes with the AXL kinase active site.