Shitao Zhang1
Lu Han
Weiwei Han- 1School of Life Science, Key Laboratory for Molecular Enzymology and Engineering of Ministry of Education, National Engineering Laboratory for AIDS Vaccine, Jilin University, Changchun, China
- 2Laboratory of Theoretical and Computational Chemistry, Institute of Theoretical Chemistry, Jilin University, Changchun, China
There are multiple drugs for the treatment of type 2 diabetes, including traditional sulfonylureas biguanides, glinides, thiazolidinediones, α-glucosidase inhibitors, glucagon-like peptide-1 (GLP-1) receptor agonists, dipeptidyl peptidase IV (DPP-4) inhibitors, and sodium-glucose cotransporter 2 (SGLT2) inhibitors. α-Glucosidase inhibitors have been used to control postprandial glucose levels caused by type 2 diabetes since 1990. α-Glucosidases are rather crucial in the human metabolic system and are principally found in families 13 and 31. Maltase-glucoamylase (MGAM) belongs to glycoside hydrolase family 31. The main function of MGAM is to digest terminal starch products left after the enzymatic action of α-amylase; hence, MGAM becomes an efficient drug target for insulin resistance. In order to explore the conformational changes in the active pocket and unbinding pathway for NtMGAM, molecular dynamics (MD) simulations and adaptive steered molecular dynamics (ASMD) simulations were performed for two NtMGAM-inhibitor [de-O-sulfonated kotalanol (DSK) and acarbose] complexes. MD simulations indicated that DSK bound to NtMGAM may influence two domains (inserted loop 1 and inserted loop 2) by interfering with the spiralization of residue 497–499. The flexibility of inserted loop 1 and inserted loop 2 can influence the volume of the active pocket of NtMGAM, which can affect the binding progress for DSK to NtMGAM. ASMD simulations showed that compared to acarbose, DSK escaped from NtMGAM easily with lower energy. Asp542 is an important residue on the bottleneck of the active pocket of NtMGAM and could generate hydrogen bonds with DSK continuously. Our theoretical results may provide some useful clues for designing new α-glucosidase inhibitors to treat type 2 diabetes.
Introduction
At present, there are multiple drugs for the treatment of type 2 diabetes, including traditional sulfonylureas (Stephen et al., 2018), biguanides (Schäfer, 1983), glinides (Chen et al., 2015), thiazolidinediones (Nanjan et al., 2018), α-glucosidase inhibitors (Kazufumi et al., 2014; Patel, 2015), glucagon-like peptide-1 (GLP-1) receptor agonists (Drucker, 2018), dipeptidyl peptidase IV (DPP-4) inhibitors, and sodium-glucose cotransporter 2 (SGLT2) inhibitors (Thornberry and Gallwitz, 2009; Kelly et al., 2019). These therapeutic drugs have been widely used in clinical trials because of their own characteristics in hypoglycemic control. For example, α-glucosidase inhibitors have been used to control postprandial glucose levels caused by type 2 diabetes since 1990 (Ríos et al., 2015; Flores-Bocanegra et al., 2017; Santos et al., 2018; Dhameja and Gupta, 2019; Usman et al., 2019; Mi et al., 2021; Tuyen et al., 2021). Acarbose (Chiasson et al., 2002) and miglitol (Satoru et al., 2015), which were clinically used for treating type 2 diabetes, may control blood glucose levels by targeting α-amylases and α-glucosidases (Lyann et al., 2010; Ren et al., 2011).
Glycoside hydrolases play significant roles in human metabolism, including digestion and decomposition of polysaccharides and biosynthesis of glycoprotein (Lovering et al., 2005). α-Glucosidases are rather crucial in the human metabolic system and are principally found in families 13 and 31 (Lovering et al., 2005). Maltase-glucoamylase (MGAM) (Sim et al., 2008) and sucrase-isomaltase (SI) (Sim et al., 2010) belong to glycoside hydrolase family 31. The main function of MGAM and SI is to digest terminal starch products left after the enzymatic action of α-amylase, which becomes an efficient drug target for insulin resistance (Van Beers et al., 1995). MGAM contains the following units: a small cytosolic domain of approximately 26 residues, a transmembrane domain (TMD) containing about 20 residues inserted into the intestinal epithelial cell membrane, an O-glycosylated linker, and two independent catalytic subunits: NtMGAM and C-terminal luminal subunit (CtMGAM) (Sim et al., 2008; Lyann et al., 2010) (Supplementary Figure S1). NtMGAM containing 864 residues (PDB ID: 3L4U) (Lyann et al., 2010) were used in this study. NtMGAM, known to be retaining α-glycosidases (Satoh et al., 2016), has received relatively little attention despite its importance and the number of different activities from a range of organisms, including mammals, plants, and microorganisms (Frandsen and Svensson, 1998; Yu et al., 1999; Lovering et al., 2005). The substrate specificities of MGAM vary and overlap from maltose (Quezada-Calvillo et al., 2008) to isomaltose (Elferink et al., 2020) and other small oligosaccharides.
Salacia reticulata (S. reticulata), a plant widely distributed in China and some other countries in Southeast Asia, is used as a traditional prescription for treating type 2 diabetes (Medagama, 2015). Sulfonium ion-containing compounds were isolated from aqueous extracts of S. reticulata by Jayakanthan and his colleagues (Jayakanthan et al., 2009), including de-O-sulfonated kotalanol (DSK). DSK comprises a 1,4-anhydro-4-thio-d-arabinitol core and a polyhydroxylated acyclic chain (Lyann et al., 2010) and can act as an inhibitor of α-glucosidase in human bodies. A previous study reported that DSK could be an efficient inhibitor of NtMGAM with a Ki of 0.03 (±0.01) μM (Lyann et al., 2010).
In order to explore the conformational changes in the active pocket and unbinding pathway for NtMGAM, molecular dynamics (MD) simulations (Zhu et al., 2018; Zhu et al., 2019; Liu et al., 2020) and adaptive steered molecular dynamics (ASMD) simulations (Zhu et al., 2018) were performed between two inhibitors (DSK and acarbose) and NtMGAM (PDB ID: 3L4U) (Lyann et al., 2010) (workflow listed in Supplementary Figure S2). Our results may provide new ideas for the further design of α-glucosidase inhibitors.
Materials and Methods
Preparation for the Structure of Protein Inhibitors
AutoDock 4.2 (Morris et al., 2009) was used for docking acarbose with NtMGAM using the Lamarckian genetic algorithm to identify a proper binding conformation with a grid box size of 66 Å × 58 Å × 66 Å points and a grid point spacing of 0.375 Å. The binding conformation with the lowest energy was chosen for simulations. The crystal structure of NtMGAM with DSK complex and the 3D structure of acarbose (PDB ID: 3JYR) (Vahedi-Faridi et al., 2010) was downloaded from Protein Data Bank (www.rcsb.org) for further studies.
MD Simulations
Simulations in our study were performed using the Amber16 package (D.A. Case et al., 2016) with the Amber ff99SB force field (Lindorff-Larsen et al., 2010). At the same time, the Gaff2 force field (Wang et al., 2004) was utilized to generate the parameterization of DSK and acarbose. It is well known that charged residues affect the environment of protein (Popović and Stuchebrukhov, 2004; Tashiro and Stuchebrukhov, 2005; Sugitani and Stuchebrukhov, 2009). H++ is an online tool that can automatically compute pKa values of ionizable groups in proteins (http://biophysics.cs.vt.edu/). We computed the ionizable groups of NtMGAM on H++ and then manually fixed the ionizable groups. All three complexes were analyzed using the MD simulations in a cubic periodic boundary box with the TIP3P water model (Bogunia and Makowski, 2020), which was prolonged to 12 Å from the protein atoms. Sodium ions were randomly added to the simulation systems for neutralization. To get the initial equilibrious structure, energy minimization was performed through the steepest descent method in 1,000 cycles. Subsequently, 50 ps of NVT (Berendsen temperature coupled with constant particle number, volume, and temperature) (Berendsen et al., 1984) and 50 ps of NPT (Parrinello–Rahman pressure coupled with constant particle number, pressure, and temperature) (Andersen, 1980) were performed to maintain the stability of the system (300 K, 1 bar). After stabilizing all thermodynamic properties, a 200 ns unconstrained MD simulation was performed with a time interval of 2 fs. The coordinates for all models were stored every 2 ps. During the simulation, the following options were specified: (I) bonds involving hydrogen are constrained and bond interactions involving H-atoms were omitted using the SHAKE algorithm (Miyamoto and Kollman, 1992), (II) the particle mesh Ewald summation algorithm (Essmann et al., 1995) was taken to calculate the long-range electrostatic interactions, (III) 1 atm constant pressure was maintained by the Langevin dynamics method (Rosenberg et al., 1986) (Guàrdia and Padró, 1985), and (IV) an optimum temperature (300 K) was maintained. MD simulations were performed three times for each system in this study (Supplementary Figures S3, S4). The root-mean-square deviation (RMSD), radius of gyration (Rg), root-mean-square fluctuation (RMSF), and solvent-accessible surface area (SASA) values were calculated using VMD (Humphrey et al., 1996).
MM-PBSA Calculations
MMPBSA.py (Miller et al., 2012) in the AmberTools17 package was employed to conduct free energy calculations for the two complexes. 200 conformations were extracted from each equilibrious trajectory (from 100 to 200 ns with an interval of 50 frames) for calculations. The binding free energies were calculated by subtracting the free energies of the receptor and the ligand DSK or acarbose from the free energy of the bound complex of two systems:
Then, the free energy change associated with each term of ΔG was calculated according to the following:
where ΔGsolvation represents true free energy. To determine the relative stability, end-state method calculations were performed to estimate the energies, according to averages from the ensemble of these snapshots:
where i is the index of a particular frame and N is the total number of frames analyzed. The gas-phase energies (Egas) can be computed from the quantum mechanical (QM) calculations and used further as a part of the force field parameterization; therefore, the Egas energies can be abstracted from the molecular mechanical (MM) energies and the corresponding force field (Miller et al., 2012).
Principal Component Analysis and Free Energy Landscapes
PCA is a common statistical multivariate method, which can select the structure of each frame in an MD trajectory as a new set of variables, called principal components (PCs), with a minimal loss of information. In our study, we employed the bio3d package in R to refine structural superposition and examine the relationship between different conformers (Grant et al., 2006). The current protocol excludes the residues displaying the largest positional differences at each round and identifies only the core residues. Following the superposition of core residues, PCA was performed to examine the conformers based on their equivalent residues. The PCs collected during MD simulation are originally the eigenvector values collected from the covariance matrix, each corresponding to the change in protein trajectory (Al-Khafaji and Taskin Tok, 2020). In order to obtain a lower-dimensional representation of the structural dataset, we project the distribution onto the subspace defined by the largest principal components. In each dimension, the corresponding eigenvalue represents the percentage of the total mean square displacement (or variance) of atom positional fluctuations. In PCA, very few dimensions are generally enough to capture about 70% of the total variance in the structures to be studied (Grant et al., 2006). Therefore, the first few eigenvectors are sufficient to provide a useful description while still retaining most of the variance in the original distribution. After PCA, the FEL (Frauenfelder et al., 1991) was obtained by calculating the joint probability distribution from the essential plane constructed from the top two eigenvectors (Singh et al., 2015).
Pathways Identified With CAVER
The software called CAVER has been widely used to analyze and visualize possible cavities and channels in protein structures. CAVER Analyst 2.0 (Jurcik et al., 2018) was employed to determine the channel position. The starting point was set at the position of the ligand for the channel computation. The minimum probe radius, clustering threshold, shell depth, and radius were set to be 0.9, 3.5, 4, and 3 Å, respectively. All the other parameters were used as default. 1,000 snapshots (we took one frame every ten frames in 10,000 from the 100–200 ns simulations) were analyzed utilizing CAVER Analyst 2.0.
Adaptive Steered Molecular Dynamics Simulations
ASMD simulations (Ozer et al., 2010; Ozer et al., 2012a; Ozer et al., 2012b; Ozer et al., 2014) were performed for NtMGAM with the two different ligands using the Amber 16 package. ASMD has been shown to alleviate the problem that many simulations must be run to converge the potential mean of force (PMF) in steered molecular dynamics (Izrailev et al., 1999) by dividing the predetermined reaction coordinate into numerous smaller stages. Each stage in the ASMD simulation contained multiple simulations that should be performed parallelly. In each stage, the trajectory with the work value closest to the Jarzynski average (JA) (Jarzynski, 1997) should be selected, and the coordinates at the end of that trajectory should be used as an initial coordinate for the next stage. Then, the JA structures were used for PMF calculation. Here, in our study, the distance between the ligand and NtMGAM in the initial conformation is 6 Å. The reaction coordinate in ASMD is predetermined to be set at 20 Å, at which the inhibitor is considered to escape from the enzyme. The stretching velocity was 10 Å/ns in this ASMD simulation, coupling a spring constant k of 40 kcal/(mol×Å2). Each simulation was split every 2 Å into seven stages; each stage contains 14 simulations to reach the final reaction coordinate of 20 Å. As the distance between the two selected atoms reached 20 Å, there were no longer any interactions between the ligand and the receptor.
Results and Discussion
Docking Pose and System Stabilize
It was reported that the structure of the NtMGAM substrate-binding site consisted of two sugar-binding sites, which in the acarbose-NtMGAM binding structure were occupied by the two nonreducing rings of acarbose (Lyann et al., 2010). To determine the docking pose, we chose to dock the DSK crystal structure to NtMGAM with AutoDock 4.2 (Morris et al., 2009). Comparing the crystal structure of DSK-NtMGAM and the docked pose (Supplementary Figure S5), it can be concluded that they were similar (RMSD 0.50), which indicated that this system may use AutoDock 4.2 software to determine the binding pose. Acarbose was docked to NtMGAM using the same method.
NtMGAM has a typical catalytic (β/α)8 barrel domain (Figure 1A). DSK and acarbose were docked in an active pocket of NtMGAM (Figure 1B). Figures 1C,D shows that the active residues around acarbose and DSK were bound around NtMGAM. It can be seen that His600, Arg526, Asp542, Asp203, Trp406, and Asp327 were anchor residues for DSK binding in Figure 1D. In contrast, in the acarbose-NtMGAM, His594, Asp321, Arg520, Asp197, and Asp536 made hydrogen bonds with acarbose, indicating that they are important residues for acarbose binding to NtMGAM.
FIGURE 1. (A) Catalytic domain of NtMGAM. The catalytic (β/α)8 barrel domain is tagged and colored differently. (B) Surface diagram of NtMGAM. (C) Binding of acarbose to NtMGAM. The active residues around acarbose binding to NtMGAM. (D) Binding of DSK to NtMGAM. The active residues around DSK binding to NtMGAM.
The MD simulations for three systems have been performed three times. The RMSD values of the Cα atom backbone of residues of three systems were calculated to evaluate the equilibrium of systems (Supplementary Figures S3, S4 and Supplementary Table S1). It can be seen in Supplementary Figure S3 that all MD simulations have got equilibrium. In Supplementary Figure S4, the RMSD of the three simulations of free-NtMGAM are slightly different from each other, with the average values of 1.67, 1.94, and 2.13 (from 100 to 200 ns), respectively. Although the RMSD values of each of the three repetitions seldom cross each other, the three simulations have all reached a state of relative equilibrium. The group with the most equilibrious and generally balanced RMSD was chosen for further study (Supplementary Table S1). The parameters of 200 ns MD simulations for three systems were listed in Supplementary Table S2. From Figure 2A, the RMSD values of NtMGAM, DSK-NtMGAM, and acarbose-NtMGAM are stabilized at about 1.67, 2.16, and 2.15 Å, respectively (Figure 2B), suggesting the structures of the three systems had reached a state of relative equilibrium. The Rg values of the three systems are shown in Figures 2C,D and are finally stabilized at 28.75 Å. In Figure 3A, the SASA values of acarbose-NtMGAM were stabilized at 34,000 Å2; meanwhile, the other two systems were stabilized at 33,000 Å2. These results indicate that three systems have attained stability and their states reached equilibrium. SASA values of single residue of the three systems were calculated for core residues. Tyr299, Phe412, Asn433, Asp542, and Phe575 had higher scores in the DSK-NtMGAM complex than the others. It was reported that Asp542 and Phe575 were important residues interacting with the hydroxyl groups of inhibitors binding to NtMGAM (Usman et al., 2019). Our results were consistent with the experimental data. In summary, all the systems were stabilized and can be used for further study.
FIGURE 2. Root-mean-square deviation (RMSD) and radius of gyration (Rg) analysis of NtMGAM with or without different ligands throughout 200 ns. (A) RMSD plot. (B) The relative frequency of the RMSD plot. (C) Gyration radius (Rg) plot. (D) The relative frequency of the Rg plot.
FIGURE 3. (A) Solvent-accessible surface area (SASA) analysis of three systems over 200 ns MD. (B) SASA of active residues in different compound structures.
Conformational Changes for Inhibitors Binding
To evaluate the deviation amount of displacement in three trajectories from MD simulations, atom positional RMSF values were calculated for the backbone atoms of three systems (Figure 4A). Figures 4B–D show that residue displacements correspond to the motions described by the first eigenvector for three complexes. These displacements represented the relative displacement of each residue caused by the motion described by a given eigenvector. It was reported that NtMGAM contains 868 residues, which can be divided into five structural domains: (I) a trefoil Type-P domain (residue No. 1–51); (II) N-terminal β-sandwich domain (residue No. 52–269); (III) catalytic (β/α)8 barrel domain (residue No. 270–651) with two inserted loops [inserted loop 1 (residue No. 367–416) and inserted loop 2 (residue No. 447–492)] protruding out between β3 and α3 and between β4 and α4, respectively; (IV) proximal C-terminal domain (residue No. 652–730); (V) distal C-terminal domain (residue No. 73–868), both with β-sandwich topologies (Sim et al., 2010) (see Figure 5A). In the DSK-NtMGAM complex, residues Ser376, Gly397, and Trp406, which are located at inserted loop 1 domain, exhibited distinct atom positional fluctuation amplitudes. This displacement widened the active site pocket, which would affect the inhibitors' binding.
FIGURE 4. (A) Average atom positional root-mean-square fluctuations (RMSF) of the backbone atoms per residue for the inhibitors bound NtMGAM. (B–D) The eigenvector components for atomic displacement along the first eigenvectors for MD-generated ensembles of free-NtMGAM, DSK-NtMGAM, and acarbose-NtMGAM, respectively.
FIGURE 5. (A) The NtMGAM domains with residue numbers labeled. (B) New cartoon diagram of NtMGAM with domains colored differently. (C) Inserted loop 1 in aligned compounds (the lowest energy conformations from FEL analysis). (D,E) NtMGAM active site pocket shown in surface representation occupied by DSK D and acarbose E. The structure of NtMGAM is shown in new cartoon and colored dissimilarly as follows: purple for catalytic domain, pink for inserted loop 1, cyan for inserted loop 2.
Subsequently, secondary structure analysis was also performed, and the corresponding average secondary structure values for each residue are shown in Figures 6A–C and Table 1. It can be seen that the proportion of α-helix in Asn498 in free-NtMGAM and acarbose-NtMGAM was about 90%, whereas in DSK-NtMGAM, it was about 59% (Table 1 and Supplementary Table S3). Supplementary Table S3 shows that the proportion of α-helix in residue 497–499 in free-NtMGAM and acarbose-NtMGAM are almost above 90%. In contrast, in the DSK-NtMGAM complex, the odds are reduced to under 70%. The results revealed that the α-helix of DSK-NtMGAM during 200 ns MD simulations disappeared partly. Residues 497–499 are located near the inserted loop 2 domain, which is quite close to the opening of the (β/α)8 barrel. Despiralization of these residues can enlarge the domain of the inserted loop 2, therefore, contributing to the architecture of the inhibitor binding site.
FIGURE 6. Conformation changes in residues Asn491-Leu493 of the three systems: (A) free-NtMGAM, (B) DSK-NtMGAM, and (C) acarbose-NtMGAM.
TABLE 1. The probability of secondary structures of residues H497 to L499.
POCASA 44 (http://altair.sci.hokudai.ac.jp/g6/service/pocasa/) (Yu et al., 2010) was utilized to predict the volume of the binding pocket. Parameters are listed as follows: the radius of probe sphere value was 1 Å, single point flag value was 10 Å, and protein depth flag value was 15Å. The active pocket conformations at 0, 100, and 200 ns were shown in Figures 7A–C. The volume of the pocket in free-NtMGAM was smaller than that of DSK-NtMGAM and acarbose-NtMGAM. Obviously, it could be considered that inhibitor of NtMGAM binding to the pocket with the nonreducing sugar ring in the −1 subsite and the reducing ring in the +1 subsite results in net retention of configuration at the anomeric center (Sim et al., 2010). Large active pockets will facilitate the inhibitor binding and entry. However, acarbose is so large that it binds to the NtMGAM active site primarily with its acarvosine unit (−1 and +1 subsite), few with its glycone rings (+2 and +3 subsite) (Sim et al., 2008).
FIGURE 7. The binding pocket conformations of free-NtMGAM (A), DSK-NtMGAM (B), and acarbose-NtMGAM (C) at 0, 100, and 200 ns with the volume of the pockets labeled.
In addition, the S group of DSK can generate charge interaction with residues (Trp400, Asp437, and Asp536) (Supplementary Figure S5), which may stabilize the DSK-NtMGAM complex.
PCA and FEL Analysis
PCA was performed to confirm whether the conformational changes of the three systems were continuous and stable, and the results are displayed through FEL (Figures 8A–C). Table 2 lists the PC1 and PC2 probabilities of the three systems. The structures of the most stable conformations of the three systems revealed that the conformational changes in the residues 497–499 domain were complex in DSK-NtMGAM (α-helix disappeared partly). In summary, we can confirm that the conformational changes in the three systems were continuous and stable, and our analyses above are reliable.
FIGURE 8. PCA based FEL analysis of NtMGAM (A). DSK-NtMGAM (B). Acarbose-NtMGAM (C) as a function of projections of the MD trajectory onto the first (PC1) and second (PC2) eigenvectors. The structures of the two most stable conformations of the three systems are presented with the conformation of Asn491 to Leu493.
TABLE 2. The probabilities of PC1 and PC2 of the three systems.
MM-PBSA Calculations
We used the end-point method to calculate the free energy between NtMGAM and the two inhibitors. The results of the Generalized Born (GB) implicit solvent method with a SASA term calculation are shown in Table 3. The binding free energies were mainly contributed by electrostatic energy, which is calculated by the molecular mechanics force field and the electrostatic contribution to the solvation free energy calculated by GB. Meanwhile, the VDW interactions are approximately consistent. As shown in Table 3, the binding free energy of the DSK-NtMGAM complex (−41.90 ± 1.14 kcal/mol) is lower than that of the acarbose-NtMGAM complex (−8.77 ± 1.08 kcal/mol). The results have the same trend as the data calculated from the Ki value from Sim’s works (Sim et al., 2008).
TABLE 3. MM-PBSA results.
ASMD Simulations
The 3D visualization of the channel of the free-NtMGAM, DSK-NtMGAM, and acarbose-NtMGAM obtained by CAVER Analyst 2.0 is listed in Figure 9A. The detailed exploration of the channel bottleneck and surrounding residues is shown in Figure 9B. The bottleneck of the channel in DSK-NtMGAM was larger than that of acarbose-NtMGAM. The surrounding residues displayed around the contour demonstrated the frequency. It can be seen that there are more residues in the DSK-NtMGAM complex comparing to that in the acarbose-NtMGAM complex.
FIGURE 9. (A) 3D visualization of candidate tunnels of (1) free-NtMGAM, (2) DSK-NtMGAM, and (3) acarbose-NtMGAM. (B) The details of the tunnel bottleneck contour over time of (1) free-NtMGAM, (2) DSK-NtMGAM, and (3) acarbose-NtMGAM.
To explore the enzyme-inhibitor interactions and the affinity of the active sites of NtMGAM via inhibitors unbinding pathway, ASMD simulations were performed on the two complexes. In Figure 10, the PMF profile displayed the energy changes during the process of pulling the ligands out of the NtMGAM channel. In the DSK-NtMGAM complex, a lower energy barrier (13 kcal/mol) should be transferred to completely dissociate DSK from the channel of NtMGAM than that in the acarbose-NtMGAM complex (22 kcal/mol).
FIGURE 10. PMF profiles along with the reaction coordinates of acarbose-NtMGAM (red) and DSK-NtMGAM (blue).
Figures 11A1–E1 show the interactional changes during the dissociation process of DSK-NtMGAM along the reaction coordinate (RC). First, for the initial coordinate (RC = 7 Å) [Figure 11A1], Asp542 and Asp443 formed a salt bridge with DSK. Then, at 8.44 Å [Figure 11B1], the free energy value dropped sharply due to the break of stronger hydrogen bonds between DSK and Asp203 and Asn449 residues. With the movement of DSK, the salt bridge between Asp443 and DSK disappeared and the salt bridge between Asp542 and DSK persisted. Nevertheless, at 9.58 Å [Figure 11C1], the residues that formed hydrogen bonds with the DSK were changed. Except for the salt bridge between Asp542 and DSK, the hydrogen bonds (salt bridge) were disappeared. Thereafter, the interconnections, including hydrogen bonds between channel residues and DSK, increased rapidly after 10.96 Å [Figure 11D1], giving rise to an increase in free energy value. Finally, at 12.72Å [Figure 11D1], DSK completely departed from the channel of NtMGAM and the curve of PMF tended to be flat. Asp542, as an important channel residue, could generate hydrogen bonds with DSK continuously, which was consistent with the results obtained in the channel analysis.
FIGURE 11. Interactions during ASMD simulations between DSK (1) or acarbose (2) with NtMGAM. (A) The initial state of the two systems before the ASMD simulations. Reaction coordinates reached (B) 8.44 Å, (C) 9.58 Å, (D) 10.96 Å, and (E) 12.72 Å.
The acarbose dissociating from NtMGAM is shown in Figures 11A2–E2. At the beginning of the ASMD simulation [Figure 11A2], acarbose was tightly fixed due to the strong hydrogen bond interactions with Tyr303, Asp327, Arg298, and His600. Subsequently, at 8.44 Å [Figure 11B2], the free energy value increased slowly because of the stronger hydrogen bond interactions between acarbose and Glu404 and Trp406. At 9.58 Å [Figure 11C2], the acarbose made hydrogen bond interactions with Asn449, Val405, and Phe450, which were stronger than DSK, resulting in the increased free energy value (Figure 11C). Thereafter, the interconnections, including the hydrogen bonds between acarbose and channel residues, disappeared at about 10.96 Å besides the only hydrogen bond with Ser448 [Figure 11D2]. Finally, at 12.72 Å [Figure 11E2], acarbose was completely departed from the unbinding pathways of NtMGAM, and the curves of PMF tended to be flat.
To sum up, compared to acarbose, DSK escaped from NtMGAM easily with lower energy. Asp542 is an important residue on the bottleneck of the active pocket of NtMGAM, which could generate hydrogen bonds with DSK continuously. Our results may provide some useful clues for designing new medicine to relieve symptoms of postprandial hyperglycemia caused by type 2 diabetes. For example, we can modify the 3D structure of acarbose to get a new compound that is suitable for the active pocket of NtMGAM.
Conclusion
At present, there are multiple drugs for the treatment of type 2 diabetes on the market, including α-glucosidase inhibitors. MGAM has become an efficient drug target for insulin resistance. In order to explore the conformational changes in the active pocket and unbinding pathway for NtMGAM, MD simulations and ASMD simulations were performed between two inhibitors (DSK and acarbose) and NtMGAM. MD simulations indicated that DSK binding to NtMGAM may lead to an enlargement of the active pocket due to the flexibility of the two domains (inserted loop 1 and inserted loop 2), which would facilitate the binding of DSK to NtMGAM. ASMD simulation results indicated that Asp542 was an important channel residue, which could continuously generate hydrogen bonds with DSK. Our results may provide some interesting thoughts for designing new medicine for the treatment of type 2 diabetes based on the molecular structure and specific intermolecular interactions between NtMGAM and DSK substrate in the binding pocket and the entrance channel.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.
Author Contributions
SZ wrote this paper and prepared the tables and figures. LH revised this paper. YW was responsible for the Supplementary Materials. XF and SW provided some revision advice. WL and WH provided the ideas and modified the papers.
Funding
This work was supported by the Overseas Cooperation Project of Jilin Province (20200801069GH).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.
Acknowledgments
Amber 16 package was purchased and authorized to use by Jilin University.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2021.711242/full#supplementary-material
Abbreviations
MGAM, maltase-glucoamylase; NtMGAM, N-terminal human maltase-glucoamylase; CtMGAM, C-terminal luminal subunit; MD, molecular dynamics; DSK, de-O-sulfonated kotalanol; PCA, Principal Component Analysis; PCs, principal components; FEL, free energy landscape; ASMD, adaptive steered molecular dynamics simulations; PMF, potential of mean force; JA, Jarzynski average; RC, reaction coordinate.
References
Al-Khafaji, K., and Taskin Tok, T. (2020). Molecular Dynamics Simulation, Free Energy Landscape and Binding Free Energy Computations in Exploration the Anti-invasive Activity of Amygdalin against Metastasis. Comput. Methods Programs Biomed. 195, 105660. doi:10.1016/j.cmpb.2020.105660
Andersen, H. C. (1980). Molecular Dynamics Simulations at Constant Pressure And/or Temperature. J. Chem. Phys. 72, 2384–2393. doi:10.1063/1.439486
Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F., Dinola, A., and Haak, J. R. (1984). Molecular Dynamics with Coupling to an External bath. J. Chem. Phys. 81, 3684–3690. doi:10.1063/1.448118
Bogunia, M., and Makowski, M. (2020). Influence of Ionic Strength on Hydrophobic Interactions in Water: Dependence on Solute Size and Shape. J. Phys. Chem. B 124 (46), 10326–10336. doi:10.1021/acs.jpcb.0c06399
Case, D. A., Betz, R. M., Cerutti, D. S., Cheatham, T. E., Darden, T. A., Duke, R. E., et al. (2016). Amber16. San Francisco, CA: University of California, San Francisco.
Chen, M., Hu, C., and Jia, W. (2015). Pharmacogenomics of Glinides. Pharmacogenomics 16 (1), 45–60. doi:10.2217/pgs.14.152
Chiasson, J.-L., Josse, R. G., Gomis, R., Hanefeld, M., Karasik, A., and Laakso, M. (2002). Acarbose for Prevention of Type 2 Diabetes Mellitus: the STOP-NIDDM Randomised Trial. The Lancet 359 (9323), 2072–2077. doi:10.1016/s0140-6736(02)08905-5
Dhameja, M., and Gupta, P. (2019). Synthetic Heterocyclic Candidates as Promising α-Glucosidase Inhibitors: An Overview. Eur. J. Med. Chem. 176, 343–377. doi:10.1016/j.ejmech.2019.04.025
Drucker, D. J. (2018). Mechanisms of Action and Therapeutic Application of Glucagon-like Peptide-1. Cel Metab. 27 (4), 740–756. doi:10.1016/j.cmet.2018.03.001
Elferink, H., Bruekers, J. P. J., Veeneman, G. H., and Boltje, T. J. (2020). A Comprehensive Overview of Substrate Specificity of Glycoside Hydrolases and Transporters in the Small Intestine. Cell. Mol. Life Sci. 77 (23), 4799–4826. doi:10.1007/s00018-020-03564-1
Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G. (1995). A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 103, 8577–8593. doi:10.1063/1.470117
Flores-Bocanegra, L., González-Andrade, M., Bye, R., Linares, E., and Mata, R. (2017). α-Glucosidase Inhibitors from Salvia Circinata. J. Nat. Prod. 80 (5), 1584–1593. doi:10.1021/acs.jnatprod.7b00155
Frandsen, T. P., and Svensson, B. (1998). Plant Alpha-Glucosidases of the Glycoside Hydrolase Family 31. Molecular Properties, Substrate Specificity, Reaction Mechanism, and Comparison with Family Members of Different Origin. Plant Mol. Biol. 37 (1), 1–13. doi:10.1023/a:1005925819741
Frauenfelder, H., Sligar, S., and Wolynes, P. (1991). The Energy Landscapes and Motions of Proteins. Science 254 (5038), 1598–1603. doi:10.1126/science.1749933
Grant, B. J., Rodrigues, A. P. C., ElSawy, K. M., McCammon, J. A., and Caves, L. S. D. (2006). Bio3d: An R Package for the Comparative Analysis of Protein Structures. Bioinformatics 22, 2695–2696. doi:10.1093/bioinformatics/btl461
Guàrdia, E., and Padró, J. A. (1985). Generalized Langevin Dynamics Simulation of Interacting Particles. J. Chem. Phys. 83, 1917–1920. doi:10.1063/1.449379
Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: Visual Molecular Dynamics. J. Mol. Graphics 14 (1), 33-38–27-38. doi:10.1016/0263-7855(96)00018-5
Izrailev, S., Stepaniants, S., Isralewitz, B., Kosztin, D., Lu, H., Molnár, F., et al. (1999). “Steered Molecular Dynamics,” in Computational Molecular Dynamics. Berlin: Springer.
Jarzynski, C. (1997). Equilibrium Free-Energy Differences from Nonequilibrium Measurements: A Master-Equation Approach. Phys. Rev. E 56, 5018–5035. doi:10.1103/physreve.56.5018
Jayakanthan, K., Mohan, S., and Pinto, B. M. (2009). Structure Proof and Synthesis of Kotalanol and De-O-sulfonated Kotalanol, Glycosidase Inhibitors Isolated from an Herbal Remedy for the Treatment of Type-2 Diabetes. J. Am. Chem. Soc. 131 (15), 5621–5626. doi:10.1021/ja900867q
Jurcik, A., Bednar, D., Byska, J., Marques, S. M., Furmanova, K., Daniel, L., et al. (2018). CAVER Analyst 2.0: Analysis and Visualization of Channels and Tunnels in Protein Structures and Molecular Dynamics Trajectories. Bioinformatics 34 (20), 3586–3588. doi:10.1093/bioinformatics/bty386
Kazufumi, N., Hiroki, O., Hajime, K., Kenei, S., Shota, F., Kyoko, W., et al. (2014). DPP-4 Inhibitor and Alpha-Glucosidase Inhibitor Equally Improve Endothelial Function in Patients with Type 2 Diabetes: EDGE Study. Cardiovasc. diabetology 13, 110. doi:10.1186/s12933-014-0110-2
Kelly, M. S., Lewis, J., Huntsberry, A. M., Dea, L., and Portillo, I. (2019). Efficacy and Renal Outcomes of SGLT2 Inhibitors in Patients with Type 2 Diabetes and Chronic Kidney Disease. Postgrad. Med. 131 (1), 31–42. doi:10.1080/00325481.2019.1549459
Lindorff-Larsen, K., Piana, S., Palmo, K., Maragakis, P., Klepeis, J. L., Dror, R. O., et al. (2010). Improved Side-Chain Torsion Potentials for the Amber ff99SB Protein Force Field. Proteins 78 (8), 1950–1958. doi:10.1002/prot.22711
Liu, Y., Zhu, J., Guo, X., Huang, T., Han, J., Gao, J., et al. (2020). How Oncogenic Mutations Activate Human MAP Kinase 1 (MEK1): A Molecular Dynamics Simulation Study. J. Biomol. Struct. Dyn. 38 (13), 3942–3958. doi:10.1080/07391102.2019.1686065
Lovering, A. L., Lee, S. S., Kim, Y.-W., Withers, S. G., and Strynadka, N. C. J. (2005). Mechanistic and Structural Analysis of a Family 31 α-Glycosidase and its Glycosyl-Enzyme Intermediate. J. Biol. Chem. 280 (3), 2105–2115. doi:10.1074/jbc.M410468200
Lyann, S., Kumarasamy, J., Sankar, M., Ravindranath, N., D, J. B., Mario, P. B., et al. (2010). New Glucosidase Inhibitors from an Ayurvedic Herbal Treatment for Type 2 Diabetes: Structures and Inhibition of Human Intestinal Maltase-Glucoamylase with Compounds from Salacia Reticulata. Biochemistry 49 (3), 443–451. doi:10.1021/bi9016457
Medagama, A. B. (2015). Salacia Reticulata (Kothala Himbutu) Revisited; a Missed Opportunity to Treat Diabetes and Obesity? Nutr. J. 14, 21. doi:10.1186/s12937-015-0013-4
Mi, C.-N., Yuan, J.-Z., Zhu, M.-M., Yang, L., Wei, Y.-M., Wang, H., et al. (2021). 2-(2-Phenylethyl)chromone Derivatives: Promising α-glucosidase Inhibitors in Agarwood from Aquilaria Filaria. Phytochemistry 181, 112578. doi:10.1016/j.phytochem.2020.112578
Miller, B. R., McGee, T. D., Swails, J. M., Homeyer, N., Gohlke, H., and Roitberg, A. E. (2012). MMPBSA.py: An Efficient Program for End-State Free Energy Calculations. J. Chem. Theor. Comput. 8 (9), 3314–3321. doi:10.1021/ct300418h
Miyamoto, S., and Kollman, P. (1992). Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 13, 952–962. doi:10.1002/jcc.540130805
Morris, G. M., Huey, R., Lindstrom, W., Sanner, M. F., Belew, R. K., Goodsell, D. S., et al. (2009). AutoDock4 and AutoDockTools4: Automated Docking with Selective Receptor Flexibility. J. Comput. Chem. 30 (16), 2785–2791. doi:10.1002/jcc.21256
Nanjan, M. J., Mohammed, M., Kumar, B. R. P., and Chandrasekar, M. J. N. (2018). Thiazolidinediones as Antidiabetic Agents: A Critical Review. Bioorg. Chem. 77, 548–567. doi:10.1016/j.bioorg.2018.02.009
Ozer, G., Keyes, T., Quirk, S., and Hernandez, R. (2014). Multiple Branched Adaptive Steered Molecular Dynamics. J. Chem. Phys. 141 (6), 064101. doi:10.1063/1.4891807
Ozer, G., Quirk, S., and Hernandez, R. (2012a). Adaptive Steered Molecular Dynamics: Validation of the Selection Criterion and Benchmarking Energetics in Vacuum. J. Chem. Phys. 136 (21), 215104. doi:10.1063/1.4725183
Ozer, G., Quirk, S., and Hernandez, R. (2012b). Thermodynamics of Decaalanine Stretching in Water Obtained by Adaptive Steered Molecular Dynamics Simulations. J. Chem. Theor. Comput. 8, 4837–4844. doi:10.1021/ct300709u
Ozer, G., Valeev, E. F., Quirk, S., and Hernandez, R. (2010). Adaptive Steered Molecular Dynamics of the Long-Distance Unfolding of Neuropeptide Y. J. Chem. Theor. Comput. 6, 3026–3038. doi:10.1021/ct100320g
Patel, S. (2015). Cerebrovascular Complications of Diabetes: Alpha Glucosidase Inhibitor as Potential Therapy. Horm. Metab. Res. 48 (2), 83–91. doi:10.1055/s-0035-1565181
Popović, D. M., and Stuchebrukhov, A. A. (2004). Electrostatic Study of the Proton Pumping Mechanism in Bovine Heart Cytochrome C Oxidase. J. Am. Chem. Soc. 126 (6), 1858–1871. doi:10.1021/ja038267w
Quezada-Calvillo, R., Sim, L., Ao, Z., Hamaker, B. R., Quaroni, A., Brayer, G. D., et al. (2008). Luminal Starch Substrate "brake" on Maltase-Glucoamylase Activity Is Located within the Glucoamylase Subunit. J. Nutr. 138 (4), 685–692. doi:10.1093/jn/138.4.685
Ren, L., Qin, X., Cao, X., Wang, L., Bai, F., Bai, G., et al. (2011). Structural Insight into Substrate Specificity of Human Intestinal Maltase-Glucoamylase. Protein & Cell 2 (10), 827–836. doi:10.1007/s13238-011-1105-3
Ríos, J., Francini, F., and Schinella, G. (2015). Natural Products for the Treatment of Type 2 Diabetes Mellitus. Planta Med. 81 (12-13), 975–994. doi:10.1055/s-0035-1546131
Rosenberg, R., Boughaleb, Y., Nitzan, A., and Ratner, M. (1986). Effective Potentials from Langevin Dynamic Simulations of Framework Solid Electrolytes. Solid State Ionics 18-19, 127–135. doi:10.1016/0167-2738(86)90099-8
Santos, C. M. M., Freitas, M., and Fernandes, E. (2018). A Comprehensive Review on Xanthone Derivatives as α-glucosidase Inhibitors. Eur. J. Med. Chem. 157, 1460–1479. doi:10.1016/j.ejmech.2018.07.073
Satoh, T., Toshimori, T., Yan, G., Yamaguchi, T., and Kato, K. (2016). Structural Basis for Two-step Glucose Trimming by Glucosidase II Involved in ER Glycoprotein Quality Control. Scientific Rep. 6, 20575. doi:10.1038/srep20575
Satoru, S., Hisakazu, N., Kitaro, K., and Hajime, H. (2015). Review: Miglitol Has Potential as a Therapeutic Drug against Obesity. Nutr. Metab. 12, 51. doi:10.1186/s12986-015-0048-8
Schäfer, G. (1983). Biguanides. A Review of History, Pharmacodynamics and Therapy. Diabete & metabolisme 9 (2), 148–163.
Sim, L., Quezada-Calvillo, R., Sterchi, E. E., Nichols, B. L., and Rose, D. R. (2008). Human Intestinal Maltase-Glucoamylase: Crystal Structure of the N-Terminal Catalytic Subunit and Basis of Inhibition and Substrate Specificity. J. Mol. Biol. 375 (3), 782–792. doi:10.1016/j.jmb.2007.10.069
Sim, L., Willemsma, C., Mohan, S., Naim, H. Y., Pinto, B. M., and Rose, D. R. (2010). Structural Basis for Substrate Selectivity in Human Maltase-Glucoamylase and Sucrase-Isomaltase N-Terminal Domains. J. Biol. Chem. 285 (23), 17763–17770. doi:10.1074/jbc.M109.078980
Singh, B., Bulusu, G., and Mitra, A. (2015). Understanding the Thermostability and Activity of Bacillus subtilisLipase Mutants: Insights from Molecular Dynamics Simulations. J. Phys. Chem. B 119 (2), 392–409. doi:10.1021/jp5079554
Stephen, C., David, M., Leiter, A. L., Siew, P. C., Giorgio, S., and Michel, M. (2018). The Place of Gliclazide MR in the Evolving Type 2 Diabetes Landscape: A Comparison with Other Sulfonylureas and Newer Oral Antihyperglycemic Agents. Diabetes Res. Clin. Pract. 143, 1–14. doi:10.1016/j.diabres.2018.05.028
Sugitani, R., and Stuchebrukhov, A. A. (2009). Molecular Dynamics Simulation of Water in Cytochrome C Oxidase Reveals Two Water Exit Pathways and the Mechanism of Transport. Biochim. Biophys. Acta (Bba) - Bioenerg. 1787 (9), 1140–1150. doi:10.1016/j.bbabio.2009.04.004
Tashiro, M., and Stuchebrukhov, A. A. (2005). Thermodynamic Properties of Internal Water Molecules in the Hydrophobic Cavity Around the Catalytic Center of CytochromecOxidase. J. Phys. Chem. B. 109 (2), 1015–1022. doi:10.1021/jp0462456
Thornberry, N. A., and Gallwitz, B. (2009). Mechanism of Action of Inhibitors of Dipeptidyl-Peptidase-4 (DPP-4). Best Pract. Res. Clin. Endocrinol. Metab. 23 (4), 479–486. doi:10.1016/j.beem.2009.03.004
Tuyen, D. T., Yew, G. Y., Cuong, N. T., Hoang, L. T., Yen, H. T., Hong Thao, P. T., et al. (2021). Selection, Purification, and Evaluation of Acarbose−an α-glucosidase Inhibitor from Actinoplanes Sp. Chemosphere 265, 129167. doi:10.1016/j.chemosphere.2020.129167
Usman, B., Sharma, N., Satija, S., Mehta, M., Vyas, M., Khatik, G. L., et al. (2019). Recent Developments in Alpha-Glucosidase Inhibitors for Management of Type-2 Diabetes: An Update. Curr. Pharm. Des. 25 (23), 2510–2525. doi:10.2174/1381612825666190717104547
Vahedi-Faridi, A., Licht, A., Bulut, H., Scheffel, F., Keller, S., Wehmeier, U. F., et al. (2010). Crystal Structures of the Solute Receptor GacH of Streptomyces Glaucescens in Complex with Acarbose and an Acarbose Homolog: Comparison with the Acarbose-Loaded Maltose-Binding Protein of Salmonella typhimurium. J. Mol. Biol. 397 (3), 709–723. doi:10.1016/j.jmb.2010.01.054
Van Beers, E. H., Büller, H. A., Grand, R. J., Einerhand, A. W. C., and Dekker, J. (1995). Intestinal brush Border Glycohydrolases: Structure, Function, and Development. Crit. Rev. Biochem. Mol. Biol. 30 (3), 197–262. doi:10.3109/10409239509085143
Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). Development and Testing of a General Amber Force Field. J. Comput. Chem. 25 (9), 1157–1174. doi:10.1002/jcc.20035
Yu, J., Zhou, Y., Tanaka, I., and Yao, M. (2010). Roll: a New Algorithm for the Detection of Protein Pockets and Cavities with a Rolling Probe Sphere. Bioinformatics 26 (1), 46–52. doi:10.1093/bioinformatics/btp599
Yu, S., Bojsen, K., Svensson, B., and Marcussen, J. (1999). alpha-1,4-glucan Lyases Producing 1,5-Anhydro-D-Fructose from Starch and Glycogen Have Sequence Similarity to Alpha-Glucosidases. Biochim. Biophys. Acta 1433 (1-2), 1–15. doi:10.1016/s0167-4838(99)00152-1
Zhu, J., Li, Y., Wang, J., Yu, Z., Liu, Y., Tong, Y., et al. (2018). Adaptive Steered Molecular Dynamics Combined with Protein Structure Networks Revealing the Mechanism of Y68I/G109P Mutations that Enhance the Catalytic Activity of D-Psicose 3-Epimerase from Clostridium Bolteae. Front. Chem. 6, 437. doi:10.3389/fchem.2018.00437
Keywords: maltase-glucoamylase, inhibitors, molecular dynamics simulations, adaptive steered molecular dynamics simulations, conformational change
Citation: Zhang S, Wang Y, Han L, Fu X, Wang S, Li W and Han W (2021) Targeting N-Terminal Human Maltase-Glucoamylase to Unravel Possible Inhibitors Using Molecular Docking, Molecular Dynamics Simulations, and Adaptive Steered Molecular Dynamics Simulations. Front. Chem. 9:711242. doi: 10.3389/fchem.2021.711242
Received: 18 May 2021; Accepted: 26 July 2021;
Published: 30 August 2021.
Edited by:
Sonja Grubisic, University of Belgrade, SerbiaReviewed by:
Dragan M. Popovic, University of Belgrade, SerbiaFrancesca Mocci, University of Cagliari, Italy
Copyright © 2021 Zhang, Wang, Han, Fu, Wang, Li and Han. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Wannan Li, bGl3YW5uYW5Aamx1LmVkdS5jbg==; Weiwei Han, d2Vpd2VpaGFuQGpsdS5lZHUuY24=