In Silico Designed Axl Receptor Blocking Drug Candidates Against Zika Virus Infection.

After a large outbreak in Brazil, novel drugs against Zika virus became extremely necessary. Evaluation of virus-based pharmacological strategies concerning essential host factors brought us to the idea that targeting the Axl receptor by blocking its dimerization function could be critical for virus entry. Starting from experimentally validated compounds, such as RU-301, RU-302, warfarin, and R428, we identified a novel compound 2' (R428 derivative) to be the most potent for this task amongst a number of alternative compounds and leads. The improved affinity of compound 2' was confirmed by molecular docking as well as molecular dynamics simulation techniques using implicit solvation models. The current study summarizes a new possibility for inhibition of the Axl function as a potential target for future antiviral therapies.

Introduction

Zika virus (ZIKV) belongs to the family of flaviviruses, which was first isolated in 1947 from rhesus monkeys in the Zika forest of Uganda.[1−4] The ZIKV epidemic has drawn global attention after a large outbreak in Brazil in 2014, which spread to over 60 countries. Transmissions to humans are usually mediated by mosquitoes. However, the virus can also be transmitted through sexual contact[5−8] as well as from mother to fetus.[9−11] Clinical observations suggest that ZIKV, transmitted from an infected woman to a fetus, mainly, targets human neural progenitor cells (hNPCs) and impairs their development.[4] There are no established drugs or vaccines against ZIKV yet. However, there are different approaches suggested regarding how to deal with ZIKV. For example, Ekins[12] and coworkers suggest as a first screening step defining a cell or target-based ZIKV in vitro assay for a high throughput compound screening and as a second step to test these compounds. The third step is to investigate the entire genome of ZIKV for already known drug targets that might show some activity against ZIKV in vitro. As a fourth step Ekins et al. suggest homology model development for ZIKV proteins similar to those that are used against dengue virus (DENV) treatment.[12] Later on, they proposed to understand the exact mechanism of action to screen the compounds with higher target specificity and, finally, to test them on animals.[12]

Another strategy to cease ZIKV replication was proposed by Mumtaz and co-workers.[13] In particular, antiviral drugs such as 7-deaza-2′-C-methyladenosine (7DMA), 2′-C-methylcytidine, ribavirin, favipiravir, and T-1105 were suggested as a potent NS5 polymerase inhibitors.[13] Although ribavirin and favipiravir have other indications, they were also tested against ZIKV, providing less antiviral effect in vitro than 7DMA. In silico drug screening against ZIKV NS2B-NS3 protease inhibitors has been performed by Yuan et al.,[14] which, ultimately, yielded two potent drugs—novobiocin and lopinavir/ritonavir. The potency of novobiocin has been confirmed by molecular dynamics (MD) simulation, as it was forming a stable complex with a protease,[14] as well as by in vivo studies, where it showed a high survival rate in novobiocin-treated mice compared to untreated ones.[14] Other medications, such as emricasan and niclosamide, were tested by Xu et al.,[4] which also aimed to inhibit ZIKV replication. We follow here a new promising and independent strategy to prevent ZIKV infection from the start, blocking the entry via the host receptor.

As mentioned above, ZIKV causes abnormalities in the developing fetal nervous system.[15] However, the exact mechanism by which the virus enters and affects neural cells is poorly investigated. Nowakowski[15,16] in his work reported that ZIKV entry into neural stem cells is mediated by the Axl receptor, one of the members of TAM family receptors.[17] The Axl receptor was found to be highly expressed in radial glia cells,[16,17] neural stem cells,[18] as well as endothelial cells,[18] and hence, it was hypothesized to be susceptible for the ZIKV entry. It is worth to note that Axl itself does not play a role of the entry factor for the virus. The virus entry into the host cells occurs only when the Axl receptor is dimerized with its ligand—the Gas6 protein.[19,20] The dimerization mechanism serves as a bridge for binding to the phosphatidylserine molecules surrounding the viral capsid. The whole process triggers the fusion of the virus inside the host cell.[19,20]

Following this path, we propose blocking the Axl dimerization with its ligand (Gas6) to prevent the virus entry into the human neural cells. This can be achieved via competitive binding of drug-like molecules to the Axl–Gas6 dimerization site. In this study, we compared already existing compounds such as warfarin,[21] R428,[22] RU-301,[23] and RU-302[23] to in silico designed novel drug-like molecules to find lead candidates with higher affinity to the Axl receptor dimerization site using molecular docking and MD approaches with implicit solvation models.

Results and Discussion

The role of the Axl receptor for ZIKV entry is still under investigation. Several labs reported that the Axl receptor may not play an important role in ZIKV infection.[24−27] Therefore, for instance, Wang et al.[24] claim that the Axl-deficient mice show the same level of replication of ZIKV as nondeficient ones. The findings of Wang were previously supported by the work of Li et al.,[25] in which they demonstrated in an Axl-knockout mouse model that the Axl receptor was not the key entry receptor for the ZIKV infection. Similar results were obtained by Hastings et al.,[26] indicating that TAM receptors (Axl, Tyro3 and Mertk) are not necessary for ZIKV entry in mice. They show this by blocking the type I interferon receptor with a MAR1-5A3 antibody, which did not have any impact on reducing the ZIKV infection in mice.[26] Wells et al.[27] reported that the Axl deletion had no effect on ZIKV entry in nondifferentiated hNPCs. In contrast, several in vitro studies in various differentiated human cell types show that blocking of the Axl receptor stops ZIKV infection.[18,28−31]

Taking into consideration all these results above, we conclude that the different observations relate to differences in the composition of entry receptors in the mouse and human. For this reason, it is well-possible that ZIKV entry is not Axl-sensitive in mice and, according to one report, may even be in human nondifferentiated neuronal progenitor cells, while in humans, the observations point to the Axl receptor as an important entry route for ZIKV, at least, for all differentiated human cell types tested.

Another caveat is that cells in vitro are not able to express the whole spectra of the entry receptors that they usually express in vivo, and different types of tissue can also differ in the composition of entry receptors.

Following our hypothesis, we used in silico methods to design new drug-like molecules for inhibition of the Axl receptor dimerization and, hence, prevent ZIKV infection in human cells. Besides the Axl protein as a drug target, there are other viral proteins and host factors as potential candidates for inhibition (Table ). For instance, vaccine strategies can target surface epitopes,[32−34] as these methods have been previously applied for flavivirus glycoproteins (glycoprotein E, for example). On the other hand, most viruses can acquire rapid drug resistance against viral polymerase or protease inhibitors.[35,36] This can be also confirmed by multiple alignments of DEAD-like helicases, glycoprotein E, and RNA-directed RNA polymerase performed for the ZIKV, West Nile (WNV), dengue (DENV), and yellow fever (YFV) viruses. As the alignments show, there are many conserved regions for RNA-directed RNA polymerase shared among all considered strains of these flaviviruses. However, less conserved regions can be observed for the DEAD-like helicases and several for the E protein C terminal domains, which indicates that the virus can undergo mutations easily. The conserved domains in the above-mentioned alignments are labeled as an asterisk sign (please refer to the separately attached .clw files for the DEAD-like helicases superfamily, E protein C terminal domain, and RNA-directed RNA polymerase in the Supporting Information). The references describing each sequence in the alignment are given in the Supporting Information.

Table 1

Possible Protein Targets and Available Drugs for Blocking ZIKV Infection Spread

	protein targets	drugs	references
ZIKV	NS5 poymerase	7DMA, 2′-C-methylcytidin, ribavirin, favipiravir, T-1105	(13)
	virus replication	niclosamide	(4)
host factor	Axl receptor	warfarin, R428	(21,22)
		BMS-777607, PF-02341066, Janssen	(65)
		RU-301, RU-302	(23)
		our new hypothesis: compounds 1–5, compounds 1′ and 2′
	caspase-3	emricasan	(4)
host factors regarding virus packaging, virus secretion, immune surveillance	secretory system, immune proteins		(66−71)

In this study, we, therefore, concentrated on the Axl receptor, a host protein required for virus entry, as a promising drug target to treat the ZIKV infection. Using the molecular docking method, we tested RU-301 and RU-302, the compounds designed by Kimani et al.[23] specifically for targeting the Axl–Gas6 dimerization site to suppress cancer disease by Gas6-induced Axl signaling, as well as warfarin and R428—known Axl blockers.[21,22] Warfarin is an indirect Axl blocker, an anticoagulant drug known for binding and inhibition of vitamin K epoxide reductase - an enzyme responsible for carboxylation of glutamic acid residues in the binding proteins, which are important for blood coagulation.[37] In this case, inhibition of epoxide reductase by warfarin would prevent carboxylation of glutamic acid residues in the epidermal growth factor domain of Gas6 protein (ligand for the Axl receptor), responsible for attachment to the membrane of the viral envelope and facilitating fusion of the virus into the host cell.[19] R428 is an anticancer drug that is intended to target the kinase domain of the Axl receptor.[22,38] However, most of the anticancer kinase inhibitors are known to be promiscuous, that is, they also target other kinases and/or have additional binding modes.[38] Therefore, the choice of the considered compounds is based also on this logic.

To evaluate the binding capabilities of these and our new derivatives, we, first, tested the abovementioned designed compounds RU-301 and RU-302 suggested by Kimani et al.[23] In Figure are shown these compounds at the Axl-Gas6 dimerization site. Both compounds fit the pocket in a similar binding mode. In particular, they are accommodated, mainly, occupying the space between E 56 and E 59 residues and exposing the pyrrolidine ring toward H 61. Slightly better binding is noticed for RU-301 compared to RU-302 (please refer to Table ), which is in agreement with experimental data. According to the results obtained by Kimani and coworkers,[23] both RU-301 and RU-302 show potent Axl–Gas6 inhibitory activity at 10 μM in oncogenic cells. RU-301 has shown slightly higher inhibitory activity than RU-302. In particular, the average percentage of Gas6 inhibition in H1299, U2-OS, and Calu-1 cells for RU-301 is 59, 52, and 44%, respectively, while the same activity for RU-302 is 45, 50, and 43%, respectively.[23] These data are in a good accordance with our docking results, showing that RU-301 binds to the dimerization site with slightly higher affinity (lower binding energy) than RU-302 (see Table for clarity).

Figure 1

Compounds RU-301 (violet red) (A) and RU-302 (medium purple) (B) at the Axl–Gas6 dimerization site. The binding pocket is shown as a secondary structure representation and is highlighted in green, the key residues are shown as a ball and stick representation and are highlighted in yellow, and the residue involved in additional contact is shown in cornflower blue.

Table 2

Lowest Binding Energies for RU-301, RU-302, Warfarin, R428, and Their Derivatives Obtained from Docking Simulations

compound	binding energies ΔGa, kcal/mol
RU-301	–6.69
RU-302	–5.72
warfarin	–5.24
1	–5.11
2	–6.2
3	–6.7
4	–8.0
5	–7.7
R428	–10.5
1′	–11.9
2′	–13.06

ΔG is Gibbs free energy = final intermolecular energy (vdW + Hbond + desolvation energy + electrostatic energy) + final total internal energy + torsional free energy – unbound systems energy. Here, vdW stands for van der Waals interactions and H-bond stands for hydrogen bond formation energy.

Next, to improve the Axl inhibition possibilities, we modified warfarin and R428, the pharmaceutically available drugs, using them as chemical “scaffolds” according to the scheme described in Figure . The summarized docking results for these novel analogues are shown in Table , where compound 2′ (ΔG = −13.06 kcal/mol) is detected as a lead molecule with the strongest affinity to Axl.

Figure 2

2D structures of warfarin and R428 (left) as chemical “scaffolds”. In silico design of novel analogues by structure modifications (right). Compound 1 has been obtained by adding the C–C=O–C group to warfarin that is indicated by X. Compound 2 was generated by the addition of one more benzene ring at the position indicated by Y. Compounds 3–5 were generated by adding methionine, phenylalanine, and tryptophan to compound 1. Compound 1′ was constructed by adding arginine to R428 at position indicated by X′ and compound 2′ was produced via the modification of compound 1′ through pyrrolidine ring extension at position indicated by Y′.

Figure shows the conformational changes of warfarin and its analogues in the Axl dimerization site, indicating their structural misfit to the binding cavity due to poor binding properties: a conformational shift inside the binding pocket toward E 56 and E 59 residues was observed for the warfarin, compounds 1 and 2 (Figure B,C). This change involves warfarin modification through the C–C=O–C extension (compound 1) and cyclohexane addition (compound 2). Eventually, these modifications contributed to the slight improvement of warfarin-based compounds via additional interactions of the C–C=O–C moiety with the Q 78 amino group of Axl (Figures B, C, S1 and S2 of the Supporting Information). On the other hand, compounds 3–5 are found to be more potent inhibitors of Axl in terms of binding energies (see Table ) compared to compounds 1 and 2 and the parental drug itself. They are more symmetrically oriented with respect to the key amino acid residues (Figures D–F and S3–S5 of the Supporting Information), and this is, mainly, due to the interactions they make with these residues. So, the amino group of compound 3 interacts with the carboxyl group of E 56, while the carbonyl group of the compound interacts with the amino group of Q 78 (Figure S3). Similarly, the amino group of the phenylalanine, which belongs to compound 4, interacts with the carboxyl group of E 56 and the oxygen atom of the compound’s aromatic ring interacts with one of the hydrogen atoms of the Q 78 amino group (Figures E and S4 of the Supporting Information). In the case of compound 5, there is a contact between the amino group of the ligand and the carboxyl group of E 56 (Figures F and S5 of the Supporting Information). There are two contacts observed with Q 78: one is between the oxygen atom of Q 78 and the hydrogen atom from the NH group of tryptophan that belongs to compound 5; the other one is between the hydrogen atom of the amino group of Q 78 and the oxygen atom of the aromatic ring of the compound. The amino group of Q 76 interacts with the carbonyl group of compound 5 (Figures F and S5 of the Supporting Information).

Figure 3

Graphical representation of warfarin (A) (purple) and its analogues (B–F) [compound 1 (B)—light sea green, compound 2 (C)—medium blue, compound 3 (D)—cyan, compound 4 (E)—orange, and compound 5 (F)—magenta] at the binding site of the Axl receptor. The receptor is represented as a secondary structure and is shown in green, the key binding residues responsible for ligand dimerization, such as E 56, E 59, T 75, T 77, and V 79, are highlighted in yellow, while the residues that make additional contacts are shown in cornflower blue. The compounds as well as the residues are shown in ball and stick representation.

R428-based analogues provided even higher binding affinity to the Axl in comparison to warfarin-based ones. The R428 drug was slightly shifted toward the E 56 and E 59 residues (Figures A and S6 of the Supporting Information). In particular, the amino group of R428 was likely to create two additional contacts: one with the carboxyl group of E 56 and the other with the carboxyl group of E 85 (Figures A and S6 of the Supporting Information). Any additional R428 modifications were introduced to enhance the structural fit to the Axl pocket. The results of docking for the compound 1′ have shown a higher affinity of binding compared to the original drug. Besides the interactions with the carboxyl groups of E 56 and E 59 inside the pocket, an additional interaction with the amino group of arginine, which is the part of the newly designed drug, with the carboxyl group of E 85 and Q 78 of the receptor resulted in better fit and stronger binding than that of the R428 (Figures B and S7 of Supporting Information). Compound 2′ has shown the best binding affinity to the Axl receptor, probably, because of more symmetrical occupancy of the receptor’s binding pocket (Figure C). The polycyclic modification of compound 2′ allowed it to be accommodated in a pocket with higher structural congruency, establishing a hydrophobic contact with E 85. The arginine part of the compound interacts with a carboxyl group of E 59. Furthermore, the CO group of Q78 makes contacts with the amino group of arginine (Figure S7 of the Supporting Information), which, in turn, strengthens the binding of compound 2′, expanding the possibilities to use the R428 derivatives as the Axl receptor inhibitors not only for ZIKV treatment, but also in cancer therapy.

Figure 4

Graphical representation of R428 (A) (red) and its derivatives (B,C) (compound 1′ (B)—magenta and compound 2′ (C)—cyan) at the binding site of the Axl receptor. The receptor is interpreted as a secondary structure representation and is shown in green. The key residues E 56, E 59, T 75, T 77, and V 79 are highlighted in yellow. The residues that are involved in additional contacts are highlighted in cornflower blue. The compounds and the residues are represented as a ball and stick representation.

Validation Considerations and Data

First of all, we started with known and validated Axl inhibitors: RU-301 and RU-302 were tested and validated in vitro to block the Axl receptor dimerization site.[23] To identify new compounds for Axl receptor inhibition, we looked at warfarin and R428, which are the Axl receptor inhibitors sufficiently validated in their action to be used in clinical trials.[21,22,39] In Table , the experimentally available IC (IC50) values are reported. As can be noticed, the calculated value of Kd and the experimentally obtained inhibitory concentrations are in the similar range, which shows the congruency of the simulation results with experimental data. The calculated Kd value for the designed compound 2′ has shown the lowest value. This means that compound 2′ would have quite a low inhibitory concentration and, hence, a very strong inhibitory activity.

Table 3

Comparison of the Calculated Dissociation Constant and Experimentally Obtained Inhibitory Concentrationsa

compound	Kd (calc.)b	IC (exp.)c
RU-301	12 μM	10 μM [22]
warfarin	139 μM
R428	18 nM	14 nM (IC50) [29]
compound 2′	0.24 nM

The dissociation constants Kd for RU-301, warfarin, R428 and the designed compound 2′ have been calculated from ΔG values.[72]

Kd (calc.)—calculated dissociation constant.

IC and IC50—inhibitory concentration and half maximal inhibitory concentration, respectively, taken from experimental data.

To further validate our results, the binding free energies (ΔGbind) based on the implicit solvation models were estimated for the best binding compound (2′) to the Axl receptor and warfarin as a reference molecule. The MM-GB/PBSA calculations using the 50 ns MD trajectories confirmed the molecular docking results: compound 2′ has much higher affinity to the Axl receptor in comparison to the reference molecule (see Table ).

Table 4

Energetic Analysis of the Receptor–Ligand Complexes Using the MM-(GB) PBSA Solvation Models

energya	Axl-warfarinb	Axl-warfarinc	Axl-2′b	Axl-2′c
ΔEvdW	–10.35 ± 1.51	–10.35 ± 1.51	–45.12 ± 2.52	–45.12 ± 2.52
ΔEEL	–1.38 ± 1.64	–1.38 ± 1.64	–481.36 ± 12.65	–481.36 ± 12.65
ΔEGB	6.59 ± 1.86		503.17 ± 11.48
ΔESURF	–1.69 ± 0.21		–5.38 ± 0.18
ΔEPB		5.22 ± 1.51		495.87 ± 11.49
ΔENPOLAR		–9.72 ± 1.11		–29.45 ± 0.97
ΔEDISPER		15.52 ± 1.14		54.45 ± 1.15
ΔGgas	–11.72 ± 2.4	–11.72 ± 2.4	–526.49 ± 12.78	–526.49 ± 12.78
ΔGsolv	4.89 ± 1.77	11.02 ± 1.68	497.78 ± 11.41	520.87 ± 11.25
ΔGbind	–6.83 ± 1.25	–0.71 ± 1.55	–28.78 ± 3.46	–5.61 ± 3.69

kcal mol–1.

MM-GBSA.

MM-PBSA.

To investigate the structural movement of the studied systems, the root mean square deviation (rmsd) of the warfarin and compound 2′, bound to the Axl receptor, with respect to the initial conformation was plotted versus simulation time (see Figure S9 in the Supporting Information). The complexes stay well-stable, supporting good binding. Moreover, the root mean square fluctuation (rmsf) parameter of Axl presented almost identical atomic movements either in the complex with warfarin or compound 2′ (Figure S9 C and D of the Supporting Information), followed by the absence of significant secondary structure changes (Figure S10A,B in the Supporting Information). In our study, we analyzed the protein–ligand interactions but not the protein folding or stability using 50 ns MD simulation. The 50 ns time scale is entirely sufficient for this investigation as an optimum time (usually 50–100 ns) and has already been confirmed by numerous publications.[40−42]

Finally, with later application and preclinical or clinical tests in mind, the toxicity potential of compound 2′ was assessed using trained Bayesian models and a sub-structure similarity search algorithm. This predicted only minor cytotoxic effects of the novel lead in various in vitro assays compared to the reference molecule (warfarin) (see Figure A and B).

Figure 5

Toxicity predictions for warfarin (A) and compound 2′ (B) using trained Bayesian models and sub-structure similarity search algorithm. Color-coded diagrams represent the signatures with low (blue), intermediate (white), or high toxic activity (red). The abbreviations of these three different biological assay systems are reported in Table .

Table 5

In Vitro Assays Used in the Potential Toxicity Prediction Based on the Trained Bayesian Models

assay		description
beta lactamase-based	HSE-blax	heat shock factor response element
	ARE-bla	antioxidant response element
	P53-bla	P53 response element
	VDR-bla	vitamin D receptor
	PPAR-gamma-bla	peroxisome proliferator-activated receptor gamma
	PPAR-delta-bla	peroxisome proliferator-activated receptor delta
	AR-bla	androgen receptor
	ER-bla	estrogen receptor
	FXR-bla	farnesoid X receptor
	GR-bla	glucocorticoid receptor
luciferase-based	TR-beta-luc	thyroid receptor beta
	AhR-luc	aryl hydrocarbon receptor
	AR-MDA-luc	androgen receptor-responsive MDA breast cancer cell line
	ER-BG1-luc	estrogen receptor-responsive ovarian adenocarcinoma BG1 cell line
various	DT40 WT	wild-type chicken DT40 bursal lymphoma cells
	DT40 REV3	DT40 cells with DNA polymerase (REV3) mutation
	ATAD5	ATPase family AAA domain-containing protein 5

Conclusions

Our hypothesis is that targeting a host factor—the Axl receptor—could stop not only ZIKV infection, but also the spread of other flaviviruses, as it uses this receptor as a medium to enter into the host cell. Among the other suggested drugs, warfarin and R428, the known Axl inhibitors, have not been considered yet as promising candidates to tackle ZIKV. We propose here novel in silico designed compounds—warfarin- and R428-based derivatives—to prevent ZIKV internalization into the neural cells. Molecular docking simulations, performed for the RU-301, RU-302, warfarin, R428, and the new derivatives, have shown improved binding properties for the designed analogues compared to the original drugs in favor of R428-based ones. In particular, the new derivatives of R428, compound 1′ and compound 2′, demonstrated the best binding capabilities to the Axl receptor amongst all the studied compounds. The results of docking have been confirmed by MD simulation analysis, showing the highest affinity for compound 2′ to the binding pocket. Additionally, the toxicity test has predicted reduced toxicity for it. Therefore, after carefully looking at a range of alternative strategies, we would suggest compound 2′ as a lead based on our computational analysis. Lead synthesis and further preclinical tests are recommended regarding cytotoxicity and teratogenicity for future biomedical applications including blocking of the ZIKV infection as well as prevention of severe complications of flaviviridae in general.

Methods

Sequence Analysis

Iterative sequence similarity searches, domain analysis, and sequence alignment followed the standard protocols as previously described.[43,44]

Structure Derivation

The 3D coordinates for the Axl–Gas6 complex have been retrieved from the protein data bank (PDB) under the reference code 2C5D.[45] In the current docking study, an Axl receptor has been considered. The 3D structures of the drugs warfarin and R428 were derived from the PubChem database under the references CID 54678486 and CID 46215462, respectively.

The 3D structures of the drugs RU-301 and RU-302 have been constructed according to the 2D scheme provided by Kimani et al.[23]

Design of the Compounds

The residue composition and the shape of the binding pocket were carefully analyzed. The design of new compounds was mainly based on their affinity. We tried several possibilities of building small molecules and filling up the binding pocket by using warfarin and R428 as starting templates. After evaluation of binding energies, we finally identified five high-affinity derivatives for warfarin and two derivatives for R428. To obtain analogues, the structures of the drugs have been modified with the help of the Avogadro software[46] (version 1.1.1, available at https://avogadro.cc/). Compound 1 has been designed by adding the C–C=O–C group to the warfarin structure at the position indicated by X (please, refer to the left side of Figure for clarity). Compound 2 has been designed by adding the cyclohexane group; the position is indicated by Y (see Figure , left side). Compounds 3–5, in turn, were designed by the addition of methionine, phenylalanine, and tryptophan, respectively, to the structure of the already existing compound 1 (see Figure ). Compound 1′ has been designed by structural modification of the existing drug R428. In particular, the structure was derived by adding arginine at position X′. The compound 1′ has been further modified by extending pyrrolidine rings three more times at position, indicated by Y′; hence, compound 2′ has been obtained (refer to Figure for clarity). All the compounds prior to docking were energy-minimized using the steepest descent algorithm to relax the structures.

Molecular Docking

Molecular docking has been performed using the AutoDock software[47] (version 4.2.6, available at http://autodock.scripps.edu). There is a lot of alternative docking software around, for instance, BINDSURF,[48] METADOCK,[49] LeadFinder,[50] the BLIND DOCKING server (available at http://bio-hpc.ucam.edu/achilles/), FlexScreen,[51,52] and Vina.[53,54] All of these have specific advantages and limitations and aficionados. However, we opted for AutoDock, as experienced AutoDock users and also because it is free of charge. In addition, it has the Lamarckian genetic algorithm implemented in a software that describes well the conformation of the ligand inside the pocket. Moreover, it is one of the most reliable and most widely used protein–ligand docking software.[55,56]

The AutoDock tools were used to generate input parameter files for docking. For the current study, the receptor was considered as a rigid molecule, while ligands contained rotatable bonds. Only protein molecules were used for the docking, while all nonprotein moieties were discarded. 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 group of charges, the so-called Gasteiger charges. The volume of the grid box was set as 60 × 60 × 60 with a 0.375 Å spacing (default value). The default values for the grid volume and spacing were used as in the AutoGrid calculation procedure, the protein is placed in a three-dimensional grid, and at each grid point a probe atom is placed. The energy of the interaction of this single atom with the protein is assigned to the grid point. Typically, grid point spacing is in the range of 0.2 to 1.0 Å, and the default value, 0.375 Å, signifies roughly a quarter of the length of the carbon–carbon single bond.

The center of the grid box was centered so that it included the important amino acids responsible for the ligand dimerization (E 56, E 59, T 75, T 77, and V 79).[45] 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.

Molecular Dynamics

All MD simulations were performed using the AMBER 16 package.[57] These simulations were performed by the AMBER 12 package using the FF99SB and general Amber force field force-fields for the Axl receptor and ligand molecules. The systems were solvated with the TIP3P water models and neutralized with Na+ ions using the tLEaP input script available from the AmberTools. Long-range electrostatic interactions were applied via the particle-mesh Ewald method.[58] The SHAKE algorithm[59] was used to constrain the length of covalent bonds, including hydrogen atoms. The Langevin thermostat was implemented to equilibrate the temperature of the systems at 300 K. A 2.0 fs time step was used for all simulations. 10 000 steps and 1 ns time period were used for minimization and equilibration with reference to all studied systems. Finally, 50 ns classical MD simulations with no constraints were performed for each of the receptor–ligand complexes using the molecular mechanics energies combined with the Poisson–Boltzmann (MM-PBSA) or generalized Born (MM-GBSA) and surface area continuum solvation approaches.[42,60,61]

According to the AMBER MM-PBSA/GBSA protocol (Miller et al., 2012; http://ambermd.org/tutorials/advanced/tutorial3/),[62] we used the explicit solvation model for all MD simulations and later the implicit solvation (MM-PBSA/GBSA) as a postprocessing end-state method to calculate free energies of molecules in solution by means of the python script (MM-PBSA.py).

Toxicity Predictions

The PolyPharma software was used to predict the ligand potential toxic effects using trained Bayesian models and sub-structure similarity search algorithm based on the experimental data from various in vitro assays (Table ). The software implements the effective Bayesian variants, which involved a Laplacian correlation to avoid any scale mismatches and floating point precision issues.[63]

Graphical Visualization

Graphical visualization and analysis have been implemented with the help of UCSF Chimera package[64] (available at http://www.rbvi.ucsf.edu/chimera).

Data Availability

All data generated or analyzed during this study are included in this published article (and its Supporting Information files).