The study was conducted on a Windows-based laptop equipped with a 1.19 GHz CPU, 8 GB RAM, and a 512 GB SSD. Throughout the research, a reliable power source and a high-speed internet connection were available. InstaDock (Mohammad et al., 2021) and Discovery Studio were utilized for virtual screening and molecular docking, while PyMOL (DeLano, 2002) and Discovery Studio Visualizer (Biovia, 2017) were used for visualization. Various resources and servers such as NCBI, UniProt (UniProt Consortium, 2019), the RCSB Protein Data Bank (Rose et al., 2017), AlphaFold (David et al., 2022), the ZINC database (Irwin & Shoichet, 2005; Irwin et al., 2020), SwissADME (Daina et al., 2017), pkCSM (Pires et al., 2015), Way2Drug PASS (Lagunin et al., 2000), and CarcinoPred-EL (Zhang et al., 2017) were utilized for data retrieval, evaluation, and analysis in this study. GROMACS was used to perform molecular dynamics simulations. The various systematic steps employed in this study are depicted in Figure.

The EphA1 structure was selected and retrieved from UniProt and AlphaFold based on structural completeness, the presence of the kinase domain, mutation-related sequence information, ligand-related information, and structural purity. The structure of EPHA1 was obtained in PDB format from the AlphaFold Protein Structure Database (David et al., 2022; UniProt Consortium, 2019) and further refined in PyMOL (DeLano, 2002) (Figure). Water molecules were removed during receptor preparation.

The ZINC database is a freely available public resource that compiles annotated, commercially available chemical structures (Irwin & Shoichet, 2005; Irwin et al., 2020). Both 2D and 3D structures of compounds can be downloaded from this database, and its interface supports rapid molecular lookup and analogue searching. ZINC maintains a catalogue of more than 230 million ready-to-dock 3D compounds for sale and also provides access to more than 750 million purchasable compounds, enabling rapid analogue searches (Irwin et al., 2020). The ZINC database has expanded substantially, and ZINC20 now contains 1.4 billion chemicals from 310 catalogues and 150 suppliers (Irwin et al., 2020). According to the 90/90/90 standard, 90% of catalogues are updated every 90 days, and 90% of compounds have been verified as available for purchase within the preceding three months (Irwin et al., 2020).

In the present study, the Lipinski rule of five was applied to select 32,901 small molecules from the ZINC database. The Lipinski rule is a cornerstone of medicinal chemistry and provides a set of threshold-based parameters to assess the drug-likeness of chemical compounds. Drug-like molecules generally satisfy the following criteria: molecular weight below 500 Da, fewer than 5 hydrogen-bond donors, fewer than 10 hydrogen-bond acceptors, and a LogP value below 5.

To investigate substrate inhibitor selectivity and understand the mechanism involved, high-throughput screening was conducted using the receptor-ligand docking method for the target AF-P21709-F1. The purpose was to examine the orientation of ligands within the protein's active site cavity. A collection of natural compounds was sourced from the ZINC database, and the ligands were acquired in pdbqt format. Molecular docking of EPHA1 was conducted to examine the bond conformations and binding affinity between the ligands and EPHA1. Regarding the docking process, a defined energy grid box with the dimensions of grid box centre (x = -1.511, y = 1.922, z = -35.295) and size parameters (X=56, Y=63, and Z=75) was used to ensure maximum binding affinity and obtain the best conformational pose for the protein-ligand interactions. The docking algorithm investigated different ligand orientations within the receptor's active site. Each ligand screening used a random seed generator and a scoring function to evaluate binding affinity.

After the virtual screening was completed, the docking results were analysed by examining the output files and log files. The most appropriate docked conformation was chosen for subsequent analysis. The PyMOL software (DeLano, 2002) and Discovery Studio Visualizer (Biovia, 2017) were employed to visualize and analyse the structure of the docked complex. Using the data as a starting point, we next went on to evaluate the compounds' potential as drugs utilising the Swiss-ADME server (http://www.swissadme.ch/) for PAINs analysis, Carcinogenicity using CarcinoPred-EL. Following this, we further evaluated the physicochemical properties, absorption, distribution, metabolism, and excretion of selected compounds using PCKCSM (http://biosig.unimelb.edu.au/pkcsm). Lastly, PASS analysis was performed using way2drug (http://way2drug.com/passonline/), where ligand activities with probabilities of activation greater than 0.71 were considered.

The shortlisted compounds were subjected to pan-assay interference compounds (PAINS) filtering using the SwissADME server to eliminate false positives. Additional physicochemical filters were applied based on total polar surface area (TPSA: 65–145 Å2) and solubility predictions. ADMET properties, including absorption, distribution, metabolism, excretion, and toxicity, were predicted using the pkCSM server. Gastrointestinal absorption, blood-brain barrier permeability, CYP2D6 inhibition, OCT2 substrate status, and AMES toxicity were evaluated. Carcinogenicity prediction was performed using CarcinoPred-EL.

Prediction of Activity Spectra for Substances (PASS) analysis was performed using the Way2Drug PASS Online server to evaluate the potential biological activities of selected compounds. Activities with a probability of activity (Pa) greater than the probability of inactivity (Pi), particularly antineoplastic activity, were considered significant. Detailed interaction analysis of selected protein–ligand complexes was performed using Discovery Studio Visualizer. Hydrogen bonds, hydrophobic interactions, van der Waals interactions, and interactions with key catalytic residues within the ATP-binding pocket were examined. Special emphasis was placed on interactions with critical residues in the kinase domain, such as Lys656 and Asp749.

MD simulations were performed for 100 nanoseconds on both free EphA1 and EphA1 with the compound ZINC12660859 using GROMACS 5.1.2 at 300 Kelvin. The simulations utilized the GROMACS G54a7 force field at the molecular mechanics level. The systems were immersed in a cubic box of water (10 Å) and solvated with the spc216 water model using the gmx solvate module. Energy minimization was performed using 1500 steepest-descent steps. During the equilibration phase, the temperature gradually increased from 0 to 300 Kelvin over 100 picoseconds, while maintaining constant volume and periodic boundary conditions. Trajectory analysis employed various GROMACS utilities, including gmx rmsd for root-mean-square deviation, gmx rmsf for root-mean-square fluctuation, gmx gyrate for radius of gyration, gmx sasa for solvent-accessible surface area, and gmx hydrogen-bonds for intra- and intermolecular hydrogen bonds. Graphs and figures were generated using QtGrace.

We conducted virtual screening of compounds to identify a competitive inhibitor for EphA1. Log files and out-files were created, containing affinity scores and docked poses for each compound in the directory. Based on binding affinities, docking scores, and binding poses, we filtered compounds using the log and output files. Our screening identified several natural compounds with favourable binding affinity scores for EphA1's binding pocket, indicating their potential as EphA1 inhibitors. Out of 32,901 compounds that were screened from the output, 320 were found to have significant EphA1 binding affinity scores ranging from -11.2 to -9.4 (Table). Among the screened compounds, ZINC03845566 exhibited the highest binding affinity (-11.2 kcal/mol), while ZINC12660859 and ZINC12661003 demonstrated binding energies of -9.9 and -9.7 kcal/mol, respectively. Although several compounds showed slightly stronger docking scores, subsequent pharmacokinetic filtering and biological activity prediction narrowed the focus to ZINC12660859 and ZINC12661003 due to their favorable drug-like properties and predicted antineoplastic potential. The observed binding energies suggest strong stabilization within the EphA1 catalytic pocket, indicating potential competitive inhibition at the ATP-binding site.

In the initial step, Discovery Studio was used to generate SMILES strings for the top hits. These SMILES IDs were subsequently subjected to PAINS analysis using SwissADME. The PAINS analysis aimed to identify pan-assay interference compounds, which are multi-target ligands that lack specificity for a single target. Eliminating such compounds is crucial because they can lead to unintended regulatory effects. Following the PAINS analysis, additional filters were applied to refine the compound selection process. Compounds violating Lipinski’s rule of five were removed. Furthermore, compounds with a total polar surface area (TPSA) outside the range of 65 to 145 were excluded. Compounds exhibiting poor solubility according to Silicos-IT Solubility (mol/l), Ali Solubility (mol/l), and ESOL Solubility (mol/l) were also eliminated from consideration. These filtering criteria helped in identifying the most promising compounds for further analysis and evaluation (Table).

After applying the filters, 24 compounds remained for further analysis. Using pkCSM, these substances were screened to determine their ADMET parameters. Remarkably, 20 compounds successfully met the ADMET parameters, indicating favourable characteristics for drug development. The physicochemical parameters evaluated for these compounds were found to fall within the range deemed suitable for drug candidacy, as outlined in Table. This observation indicates that the selected compounds possess desirable properties, making them promising candidates for further investigation and drug development.

In addition, a PASS analysis was performed on the 20 selected compounds using SMILES strings via the way-2-drugs-PASS Online tool. This analysis aimed to identify compounds with anti-cancer activity, thereby offering therapeutic leads for EphA1 targeting. After careful analysis, two compounds were chosen based on their demonstrated anti-cancer properties. The activities and corresponding Pa (probability of activity) and Pi (probability of inactivity) values for these selected compounds are presented in Table. These compounds exhibit promising characteristics and warrant further exploration as potential anti-cancer agents targeting EphA1.

To analyse the binding modes and interaction patterns of the two selected compounds, Discovery Studio was used. A set of 18 docked conformers was produced for each of the two natural compounds using the output files. The analysis focused on examining the interacting residues within the compounds. The Discovery Studio software was used to identify and visualize hydrogen bonding and other interactions between the compounds and EphA1. Notably, residues within the kinase domain of EphA1 were observed to play a significant role in mediating numerous interactions with the compounds. This analysis provides valuable insights into the binding mechanism and potential interactions between the selected compounds and EphA1. Residues of EphA1's kinase domain have been shown to offer a considerable number of interactions, including Asp749, Gly769, Gly633, Asp767, Glu673, Lys656, Arg753, Asn754, Val638, Ile630, Leu,770, Leu756, Ala709, Gly631, Ala654, Phe704, Met705, Gly708, Phe635, Glu634, Ser766, Leu686, Gly708 (Table). Each compound interacts with ATP binding site or is closely associated with ATP binding sites. ATP binding sites, Lys656, as of EphA1, are inhibited either by a hydrogen bond or a strong bond in each compound, and active site Asp749 of EphA1 closely forms bonds with compounds.

Compound ZINC12660859 is forming hydrogen bonds to Asp749, Gly769, Gly633 and several other interactions to656, Arg753, Asn754, Asp767, Val638, Ile630, Leu,770, Leu756, Ala709, Gly631, Ala654, Phe704, Met705, Gly708, Phe635, Glu634. Compound ZINC12661003 is interacting to GlAsp767, Glu673 via Hydrogen bonding, and Lys656, Ser766, Leu686, Ala709, Gly708, Val638, Met705, Phe704, Ala654, Leu756, Ile630 are participating in other interactions (Figure). A significant number of interactions, including Hydrogen bonding with the active site Asp 749 and Vander wall’s interaction with ATP binding site Lys656, formed between EphA1 and ZINC12660859, demonstrating a significant affinity for binding and its further implications as a medicinal compound targeting EphA1 (Figure).

An effective way to examine the atomic-level structural specifics and dynamic behaviour of protein-ligand complexes is to use the MD simulation approach. To acquire a thorough understanding of the stability and dynamics of the EphA1-ZINC12660859 complex, we conducted all-atom MD simulations for a duration of 100 ns on both the EphA1-ZINC12660859 complex and the apo EphA1 protein. We aimed to uncover numerous systematic and structural factors that affect the complex's behaviour through rigorous research into intermolecular interactions, conformational changes, hydrogen-bonding patterns, and solvent effects. Our research contributes to the field of drug development and targeted therapies by providing a deeper understanding of atomic-level protein-ligand interactions.

Root-Mean-Square Deviation

The root-mean-square deviation (RMSD) is a well-known metric frequently used to analyse structural heterogeneity in proteins. It is crucial as a tool for assessing and analysing the dynamic behaviour and conformational changes of protein structures. RMSD provides vital insights into protein structural deviation and dynamics by quantifying the discrepancy between the atomic positions of a target structure and those of a reference structure. The structural dynamics of the protein EphA1 before and after ligand binding revealed low mean RMSD values of 0.20 nm and 0.26 nm for EphA1 and EphA1-ZINC12660859, respectively (Figure A). During the 100 ns simulation, both systems exhibited stable trajectories, with RMSD values approaching equilibrium. Analysis of the probability distribution function (PDF) revealed a significant increase in EphA1 stability, as evidenced by a pronounced peak in the distribution during compound binding. This finding indicates that the binding of ZINC12660859 positively influenced the overall stability of the complex.

Root-Mean Square Fluctuation

The root mean square fluctuation (RMSF) analysis has been widely employed to investigate protein dynamics during MD simulations. It offers valuable information on the influence of ligand binding on protein vibrations. Through the examination of the RMSF plot, we examined the residual dynamics of EphA1 prior to and following ligand binding. Notably, the protein-ligand complex exhibited a significant decrease in RMSF fluctuations, indicating increased stability upon binding of ZINC12660859 (Figure B). The RMSF analysis revealed a significant decrease in residual fluctuation within the POT1 protein. These findings are further corroborated by PDF analysis, which clearly indicates a pronounced decrease in residual flexibility following complex formation.

Radius of gyration

The radius of gyration (Rg) is a direct indicator of a protein's tertiary structure and overall conformation. It is commonly used to assess the compactness and folding behaviour of proteins. In our study, we investigated the stability of both EphA1 and the EphA1-ZINC12660859 complex by analysing their respective Rg values. The average Rg for EphA1 in its apo form was 1.92 nm, whereas for the EphA1-ZINC12660859 complex, it was found to be 1.93 nm (Figure A). The Rg plot revealed a slight increase of approximately 0.01 nm in Rg upon binding of ZINC12660859 to EphA1. This minor increase is likely due to the compound's packing effect. Remarkably, the introduction of ZINC12660859 did not result in any notable alterations to the structure of EphA1. Throughout the simulation trajectory, EphA1 maintained a constant equilibrium of Rg, indicating the enduring stability of the complex.

Solvent Accessible Surface Area

The solvent-accessible surface area (SASA) of a protein refers to the region on its surface that is accessible to the surrounding solvent. Analysing SASA is an important method for examining protein folding and stability. Upon analysis of the SASA plot, no notable alterations in SASA values were observed, indicating the stability of EphA1 interactions with ZINC12660859 (Figure B). The SASA values for the EphA1 and EphA1-ZINC12660859 complexes were 140.42 nm2 and 139.62 nm2, respectively. The distribution of surface-area values exhibited a consistent pattern across both systems. Notably, there was a marginal reduction in the average surface area upon binding of the compound to EphA1. This reduction suggests the formation of stable interactions or binding between the protein and ligand.

Hydrogen Bond Analysis

The creation of hydrogen bonds (H-bonds) is critical in the kinetics of protein folding. Protein conformational alterations are primarily governed by the disruption and formation of hydrogen bonds. In this study, we analysed the temporal progression of hydrogen bonds to assess the integrity of intramolecular bonding within the EphA1-ZINC12660859 complex. The plot we obtained indicated negligible alterations in the number of hydrogen bonds in EphA1 upon complex formation with ZINC12660859 (Figure). In the context of EphA1 intramolecular interactions, we observed that the average number of hydrogen bonds formed was 189 before binding of ZINC12660859, and this count increased to 190 after binding. A minor rise in the number of hydrogen bonds was noted, which can be attributed to the heightened compactness resulting from the binding of the ligand (Figure A). Analysis of the PDFs for both systems confirmed the consistent stability of intramolecular hydrogen bonds. From the generated plots, it can be inferred that the intramolecular hydrogen bonds in EphA1 remained stable throughout the simulation, even after compound binding. Throughout the simulation, ZINC12660859 was discovered to bind into the specified pocket of EphA1, creating 1–2 stable hydrogen bonds with little variation (Figure B). There were further occasions where ZINC12660859 showed 3–4 hydrogen bonds with larger variations. ZINC12660859 remained at its original docking site on EphA1, consistent with our predictions based on the observed intermolecular hydrogen bonding. These hydrogen bonds play an important role in the stability of the complex structure, preventing the compound from migrating or dissociating from its binding site on EphA1 during the simulation.