Prediction of hyaluronic acid target on sucrase-isomaltase (SI) with reverse docking and molecular dynamics simulations for inhibitors binding to SI

Auricularia cornea (E.) polysaccharide is an important component of A. cornea Ehrenb, a white mutant strain of Auricularia with biological activities, such as enhancement of human immune function and cancer prevention. The hyaluronic acids (HAs) are important components of the A. cornea polysaccharide and have extremely high medicinal value. In this study, we used HA to search the target protein sucrase-isomaltase (SI). In addition, we also performed molecular dynamics (MD) simulations to explore the binding of three inhibitors (HA, acarbose and kotalanol) to SI. The MD simulations indicated that the binding of the three inhibitors may induce the partial disappearance of α helix in residues 530–580. Hence, the hydrogen bond for Gly570-Asn572, which was near the catalytic base Asp471 in SI, was broken during the binding of the three inhibitors. We reveal a new inhibitor for SI and provide reasonable theoretical clues for inhibitor binding to SI.


Introduction
Auricularia auricula is an edible and medicinal fungus ranking fourth in production among the worldwide [1][2][3][4]. Scientific research is extensive on this species because of the abundant resources of A. auricula in China [5,6]. cornea Ehrenb [7]; the variant is highly nutritious and has high economic value [8].
In 2018, Li Yu et al. pointed out that A. cornea Ehrenb suppressed the levels of total cholesterol and triglyceride and enhanced levels of hepatic glycogen and high-density lipoprotein cholesterol, which may be involved in diabetes mellitus (DM) [7], which is a progressive metabolic disease [9].
The polysaccharide of A. cornea Ehrenb has been widely studied; it is an important chemical component for regulating human life activities. The polysaccharide content of A. cornea Ehrenb is reportedly better than other Auricularia species [7] is widely used in medicine, food science and other fields. However, we still do not know which components in the polysaccharide of A. cornea Ehrenb are involved in metabolic disease. Hyaluronic acid (HA) is the important compound of polysaccharide of A. cornea Ehrenb. HA belongs to the glycosamino glycans connected with β 1,4 glycosidic bond. Considering its unique physical and chemical properties, HA is widely used in medicine [10][11][12]. So, in this study, HA was used to search the target protein in metabolic disease. Sucrase-isomaltase (SI) functioned as an attractive target for inhibition by α-glucosidase inhibitors as a means of controlling blood glucose levels in individuals with type 2 diabetes [13,14]. In addition, we also used molecular dynamics (MD) simulations to explore the binding mode of the three inhibitors (HA, acarbose and kotalanol) to SI. Our results will provide useful clues for further A. cornea Ehrenb study.

Reverse docking
The 2D structure of HA was used for SwissTargetPrediction [15] to search for the target protein.

Protein preparation
The 3D structure of the SI was obtained from the RCSB Protein Data Bank (www.rcsb.org) (PDB ID: 3LPO) [13]. The 3D structures of HA, acarbose and kotalanol were downloaded from the Chemspider database (www.chemspider.com). In this study, AutoDock Vina [16] was used to construct protein-inhibitor complexes. In the AutoDock Vina configuration files, the parameter num_modes was set to 9 Å. We identified the receptor binding pocket based on the point of the substrate binding to the SI. Hence, we kept all the rotatable bonds in ligands flexible during the docking procedure, and we kept all the protein residues inside the binding pockets rigid. The Kollman charges were used to convert all receptors and ligands to the PDBQT format using the AutoDockTools package [17]. The search spaces were 30 × 30 × 30 Å. The docking results were clustered automatically.

Molecular dynamics simulations
All MD simulation courses were conducted by using the Amber 16 package [18]. The AMBER ff99SB forcefield [19] was applied to the SI-inhibitor complexes. Then, they were solvated using the TIP3P water model [20] with the box at 12 Å. The Counter ions (Na+ and Cl-) were assigned to neutralise the three systems. Based on each of the prepared systems, energy minimization was used for the solvent complexes. Then, the constraints were used on protein backbone atoms. During the energy minimisation, the steepest descent algorithm [20] and conjugate gradient algorithm [21] were used successively. Subsequently, each system underwent a gradual heating process for 500 ps from 0 to 300 K and then was equilibrated for at 300 K for another 500 ps. Finally, the whole system was performed 100 ns MD simulation at 1 atm constant pressure and 300K constant temperature conditions. We used the Langevin Nosé-Hoover thermostat [22,23] and the Parrinello-Rahman method [24,25] to maintain a constant temperature and a constant pressure of 1 atm in each system. The LINCS algorithm were constrained by all bonds [26]. The long-range electrostatic interactions were performed with the particle mesh Ewald method [27] with a grid size of 1.2 Å. The periodic boundary conditions were implemented in all directions along the simulation box.

Data analyses
PCA was performed using Bio3D version 2.3.0 to study the collective motions in 100 ns of protein-ligand complex [28,29]. This method used the calculation and diagonalisation of the covariance matrix. The covariance matrix was calculated as follows: where xi/xj is the coordinate of the ith/jth atom of the systems. Free energy landscape (FEL) is a map of all possible conformations of molecular entities [30] and can be used to understand the stability, folding and function of the protein. The FEL can be constructed as follows: where KB is Boltzmann constant, T is the temperature of simulation systems, and 300 K is set in the current calculations. P(X) is the probability distribution of the molecular system along the PCs.

MM-PBSA calculations
Molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) is a popular method to calculate the binding free energy between protein and ligands. It is more accurate than most scoring functions of molecular docking and less computationally demanding than alchemical free energy methods with Amber 16 package [18]. The Molecular Mechanics-Generalised Born Surface Area (MM/GBSA) method in AMBER16 package [18] was used to calculate the binding free energies. A total of 200 snapshots were chosen evenly from the MD trajectory.

Reverse docking
Firstly, HA was used to search for target protein with SwissTargetPrediction [15]. Among these results, the possibility of binding score between the target proteins and HA is not very different. SI is an attractive target for inhibition by α-glucosidase inhibitors to control blood glucose levels related to type 2 diabetes [14]. In 2018, Li Yu et al. pointed out that A. cornea Ehrenb was involved in DM [7], a progressive metabolic disease [9]. In addition, HA is a disaccharide; so, SI was selected for further study.

Binding pose of the inhibitors to SI
SI is composed of the N-and C-terminal duplicated catalytic domains ( Fig 1A). The N-terminal catalytic domain of N-terminal SI has a broader specificity for both 1,4-and 1,6-oligosaccharides [13]. Acarbose and kotalanol (the known α-glucosidase inhibitors) [13] and HA in Fig 1B-1D are docked with AutoDock Vina. The complex with SI bound with kotalanol were upload from PDB (PDB code 3LPP). The SI-kotalanol interaction was shown in S1A Fig. We also docked kotalanol to SI (S1B Fig). It can been seen that the docked kotalanol was similar to the crystal kotalanol, which indicated the docking software, AutoDock Vina, was reliable.

Conformational changes for the inhibitors binding to SI
To explore the conformational changes for the inhibitors binding to SI, MD simulations for the four systems (SI, SI-acarbose, SI-kotalanol and SI-HA) were performed with Amber 16 software [18]. The parameters of MD simulations are listed in Table 1. According to Fig 3A  (RMSD plot), all four complexes were stable. MD simulation for free SI at a trajectory of 100 ns and with RMSD value of 2.4 Å was used as a reference. According to the curves where SI combined with any inhibitor, the value of RMSD became lower than that in free protein in the catalytic domain (residues 297-681). SI with kotalanol had a more drastic score that other inhibitor-enzyme complexes either in the full protein or in the catalytic domain (residues 297-681). Our results indicated that the conformational changes in a protein were more than those that occurred in a protein with inhibitors. When SI combined with an inhibitor, the conformation became stable. Kotalanol is a better inhibitor with lower Ki compared with the others, and so, it may induce larger conformational changes in SI [13].
The rigidity of the protein system was examined using R g values. The R g plot of the α-carbon atoms versus time was obtained and presented in Fig 3B. The R g values retained their stability throughout the 100 ns time of the MD simulation, which corresponded to the  simulation. These values were used to facilitate the interpretation of the secondary structure. The R g value for four systems (in full protein) was stable at about 28 Å, whereas the four systems, i.e., the catalytic domain (residues 297-681), were stable at about 21 Å. The three inhibitors made the secondary structure of the catalytic domain (residues 297-681) more stable and compact.
The overall conformational changes were further validated by the SASA graph, which was plotted against the MD simulation time, as shown in Fig 3C. The probabilities based on the SASA plots indicated that the four complexes had similar values. However, in Fig 3D, the SASA score of Tyr299, Trp327 and Lys480 were shown to be lower than that in the free SI. At the same time, Asp443, Asp542 and His600 had higher SASA score than SI. All six residues were located at the active site. These changes may cause some movement for inhibitor binding.  To explore the secondary structural changes related to the binding of three inhibitors, secondary structures during MD simulations were calculated, as shown in Fig 5A-5D. Fig 5A-5D shows the location of residue Gly530 to Glu580. The conformational changes caused by the three inhibitors were investigated and compared with those caused by the free SI. To investigate the conformational changes, we obtained the difference in the secondary structure (DSSP). The DSSP of residues Gly530 to Glu580 differed from that in the SI-inhibitor. The α helix of residues Gly530 to Glu580 in SI-inhibitors partly disappeared.
Among them, we truncated the most different parts, as presented in Table 2. In residues Gly530 to Glu580, free SI had a probability of 86.50%-100.00%, whereas the probabilities of the other systems were almost decreased. As Gly530 to Glu580 belongs to the catalytic domain (residues 297-681), so the conformation change on the 530-580 domain induced by inhibitors binding may be the main factor affecting the efficiency of enzyme catalysis.
According to Fig 6A, 6C and 6E, free SI stayed stable within the 1.5 Å RMSD value. However, the RMSD value of the three inhibitor-enzymes were kept at a high level, showing that  the three inhibitors binding to SI induced Gly530 to Glu580 domain flexibly. Fig 6B, 6D and 6F show the relative frequency of RMSD plot of Gly530 to Glu580 domain for the free SI with the SI-inhibitor complex. Thus, the three inhibitors binding to SI caused drastic changes in the Gly530 to Glu580 domain. Residues Gly530 to Glu580 contained the acid/base catalyst, Asp571. Seen from Fig 7A, Asn572 made hydrogen bonds with Gly570. The relative frequency of distance between Asn572 and the O atom of Gly570 of the four complexes are shown in Fig 7B. The hydrogen bond between Asn572 and Gly570 contained all the MD simulations (about 2 Å), whereas this may become weak in the inhibitor-SI complexes (SI-HA, about 3 Å), or it disappeared. Residues Gly530 to Glu580 contained acid/base catalyst, Asp571. The hydrogen bond weakened and disappeared, which may be useful to inhibitor binding.

Principal component and free energy landscape analysis
We probed the internal dynamics of different system, and the results were depicted in the figure. The four systems exhibited obvious differences in the correlated extents of protein motion. The Gibbs free energy landscape (FEL) was calculated using the first two principal components as reaction coordinates. Using principal component analysis (PCA), Helmoltz free energy change was calculated. The sum percentage of PC1 and PC2 for four systems were shown in Table 3, and the FELs obtained from the simulations were plotted, as shown in Fig  8A-8D. The FEL can provide remarkable information about the different conformational states accessible to the protein in the simulation. As shown in the figure, the inhibitor-protein was significantly different from the free protein. This energy minimum corresponded to a structure with some loss of irregular secondary structures, such as coils and turns. The structures of the two most stable conformations of SI revealed that the conformational changes in the α helix in residues 530-580 domain (partly disappeared) existed all SI-inhibitor complexes. This finding was consistent with the previous analysis.

MM-PBSA calculations
By calculating the potential energy in the vacuum, van der Waals, electrostatic interactions and net non-bonded potential energy between the protein and ligands were calculated, as

PLOS ONE
Prediction of hyaluronic acid target on sucrase-isomaltase (SI) shown in Table 4. An average binding energy equal to -66.47 kcal /mol was achieved for SIkotalanol, -54.39 kcal /mol for SI-HA and -49.58 kcal /mol for SI-acarbose. This finding was consistent with the previous analysis, in which HA had a higher Ki than kotalanol.

Conclusion
The HA is an important component of the A. cornea (E.) polysaccharide and has extremely high medicinal value. We used HA as the lead compound to search the target protein sucrase-