Structure alignment and pharmacophoric models
In the flexible structure alignment of the seven selected GSK-3β inhibitors (I–VII), the optimal scoring pose yielded U = 41.8961, F = -108.0886, and S = -66.1924, as illustrated in Fig. 2. Our analysis indicates that these metrics encapsulate complementary aspects of the superposition. U represents the alignment and uniformity across the ligand set, with higher values indicating a more consistent standard frame of reference37. F is a field score that measures the electrostatic and van der Waals field overlap between the superposed molecules, where more negative values denote better field agreement38. S quantifies steric complementarity or shape overlap, with more negative values indicating an improved steric fit39. The combination of a high U with strongly negative F and S values suggests a well-conserved arrangement of key pharmacophoric elements among the seven inhibitors (I-VII), establishing a solid foundation for ligand-based pharmacophore extraction.
Fig. 2
Structure alignment of seven representative GSK-3β inhibitors (I-VII) and derived pharmacophore model.
A pharmacophore based on ligands was constructed from the alignment by identifying recurring aromatic, hydrophobic, and hydrogen-bonding features in the superposed compounds. Several feature types were co-localized: an aromatic and hydrophobic pair at (0.62, 1.85, -0.75), another aromatic/hydrophobic pair at (-3.18, -1.04, 0.98), and a third aromatic/hydrophobic pair at (-0.66, -1.04, -0.05). Two distinct pairs of hydrogen-bond acceptor/donor features were identified at (2.61, 3.60, -0.91) and (-1.67, 2.48, -0.79), illustrating their 3D spatial arrangement as shown in Fig. 3A. The .ph4 pharmacophore query was applied to the ZINCPharmer collection, yielding 1,085 hits that met our fitting criteria.
Fig. 3
Visual representation and 3D spatial arrangement of the generated pharmacophore features used for virtual screening. (A) The ligand-based consensus model and (B) the structure-based model derived from PDB 3GB2. The feature colors indicate the corresponding pharmacophoric elements in the displayed model [green = hydrophobic, orange = hydrogen bond acceptor, gray or white = overlapping aromatic or excluded volume features].
A structure-based pharmacophore was generated from the GSK-3β crystal structure (PDB ID: 3GB2) using the Pharmit webserver, deriving features from the co-crystallized ligand. The model consisted of three hydrogen-bond acceptor features positioned at approximately (36.17, 14.51, -0.11), (35.74, 15.83, 0.02), and (25.26, 15.31, 4.45), along with five hydrophobic points clustered around (28.76, 16.72, 3.67), (33.28, 18.26, 1.46), (31.55, 19.32, 2.23), (35.27, 14.92, 0.54), and (35.23, 12.32, 0.95), as shown in Fig. 3B. The pharmacophore was exported as a .json query and screened in ZINCPharmer, resulting in 36 candidate hits. Only compounds aligned with the pharmacophore with an RMSD ≤ 1.0 Å were selected for further docking and ADME evaluation in ligand- and structure-based searches.
Initial molecular docking filtering
Molecular docking is essential in structure-based drug discovery, as it predicts the orientation and interaction of small molecules within a protein’s binding site40,41. In our hierarchical virtual screening pipeline, hits were sequentially filtered based on geometric pharmacophore mapping (RMSD ≤ 1.0 Å), predicted binding affinity, and strict ADME/developability thresholds. It offers a rapid and cost-effective approach to prioritizing extensive compound libraries for experimental evaluation. Docking algorithms sample ligand conformations and poses, evaluating them with scoring functions that approximate binding affinity. This process enables researchers to rank hits from virtual screens, identify key protein-ligand contacts such as hydrogen bonds, hydrophobic pockets, and π-π interactions, and generate hypotheses for structure–activity relationships and lead optimization42,43,44. Molecular docking was conducted for all compounds obtained from the ligand- and structure-based pharmacophore screens, using the energy of the co-crystallized ligand (-6.8 kcal/mol) as a benchmark for filtering. Of the 1,085 hits identified by the ligand-based pharmacophore, 957 compounds (approximately 88.2%) demonstrated predicted binding affinities exceeding − 6.8 kcal/mol and were consequently retained. Of the 36 candidates identified through the structure-based pharmacophore, 22 compounds (approximately 61.1%) exhibited superior scores to the co-crystallized ligand and were subsequently advanced.
ADME filtering
In-silico ADME prediction was conducted to efficiently evaluate the drug-likeness and developability of virtual hits before experimental testing. Computational ADME tools assess essential properties that affect absorption, distribution, metabolism, and excretion45. We prioritized BBB permeability and Lipinski compliance due to compounds needing to achieve therapeutically relevant concentrations in the central nervous system for our target indication, Alzheimer’s disease46. As an explicit exclusion criterion, compounds presenting more than one Lipinski violation were strictly rejected to minimize developability concerns. SwissADME was utilized to predict BBB permeation and associated properties, including P-glycoprotein efflux propensity47 and polar surface area, allowing for the deprioritization of potent docking hits that do not possess the necessary physicochemical profile for brain exposure48. For the ligand-based pharmacophore hits, we established an additional criterion of at least “moderate” aqueous solubility49. This requirement is essential as sufficient solubility facilitates systemic exposure, allows for reliable in-vitro testing and formulation, and ensures that a compound’s lipophilicity-driven blood-brain barrier potential is not compromised by inadequate dissolution. Collectively, these ADME filters aid in balancing predicted potency with a feasible opportunity for advancing a CNS-active candidate.
All ligand-based hits (LB1-LB3) satisfy the “moderately soluble” ESOL classification, serving as an additional filter for the ligand-based pharmacophore series (Table 1). In contrast, the structure-based hit (SB1) is predicted to be poorly soluble by Silicos-IT, despite being classified as “moderately soluble” by ESOL, indicating a potential concern for the developability of (SB1). All compounds show high predicted gastrointestinal absorption and are predicted to be BBB-permeant, which is desirable for a CNS target such as GSK-3β in Alzheimer’s disease. Two compounds are predicted to be P-glycoprotein substrates (SB1 and LB2), which could limit brain exposure despite predicted BBB permeation and thus should be considered when prioritizing compounds (Fig. 4). All compounds adhere to Lipinski’s Rule-of-Five with no violations and exhibit identical predicted oral bioavailability scores of 0.55. Synthetic accessibility scores vary from approximately 2.2, indicating ease of synthesis, to 3.7, representing greater difficulty, with PAINS alerts recorded as zero throughout the series.
Table 1 Predicted physicochemical and ADME properties of prioritized hits from the structure- and ligand-based pharmacophore screens, as computed by SwissADME.Fig. 4
SwissADME BOILED-Egg model predictions for prioritized GSK-3β inhibitor hits.
In-silico toxicity profiling
Following the ADME evaluation, the safety profiles of the top candidates, SB1 and LB1, were assessed using the ProTox-3.0 webserver. The predicted acute oral toxicity and specific endpoint liabilities provided critical insights into the developability of these compounds. SB1 demonstrated a highly favorable acute toxicity profile, with a predicted LD50 of 2500 mg/kg, placing it in Toxicity Class 5 (indicating it is generally considered to have low acute toxicity). However, specific endpoint predictions highlighted potential risks for respiratory toxicity (probability: 0.82) alongside mild alerts for neurotoxicity and mutagenicity, as illustrated in Fig. 5.
In contrast, LB1 exhibited a predicted LD50 of 521 mg/kg, classifying it into Toxicity Class 4 (harmful if swallowed). This indicates a moderate acute toxicity profile compared to SB1. The endpoint analysis for LB1 showed potential liabilities for respiratory toxicity (probability: 0.69) and nephrotoxicity (probability: 0.59).
Fig. 5
In-silico toxicity profiling and endpoint prediction of the top prioritized GSK-3β candidates. The figure illustrates the predicted acute oral toxicity (LD50), designated toxicity class, prediction accuracy, and radar plots detailing specific organ and endpoint liabilities for the structure-based hit SB1 (left) and the ligand-based hit LB1 (right).
Notably, both compounds yielded “Active” predictions for blood-brain barrier (BBB) permeability within the ProTox-3.0 screening. Rather than a liability, this corroborates our SwissADME findings and confirms a crucial pharmacokinetic prerequisite for Alzheimer’s disease therapeutics. While these in-silico toxicity alerts serve as a valuable early-warning system to guide future structural optimization, they are strictly predictive and necessitate rigorous in vitro and in vivo toxicological validation during subsequent lead optimization phases.
Molecular docking of final hits for interaction analysis
The docking protocol was validated by re-docking the co-crystallized ligand into the binding pocket of GSK-3β (PDB ID: 3GB2). The predicted binding pose achieved a root-mean-square deviation (RMSD) of 1.45 Å compared with the experimental crystallographic conformation, indicating good reliability of the docking setup and justifying its use for subsequent analyses (Fig. 6).
Fig. 6
Validation of the docking protocol.
The docking scores of the six prioritized hits were compared with that of the co-crystallized ligand (-6.80 kcal/mol) to evaluate their predicted binding poses and relative scoring metrics toward GSK-3β. All candidate molecules demonstrated more favorable predicted binding affinities than the reference ligand. However, it is important to note that these docking scores are semi-quantitative estimates used for compound prioritization; more negative scores do not inherently guarantee proportionally higher in-vitro inhibitory potency. The structure-based hit SB1 yielded a binding affinity of -7.16 kcal/mol, while the ligand-based hits LB1, LB2, and LB3 showed values of -7.28, -7.00, and − 7.39 kcal/mol, respectively (Table 2). These results suggest that all four hits bind more tightly than the crystallographic reference, with LB1 and LB3 emerging as the most promising candidates for further structural and dynamic investigation.
Table 2 Binding affinities of the structure-based hit (SB1), ligand-based hits (LB1-LB3), and the co-crystallized ligand within the active site of GSK-3β.
The co-crystallized ligand of GSK-3β formed a single hydrogen bond between its oxadiazole scaffold and Lys85, along with multiple hydrophobic contacts involving Ile62, Phe67, Val70, Ala83, Lys85, Val110, Leu188, and Cys199 (Fig. 7). In contrast, the newly identified hits established richer interaction profiles. The structure-based hit SB1 formed two hydrogen bonds, one between its sulfur atom and Lys183 and another between its oxygen atom and Ser203, in addition to π-anion interactions between its naphthyl group and Asp181/Asp200. Similarly, LB1 reproduced the Lys183 and Ser203 hydrogen bonds via its dihydroxy-substituted phenyl ring, accompanied by π-anion interactions with Asp181 and Asp200. LB2 also maintained the Lys183 hydrogen bond and introduced an additional interaction with Asn186; its extended conjugated system enabled both π-anion and π-cation interactions. LB3 showed an interaction profile comparable to LB2 but lacked the Asn186 hydrogen bond.
These additional hydrogen bonds and π-based interactions are predicted to enhance affinity and specificity toward GSK-3β. Hydrogen bonding with residues such as Lys183, Asn186, and Ser203 likely improves binding stability by anchoring the ligands within the active site. In contrast, π-anion and π-cation interactions with Asp181, Asp200, and positively charged residues further reinforce electrostatic complementarity. These enriched interaction profiles suggest that the identified hits establish robust structural contacts with GSK-3β, yielding highly favorable computational scoring metrics. While these predictive profiles are encouraging, subsequent in-vitro enzymatic assays are strictly required to translate these in-silico docking scores into confirmed inhibitory potency.
Fig. 7
Two-dimensional interaction diagrams of the co-crystallized ligand (A), SB1 (B), LB1 (C), LB2 (D), and LB3 (E) docked within the active site of GSK-3β (PDB ID: 3GB2).
MD stability and interactions
Molecular dynamics (MD) gives biomolecular systems a time-resolved, physics-based picture that static crystal structures cannot. MD shows whether a protein retains its folded architecture or undergoes functionally important rearrangements, whether ligands hold their binding posture or diffuse and reorient, and on what time frames these events occur by propagating atom positions under realistic force-field stresses50. Before doing more expensive studies, this temporal information connects structural snapshots to thermodynamic and kinetic behavior to inform theories about binding affinity, selectivity, and hit stability51,52.
In MD analysis, root-mean-square deviation (RMSD) is basic but powerful. After aligning the protein, ligand RMSD measures ligand binding strength, while backbone RMSD records global protein structural drift and checks system equilibration immediately53. Low and tightly distributed backbone RMSD values indicate the protein retains its crystallographic fold and the complex does not experience significant, artifactual rearrangements; low ligand RMSD values indicate pose retention and fewer major ligand translations or rotations inside the binding pocket54. Two RMSD measurements assess simulation fidelity and binding stability across ligands or pharmacophore-derived models55,56.
The 100-ns simulations indicate notable variations between the co-crystallized ligand (Co), the structure-based pharmacophore model SB1, and the ligand-based model LB1 (Fig. 8). SB1 has the lowest protein perturbation (backbone RMSD ≈ 0.180 Å), compared to the co-crystallized ligand (Co: ≈ 0.202 Å) and slightly lower than LB1 (≈ 0.197 Å) in the full 0-100 ns traces. This suggests that the protein accommodates SB1 with less global drift. The ligand RMSD of LB1 is ≈ 0.815 Å, considerably lower than the co-crystallized ligand (Co: ≈ 1.000 Å) and slightly lower than SB1 (∼ 0.997 Å) along the trajectory, indicating a tighter binding site retention.
Fig. 8
RMSD profiles for five 100-ns MD runs on GSK-3β, showing backbone stability (left panel) and ligand positional fluctuations (right panel).
In the final 20 ns of the production-like interval (80–100 ns), SB1 has a lower and more tightly distributed backbone RMSD than the co-crystallized ligand (Co: ≈ 0.216 Å) and LB1 (≈ 0.231 Å), indicating less global structural perturbation in the equilibrated complex. During the production phase, LB1 remains the most conformationally constrained ligand in the binding site, with a ligand RMSD of 0.721 Å, compared to the co-crystallized ligand at 1.081 Å and SB1 at 1.115 Å. More insight comes from predicted stabilization timeframes. The rolling-mean stability test shows that SB1’s backbone stabilizes early (28.7 ns) and remains narrowly distributed throughout the simulation, while the co-crystallized ligand and LB1 stabilize later (Co ~ 74.7 ns and LB1 ~ 76.6 ns). For ligand RMSD, the co-crystallized ligand stabilizes at 73.7 ns, LB1 at 78.6 ns, and SB1 at ≈ 86.0 ns, indicating that SB1 has a rapidly equilibrated protein backbone. Still, its ligand requires longer rearrangements to stabilize.
The root-mean-square fluctuation (RMSF) is a key indicator for ligand-induced stabilization because it gives residue- or frame-resolved information about a protein’s local flexibility during MD simulations57. RMSF complements RMSD by showing which protein regions (or the whole-protein average, when reported as a single curve) become more rigid or flexible with different ligands. Reduced RMSF in ligand-bound simulations often indicates stronger, more persistent contacts and/or a reduction in local conformational sampling, which can favor high-affinity binding. RMSF comparisons in GSK-3β reveal how new pharmacophore-derived molecules impact backbone and side-chain mobility compared to the crystallographic reference58,59.
The average RMSF values for all simulated frames were: Co = 0.1051 Å, SB1 = 0.0927 Å, and LB1 = 0.1022 Å (Fig. 9). Compared to the co-crystallized ligand, SB1 decreases GSK-3β average fluctuation by ≈ 11.8%, while LB1 only reduces the mean by ≈ 2.8%. Comparing solely equilibrated segments (frames from stable plateau start to end), the mean RMSF values are Co = 0.0948 Å, SB1 = 0.0867 Å, and LB1 = 0.1003 Å. SB1 has a lower RMSF than Co (≈ 8.6% lower) in the equilibrated regime, but LB1 has a slightly higher mean (≈ 5.8% higher). Maximum instantaneous fluctuations (peaks) are Comax = 0.4371, SB1max = 0.4020, and LB1max = 0.3491. Co showed 22 large spikes (RMSF > 0.20 Å), SB1 14, and LB1 25, showing that SB1 exhibits lower average fluctuations and fewer large intermittent excursions than the other systems.
Fig. 9
RMSF values for the five trajectories for the 100-ns production run.
Many kinase inhibitors use hydrogen bonds to anchor ligands in the ATP-binding pocket and determine specificity. For GSK-3β, a strong hydrogen bond network helps stabilize the inhibitor’s conformation, reducing conformational entropy and increasing apparent affinity. Persistent hydrogen bonds can also enable crucial polar interactions that withstand solvent competition and keep the ligand catalytically relevant over inhibitory periods. Counting and characterizing protein-ligand hydrogen bonds through MD trajectories is not just descriptive, but also links microscopic contact patterns to thermodynamic and structural signatures (fold stability, binding-pose retention, and selectivity) that determine GSK-3β inhibitor efficacy60,61,62.
H-bond time series indicate model-dependent changes that explain RMSD behavior (Fig. 10). The co-crystallized ligand has sparse hydrogen bonding (usually 0–2 H-bonds) while SB1 generates a more persistent collection of polar contacts (3–6 H-bonds) during the run, indicating prolonged interactions. In contrast, LB1 forms fewer H-bonds (usually 0–3). Still, these contacts appear more consistent during the production window, which matches its superior ligand RMSD (i.e., tighter pose retention) even when its backbone stabilization occurs later. This suggests LB1 achieves pose stability through fewer well-maintained interactions than many transient ones. The H-bond counts fluctuate across all traces, with transient spikes and short-lived losses. SB1’s abundant H-bonding correlates with its early, low backbone RMSD, while LB1’s steadier (though fewer) H-bonds correlate with its lower ligand RMSD in the production period.
Fig. 10
Number of hydrogen bonds generated by the two pharmacophore models and the co-crystallized ligand over the 100-ns production run.
MD trajectories report radius of gyration (Rg) and solvent-accessible surface area (SASA) to quantify macromolecular compactness and solvent exposure. Small, low-variance Rg indicates that the system remains compact and that no large unfolding or domain-separation events occur during the simulation. In contrast, sustained increases or large fluctuations in Rg indicate the structure’s transient expansion or breathing motions. SASA measures the protein surface accessible to solvent and is sensitive to side-chain rearrangements, loop motions, and ligand repositioning. Increases in SASA indicate pocket opening or ligand egress, while decreases indicate burial or more intimate residue-ligand packing63,64.
Throughout all systems (Co, SB1, and LB1-LB3), the protein maintains a compact structure for most of the 100-ns trajectories (Fig. 11). The Rg values begin at approximately 2.08–2.09 Å and predominantly remain within a narrow range of 2.10–2.14 Å. The raw Rg traces exhibit only minor, rapid fluctuations on the order of 0.01–0.04 Å, suggesting the absence of global unfolding. Despite the overall stability, spikes are observed at several intervals, particularly around the ~ 28 ns region and in the ~ 54–58 ns window. These occurrences are consistent with short-lived breathing motions or local rearrangements, rather than a persistent loss of compactness. The SASA traces provide additional insights: initial SASA values exhibit minor variations across systems (e.g., Co ≈ 166.8 nm2 and SB1 ≈ 170.1 nm2), while most trajectories fluctuate around ~ 170–179 nm2, with occasional peaks reaching ~ 181–184 nm2. This suggests intermittent increases in solvent exposure, such as pocket opening or side-chain reorientation, which temporally align with certain Rg excursions.
Fig. 11
Rg (left panel) and SASA (right panel) values for the five trajectories.
Binding free energy calculations (MM/PBSA and MM/GBSA)
To provide a more rigorous thermodynamic validation of the binding affinities beyond the static docking scores, the Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) and Generalized Born Surface Area (MM/GBSA) methods were employed. Binding free energies (ΔGbind) and their corresponding energy components were calculated over the stabilized production trajectories for the co-crystallized ligand (Co), the ligand-based hit (LB1), and the structure-based hit (SB1). The energy decomposition analysis is summarized in Table 3. The calculations conclusively support the earlier docking and structural stability predictions. Both prioritized hits demonstrated notably stronger binding free energies compared to the co-crystallized reference. Under the MM/GBSA model, the co-crystallized ligand exhibited a ΔGbind of -8.75 ± 6.93 kcal/mol. In contrast, LB1 and SB1 displayed highly favorable binding energies of -25.74 ± 3.41 kcal/mol and − 27.68 ± 7.20 kcal/mol, respectively. A parallel trend was observed using the MM/PBSA method, confirming the robustness of these thermodynamic profiles.
Furthermore, the energy decomposition provides critical insights into the distinct binding mechanisms of the two candidates. The binding of LB1 is predominantly driven by highly favorable van der Waals interactions (ΔEvdW = -36.79 ± 3.61 kcal/mol) and non-polar solvation energy (-4.74 ± 0.44 kcal/mol). Conversely, the exceptionally strong affinity of SB1 is overwhelmingly dictated by massive electrostatic contributions (ΔEele = -214.36 ± 33.80 kcal/mol). This profound electrostatic stabilization of SB1 perfectly corroborates our previous MD observations, specifically the formation of a robust and persistent network of 3–6 hydrogen bonds (Fig. 9) that tightly anchors the ligand within the ATP-binding cleft despite heavy polar solvation penalties. Together, these thermodynamic calculations validate the structural integrity observed during the MD simulations and strongly support the advancement of both LB1 and SB1 as potent GSK-3β inhibitors.
Table 3 Energy decomposition of MM/GBSA and MM/PBSA binding free energies. Values are in kcal/mol and reported as mean ± SD.Validation of ligand-based and structure-based pharmacophore models
The validation of ligand-based and structure-based pharmacophore models was conducted using active-decoy datasets to assess their efficacy in distinguishing known GSK-3β active compounds from LUDe-generated decoys. ROC curve analysis utilized negative RMSD values as the ranking score, as lower pharmacophore RMSD values indicate superior alignment with pharmacophoric features. In the ligand-based pharmacophore model, the validation dataset included three active compounds and 150 decoy compounds. The model effectively identified 2 of the 3 active compounds and 7 of the 150 decoys. This resulted in a sensitivity of 66.67% and a specificity of 95.33%. ROC curve analysis yielded an AUC value of 0.8233, suggesting a strong preliminary discriminatory capability. The standard error was 0.1692, accompanied by a 95% confidence interval of 0.4917-1.000 and a p-value of 0.0555, which is considered borderline. These results indicate that the ligand-based pharmacophore model effectively identified active compounds while largely excluding decoys, demonstrating strong selectivity.
In the structure-based pharmacophore model, the validation dataset included one active compound and fifty decoys. The model identified the active compound along with 11 of the 50 decoys. This resulted in a sensitivity of 100% and a specificity of 78.00%. ROC curve analysis indicated an AUC value of 0.8200, suggesting satisfactory preliminary discriminatory performance. The standard error was 0.05433, with a 95% confidence interval ranging from 0.7135 to 0.9265, and a p-value of 0.2770. While the structure-based model effectively identified the active compound, the retrieval of 11 decoys suggests reduced selectivity in comparison to the ligand-based model.
Both pharmacophore models demonstrated AUC values exceeding 0.8, indicating their efficacy in differentiating active compounds from decoys. In contrast, the ligand-based model exhibited enhanced decoy rejection, retrieving merely 4.67% of the decoys, whereas the structure-based model retrieved 22.00%. Consequently, the ligand-based pharmacophore model demonstrated greater selectivity, whereas the structure-based pharmacophore model exhibited acceptable yet comparatively lower selectivity in preliminary performance (Fig. 12). Given the restricted availability of experimentally validated active compounds for validation, especially concerning the structure-based model, the ROC results must be regarded as preliminary retrospective validation. The acceptable AUC values obtained for both models validate their application in subsequent virtual screening, molecular docking, ADME/BBB filtering, toxicity prediction, and molecular dynamics simulation workflows.
Fig. 12
ROC curve analysis of the pharmacophore models.