Comprehensive in silico discovery of c-Src tyrosine kinase inhibitors in cancer treatment: A unified approach combining pharmacophore modeling, 3D QSAR, DFT, and molecular dynamics simulation

2.1 Computational studies

Computational studies were performed using HP Intel(R) i5 CPU with 6.0 GB RAM. Schrodinger Suite 2018–4 was employed for pharmacophore modeling and virtual screening. MD simulations was performed using Gromacs 2020.4 on Intel Xenon E3-1245 processor 32 GB RAM, and NVIDIA Quadro P5000 GPU card (Iqbal et al., 2021).

2.2 Dataset generation and alignment

The dataset comprised 34 purine derivatives, recognized as c-Src tyrosine kinase inhibitors, sourced from literature (Sharma, 2016). ChemSketch was employed to generate 2D molecular structures (Table S1 in supplementary material), which were later converted to 3D using Schrodinger's Ligprep module. Energy minimization for low-energy conformers was carried out by OPLS_2005 forcefield (Dixon et al., 2006). Each ligand generated up to 32 stereoisomers, considering all ionization states at pH 7.0 with Epik. Ligands underwent energy minimization until reaching an RMSD of 0.01 Å, forming the basis for modeling studies.

2.3 Pharmacophore model generation

The pharmacophore model was constructed using the Phase software's “Develop Pharmacophore Hypothesis” protocol, utilizing aligned conformations of purine derivatives. (Mohan et al., 2020). The dataset comprised c-Src tyrosine kinase inhibitors categorized into active (pIC50 > 6.40) and inactive (pIC50 > 5.80) sets based on pIC50 values. The pharmacophore model, featuring 4 to 5 sites, was developed with a maximum of 5 sites and a minimum of 4 sites. Phase (v4.0) performed flexible ligand superposition, with the most active compound as the template, considering default settings of 10 conformations per rotatable bond and up to 100 conformers. Conformations of the compounds were varied, and only one alignment per ligand was retained.

Pharmacophore feature like hydrogen-bond acceptor (A), hydrogen-bond donor (D), hydrophobic group (H), negatively charged group (N), positively charged group (P), and aromatic ring (R) were assigned to the molecules using Phase's predefined features. Multiple common pharmacophore hypotheses were generated, scored, and ranked based on active and inactive survival scores.

2.4 Validation of pharmacophore model

The highly-ranked pharmacophore hypothesis underwent validation via PLS (Partial Least Square) analysis. The Phase module was employed to develop an atom-based 3D-QSAR model for predicting the potential c-Src tyrosine kinase inhibitory activity (Rücker et al., 2007). Atom-based QSAR models consider occupancy properties such as hydrogen bond donor, hydrogen bond acceptor, hydrophobic, negative ionic, and positive ionic features. Molecule alignment for the 3D QSAR model utilized Phase (v4.0) shape screening, aligning c-Src tyrosine kinase inhibitors with DDRRR_1 pharmacophore model. The dataset was split into a 70 % training set and a 30 % test set with default parameter settings. The generated QSAR models were ranked based on statistical parameters including R2 (correlation coefficient of the training set), Q2 (correlation coefficient of the test set), SD (Standard Deviation), Pearson-r values, and Y-randomization (Khamouli et al., 2022a; Khamouli et al., 2018). Additional external validation tests included Tropsha and Golbraikh criteria (Golbraikh and Tropsha, 2002; Tropsha, 2010), rm2 metric analysis (Roy et al., 2013), and PLS factor 5 to establish QSAR model robustness and predictiveness. Contour maps facilitated model visualization.

The pharmacophore model's accuracy and specificity in virtual screening were validated using 1000 decoy molecules enriched with 10 active molecules from the Schrödinger database. The Phase module's “Hypothesis Validation Tool” utilized the hypothesis file, decoy, and active datasets to calculate performance parameters. Statistical parameters included RF (Enrichment Factors), RIE (Robust Initial Improvement), BEDROC (Boltzmann Enhanced Discrimination of Receiver Operating Characteristic), AUC (Area Under Curve), and ROC (Receiver Operating Characteristics) were computed to assess and validate the accuracy of the pharmacophome model in virtual screening (Rizzi and Fioni, 2008).

2.5 Database creation and screening

The 3D coordinates of 111,415 compounds were obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/) and processed using Ligprep. Tautomers, ionization states at pH 7.0, and partial atomic charges were computed using OPLS-2005 forcefield. Virtual screening, employing DDRRR_1 pharmacophore and “PhaseLigand Screening” efficiently pinpointed potential c-Src tyrosine kinase binders. Maestro's “Virtual Screening Workflow” was applied with prefiltering based on Lipinski's Rule of 5 (such as Molecular mass < 500 Da, Number of H-bond donor < 5, Number of H-bond acceptor < 10, Octanol-Water partition coefficient (log P) < 5) to eliminate unfit ligands (Friesner et al., 2006). The filtered database underwent a three-step docking process: Glide HTVS (High Throughput Virtual Screening), Glide SP (Standard Precision) of top 10 % hits in HTVS, and Glide XP (Extra Precision) of top 10 % hits in SP. Only compounds with docking scores below −10.5 kcal mol⁻¹ were retained. ADMET properties of the top selected compounds were predicted using Maestro's QikProp and the ProTox-II web server.

2.6 Protein preparation and receptor grid generation

The c-Src tyrosine kinase crystal structure (PDB ID: 1Y57) (Cowan-Jacob et al., 2005) was downloaded from Protein Data Bank (http://www.rcsb.org) having 1.91 Å resolution. “Protein Preparation Wizard” of Schrodinger removed water molecules, added hydrogen atoms, assigned atom charges, and energy-minimized using OPLS-2005 force field keeping RMSD (Root Mean Square Deviation) ≤ 0.3 Å. The glide grid file was generated using Schrödinger, with coordinates extracted from the crystallized ligand MPZ at the active site (X: 18.69 Å, Y: 13.24 Å, Z: −19.81 Å). Docking with the co-crystal ligand in XP mode was performed using Glide, and RMSD was validated (Supplementary Figure S1).

2.7 Molecular docking

Maestro's “Virtual Screening Workflow” was employed to dock and score compounds based on best fitness scores. The workflow offers different levels of docking precision, including HTVS, SP, and XP (Samad et al., 2022). Initially, HTVS calculations were performed, followed by SP and finally XP mode to refine the ligand poses. The ligand-receptor interactions were visualized using Discovery Studio 2020 (Selvaraj et al., 2012).

2.8 DFT calculation

Top screened compounds underwent DFT (Density Functional Theory) calculations using Gaussian 09 software to determine LUMO (Lowest Unoccupied Molecular Orbital), and the HOMO (Highest Occupied Molecular Orbital) quantum properties. B3LYP functional theory with a 6-31G* basis set and hybrid DFT with Becke's 3-parameter exchange potential were applied for HOMO and LUMO calculations.

2.9 In silico ADMET prediction

The top five compounds from virtual screening and docking underwent ADME (Absorption, Distribution, Metabolism, and Excretion) analysis with Maestro's QikProp module. This analysis predicted various ADME properties, including solubility (LogS), partition coefficient (LogPo/w), permeability (PCaco), skin permeability (LogKp), blood–brain partition coefficient (LogBB), human serum albumin binding (LogKhsa), hERG inhibition (loghERG), oral absorption, and Rule of Five and Rule of Three violations. Additionally, toxicity assessment was performed using the ProTox-II web server (https://tox-new.charite.de/protox_II/), predicting mutagenicity, carcinogenicity, immunotoxicity, and hepatic toxicity.

2.10 Molecular dynamics simulation (MDS)

MDS was performed to assess protein–ligand complex stability (Khan et al., 2022). Gromacs 2020.4 was employed to conduct MDS using CHARMM36-all atom force field and TIP3P water. CGenFF (CHARMM General Force Field) was used to generate ligand topology. The complex was placed in a dodecahedron box with 10 Å minimum distance from box boundaries. Neutralization included 0.15 M NaCl and counterions (Na+/Cl-). The system's energy underwent steepest descent minimization (max. 50,000 steps). NVT (Isothermal-Isochoric) and NPT (Isothermal-Isobaric) ensembles were run at 300 K and 1.0 bar, using Berendsen thermostat (Berendsen et al., 1984) and Parrinello-Rahman barostat (Parrinello and Rahman, 1981). A 100 ns production run with 2 fs time-step used the Leap-frog integrator and applied LINCS constraints. RMSD, RMSF (Root Mean Square Fluctuation), and number of hydrogen bonds were analyzed during MDS for system dynamics and stability.

2.11 Principal component analysis (PCA)

Conformational flexibility of c-Src, with and without CID_70144047 was analyzed using Bio3D's PCA (Grant et al., 2021). This involved removing translational and rotational motions, superimposing atomic coordinates onto a reference structure, and calculating a positional covariance matrix to represent atomic coordinate relationships and fluctuations. The diagonalization of this matrix provided eigenvalues and eigenvectors, indicating mean-square fluctuation along each eigenvector. The covariance matrix (C) was calculated to evaluate c-Src's conformational flexibility and its interactions with ligands as follows: $$C_{\mathit{ij}}=\langle \left(x_i-\langle x_i\rangle \right)\left(x_j-\langle x_j\rangle \right)\rangle \left(i,j=1,2,3,\cdots .,3N\right)$$ where, N, x_i/j and < x_i/j > represent the number of Cα-atom, the Cartesian coordinate of the i^th/j^th Cα-atom, and time average of all the conformations respectively.

3.1 Pharmacophore modelling

A dataset of 34 compounds, categorized into active and inactive sets based on experimental data, was used to generate 22 pharmacophore variants. The top five hypotheses, detailed in Table 1, were scored, and pharmacophore hypothesis DDRRR_1 emerged as the best with favorable scores in survival activity (5.723), survival-inactive (2.551), vector (1), volume (0.882), selectivity (2.647), and site matches (10). The alignment of active and inactive compounds with their pharmacophoric features is depicted in Fig. 1A,B.

Close

Fig. 1

(A) Mapping of the active compounds onto the pharmacophore, and (B) Mapping of inactive compounds onto the pharmacphore. Represents pharmacophore model DDRRR_1 with respect to (C) intersite distances in Å, and (D) angles between the pharmacophoric points.

Export to PPT

The identified five-point pharmacophore (DDRRR_1) includes two hydrogen bond donor (D) features related to the amine group at the 9′ position of the purine ring and the amine group between the purine and phenyl rings. Additionally, three aromatic ring (R) features are associated with the purine rings and the phenyl ring on both adjacent sides of the compounds. These chemical features are deemed essential for the biological activity within the dataset, emphasizing their significance for screening and designing potential compounds.

The model's discriminatory ability was evaluated using a test against 1000 decoy molecules from the Schrodinger database. It successfully identified all active compounds in the hit list, showcasing strong performance. The pharmacophore model DDRRR_1 exhibited an RIE value of 18.14, and an EF 1 % value of 101, indicating its superior ranking compared to random distribution. The AUC of the ROC curve, a reliable metric for model evaluation, yielded a favorable value of 1 for DDRRR_1. Additional details on angles and distances among various pharmacophoric features of the selected model are provided in Fig. 1C,D and Table S2 in the supplementary material.

3.2 Atom based 3D- QSAR and validation

The top-scoring hypothesis, DDRRR_1, was subjected to a QSAR study using PLS regression, aligning pharmacophoric features with optimized bioactive ligands. PLS analysis, a chemometric technique, correlated activity (dependent variable) with binary values (independent variable) to develop the QSAR model. The dataset, divided into training (24 molecules) and test (9 molecules) sets, yielded significant statistical outcomes with five factors. The resulting 3D QSAR model was validated on both sets to assess its reliability.

The selection of the best QSAR model was guided by multiple criteria such as R2 and Q2 values, which represent how well the model fits the training data and the predictive power of the generated model respectively. It displayed a robust cross-validated R2 value (>0.6) and a high Q2 value (>0.5). Moreover, the SD values were below 0.3, the root mean square error (RMSE) was minimal, and the F-test value (indicating the overall statistical significance of the model) was substantial. These factors collectively informed the choice of the best QSAR model. The hypothesis exhibited a strong correlation with an R2 value of 0.926 for the training set, accompanied by a high F-test value of 47.9. The Pearson-R value of 0.950 indicated a significant correlation, further reinforced by the corresponding low P-value (4.03 e-10). The SD was 0.085, and the RMSE value was 0.120, as detailed in Table 2. Visual representation of the predicted versus actual pIC50 values for both the training and test sets is presented in Figure S3 (supplementary material), while Table S1 (supplementary material) provides the actual pIC50 and residual pIC50 values for both sets.

According to Galbraith and Tropsha criteria along with r2m, matrix the proposed QSAR model exhibits strong predictive abilities (Table 3) (Golbraikh and Tropsha, 2002; Tropsha, 2010). With a robust R2 and Q2 values of 0.903 and 0.895 for the test set, the model surpasses the established thresholds of 0.5. The R2pred exceeds 0.5, and the Degree of Lack of Relationship falls within the satisfactory range of 0.85 to 1.15 for both K and K', indicating a robust relationship between variables. The r2m value, surpassing 0.5, further enhances the model's credibility. External validation results affirm the high agreement for the atom-based 3D-QSAR models, highlighting their potential in predicting the activities of new derivatives.

Furthermore, the Y-randomization test was conducted through 10 random shuffles of the dependent variable i.e. biological activity (supplementary Table S3). These shuffles, preserving the original independent variable matrix, produced new QSAR models with notably lower Q2 and R2 values, highlighting the importance of the initial dependent variable arrangement for meaningful correlations.

3.3 Contour plot analysis

Contour plot analysis revealed key pharmacophoric requirements in the molecular structure, where blue cubes indicated positive contributions, and red cubes represented negative contributions. For highly active compound 20 (Fig. 2A), favorable hydrogen bond donor groups were indicated by blue cubes at specific positions, such as the NH groups of the purine ring (R8 feature) and the NH between the phenyl ring (R9 feature) and pyrimidine (R10 feature) rings, while red cubes around the bridged sulfone group at the 4″ position suggested potential decreased activity. The withdrawing features (Fig. 2B) highlighted by blue cubes were mainly located around the purine structure due to the amine (NH) group at the 6″ position of the purine ring. Numerous blue cubes around the purine and phenyl rings (Fig. 2C) indicated the dominance of hydrophobic moieties in activity. Enhanced activity was associated with positive ionizable groups at the 3″ and 5″ positions of the benzene structure, while the 4″ position of the benzene ring, indicated by red cubes, hinted at potential decreased activity (Fig. 2D).

Close

Fig. 2

Pictorial representation of the contour plots generated of the active compound (20) in the context of (A) hydrogen bond donor effects, (B) electron withdrawing groups, (C) hydrophobic interactions, and (D) negative ionizable groups.

Export to PPT

3.4 Virtual screening and molecular docking

The pharmacophoric model DDRRR_1 guided a 3D database search for compounds meeting its criteria. Initially, compounds with pharmacophore fitness scores ranging between 2 and 1.3 were selected. Further refinement through Glide’s “Virtual Screening Workflow”, including Lipinski's rule of five and molecular docking, resulted in the identification of five hits (CID_70144047, CID_11291395, CID_22844919, CID_134165918, and CID_60151031) characterized by favorable interactions with the active site of Src tyrosine kinase.

Molecular docking explored the binding modes of the top five hits with Src tyrosine kinase. The binding interaction poses were carefully examined and visualized using the Maestro interface and Discovery Studio Visualizer software (Fig. 3). Analysis of docking scores, and binding free energy highlighted substantial contributions from hydrogen bonding, hydrophobic forces, and π-π stacking interactions in the binding process with nearby residues of Src tyrosine kinase (Table 4).

Close

Fig. 3

2D and 3D interactions of the top hits with the active sites of c-Src tyrosine kinase. (A) CID_70144047, (B) CID_11291395, and (C) CID_22844919, (D) CID_134165918, (E) CID_60151031, and (F) Co-crystalized ligand respectively.

Export to PPT

The top five hits exhibited high scores toward the target, ranging from −10.502 to −11.401 kcal/mol, surpassing the reference ligand (MPZ) with a docking score of −9.488 kcal/mol, indicating stronger binding affinity. Additionally, these hits demonstrated better binding energies (-46.306 to −53.731 kcal/mol) compared to the reference compound, reflecting robust binding affinity towards c-Src tyrosine kinase. For instance, CID_70144047 formed four conventional hydrogen bonds with ASP404, THR338, and MET341 (two bonds), and engaged in hydrophobic interactions with VAL281, LEU393, VAL173, ALA393, and LEU273. Interactions with key residues like SER345, ASP348, ASP404, THR338, MET341, GLU310, and LYS295 confirmed the placement of hits within c-Src tyrosine kinase's binding pocket (Table 4). All the hits, like reference ligand, formed hydrogen bonds with Met341, indicating their potential as competitive inhibitors (Grant et al., 2006).

3.5 Frontier molecular orbitals (FMO)

FMOs (Frontier Molecular Orbitals) encompass a molecule’s HOMO and the LUMO, representing its electron-donating capacity and electron-accepting capacities respectively. The energy gap (ΔE) between HOMO and LUMO, crucial for understanding reactivity, kinetic stability, and chemical properties, was calculated in this study (Karabacak et al., 2012). The corresponding energy values and contour diagrams of these molecular orbitals can be found in Table S4 (supplementary materials) and Fig. 4, respectively. The DFT analysis revealed HOMO and LUMO values from −4.136 to −7.319 and −1.225 to −2.530, respectively, with an energy gap (ΔE) ranging from 1.877 to 4.979 (Table S4), providing insights into electron-binding properties and indicating kinetic stability and chemical reactivity, with smaller gaps suggesting higher reactivity.

Close

Fig. 4

The predicted DFT plots of the HOMO, and LUMO. (A) CID_70144047, (B) CID_11291395, (C) CID_22844919, (D) CID_134165918, and (E) CID_60151031.

Export to PPT

3.6 In silico predicted ADMET

We employed in silico tools, including ADME prediction and QikProp (Ioakimidis et al., 2008), and ProTox-II, to access the ADMET features of the top five compounds. All five compounds adhered to Lipinski's rule of five and the rule of three, indicating favorable drug-like properties. QikProp analysis revealed promising pharmacokinetic parameters, with partition coefficient (QPlogPo/w; in 1.085 to 2.039 range), aqueous solubility (QPlogS; in −4.439 to −3.580 range), human serum albumin (QPksha) binding (in −0.541 to 0.238 range), brain/blood barrier values (in −1.890 to −0.997 range), percentage of human oral absorption (in 74.568 to 82.345 range), and permeability across the gut-blood barrier (QPPCaco; in 179.402 to 410.245 range) were within acceptable limits Table S5 (supplementary materials) (Artursson and Karlsson, 1991). Additionally, compound permeability through Madin-Darby Canine Kidney Cells (QPlogMDCK), indicating transport across distal renal epithelia, was within the desired range. Furthermore, predicted toxicity evaluations suggested non-mutagenic, non-carcinogenic, non-immunotoxic, and non-cytotoxic properties for compounds CID_70144047 and CID_22844919, while the reference compound MPZ showed positive immunotoxicity and hepatoxicity (Table S6 supplementary materials).

3.7 Analysis of molecular dynamic simulation

Based-on the findings of multistep molecular docking, DFT analysis, and ADMET assessment, the compound CID_70144047 was found to be the most promising inhibitor of c-Src. Hence, it was selected for further analysis by molecular dynamics simulation.

3.7.1

3.7.1 Root mean square deviation (RMSD)

The stability of a protein–ligand complex can be assessed through RMSD, which measures the deviation of the complex's structure from its initial conformation (AlAjmi et al., 2018; Jairajpuri et al., 2021). For c-Src, Cα-atoms RMSD remained relatively stable during the 20–100 ns simulation, averaging 4.2 Å (Fig. 5A). Similarly, the RMSD of c-Src-CID_70144047 complex fluctuated initially but reached a stable state as favorable interactions were formed. The average RMSD for the c-Src-CID_70144047 complex was 4.8 Å, which was comparable to c-Src alone, confirming the formation of a stable complex within the acceptable threshold of 2.5 Å.

Close

Fig. 5

Molecular dynamics simulation of c-Src tyrosine kinase in the absence and presence of CID_70144047. (A) Root mean square deviation (RMSD), (B) Root mean square fluctuation (RMSF), (C) Intramolecular hydrogen bonds, and (D) Intermolecular hydrogen bonds.

Export to PPT

3.7.4

3.7.4 Analysis of principal component analysis (PCA)

PCA, a dimensionality reduction method, was employed to analyze the overall motion of c-Src with and without CID_70144047 during simulations. Conformational sampling of c-Src along PC1 (Principal Component 1) and PC2 (Principal Component 2) axes was accessed, red and blue dots representing different conformational states and clusters indicating distinct and energetically favorable regions (Fig. 6). In the absence of any ligand, c-Src occupied a conformational space ranging from −70 to + 50 along PC1 (47.2 %) and −60 to + 40 along PC2 (19.75 %) (Fig. 6A). Conversely, in the presence of CID_70144047, the conformational space covered −70 to + 50 along PC1 (47.84 %) and −50 to + 50 along PC2 (21.75 %) (Fig. 6B). Notably, the first three eigenvalues for c-Src alone and c-Src with CID_70144047 encompassed 80.4 % of the conformational variances (Fig. 6C,D), indicating that the flexibility of c-Src in the presence of CID_70144047 resembles that observed when c-Src is unbound.

Close

Fig. 6

Principal component analysis (PCA) of c-Src tyrosine kinase in the absence and presence of CID_70144047. (A,C) c-Src tyrosine kinase alone, and (B,D) c-Src tyrosine kinase and CID_70144047 complex.

Export to PPT