Exploring of hypoglycemic mechanism of a Chinese Medicine Xiao-Ke-An based on target dipeptidyl peptidase-4 ： a molecular docking and molecular dynamics study

: Xiao-Ke-An, a traditional Chinese medicine formula against diabetes in China, is reported inhibitory activity on Dipeptidyl peptidase-4 (DPP-4). However, its molecular mechanism is still unclear due to the complexity of its components. In this study, a combination of molecular docking, molecular dynamics (MD) and molecular mechanics/Poisson Boltzmann surface area (MM-PBSA) method was conducted. The results show that multiple components of Xiao-Ke-An can take effect on DPP-4 together, which could be the part of the hypoglycemic mechanism of Xiao-Ke-An. Of the eight selected compounds, MOL003714 shows the highest activity for inhibiting DPP-4 with the binding free energy of -44.14 kcal·mol -1 via binding to the S 2 and S' 1 subsites of DPP-4, which could be a potent inhibitor for DPP-4. This study provided theoretical basis for the discovery and modification of natural DPP-4 inhibitors and the deeply insight into Xiao-Ke-An. In this study, DPP-4 was used as the target protein, eight representative compounds, MOL003959 and MOL007062 were screened to explore the possible hypoglycemic mechanism of Xiao-Ke-An. It is concluded that (1) The spacious active site of DPP-4 makes it possible that multiple components in Xiao-Ke-An could inhibitor DPP-4 synergistically with different binding modes; (2) The hydrophilic S 2 subsite is a key subsite of DPP-4, so polyhydroxy components prefer to inhibit the activity of DPP-4; (3) MOL003714, which has several phenolic hydroxyl groups and glycocyclic hydroxyl groups and the lowest binding free energy, could be developed into a potent inhibitor for DPP-4. Binding modes of the eight representative compounds of Xiao-Ke-An with DPP-4 theoretically explained hypoglycemic mechanism in molecular level. We expected that Xiao-Ke-An will get more attention in the future. results


Introduction
Glucagon-like peptidase-1 (GLP-1) is a 30-amino-acid incretin secreted by L-cells in the lower digestive tract. It can promote the transcription of insulin genes and proliferation of β-cell and increase the biosynthesis and secretion of insulin by acting on pancreatic β-cells in a glucose-dependent manner, at the same time, it can inhibit the secretion of glucagon of α-cell.
Thereby, GLP-1 plays an important role in maintaining blood glucose level. Dipeptidyl peptidase-4 (DPP-4) is a kind of serine proteinase mainly expressed in small intestine, liver, kidney, pancreas and spleen. In vivo, it can resolve GLP-1, resulting in the decreasing of insulin secretion and the occurrence of diabetes [1][2][3]. DPP-4 has become a hot target in the development of hypoglycemic drugs recently. Its active pocket is composed of five subsites, S 1, S 2, S' 1 , S' 2 and S 2 -extensive subsites ( Figure 1). It recognizes N-terminal of GLP-1 from amino terminus preferentially having proline or alanine residues at penultimate position through the S 2 subsite, then the dipeptide was cleaved by Ser630 of the S 1 subsite [4]. Based on these subsites, especially the S 1 and S 2 subsites, several DPP-4 inhibitors have been successfully developed and available on market, such as sitagliptin, linagliptin, vildagliptin, saxagliptin, etc [3,5]. These drugs have different binding modes in DPP-4 active pocket [6][7][8]. Vildagliptin and saxagliptin were designed as peptidemimetics inhibitors, whose cyanopyrrolidine moieties bind to the S 1 subsite and form a covalent bond between nitrile group and Ser630. The hydroxy adamantyl groups bind to the S 2 subsite, and the nitrogen atom of amino group forms a salt bridge that connected Glu205 and Glu206 [9,10]. In addition to the S 1 and S 2 subsites, alogliptin and linagliptin bind to the S' 1 subsite, forming π-π interaction between Tyr547 and their uracil rings, and the phenyl of the quinazoline substituent in linagliptin forms a π-π interaction with Trp629 in the S' 2 subsite [6,11]. The triazolopyrazine moiety of sitagliptin, 1-phenylpyrazol and piperazine moiety of teneligliptin bind to the S 2 -extensive subsite (Arg358, Phe357) mainly by forming hydrogen bond or π-π interaction with Phe357 [12,13]. The less cardioprotective effect of DPP-4 inhibitors and their low risk of causing hypoglycaemia makes them become widely used as second-line agents [5,14,15] . Other excellent natural DPP-4 inhibitors, such as chrysin [16], galangin [17], plectranthoic acid [18], iso-daphnetin [19,20], polyphenols, flavonoids [21] and cyanidin 3,5-diglucoside [22] [26], is widely used for the treatment of type II diabetes mellitus in China with low cost and minimal side effect [27]. Wu Xue and his colleagues verified that some components in Xiao-Ke-An have DPP-4 activity, and screened out some salvianolic acids and saponins with strong inhibitory effect on DPP-4 [28]. The subsequent studies proved that many components in Xiao-Ke-An could activate AKT/GSK-3β pathway by inhibiting DPP-4 and improving glucose and lipid metabolism [27]. However, systematical screening of all possible components with DPP-4 activity in Xiao-Ke-An has not been studied and their binding modes between the active components and DPP-4 are still unclear.
In this study, a database containing 819 components of Xiao-Ke-An for virtual screening was obtained by molecular docking method, and the binding modes between active components and DPP-4 were furtherly explored through molecular dynamics simulation. It was expected the results could explain the hypoglycemic molecular mechanism of Xiao-Ke-An and promote the discovery of new DPP-4 inhibitors and the deeply insight into Xiao-Ke-An. The 3D structure of the components of Xiao-Ke-An were downloaded from the traditional Chinese medicine systems pharmacology database and analysis platform(TCMSP, https://tcmspw.com/tcmsp.php) [29], and optimized by using MMF94 force.

Preparation of DPP-4
The crystal structure of DPP-4 was derived from RCSB protein data bank database(PDB ID: 3VJK) [30]. 3VJK is composed of DPP-4 and its ligand, an inhibitor named tenegliptin. It has four same chains, A, B, C and D, and the active site is in each chain, so only chain A was used for molecular docking and the following molecular dynamics simulation in this study. Prior to the molecular docking, crystal water molecules and the ligand were removed, and hydrogen atoms were added with Autodock Tools 1.5.6 package[31].

Molecular docking
Autodock Vina 1.1.2 [32] was used for molecular docking. All small molecules were kept flexible and the receptor DPP-4 was kept rigid. The exhaustiveness parameter was set to 10, the center coordinates of the box was ( 48.68, 62.15, 32.28), the number of grid points of the box was (30,26,30), the grid spacing was set to 1 Å, and the number of generated modes was set to 10. All other parameters were set as default values.

Molecular dynamics simulation
Molecular dynamics simulations were conducted with Amber18 package [33,34]. The optimal conformation of each docking result was used for the initial structure of dynamic simulation. Prior to the dynamic simulation, Gaussian 09 [35] was used to construct small molecule force field under the level of HF/6-31G*, and RESP (Restrained ElectroStatic Potential) was used to generate atom charge. Topology and coordinate files were generated through tleap. All parameters were established under ff14SB and GAFF force field [36]. Each system was solvated in an octahedral water box with TIP3P water model and periodic boundary of 10 Å. Na + was added as the counter-ion to neutralize the system.
The energy minimization was divided into three steps: Firstly, energy minimization with the complex constrained. The maximum number of cycles was set to 4000 steps. Steepest descent method was used for the first 2000 steps, and the conjugate gradient method for the next 2000 steps. Secondly, energy minimization with the protein skeleton constrained. The maximum number of cycles was 5000 with the steepest descent method for the first 2500 steps and the conjugate gradient method for the 2501-5000 steps. Finally, energy minimization without any constraints. The maximum number of cycles was 10000. The steepest descent method was used in the first 5000 steps, and the conjugate gradient method was used in 5001 to 10000 steps.
The system was firstly heated for 40 ps from 0 K to 300 K in a NVT ensemble by Berendsen control method, then equilibrated at 300 K for another 20 ps. 100 ns production simulation was finally performed at constant pressure and constant temperature under periodic boundary conditions by Berendsen control method.

Binding free energy calculation
The final 5 ns stable MD trajectory was selected to calculate the binding free energy with MM-PBSA [37] method in Amber18. The binding free energy (∆G binding )is calculated by the following equations [38]: ∆G complex , ∆G protein and ∆G ligand are the free energies of complex, protein and ligand respectively. ∆E gas and ∆G sol are the gas phase interaction energy and the solvation free energy respectively; ∆E ele and ∆G vdW are the electrostatic energy and van der Waals energy respectively; ∆G polar and ∆G nopolar are the polar contribution and van der nopolar contribution respectively. The entropy contribution (∆S) was ignored in this study, considering that the computation of entropy contribution by normal mode analysis (NMA) approach is computationally expensive and tends to have a large error that introduces significant uncertainty in the result [39].

Results and discussion
3.1 Molecular Docking 45 compounds with higher docking score than tenegliptin (ΔG score = -8.4 kcal·mol -1 ) were initially screened from molecular docking and furtherly screened according to drug-likeness (DL) parameter (DL>0.18) in TCMSP database. 17 compounds were finally obtained. These compounds were divided into three groups, phenols group, steroids group and other group based on their structure characteristic. Finally, eight compounds with the lowest binding energies of each sub-group were selected for the following research, which are shown in Table 1. Phenols group includes dimetbyl lithosper-mate B, jionoside A and physcion-8-O-beta-D-gentiobioside; Steroidal group includes polygosides A and ginsenoside Rg5, other group contains three compounds, neoprzewaquinone A, dauricine and limonin.   (Figure 2g, Figure 3g). In the eight systems, MOL007062 and MOL003959 are the most stable, which could be seen from both of their RMSDs and the superposition results (Figure 2f, 3f, 2h, 3h).

Binding free energy
The binding free energies calculation for 2500 snapshots of the eight complexes were performed from the final 5 ns simulation by MM-PBSA method, and the results were shown in Table 2. As shown in Table 2

Decomposition free energy and binding mode analysis
To explore the interaction modes of each molecule that binds to DPP-4, decomposition free energy of the selected residues were calculated by using MM-PBSA method. The decomposition free energy of the residues are shown in Figure 5. kcal·mol -1 , -1.90 kcal·mol -1 , -3.10 kcal·mol -1 and -1.09 kcal·mol -1 , respectively. Other residues, such as Tyr666 and Arg669, maybe also form van der Waals force with MOL005400.
MOL003959, a terpenoid compound, is composed of alkylene oxide and cycloalkanone. In

Comparison of Binding mode
The binding sites of each compounds are summarized in Table 3  ----+ Note: "+" represents the subsite is occupied "-" represents the subsite is not occupied subsites gradually in the process of dynamic simulation, also not be the potential DPP-4 inhibitors.

Conclusion
In this study, DPP-4 was used as the target protein, eight representative compounds,