Abstract
Objectives
HER2 overexpression is almost invariably associated with advanced breast cancer disease and poor prognosis, hence its extensive review. This study therefore aims to discover and analyze potential HER2 inhibitors through computational methods to advance drug discovery and optimization.
Methodology
A ligand-based virtual screening (LBVS) approach was employed to screen compounds from the ChEMBL database. From 8900 initial matches, 39 candidate compounds were selected based on structural similarity and ADME properties. Molecular docking was performed to assess binding affinity with HER2, followed by molecular dynamics (MD) simulations to evaluate complex stability. Additionally, a QSAR (quantitative structure-activity relationship) model was established to elucidate key structural features influencing inhibitory activity.
Results
Five lead compounds were prioritized based on strong docking scores (<−8.4 kcal/mol). Among them, compound 2048788 (−11.0 kcal/mol, predicted pIC50 ≈ 8.6) and compound 3956509 (pIC50 ≈ 8.4) showed superior binding affinity and pharmacokinetic properties compared to FDA-approved drugs (doxorubicin, letrozole, lanatuzumab). MD simulations confirmed complex stability. The initial QSAR model showed low predictive power (R2 = 0.18, RMSE = 1.19), but after feature selection, performance improved significantly (RMSE = 0.57). Key positive contributors included hydrogen bond donor count (r = 0.63), lipophilicity (LogP, r = 0.60), and sp3 carbon fraction (r =0.60), while excessive polarity and aromaticity reduced activity. Compounds within the 450-500 Da molecular weight range exhibited the highest activity (pIC50 = 8.0-8.6).
Conclusion
This study integrated virtual screening, docking, MD simulations, and QSAR modeling to identify compound 2048788 as a highly promising HER2 inhibitor. These findings provide a strong foundation for further optimization and the preclinical development of targeted HER2 therapies.
Introduction
Breast cancer affects millions worldwide and continues to impose a significant burden on global health. Human epidermal growth factor receptor 2 (HER2)-positive breast cancer represents one of the most aggressive molecular subtypes, characterized by high recurrence rates, rapid disease progression, and poor prognosis.
1
HER2 is a transmembrane tyrosine kinase receptor within the epidermal growth factor receptor (EGFR/HER) family, acting as a key regulator of cellular processes including differentiation, proliferation, apoptosis, and survival.
2
Approximately 15%–20% of breast cancer cases are associated with HER2 gene amplification or overexpression, leading to increased activation of oncogenic pathways such as PI3K/AKT and MAPK, which in turn drive malignant cell proliferation and metastatic spread.3,4 HER2 is widely recognized as a viable therapeutic target, consistent with its critical role in tumor progression.
5
Multiple HER2-targeted therapies have demonstrated clinical success, including small-molecule tyrosine kinase inhibitors, antibody-drug conjugates, and monoclonal antibodies (eg, trastuzumab).
6
However, resistance frequently arises from secondary mutations, compensatory signaling, or conformational changes in the receptor. An example of functional group differentiation reveals that these substructures exhibit receptor selectivity through halide, hydroxyl, and methoxy groups. By modifying hydrophobic and electron interactions, quinoline and benzimidazole skeletons have demonstrated potent binding capabilities.
7
Figure 1 schematically depicts these two skeletons binding to HER2 as inhibitors and substrates; skeletal optimization has recently been further investigated computationally. One study demonstrated the inhibitory activity of novel heterocyclic compounds via molecular docking and molecular dynamics studies, revealing stable binding interactions with hinge residues.
8
The addition of hydrophobic substituents to pyrimidine scaffolds and fused heterocyclic derivatives enhances docking scores.
9
Structurally optimized nitrogen-containing heterocyclic compounds can yield leukotriene receptor antagonists and drug analogs.
10
Quinazoline analogues serve as scaffold-optimized HER2 inhibitors exhibiting superior performance compared to clinical reference drugs, including enhanced efficacy.
8
Collectively, these findings underscore the critical role of scaffold design, functional groups, and enzymatic interactions on stability, binding strength (force constants), and inhibitory properties. The Key Lead Scaffolds Along With Their Mechanistic Roles of Major EGFR/HER2 Inhibitor Classes Used in Drug Discovery and Development
Ligand-based virtual screening (LBVS), molecular docking, and molecular dynamics simulations have significantly accelerated the development of effective HER2 inhibitors. 11 These computational methods aid in predicting the binding affinity of potential drugs toward their targets, the stability of interactions, and molecular changes over time. For instance, molecular dynamics simulations provide insights into the stability and flexibility of drug-receptor complexes, while molecular docking estimates binding energies and visualizes how drugs interact with receptors. Conversely, LBVS leverages the known structures and properties of existing ligands to identify novel candidate drugs with similar characteristics. 12
By combining these computational techniques, scientists can rapidly identify promising compounds and ascertain which chemical scaffolds are most likely to yield drug-like molecules. In our study, we employed an integrated computational approach to screen the ChEMBL compound library, focusing on interactions between diverse molecules and HER2 protein. We focused specifically on two aspects: the core structure of the molecules and the functionality of their functional groups, including π-π stacking interactions, hydrophobic regions, hydrogen bond donors or acceptors, and electronic effects. 13 Through methods such as ligand-based virtual screening, molecular docking, and molecular dynamics simulations, we were able to comprehensively evaluate the binding strength and stability of compounds to the HER2 protein. Upon identifying the most promising compounds, we compared them to currently clinically approved inhibitors. This enables us to apply rational drug design principles, transforming promising molecular frameworks into potent lead candidates. Ultimately, our research not only uncovers novel substances with potential anti-breast cancer properties but also charts a promising pathway for developing more effective breast cancer therapies.
Materials and Methods
Software and Hardware
This study was conducted over a period of 8 months, from August 2024 to March 2025. This study employed a suite of straightforward software and databases. BIOVIA Discovery Studio 2.5.5 (Accelrys Inc, USA, https://www.accelrys.com/) was utilized for molecular modelling and simulation. The Protein Data Bank (https://www.rcsb.org/) contained the 3D crystal structure of the human HER2 target. The cheminformatics component of this study employed virtual ligand screening via the ChEMBL database (https://www.ebi.ac.uk/chembl/) to identify analogous compounds. SwissADME (SIB Swiss Institute of Bioinformatics (https://www.swissadme.ch/), was employed for ADME analysis. Autodock Vina software, running on the Windows 10 operating system, was utilized for docking studies. DRUGBANK (https://go.drugbank.com/) provided information on reference ligands representing three FDA-approved anticancer drugs: doxorubicin, letrozole, and lenalidomide. Molecular properties of the compounds were analyzed using the Molinspiration cheminformatics online tool (https://www.molinspiration.com/).
Human Epidermal Growth Factor Receptor2 (HER2) Preparation
Since HER2 is believed to promote cancer, we selected it as our target protein. Clinically, approximately 15% to 30% of breast tumors exhibit elevated HER2 levels, which are closely associated with tumor aggressiveness, poor patient prognosis, and treatment resistance. Targeting this protein is crucial, as its kinase domain regulates the signaling pathways responsible for uncontrolled cell proliferation. The HER2 crystal structure (PDB ID: 3PP0) clearly reveals the kinase domain with its co-crystallized ligand. This enables us to further elucidate its binding site. To validate the molecular docking approach, we re-docked the co-crystallized ligand to the HER2 active site. The predicted binding sites were compared with experimentally determined conformations. During receptor preparation, all non-essential molecules, including water, ions, and crystal ligands, were removed. The cleaned structures were saved in PDBQT format for use with AutoDock Vina (version 1.5.7) for molecular docking. The amino acid sequence for HER2 was obtained from BIOVIA Discovery Studio 2.5.5. Receptor refinement was performed using SWISS-PdbViewer 4.1.0 (https://spdbv.vital-it.ch/) via energy minimization to relax the receptor structure and resolve spatial conflicts. The coordinates of the active site grid box were meticulously defined to encompass the binding pocket. We calculated the root mean square deviation (RMSD), with results indicating RMSD values ≤ 2.0 Å. This validated the reliability and consistency of the docking process. Despite the availability of multiple HER2-targeted therapies, challenges such as acquired resistance and limited efficacy persist. This underscores the ongoing need for novel inhibitors and supports the focus of this computational study on HER2.
Ligand-Based Virtual Screening
This study employed LBVS (a sophisticated computational strategy) to identify potential inhibitors of HER2. Researchers extracted approximately three million small molecules from the chemical database ChEMBL (https://www.ebi.ac.uk/chembl/) to construct a ligand library. This database was selected due to its use of robust algorithms and support for ligand-based virtual screening. This screening yielded diverse bioactive compounds, ensuring coverage of a broad spectrum of potential inhibitors and drug-like molecules suitable for HER2-targeted inhibition. 14 The screening commenced with the established HER2-targeting co-crystalline natural ligand 2-{2-[4-({5-chloro-6-[3-(trifluoromethyl)phenoxy] pyridin-3-yl}amino)-5H-pyrrolo[3,2-d] pyrimidin-5-yl]ethoxy}ethanol, present in the crystal structure of the HER2 receptor. Researchers collected Simplified Molecular Input Line Entry System (SMILES) data for three FDA-approved anticancer drugs (doxorubicin, letrozole, and lanatuzumab) from the drug library database (https://go.drugbank.com/). 15 These comparator drugs have been extensively studied for their interactions with HER2, providing reference points for binding affinity and interactions. To ensure uniqueness, Open Babel employed InChIKey matching to remove duplicate structures, eliminating any redundant information from the compound pool prior to further analysis.
Drug-likeness Studies
After candidate molecules were identified, Lipinski’s five rules (Ro5) were applied to screen these candidates for similarity to known drugs, thereby predicting oral bioavailability and various properties of the candidate compounds. This process assessed the suitability of candidate drug selection and eliminated compounds exhibiting poor pharmacokinetic properties. These variables include partition coefficient, lipophilicity (log P ≤ 5), molecular weight (MW ≤500 g/mol), hydrogen bond acceptors (HBA < 10), hydrogen bond donors (HBD < 5), topologically polar surface area (TPSA ≤ 140), and rotatable bond count (RB ≤ 10). This systematic screening approach ensures candidate compounds possess desirable pharmacokinetic properties and physicochemical characteristics. Furthermore, the Veber rule, Egan rule, and Ghose filter were applied to evaluate rotatable bonds and topological polar surface area, confirming compounds meet the parameters required for favorable permeability and solubility. Compounds were evaluated for gastrointestinal absorption, lipophilicity, impact on blood-brain barrier permeability, water solubility, skin permeability, CYP450 inhibition, lead similarity, synthetic accessibility, and physicochemical properties. Gastrointestinal and brain permeability was predicted using the boiled egg model. 16 LogS (water solubility) was predicted to forecast how solubility impacts oral bioavailability. The Molinspiration cheminformatics tool utilized the SMILES model.
Molecular Docking Study
We conducted molecular docking studies to predict the binding efficacy of screened compounds to the HER-2 receptor. Using Discovery Studio 2.5.5 software, we prepared the three-dimensional crystal structure of the receptor for docking by removing water molecules, adding missing atoms, and optimizing the shape. Molecular docking was performed using Autodock Vina software. Polar hydrogen atoms were added to the protein, and charges were assigned. The grid x, y, and z coordinates were set to 26 × 26 × 26 Å), with the ligand screening level set to 100. Within the active site of the HER2 receptor, co-crystallized ligands were re-docked to validate the docking methodology. The predicted conformations were subsequently compared with experimentally determined crystal structures. As the re-docked ligands aligned with the experimental binding orientation and exhibited a root mean square deviation (RMSD) below 2.0 Å within the acceptable threshold range for docking validation. For the comparison, the binding affinities of reference anti-cancer drugs (doxorubicin, letrozole, and lanatitinib) were calculated. Among compounds with documented binding affinities to screened ligands, those exhibiting strong binding affinity to HER2 were identified using threshold energy filtering.
Molecular Dynamics Simulation
In order to understand the behaviors and stability of protein-ligand complexes in binding, we employed molecular dynamics simulations. These simulations enabled us to observe the natural motions and flexibility of proteins, particularly within the most strongly bound complexes. For this purpose, we utilized the iMODS server (https://imods.iqf.csic.es/), a tool designed to explore how proteins move and transform between different conformations. iMODS provided an intuitive understanding of the diverse ways proteins alter their structure, even plotting actual pathways between similar protein conformations. To assess protein stability, we examined several internal characteristics: the degree of deformation in different protein segments, eigenvalues (correlated with structural rigidity), B-factors (indicating atomic mobility), and the covariance matrix (revealing how different protein segments move in concert). All these measurements stem from a technique called Normal Mode Analysis (NMA). We summarized the protein’s stability using graphs and models illustrating the flexibility of the main chain, B-factor values, eigenvalues, the covariance matrix, and the elastic network model. In summary, these simulations provide valuable insights into how HER-2 changes when ligands bind to it. This information helps us identify ligands that not only bind well but also promote protein stability and motion, making them promising candidates for future therapeutic development.
QSAR Analysis
The dataset of HER2 inhibitors and their biological activity values (IC50) was sourced from the ChEMBL database. 17 The original dataset was curated by removing duplicate entries, inorganic molecules, and compounds lacking experimental activity data. The final dataset was standardized into SMILES format for subsequent analysis. Molecular descriptors and fingerprints, encompassing physicochemical, topological, and structural parameters, were computed using PaDEL-Descriptor 18 and RDKit 19 within an Ubuntu environment. Prior to modelling, descriptors with constant values or high collinearity (r > 0.9) were excluded, and the dataset was normalized using Min-Max scaling. The curated data were partitioned into a training set (80%) and a test set (20%). QSAR models were constructed using machine learning algorithms implemented in Scikit-learn (Pedregosa et al, 2011), including multiple linear regression (MLR), random forests (RF), and support vector regression (SVR). Model hyperparameters were optimized via grid search combined with 10-fold cross-validation. Model performance was evaluated based on statistical metrics including R2, Q2, RMSE, and MAE, whilst external validation was conducted using the test set according to OECD guidelines (OECD, 2007). 20 The optimized QSAR models were subsequently applied to the ChEMBL virtual screening dataset to identify compounds predicted to possess high inhibitory potential against HER2 for subsequent molecular docking and molecular dynamics simulations.
Results
The Energy Minimization of the HER-2 Model in kJ/moL Before Performing Molecular Docking
The Crystal Structure of HER2 Receptor (PDB ID: 3PP0) Displaying the Kinase Domain With Its Native Ligand Bound to the Active Site. The Co-crystallized Ligand Shown in Stick Representation and Highlighting Its Orientation and Interactions Within the Binding Pocket which Served as a Reference Structure for Docking Validation and Characterization of the Architecture of HER2 Active Site
The Molecular Properties, Medicinal Significance and Interactions of Selected HER2 Inhibitors Along With Their ChEMBL IDs, Similarity Scores and Binding Affinities
The Drug Likeness Filters, Lipophilicity, Pharmacokinetics and Physiochemical Properties of Top Five Virtually Screened Hits
The compound (CHEMBL ID: 1923019) exhibited a stable docking configuration and achieved a notable binding affinity of −8.5 kcal/mol. There were several critical residues involved in hydrophobic interactions, including VAL773, VAL839, ARG840, ASP871, ILE872, ASP950, MET953, PRO967, MET979, and PHE986. Conventional hydrogen bonding was also observed with ARG985, ARG978, and LYS957, with bond lengths ranging from 2Å. Moreover, LEU869 formed a sigma interaction (3.98 Å), GLU975 was involved in halogen bonding (3.24 Å), and TYR772 and ILE954 contributed to π-π stacking at distances of 4.30 Å and 5.87 Å, respectively (Figure 3). The ADME profile of this chemical compound revealed moderate gastrointestinal absorption and a high kinase inhibition rate of 86.7%. Structurally, it consists of 35 heavy atoms, eight hydrogen bond acceptors, and three hydrogen bond donors. With an iLOGP value of 3.11, the compound exhibited favorable lipophilicity, which supports effective systemic absorption. Although it exhibited moderate solubility and a bioavailability score of 0.55, its negative Log Kp value indicates low skin permeability for transdermal applications (see Table 3). The compound (CHEMBL ID: 2048788) demonstrated reliable docking interactions and a binding affinity of −8.6 kcal/mol. It was also scored at 0.55 for bioavailability, indicating favorable oral drug-likeness. The compound showed van der Waals interactions with LYS765, VAL773, VAL839, ARG840, LEU869, ASP950, PRO967, and GLU971 (2.08-3.35 Å) and formed hydrogen bonds with ASN708, TYR772, ALA756, and LYS957. π-Cation interactions were established with ARG978 and ARG985 (4.17-4.89 Å), and a π-π stacking interaction was observed with ILE954 (4.83 Å) (Figure 3). Functionally, the molecule inhibits kinase activity by 73.3%, confirming its therapeutic relevance as a HER2 inhibitor. It consisted of 34 heavy atoms, five hydrogen bond acceptors, and three donors. The ADME results indicated a high GI tract absorption, with a synthetic accessibility score of 3.53. With an iLOGP of 3.86, the compound was highly lipophilic, favoring better permeability and absorption (see Table 3). Molecular Docking Analysis Interactions of Top Selected Ligands With the HER2 Receptor. The Docking Poses are Illustrated for (A) ChEMBL1923019, (B) ChEMBL2048788, (C) ChEMBL3956509, (D) ChEMBL3903447, and (E) ChEMBL3327013, Depicting Their Orientations and Key Interactions Within the Active Site of the HER2 Kinase Domain. Additionally, Bioavailability Radar Plots, the BOILED-Egg Models and Enzyme Inhibition Profiles of These Compounds are Presented to Highlight Their Pharmacokinetic Properties and Drug-likeness
Similarly, compound CHEMBL3956509 showed potent binding affinity (−9.0 kcal/mol) and bioavailability (0.55), reflecting excellent drug-likeness potential. Noteworthy molecular interactions included a vander Waals bond with ASP769 (2.75 Å), hydrogen bonding with ASP950 and LYS957 (4.61 Å and 2.37 Å, respectively) and π-cation interactions with LYS765, ASP871, GLU975 and ARG978 (ranging from 3.88 to 4.92 Å). Additionally, a π-π stacking interaction was observed with TYR772 (4.23 Å), and ILE954 was involved in a π-alkyl interaction (5.15 Å) (Figure 3). This molecule is comprised of 34 heavy atoms, six acceptors and two hydrogen bond donors. ADME profiling revealed high GI absorption and an iLOGP of 3.33, indicating strong lipophilic behavior favorable to membrane permeability and drug distribution (Table 3). The compound (CHEMBL ID: 3903447) also achieved the highest binding affinity of −9.1 kcal/mol, accompanied by a bioavailability score of 0.55, which reinforces its potential for oral administration. The compound formed hydrogen bonds with LYS957 (3.04 Å), π-π stacking with TYR772 (4.23 Å), and π-cation interactions with multiple residues (LYS769, ASP950, GLU975, ASP871, and ARG978) (3.47-4.92 Å) (Figure 3). Functionally, this molecule demonstrated kinase inhibition (40.0%) and was comprised of 35 heavy atoms. It contained nine hydrogen bond acceptors and two donors, and it displayed high GI absorption and synthetic accessibility of 3.36. Its iLOGP value of 2.95 confirmed lipophilicity (see Table 3). The compound (CHEMBL ID: 3327013) exhibited a binding affinity of −8.7 kcal/mol and a bioavailability score of 0.50, indicating its potential as a drug candidate. The key binding features included hydrogen bonds with LYS765, LYS957, and ARG985 (2.07 - 2.36 Å), halogen bonding with GLU971 (3.17 Å), and a π-donor hydrogen bond with TYR772 (2.21 Å). Additionally, alkyl interactions were observed with ARG840, LEU869, and ILE954 (3.81 - 5.19 Å), and π–anion interactions were noted with ASP871 and GLU975 (3.97 - 4.15 Å) (Figure 3). Pharmacologically, it inhibited kinase function by 53.3% and other enzymatic activities by 26.7%. With 31 heavy atoms, the compound featured eight acceptors and two hydrogen bond donors. ADME prediction indicated efficient gastrointestinal absorption and strong lipophilicity (iLOGP = 2.95) (Table 3).
The interactions between the targeted ligands and the HER2 active site were comparable to those of the reference compounds, which were FDA-approved anti-cancer drugs. Molecular interaction analysis showed that doxorubicin formed conventional hydrogen bonds with ASP871 (3.18 Å) and LYS957 (1.99 Å and 2.48 Å). ARG978 exhibited an unfavorable donor–donor interaction (2.71 Å), while ARG840 and LEU869 engaged in alkyl interactions (4.42 Å and 5.10 Å, respectively). LYS765 formed π–cation interactions (4.36 Å and 4.72 Å), and Gln975 participated in a π–anion interaction (3.71 Å). Letrozole was bound to the receptor via van der Waals interactions with ASN708, LYS765, GLY776, ILE872, and ARG840. TYR772 engaged in a π-π T-shaped interaction (4.69 Å), while VAL773 and VAL839 exhibited π-alkyl interactions (5.40 Å and 4.90 Å, respectively). TYR835 formed a conventional hydrogen bond (2.52 Å), and LEU869 was involved in a π-σ interaction (Figure 4). The Docked Poses of Three FDA-Approved Drugs and Top Five Screened Hits Within the HER2 Active Site. (A) Zoomed View of the Binding Orientations of the Control Drugs (Doxorubicin, Letrozole, and Neratinib) Within the HER2 Active Pocket. (B) Zoomed View of the Docking Conformations of the Five Screened Compounds, Illustrating Their Interactions and Orientations in Comparison to the Reference Controls
Similarly, neratinib showed binding through van der Waals interactions with ASN708, VAL839, VAL773, LEU869, GLY776, GLU971, SER974, ARG878, ARG840, and PRO967. ARG985 formed three conventional hydrogen bonds (3.85 Å, 4.86 Å, and 6.30 Å), while ALA706 and TYR772 formed carbon-hydrogen bonds (3.54 Å and 3.57 Å). LYS765 was involved in a π-cation interaction, and ASP871 formed two π-anion interactions (3.87 Å and 5.10 Å). GLU975 was also engaged in a π-anion interaction (5.20 Å) (Figure 4). Compared to the reference anticancer drugs doxorubicin, letrozole, and neratinib, the five selected ligands demonstrated significantly stronger binding affinities (−8.5 to −9.1 kcal/mol) and more favorable interaction profiles with the HER2 active site. Doxorubicin exhibited limited hydrogen bonding and even unfavorable donor–donor interactions (eg, with ARG978), letrozole primarily relied on weak van der Waals and π–alkyl interactions, and neratinib formed conventional interactions but lacked structural diversity and strength in its contacts. By contrast, the selected ligands, particularly CHEMBL1923019, CHEMBL2048788, and CHEMBL3956509, engage in multiple, diverse non-covalent interactions. These include strong hydrogen bonds with LYS957 and ARG978, π-π stacking with TYR772, halogen bonding with GLU975, and π-cation interactions with key residues such as ARG985. These ligands also outperformed the reference drugs in terms of pharmacokinetic properties, demonstrating higher gastrointestinal absorption and favorable iLOGP values (2.95-3.86). They also exhibited consistent bioavailability scores (0.50 - 0.55). Their kinase inhibition rates (up to 86.7%) also underscore their therapeutic superiority, indicating enhanced stability, binding efficiency, and drug-likeness compared to standard HER2-targeted treatments.
The conformation with the best docking energy was selected for further analysis, and the corresponding PDB file was saved to facilitate the subsequent molecular dynamics simulation of the HER2 receptor and five ligand complexes. The analysis of these complexes was conducted utilizing the iMODS approach. The deformability data revealed that high flexibility was primarily present in the hinge regions of the receptor, which are known to be flexible. The elevated points on the deformability graph are indicative of regions such as loops, terminal residues, and solvent-exposed areas, which have the capacity to undergo substantial conformational changes. These flexible regions were typically implicated in ligand binding or functional motions. The regions in the plot exhibited a flat profile corresponding to rigid areas, such as alpha-helices, beta-sheets, or tightly packed core domains. It is noteworthy that sharp peaks emerged at specific indices, with considerable flexibility observed around residues 300 and 350 in all complexes. The overall deformability remains moderate, suggesting that the receptor maintains a generally stable structure, with certain regions exhibiting functional flexibility. The comparable outcomes witnessed across the five receptor-ligand complexes can be ascribed to the conserved binding site, the chemical similarity of the ligands, and their capacity to stabilize the receptor in a comparable conformation, inducing analogous conformational changes and dynamic behavior in all complexes (Figure 5). The Normal Mode Analysis (NMA) B-factor profiles, which represent the mobility and amplitude of atomic displacements, demonstrated that receptor-ligand complexes exhibited analogous vibrational patterns. The complex eigenvalues were consistently low, with a mean value of 8.218923e-05, indicative of low-frequency vibrational modes. Low eigenvalues indicated the occurrence of global or collective motions, such as hinge bending or breathing at the active site, with flexibility in specific domains. The observed low-energy conformational modes imply that ligand binding stabilizes these motions, enhancing the functionality of the receptor and possibly facilitating favorable allosteric inhibition or selective conformational targeting (Figure 5). Molecular Dynamics Versus Normal Mode Analysis (NMA) of HER2-Ligand Complexes. Rows Denote Specific Compounds; Columns Represent Particular Analyses. (A) Main Chain Deformability Reveals Flexible Loops and Hinge Regions. (B) B-Factor Plot Comparing Crystallographic Data With NMA-Predicted Mobility. (C) Eigenvalue Plot Depicts Overall Structural Rigidity: Lower Eigenvalues Indicate Greater Flexibility. (D) Variance Plot Showing Individual (Red) and Cumulative (Green) Contributions per Mode. (E) Covariance Plot Displaying Correlated (Red), Uncorrelated (White), and Anticorrelated (Blue) Motions of Residues, Calculated From Atomic Fluctuations. (F) Elastic Network Model Illustrating Connections Between Residues, With Darker Nodes Indicating Stronger Interactions
Principal Component Analysis (PCA) was utilized to provide a visual representation of the conformational motions of the receptor-ligand complexes. A consistent pattern in the results was observed across all complexes. The Mode Index vs Variance (%) plot demonstrated that the first four principal components accounted for a significant portion of the overall conformational motion, with the first mode contributing 20%, the second 18%, the third 16%, and the fourth 10%, collectively accounting for 64% of the total variance. This finding suggests that receptor-ligand complexes exhibit dominant, collective motions in specific low-frequency directions. It has been established that such motions are typically associated with functionally relevant conformational changes, including domain rearrangements, loop flexibility, and active site breathing. Modulation or stabilization of these changes has been shown to be possible through ligand interactions (Figure 5). The correlation matrix from MD simulation demonstrated analogous dynamic interactions between residues in all receptor-ligand complexes. The red regions indicate positive correlations, thus highlighting residues that have moved in a coordinated manner, particularly in the vicinity of the ligand-binding site. This finding provides support for the induced-fit mechanism, in which ligand binding induces conformational changes in the receptor, thereby stabilizing its structure and promoting cooperative motions between residues. It was evident that the blue regions in the matrix indicated negative correlations, thus suggesting that certain residues moved in opposite directions. This phenomenon may be indicative of compensatory motion or of allosteric effects. The existence of such negative correlations indicates the potential for allosteric communication, whereby the binding of a ligand at one site could trigger structural changes in distant regions of the receptor. The main diagonal of the correlation matrix confirmed that each residue maintains independent, harmonious motion, as expected in dynamic systems. The results of the receptor-ligand MD simulations suggested that ligand binding induced coordinated movements within the receptor and stabilized its conformation. The simulations also demonstrated that ligand binding caused structural shifts in distant regions, potentially via allosteric communication. This analysis demonstrated consistency across all five receptor-ligand complexes, thereby underscoring the observation of analogous dynamic behavior and interactions (Figure 5).
When we initially established the QSAR model for the five selected HER2-targeting ligands, their predictive accuracy proved highly limited, with an R2 value of merely 0.18 and an RMSE of 1.19. However, after optimization of the model by selecting the most relevant features, its performance improved significantly; the RMSE was decreased to 0.57, which indicates a closer match between predicted and actual pIC50 values. Upon further investigation of molecular descriptors, we identified certain properties that play a crucial role in HER2 inhibition. Specifically, possessing more hydrogen bond donors (correlation coefficient r = 0.63), higher lipophilicity (LogP, r = 0.60), and stronger sp3 hybridization (Fraction CSP3, r = 0.60) were all closely associated with enhanced activity. This suggested that augmenting these properties could yield more potent future compounds.
QSAR Model Performance and Correlations Between Descriptors With HER2 Inhibitory Activity (pIC50)
QSAR Analysis of HER2 Inhibitors Illustrating the Alignment Between Predicted vs Experimental pIC50 Values, Correlations Between Descriptors and Activity, and Molecular Weight-Potency Trends, Emphasizing the Key Structural Features Driving Inhibitory Activity
The experimental pIC50 values for the selected HER2 ligands exhibited a clear activity gradient, with two compounds demonstrating potent inhibitory effects. Compound 2048788 (pIC50 ≈ 8.6) and compound 3956509 (pIC50 ≈ 8.4) exhibited the highest activity, indicating their potential for further development as HER2-targeted compounds. Compound 1923019 (pIC50 ≈ 8.0) also exhibited significant activity, suggesting that it is a HER2-associated ligand with moderate therapeutic potential. In contrast, compounds 3327013 (pIC50 ≈ 6.0) and 3903447 (pIC50 ≈ 5.0) exhibited slightly weak inhibitory effects. Importantly, the dynamic range exceeded three logarithmic units (∼3.6 log units pIC50 range, 5.0-8.6), which is ideal for QSAR methods, as the differential range between weak and strong binders was sufficiently large to enhance the robustness or indicative power of large-scale structure-activity relationship analyses (Figure 6). The efficacy of the random forest model was assessed by comparing experimental values with predicted values for five HER2 ligands (Figure 6). The scatter plot exhibits a distinct near-diagonal trajectory, wherein higher experimental activity correlates with elevated predicted values, which indicates the model accounts for overall structure-activity relationships. Some minor biases may exist within the model (eg, overestimation of weak inhibitors (pIC50 ≈ 5-6) and underestimation of strong inhibitors (pIC50 ≈ 8-8.7)), but in terms of predictive range, the model was reasonably accurate, yielding no extreme outliers above or below experimental values; this indicates regression, suggesting the model is expected to correctly place compounds within the correct model in a small dataset. The random forest model demonstrated promising predictive capability and, more importantly, exhibited meaningful relative consistency with experimental results, indicating its suitability for preliminary QSAR analysis of HER2 inhibitors. Analysis of molecular weight against experimental pIC50 values revealed a distinct favorable potency window within the 450-500 Da range (Figure 6). Ligands with molecular weights of 458 Da, 480 Da, and 506 Da consistently exhibited high pIC50 values (eg, pIC50 = 8.0 - 8.6), indicating that ligands within the medium-to-high molecular weight range are potent HER2 inhibitors. While this trend does not exhibit a strictly linear relationship, the overall pattern indicates that only compounds with slightly larger molecular weights possess higher activity.
Discussion
HER2 overexpression is closely associated with an aggressive disease phenotype, and breast cancer remains a significant clinical challenge. Although HER2 inhibitors have markedly improved patient outcomes, issues such as drug resistance and off-target effects necessitate the development of novel scaffolds with superior pharmacological properties. Beyond conventional drug discovery approaches, ligand-based virtual screening represents an effective strategy for identifying novel HER2 inhibitors from compound libraries.21,22 Structural diversity in compounds may assist researchers in developing inhibitors capable of circumventing potential resistance mechanisms observed with currently marketed HER2-targeted therapeutics. To avoid ADME-related issues in later development stages, establishing suitable ADME properties during drug discovery is paramount. 23
Virtual screening of the ChEMBL database has yielded 39 compounds exhibiting favorable ADME characteristics; among these, five compounds demonstrated superior binding affinity compared to reference drugs: doxorubicin, letrozole, and lanatitinib. Further investigation into the contribution of their structural and functional groups will provide a mechanistic basis for their binding behavior. Results indicate that the binding behavior of these compounds closely correlates with that of the reference compounds, exhibiting no anomalous behavior during the docking process. Interactions indicate enhanced hydrogen bonding compared to reference compounds. Several identified compounds contain halogen substituents (fluorine and chlorine), enabling halogen bonds with residues GLU975 and GLU971. Halogen interactions are widely recognized for enhancing target affinity and specificity in kinase inhibitors. 24 The Michael receptor moiety of lanatitinib is formed via covalent bonds, whereas halogenated quinazoline derivatives have been demonstrated to enhance the reversible interactions within HER family kinases. 25 The presence of halogen bonds in our most promising compounds suggests a similarly potent anchoring mechanism, enabling robust fixation within the HER2 pocket. Furthermore, the aromatic heterocycle (quinoline/quinazoline skeleton) in our lead compound engages in π-π stacking interactions with TYR772 and ILE954. Earlier studies reported similar aromatic rings stabilizing EGFR/HER2 inhibitor orientation via π-π stacking and hydrophobic effects. 26 Although Neratinib relies on the covalent modification of CYS805 in EGFR/HER2, the results of our study are matched (with non-covalent interactions), suggesting that reversible binding may reduce off-target reactivity.
Among several matched results screened, the presence of basic amino groups engaged π-cation interactions with ARG978 and LYS765, a motif recurrently observed in kinase domain recognition. 27 Such cationic interactions may be crucial for stabilizing the ligand orientation near the ATP-binding pocket. Previous computational reports on curcumin derivatives indicated that these cationic interactions are also paramount to HER2 inhibition. 25 Among the top five compounds, distinctive hydrophobic interactions at residues VAL773, VAL839, and MET953 represent consistent findings. These also form the core of action for FDA-approved HER2 inhibitors (eg, lapatinib and neratinib), which fundamentally exploit hydrophobic anchoring within the kinase hinge region. 28 This indicates that hydrophobic substitutes play a considerably important role in the stable occupation of the active site. Molecular dynamics simulations reinforced the stability of these interactions, showing consistent flexibility in hinge regions and low-frequency vibrational modes, in agreement with previous reports. 24 Importantly, compounds I and V bind with hydrogen bonding, π-π stacking, and halogen interactions, exhibiting the strongest and most stable structural features, suggesting that multifunctional interactions may support their superior inhibitory activity. Collectively, these results indicate that structural features such as halogen substituents, aromatic skeletons, and cationic functional groups are key determinants of HER2 inhibitory potential. In comparison to existing HER2 inhibitors, these identified compounds offer alternative scaffolds with diverse interaction profiles, potentially advantageous for combating resistance. Future optimization of these scaffolds alongside in vitro and in vivo validation may facilitate the emergence of next-generation HER2-targeted therapies. Though these computational findings represent valuable preliminary insights, certain limitations of the study should be acknowledged. This work describes a computational drug discovery approach, rendering its findings predictive rather than conclusive. Computational models cannot capture the full complexity of biological systems, and their outcomes depend on the accuracy of selected protein structures, docking algorithms, and parameter settings, potentially introducing bias. Consequently, experimental validation using in vitro enzyme inhibition assays, cell-based approaches, and in vivo pharmacological assessments is essential to confirm the therapeutic relevance of the identified compounds.
Conclusion
Novel and cost-effective therapeutic agents are being identified using computational methods for drug discovery. In-silico approaches reduce time and resources by predicting drug candidates and evaluating ADME properties before experiments. Targeting overexpressed HER-2 receptors, which are critical drivers of breast cancer progression, remains a promising therapeutic approach. Five potential HER-2 inhibitors were identified in this study from a ChEMBL compound library (IDs: 1923019, 2048788, 3956509, 3903447, 3327013) with binding affinities ranging from −9.2 to −10.6 kcal/mol, outperforming the reference drugs Doxorubicin (−7.8), Letrozole (−7.5%), and Neratinib (−9.0%). The candidates also demonstrated excellent absorption and permeability (Lipinski compliance, TPSA 120, LogP 2.5-4.0). Hydrophobic interactions, hydrogen bonding, and stacking with major HER-2 residues were observed in docking, while MD simulations demonstrated structural stability (RMSD <2.5 Å, consistent RMSF profiles). QSAR analysis further validated key structural descriptors influencing their activity, reinforcing the reliability of these hits. Collectively, these compounds demonstrate superior binding efficiency, stability, pharmacokinetics, and QSAR-predicted activity, highlighting their potential as strong candidates for early-stage HER-2–targeted drug development.
Footnotes
Acknowledgements
The authors gratefully acknowledge the support of the Faculty of Pharmacy, Gomal University, for providing computational facilities during this research. We also extend our thanks to colleagues and collaborators who contributed their insights and assistance throughout the study. The authors did not receive any specific grants from funding agencies in the public, commercial, or not-for-profit sectors.
Funding
The authors received no financial support for the research, authorship, and/or publication of this article.
Declaration of Conflicting Interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Ethics Statement
This research was conducted entirely through computational and in-silico approaches. No human participants, clinical data, or animal models were involved in the study; therefore, ethical approval and informed consent were not required.