The diabetic-related receptor proteins identified in Section 2.1 were preprocessed using Discovery Studio 2019 win, which primarily involved repairing missing residues, optimizing protonation states, and removing crystallization water molecules to ensure the structural integrity of the receptor. Then, using the semi-flexible docking method in the CDOCKER module (Tessaro and Scapozza, 2020), the residues within a 10 Å range of the eutectic ligand were defined as the active pocket, and the binding process between the ligand and receptor protein was simulated. In the end, according to the energy score and screening results, the best binding mode was selected as the initial conformation of the subsequent molecular dynamics simulation.

Based on the results from Section 2.2, this study performed molecular dynamics simulations to investigate the binding of rhamnocitrin and folic acid with the target protein GSK3β. Molecular dynamics simulation can simulate the movement and interaction of molecules at the atomic level, and analyze the dynamic changes of proteins, ligands and the environment, so as to provide information on the behavior of molecular conformational change, binding stability and dynamics between protein and ligand. In this study, Schrodinger 2024 Linux program was used to construct a molecular dynamics simulation system. In the simulation system, the receptor-ligand complex was placed in a TIP3P water solvent model, with the box boundary set at least 12 Å away from the complex to ensure adequate solvation of water molecules and avoid boundary effects. The ion concentration in the system was set to 0.154 M, with Na⁺ and Cl⁻ used to neutralize the charge of the system in the simulation process. To enhance simulation accuracy, the OPLS3e force field (Tirado-Rives and Jorgensen, 2008) was employed in the system, which effectively describes non-bonded interactions and binding modes between molecules, particularly suited for the study of protein-ligand complexes (Wu et al., 2022; Collier et al., 2020).

In the initial energy optimization, the steepest descent method was adopted for 5,000 steps to minimize unreasonable contact and high-energy conformations, with a convergence threshold set at 10 kJ/mol/nm to ensure effective adjustment of the force field parameters. Then the conjugate gradient method was applied for 2000 steps to optimize the thermodynamic state of the system in an effort to guarantee its stability (Levitt and Lifson, 1969; Powell, 1977; Fletcher, 2013).

When entering the equilibrium phase, the NVT ensemble simulation of 100 ps with a time step of 2 fs was first performed, and the system was gradually heated to 300 K so as to eliminate the influence of the initial structure and make sure thermodynamic equilibrium. This was followed by switching to the NPT ensemble for further equilibration at a constant pressure of 1 bar, with a simulation time of 100 ps to stabilize the density and pressure. The primary purpose of this stage was to allow the system to reach a stable thermodynamic state, preparing it for subsequent production simulations (Zhao et al., 2007; Wen et al., 2003).

The formal molecular dynamics simulations ran for 100 ns, during which the temperature was maintained at 300 K, the pressure at 1 bar, with the time step of 2 fs and trajectory sampling interval of 10 ps. In order to ensure the stability and accuracy of the simulation, physical properties such as temperature, pressure and volume were regularly monitored during simulation process to ensure they remained within the expected range.

On the basis of molecular dynamics simulation, in order to further analyze the interactions between ligands and key residues of diabetes-related target proteins, the Discovery Studio 2021 program was utilized in this study to carry out alanine flexible scanning. Alanine scanning replaces critical residues with alanine to identify residues that play a crucial role in ligand binding, and to assess their impact on the binding interface and affinity. The scanning targets are all amino acid residues within a 3 Å radius of the ligand-binding interface.

In this study, density functional theory (DFT) (Zhao et al., 2004; Qin et al., 2018) was used to evaluate the chemical activity of Astragalus-Codonopsis drugs on core components. Firstly, the geometric structures were optimized at the B3LYP/6-311G D3 (d, p) level (Harder et al., 2016)by the using of Gaussian 09w Linux program. After optimization, single-point energy calculations were performed under the same basis set, and energy distribution data for the Highest Occupied Molecular Orbital (HOMO) and the Lowest Unoccupied Molecular Orbital (LUMO) were extracted (Zhang and Musgrave, 2007). By analyzing the HOMO-LUMO gap and orbital spatial distribution, the electronic excitation properties of the core components and their electron transfer activity at the GSK3β protein binding interface were evaluated. In addition, the HOMO-LUMO orbital distribution and electron density plots were analyzed and visualized using Multiwfn_3.8 (Lu and Chen, 2012; Lu, 2024).

In this study, the binding pattern and affinity of 13 core components from Astragalus and Codonopsis Radix with their core targets were evaluated by molecular docking technology. Ultimately, two ligand-protein complexes with stronger binding effects were selected and assessed. As shown in Figure 6, among the 13 core compounds, compound Rh (rhamnocitrin) exhibited relatively stable binding with the GSK3β protein, especially the formation of hydrogen bonds with residues Lys205 and Asn213, and interacting with other amino acid residues through van der Waals forces, resulting in a docking score of 119.1087. This indicates that it has strong binding affinity. In contrast, the binding of compound Fa (folic acid) with protein GSK3β involved more amino acid residues, forming multiple hydrogen bonds between Ser66, Gly63, Phe67, Gly65, and Asp133, with a docking score of 130.2159, indicating higher binding stability.

The alanine scanning results indicate that the binding energies of mutations at residues Arg180, Lys165 and Val214 on chain A, and Lys292 and Phe293 on chain B, were all greater than 0.5 after alanine substitution for the Rh-GSK3β complex (Table 4). This suggests that the binding interaction between GSK3β protein and rhamnocitrin was weakened, leading to a disruption in the stability of the complex after the mutation. Similarly, in the Fa-GSK3β complex (Table 5), when the residues Gly65, Val70, Lys85, and Leu132 were all replaced with alanine on protein chain A, the binding energies increased to varying degrees, which meant that the stability of the complex was more destroyed. These results suggest that the aforementioned amino acids play a magnificent role in the binding process between folate and the key protein GSK3β, which is consistent with the previous molecular docking findings.

Based on the results of molecular docking, we simulated the molecular dynamics of two groups of receptor-ligand complexes (Rh-GSK3β, Fa-GSK3β) over a period of 100 ns The complexes formed by rhamnocitrin (Rh) and folic acid (Fa) with GSK3β exhibited structural stability, as demonstrated by molecular dynamics simulations. Figures 7A,B show the changes of RMSD (Root Mean Square Deviation) of the complex formed by the combination of active ingredient Rh and Fa with protein GSK3β in a vivo environment, respectively. In the FA-GSK3β complex, the binding of compound Fa with protein GSK3β shows a significant conformational reversal at the 20th ns, after which the RMSD fluctuation maintained around 3.6 Å steadily, exhibiting good uniformity, and the conformational changes of receptor and ligand were synchronized. The RMSD fluctuation of complex Rh-GSK3β began to rise around the 30th ns, and the fluctuation became stable at the 40th ns. Finally, the conformation changes of protein stabilized around 4 Å, while the conformation change of ligand stabilized around 3.2 Å. Subsequently, we analyzed the changes in RMSF (Root Mean Square Fluctuation) of the two groups of complex systems (Figures 7C,D), the value of which represents the flexibility of the corresponding amino acid residues. The results show that amino acid residues in contact with the core components in contact with the core components exhibit minimal fluctuation, suggesting lower flexibility and higher rigidity, that is, the fluctuation of the residues of the receptor remain relatively stable, and the binding of the core components Rh and Fa to the target protein GSK3β maintain stable, reducing the possibility of deviation due to large fluctuations.

In order to further investigate the physical and chemical pathways of selected ligand compounds binding to the receptor active pockets bind, the interaction analysis of the ligand-receptor binding interface was conducted (Figures 8A,B). In the Rh-GSK3β system, the hydrogen bond interaction between the compound Rh and residues of Asn95 and Gln99 is formed, with the duration of this interaction exceeding 100% and 50%, respectively. In addition, the ligand Rh exhibits hydrophobic interactions with residues Val214, Tyr216, Ser261, Val263 and Asn285. In FA-GSK3β system, ligand Fa forms continuous hydrogen bonds with Ile62, Ser66 and Val135 residues, as well as with Asn64, Gln185 and Asp200 residues through water bridges. The hydrogen bond occupancy timeline diagram of the complex also confirmed this point (Figures 8C,D).

Next, we specifically focused on investigating the impact of the ligands on the interactions within the complex systems (Figure 9). The results indicate that the carbonyl, hydroxyl, and oxygen ions of the Rh ligand contribute most of the interaction of the complex interface. For example, the carbonyl group interacts with Arg180 residue of the A chain by dispersion-induced interaction in 25% of the total simulation time, while the hydroxyl group forms a water-bridge interaction with the Asp264 residue of chain B, with the interaction duration accounting for 28%. In the FA-GSK3β system, the N and O atoms of the ligand Fa play a crucial role in the binding process of the interface, forming persistent chemical interactions with several amino acid residues (such as water-bridge and salt-bridge interactions).

Finally, the torsion angle distribution histogram was used to quantitatively analyze the conformational changes of the ligands (Figure 10). The study revealed that the C2-C1-O1-C3 dihedral angle of the Rh ligand exhibits a bimodal distribution (peaks at ± 45°), indicating that the ligand is embedded in small cavities within the receptor pocket, demonstrating significant stability. In the Fa-GSK3β system, the rotation of the C-C bonds at positions 9 and 12 of the Fa ligand is significantly inhibited, suggesting that these portions are inserted into the binding pocket of the GSK3β protein, preventing bond rotation and indicating a stable binding interaction with the protein.

The computational results (Figure 11) show that the low HOMO-LUMO gaps (0.12 eV for Rh, 0.21 eV for Fa) suggest enhanced electron transfer activity, facilitating interactions with polar residues of GSK3β. The HOMO orbital of rhamnolitrin is predominantly localized on the double bond and hydroxyl group of the pyran ring, as well as the phenolic hydroxyl group of the attached benzene ring, while the LUMO orbital is distributed over the carbon-oxygen double bond of the pyran ring and the carbon-carbon single bond connecting to the benzene ring. This configuration suggests that the electron tends to transfer from the benzene ring to the pyran ring. Furthermore, its low energy gap facilitates electron transfer within the GSK3β protein’s binding pocket, leading to the rearrangement of electron clouds in the polar amino acid residues, thereby strengthening the hydrogen bonding and enhancing the binding stability between rhamnetin and GSK3β. In addition, the delocalized π-electron system in the structure of rhamnolitrin can engage in π-π stacking or cation-π interactions with aromatic residues (such as Phe and Tyr), synergistically enhancing the binding stability. The HOMO orbital of the folic acid molecule is primarily distributed on the amine group and the phenyl ring of the para-aminobenzoic acid, while the LUMO orbital is concentrated on the adjacent pteridine structure. What’s more, the low energy gap also enables electronic transfer potential. It is similar to the results observed for rhamnetin. Both molecules enhance their binding stability to the GSK3β protein through electronic transfer and π-π stacking effects.

In addition, the electron density map in Figure 12 further validates the binging mechanism of rhamnolitrin with target protein GSK3β. The electron density of folic acid is primarily concentrated on the pteridine structure and the carbonyl group of the amide bond, which may serve as potential binding sites for stable interaction with the target protein. The electron density distribution of rhamnolitrin is uniform, suggesting that its binding to GSK3β primarily relies on the delocalized π-electrons and the π-π accumulation with amino acid residues.

Through electrostatic potential analysis (Figure 13), it was found that the electron density of rhamnolitrin is primarily concentrated on the carbonyl oxygen, hydroxyl and phenolic hydroxyl groups of the pyrene ring, which facilitates stable binding to the protein binding pocket through hydrogen bonding. The electron density of folic acid is concentrated on the nitrogen atom of the pyrazine ring in the pteridine structure and the carbonyl group of the benzamide. The results of the electrostatic potential analysis are consistent with the aforementioned computational findings, indicating that rhamnolitrin and folic acid achieve stable binding to the GSK3β protein binding pocket through hydrogen bonding and π-π stacking interactions at their active sites.

This study systematically elucidates the mechanism by which the “Astragalus-Codonopsis” herbal pair regulates GSK3β-mediated insulin signaling through multi-component synergy by the use of multi-level computational biology validation. The network pharmacology revealed that GSK3β is at the core of PPI network, which is highly consistent with the clinical studies linking GSK3β overexpression to insulin resistance. Molecular docking showed that rhamnolitrin and folic acid can bind stably to the binding pocket of GSK3β protein through various interactions such as hydrogen bonds. Alanine virtual scanning further validated the key amino acid residues of protein GSK3β that can bind with these core active ingredients.

The subsequent kinetic simulation revealed the dynamic interaction patterns of rhamnolitrin and folic acid with GSK3β protein in physiological environment. The RMSF and RMSD plots displayed the conformational changes of the complex formed between the compounds and GSK3β protein after the simulation started, with the receptor-ligand complex maintaining stability at the end of the simulation. Furthermore, the hydrogen bonding time diagram and key residue interaction analysis both elucidated the primary chemical interactions between two core components and the GSK3β protein binding pocket. The subsequent torsion distribution histogram visually demonstrated that rhamnolitrin and folic acid were stably inserted into the active cavity of the protein GSK3β, forming a stable “lock-key chelation structure”. Finally, the DFT theoretical calculation analyzed the quantum mechanical properties of rhamnolitrin and folic acid, further clarifying the binding mechanism: the phenolic hydroxyl group and carbonyl oxygen in the structure of rhamnolitrin serve as binding sites, potentially forming hydrogen bonds with the polar amino acid residues of the GSK3β protein. Meanwhile, the delocalized π electrons will also generate π-π accumulation to further stabilize the binding process. The folic acid achieves stable binding with the GSK3β protein through the pteridine part of its structure and the carbonyl oxygen of its amide bond. This may be the key mechanism by which the “Astragalus- Codonopsis” herbal pair reduces (QIAN, X, et al., 2019) GSK3β kinase activity and regulates the insulin receptor signaling pathway, thereby exerting its hypoglycemic effect.

This research has carried out methodological innovation and technical integration based on similar studies in the existing literature. Compared with the conventional research paradigm that combines network pharmacology with molecular docking, this project has constructed a multi-dimensional research system: By integrating technical means such as network pharmacology prediction, molecular docking verification, dynamic analysis of molecular dynamics simulation, key site identification of alanine scanning mutation technology, and DFT theoretical calculation, the multi-target action network and molecular mechanism of drug synergistic therapy for type 2 diabetes mellitus (T2DM) were systematically revealed. This multi-dimensional research strategy not only provides a theoretical basis for clinical transformation but also points out multiple potential targets for the development of new anti-diabetic drugs.

In this study, the molecular mechanism model of “Astragalus-Codonopsis” on the treatment of type 2 diabetes was established: rhamnolitrin and folic acid can bind to GSK3β protein through hydrogen bonding, van der Waals force and hydrophobic interaction, and regulate the insulin receptor signaling pathway to play a pharmacological role in its hypoglycemic effect. The results of this study not only provide a new method for the modernization of TCM compounds, but also lay a theoretical foundation for the development of multi-target hypoglycemic preparations based on natural products.