Background

GSK3B (Glycogen Synthase Kinase 3 Beta) is a crucial serine-threonine kinase enzyme, encoded by the GSK3B gene, that acts as a central regulator in numerous cellular processes, including metabolism (like glucose), cell development, Wnt signaling (affecting beta-catenin), and apoptosis (programmed cell death). It's known for its role in negatively regulating pathways and its involvement in diseases like Alzheimer's, Parkinson's, bipolar disorder, and certain cancers, making it a significant target for therapeutic research. Explore CD ComputaBio's GSK3B inhibitors discovery services, offering innovative solutions for drug development and precision medicine. Our expertise accelerates the identification and optimization of GSK3B inhibitors to support your research in diverse areas.

Computational Methods

Molecular Docking

The crystal structure of the GSK3B protein was obtained from the AlphaFold database, and the 3D structure of the small molecule Artemisinin was downloaded from the PubChem database, followed by energy minimization under the MMFF94 force field. Molecular docking was performed using AutoDock Vina 1.2.3. Prior to docking, the receptor protein was processed in PyMol 2.5.2 to remove water molecules, salt ions, and any co-crystallized small molecules. A docking box was then set to enclose the entire protein. All processed small molecules and receptor structures were converted to the PDBQT format required by AutoDock Vina using ADFRsuite 1.0. During docking, the global search exhaustiveness was set to 32 while other parameters remained at their defaults. The highest-scoring pose was taken as the binding conformation and visualized in PyMol 2.5.2.

Molecular Dynamics (MD) Simulation

Full-atom MD simulations were performed in AMBER 24 starting from the docked protein–ligand complex. Ligand charges were obtained via the antechamber module and Gaussian 09 using Hartree–Fock (HF) SCF/6-31G* calculations. The ligand and protein were described by the GAFF2 small-molecule force field and the ff14SB protein force field, respectively. All systems were prepared using LEaP by adding hydrogen atoms, enclosing the complex in a truncated octahedral TIP3P water box with a 10 Å buffer, and adding Na+/Cl− ions to neutralize the system, after which topology and parameter files were generated.

Energy minimization comprised 2,500 steps of steepest descent followed by 2,500 steps of conjugate gradient. Subsequently, the system was heated at constant volume for 200 ps from 0 K to 298.15 K. An NVT ensemble run of 500 ps at 298.15 K allowed solvent molecules to distribute uniformly, followed by a 500 ps NPT equilibration. A 100 ns NPT production run was then carried out under periodic boundary conditions. A 10 Å nonbonded cutoff was used; long-range electrostatics were treated by the Particle Mesh Ewald (PME) method; the SHAKE algorithm constrained bonds involving hydrogens; and temperature was controlled with a Langevin thermostat with collision frequency γ = 2 ps⁻¹. The pressure was set to 1 atm, the integration time step to 2 fs, and snapshots were saved every 10 ps for analysis.

MM/GBSA Binding Free Energy Calculation

Binding free energies between protein and ligand were computed using the MM/GBSA method based on 90–100 ns of MD trajectories, according to: ΔG_bind = ΔG_complex − (ΔG_receptor + ΔG_ligand) = ΔE_internal + ΔE_vdW + ΔE_elec + ΔG_GB + ΔG_SA. Here, ΔE_internal includes bond (E_bond), angle (E_angle), and torsion (E_torsion) terms; ΔG_GB and ΔG_SA are the solvation free energy contributions, with ΔG_GB the polar term (GB model with igb = 2) and ΔG_SA the nonpolar term, estimated as ΔG_SA = 0.0072 × ΔSASA. Configurational entropy was neglected due to its high computational cost and limited accuracy.

Computational Results Analysis

Docking Results Analysis

Docking is an efficient approach to explore small-molecule interactions with a target protein. Using Vina 1.2.3, Artemisinin was docked into GSK3B. The binding mode shows Artemisinin forming a hydrogen bond with LYS-104 of GSK3B, strengthening the complex. The ligand also engages in π–π interactions with PHE-109 and hydrophobic contacts with LEU-5, LEU-4, LYS-104, PHE-109, and LEU-8, providing substantial van der Waals stabilization.

Figure 1. Binding mode of GSK3B with Artemisinin obtained from docking (left: overall view; right: close-up).

A negative docking score indicates potential binding—the lower (more negative), the stronger the predicted affinity. For this complex, the docking score for Artemisinin bound to GSK3B was −7.776 kcal/mol, suggesting favorable binding and motivating subsequent MD analysis.

MD Results Analysis

Panels (A–F) track (A) Ligand RMSD, (B) Complex RMSD, (C) Protein RMSF, (D) Complex radius of gyration (Rg), (E) number of hydrogen bonds, and (F) complex solvent-accessible surface area (SASA) over simulation time.

Figure 2. Time series of MD observables (A–F).

Ligand RMSD (A): Throughout the 50 ns trajectory, the ligand RMSD remained low (~0.01–0.05 nm) with minor fluctuations and no upward trend, indicating conformational stability within the binding pocket without large-scale displacement—consistent with stable receptor–ligand interactions.

Complex RMSD (B): The complex RMSD rose rapidly during 0–10 ns and then plateaued around 0.20–0.22 nm, suggesting early equilibration. The modest fluctuations imply a structurally stable complex without substantial conformational rearrangements.

Protein RMSF (C): Residue-wise RMSF values were generally low, with most residues fluctuating under 0.1 nm. Low flexibility in the binding-site region further supports complex stability.

Complex Rg (D): The radius of gyration stayed near 2.13–2.16 nm with a narrow fluctuation range and no trend, indicating a compact complex without expansion or collapse—evidence of consistent global conformation.

Number of H-bonds (E): The protein–ligand hydrogen bond count fluctuated between 0 and 3, frequently showing 1–2 H-bonds. This points to persistent yet dynamic hydrogen-bonding that helps stabilize binding while reflecting some flexibility at the site.

Complex SASA (F): SASA varied within ~170–190 nm², initially somewhat higher and then stabilizing with small fluctuations. This suggests minor early adjustments followed by convergence and stable surface exposure.

Overall, the protein–ligand complex exhibited robust structural and binding stability during MD.

MM-GBSA Results Analysis

Table 1. Binding free energies and energy components predicted by MM/GBSA (kcal/mol).

ΔE_vdW: van der Waals energy; ΔE_elec: electrostatic energy; ΔG_GB: electrostatic contribution to solvation; ΔG_SA: nonpolar contribution to solvation; ΔG_bind: binding free energy.

Based on MD trajectories, MM-GBSA calculations yielded a binding free energy of −23.88 ± 0.63 kcal/mol for GSK3B_Artemisinin. The negative value supports appreciable affinity; more negative indicates stronger binding. Energy decomposition points to van der Waals interactions as the dominant favorable term, followed by the nonpolar solvation contribution.

Figure 3. Top-10 residue contributions to total binding free energy (more negative indicates greater stabilization).

The plot highlights the top 10 residues contributing most to binding (all negative contributions). ALA123 (≈ −2.0 kcal/mol) and LEU387 (≈ −1.8 kcal/mol) are primary hotspots, providing the strongest driving forces. MET123 and PHE456 also furnish notable hydrophobic/van der Waals support. Hydrophobic residues (LEU, MET, ALA) dominate, suggesting a hydrophobic pocket; the aromatic ring of PHE753 may further engage in π-stacking or hydrophobic insertion. These residues collectively lower the binding free energy and represent important targets for structure-based optimization.

Free Energy Landscape (FEL) Analysis

Figure 4. FEL of the GSK3B–Artemisinin complex as a function of RMSD (nm) and Rg (nm).

The FEL map uses color to denote free-energy levels (red: higher; blue: lower). The complex exhibits a low-energy basin near RMSD ≈ 2.3 nm and Rg ≈ 2.135 nm, indicating particularly stable conformational states.

Related Services

CD ComputaBio provides you with comprehensive services, including but not limited to:

• Protein-Small Molecule Docking Service
• Protein Molecular Dynamics Simulation
• Binding Free Energy Analysis Service

Connect with Us Anytime!

Studies on the interactions between proteins and inhibitors are crucial for understanding the mechanisms of biological regulation, and they offer a theoretical foundation for the design and discovery of new drug targets. If you require services for inhibitor discovery, please don't hesitate to get in touch with us.

References:

1. Eberhardt, J.; Santos-Martins, D.; Tillack, A. F.; Forli, S., AutoDock Vina 1.2.0: New docking methods, expanded force field, and python bindings. J. Chem. Inf. Model. 2021, 61(8), 3891–3898.
2. DeLano, W. L., PyMol: An open-source molecular graphics tool. CCP4 Newsletter on Protein Crystallography 2002, 40(1), 82–92.
3. Ravindranath, P. A.; Forli, S.; Goodsell, D. S.; Olson, A. J.; Sanner, M. F., AutoDockFR: Advances in Protein–Ligand Docking with Explicitly Specified Binding Site Flexibility. PLoS Comput. Biol. 2015, 11(12), e1004586.
4. Salomon-Ferrer, R.; Case, D. A.; Walker, R. C., An overview of the AMBER biomolecular simulation package. WIREs Comput. Mol. Sci. 2013, 3(2), 198–210.
5. Frisch, M.; Trucks, G.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G., Gaussian 09, Revision D.01, Gaussian, Inc., Wallingford CT, 2009.
6. Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A., Antechamber: an accessory software package for molecular mechanical calculations. J. Am. Chem. Soc. 2001, 222, U403.
7. Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A., Development and testing of a general AMBER force field. J. Comput. Chem. 2004, 25(9), 1157–1174.
8. Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C., ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 2015, 11(8), 3696–3713.
9. Mark, P.; Nilsson, L., Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K. J. Phys. Chem. A 2001, 105(43), 9954–9960.
10. Sagui, C.; Darden, T. A., Molecular dynamics simulations of biomolecules: long-range electrostatic effects. Annu. Rev. Biophys. Biomol. Struct. 1999, 28(1), 155–179.
11. Kräutler, V.; Van Gunsteren, W. F.; Hünenberger, P. H., A fast SHAKE algorithm to solve distance constraint equations for small molecules in MD simulations. J. Comput. Chem. 2001, 22(5), 501–508.
12. Larini, L.; Mannella, R.; Leporini, D., Langevin stabilization of molecular-dynamics simulations of polymers by means of quasi-symplectic algorithms. J. Chem. Phys. 2007, 126(10), 104101.
13. Hou, T.; Wang, J.; Li, Y.; Wang, W., Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. Accuracy of binding free energy calculations based on MD simulations. J. Chem. Inf. Model. 2011, 51(1), 69–82.
14. Genheden, S.; Ryde, U., The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin. Drug Discov. 2015, 10(5), 449–461.
15. Rastelli, G.; Di Rienzo, A.; Degliesposti, G.; Sgobba, M., Fast and accurate predictions of binding free energies using MM-PBSA and MM-GBSA. J. Comput. Chem. 2010, 31(4), 797–810.
16. Nguyen, H.; Roe, D. R.; Simmerling, C., Improved Generalized Born solvent model parameters for protein simulations. J. Chem. Theory Comput. 2013, 9(4), 2020–2034.
17. Weiser, J.; Shenkin, P. S.; Still, W. C., Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO). J. Comput. Chem. 1999, 20(2), 217–230.

• Tubulin Polymerization Inhibitor Design
• Protein-ligand Interactions Case Study