ORIGINAL RESEARCH article

Front. Chem., 30 August 2021

Sec. Theoretical and Computational Chemistry

Volume 9 - 2021 | https://doi.org/10.3389/fchem.2021.711242

Targeting N-Terminal Human Maltase-Glucoamylase to Unravel Possible Inhibitors Using Molecular Docking, Molecular Dynamics Simulations, and Adaptive Steered Molecular Dynamics Simulations

  • 1. School of Life Science, Key Laboratory for Molecular Enzymology and Engineering of Ministry of Education, National Engineering Laboratory for AIDS Vaccine, Jilin University, Changchun, China

  • 2. Laboratory of Theoretical and Computational Chemistry, Institute of Theoretical Chemistry, Jilin University, Changchun, China

Abstract

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

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

FIGURE 3

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

FIGURE 5

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

TABLE 1

ResidueFree-NtMGAMDSK-NtMGAMAcarbose-NtMGAM
α-HelixLoopα-HelixLoopα-HelixLoop
H4970.940.060.590.410.890.11
N4980.940.060.590.410.920.08
L4990.940.050.60.400.920.08

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

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

TABLE 2

PC 1 (%)PC 2 (%)
Free-NtMGAM16.466.28
DSK-NtMGAM17.036.44
Acarbose-NtMGAM19.615.46

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

DSKAcarbose
VDWAALS−15.84 ± 0.39−29.86 ± 1.01
EEL−586.14 ± 2.32−81.12 ± 4.95
EGB565.10 ± 1.55108.44 ± 4.87
ESURF−5.02 ± 0.06−6.22 ± 0.19
ΔG gas−601.98 ± 2.49−110.99 ± 4.99
ΔG solv560.08 ± 1.50102.21 ± 4.74
ΔTOTAL−41.90 ± 1.14−8.77 ± 1.08
Ki (μM)0.03 ± 0.0162 ± 13
ΔGexp−10.32−5.76

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

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

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

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.

Statements

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).

Acknowledgments

Amber 16 package was purchased and authorized to use by Jilin University.

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.

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

  • 1

    Al-KhafajiK.Taskin TokT. (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. 10.1016/j.cmpb.2020.105660

  • 2

    AndersenH. C. (1980). Molecular Dynamics Simulations at Constant Pressure And/or Temperature. J. Chem. Phys.72, 23842393. 10.1063/1.439486

  • 3

    BerendsenH. J. C.PostmaJ. P. M.van GunsterenW. F.DinolaA.HaakJ. R. (1984). Molecular Dynamics with Coupling to an External bath. J. Chem. Phys.81, 36843690. 10.1063/1.448118

  • 4

    BoguniaM.MakowskiM. (2020). Influence of Ionic Strength on Hydrophobic Interactions in Water: Dependence on Solute Size and Shape. J. Phys. Chem. B124 (46), 1032610336. 10.1021/acs.jpcb.0c06399

  • 5

    CaseD. A.BetzR. M.CeruttiD. S.CheathamT. E.IIIDardenT. A.DukeR. E.et al (2016). Amber16. San Francisco, CA: University of California, San Francisco.

  • 6

    ChenM.HuC.JiaW. (2015). Pharmacogenomics of Glinides. Pharmacogenomics16 (1), 4560. 10.2217/pgs.14.152

  • 7

    ChiassonJ.-L.JosseR. G.GomisR.HanefeldM.KarasikA.LaaksoM. (2002). Acarbose for Prevention of Type 2 Diabetes Mellitus: the STOP-NIDDM Randomised Trial. The Lancet359 (9323), 20722077. 10.1016/s0140-6736(02)08905-5

  • 8

    DhamejaM.GuptaP. (2019). Synthetic Heterocyclic Candidates as Promising α-Glucosidase Inhibitors: An Overview. Eur. J. Med. Chem.176, 343377. 10.1016/j.ejmech.2019.04.025

  • 9

    DruckerD. J. (2018). Mechanisms of Action and Therapeutic Application of Glucagon-like Peptide-1. Cel Metab.27 (4), 740756. 10.1016/j.cmet.2018.03.001

  • 10

    ElferinkH.BruekersJ. P. J.VeenemanG. H.BoltjeT. J. (2020). A Comprehensive Overview of Substrate Specificity of Glycoside Hydrolases and Transporters in the Small Intestine. Cell. Mol. Life Sci.77 (23), 47994826. 10.1007/s00018-020-03564-1

  • 11

    EssmannU.PereraL.BerkowitzM. L.DardenT.LeeH.PedersenL. G. (1995). A Smooth Particle Mesh Ewald Method. J. Chem. Phys.103, 85778593. 10.1063/1.470117

  • 12

    Flores-BocanegraL.González-AndradeM.ByeR.LinaresE.MataR. (2017). α-Glucosidase Inhibitors from Salvia Circinata. J. Nat. Prod.80 (5), 15841593. 10.1021/acs.jnatprod.7b00155

  • 13

    FrandsenT. P.SvenssonB. (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), 113. 10.1023/a:1005925819741

  • 14

    FrauenfelderH.SligarS.WolynesP. (1991). The Energy Landscapes and Motions of Proteins. Science254 (5038), 15981603. 10.1126/science.1749933

  • 15

    GrantB. J.RodriguesA. P. C.ElSawyK. M.McCammonJ. A.CavesL. S. D. (2006). Bio3d: An R Package for the Comparative Analysis of Protein Structures. Bioinformatics22, 26952696. 10.1093/bioinformatics/btl461

  • 16

    GuàrdiaE.PadróJ. A. (1985). Generalized Langevin Dynamics Simulation of Interacting Particles. J. Chem. Phys.83, 19171920. 10.1063/1.449379

  • 17

    HumphreyW.DalkeA.SchultenK. (1996). VMD: Visual Molecular Dynamics. J. Mol. Graphics14 (1), 33-3827-38. 10.1016/0263-7855(96)00018-5

  • 18

    IzrailevS.StepaniantsS.IsralewitzB.KosztinD.LuH.MolnárF.et al (1999). “Steered Molecular Dynamics,” in Computational Molecular Dynamics. Berlin: Springer.

  • 19

    JarzynskiC. (1997). Equilibrium Free-Energy Differences from Nonequilibrium Measurements: A Master-Equation Approach. Phys. Rev. E56, 50185035. 10.1103/physreve.56.5018

  • 20

    JayakanthanK.MohanS.PintoB. 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), 56215626. 10.1021/ja900867q

  • 21

    JurcikA.BednarD.ByskaJ.MarquesS. M.FurmanovaK.DanielL.et al (2018). CAVER Analyst 2.0: Analysis and Visualization of Channels and Tunnels in Protein Structures and Molecular Dynamics Trajectories. Bioinformatics34 (20), 35863588. 10.1093/bioinformatics/bty386

  • 22

    KazufumiN.HirokiO.HajimeK.KeneiS.ShotaF.KyokoW.et al (2014). DPP-4 Inhibitor and Alpha-Glucosidase Inhibitor Equally Improve Endothelial Function in Patients with Type 2 Diabetes: EDGE Study. Cardiovasc. diabetology13, 110. 10.1186/s12933-014-0110-2

  • 23

    KellyM. S.LewisJ.HuntsberryA. M.DeaL.PortilloI. (2019). Efficacy and Renal Outcomes of SGLT2 Inhibitors in Patients with Type 2 Diabetes and Chronic Kidney Disease. Postgrad. Med.131 (1), 3142. 10.1080/00325481.2019.1549459

  • 24

    Lindorff-LarsenK.PianaS.PalmoK.MaragakisP.KlepeisJ. L.DrorR. O.et al (2010). Improved Side-Chain Torsion Potentials for the Amber ff99SB Protein Force Field. Proteins78 (8), 19501958. 10.1002/prot.22711

  • 25

    LiuY.ZhuJ.GuoX.HuangT.HanJ.GaoJ.et al (2020). How Oncogenic Mutations Activate Human MAP Kinase 1 (MEK1): A Molecular Dynamics Simulation Study. J. Biomol. Struct. Dyn.38 (13), 39423958. 10.1080/07391102.2019.1686065

  • 26

    LoveringA. L.LeeS. S.KimY.-W.WithersS. G.StrynadkaN. C. J. (2005). Mechanistic and Structural Analysis of a Family 31 α-Glycosidase and its Glycosyl-Enzyme Intermediate. J. Biol. Chem.280 (3), 21052115. 10.1074/jbc.M410468200

  • 27

    LyannS.KumarasamyJ.SankarM.RavindranathN.DJ. B.MarioP. 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. Biochemistry49 (3), 443451. 10.1021/bi9016457

  • 28

    MedagamaA. B. (2015). Salacia Reticulata (Kothala Himbutu) Revisited; a Missed Opportunity to Treat Diabetes and Obesity?Nutr. J.14, 21. 10.1186/s12937-015-0013-4

  • 29

    MiC.-N.YuanJ.-Z.ZhuM.-M.YangL.WeiY.-M.WangH.et al (2021). 2-(2-Phenylethyl)chromone Derivatives: Promising α-glucosidase Inhibitors in Agarwood from Aquilaria Filaria. Phytochemistry181, 112578. 10.1016/j.phytochem.2020.112578

  • 30

    MillerB. R.3rdMcGeeT. D.Jr.SwailsJ. M.HomeyerN.GohlkeH.RoitbergA. E. (2012). MMPBSA.py: An Efficient Program for End-State Free Energy Calculations. J. Chem. Theor. Comput.8 (9), 33143321. 10.1021/ct300418h

  • 31

    MiyamotoS.KollmanP. (1992). Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem.13, 952962. 10.1002/jcc.540130805

  • 32

    MorrisG. M.HueyR.LindstromW.SannerM. F.BelewR. K.GoodsellD. S.et al (2009). AutoDock4 and AutoDockTools4: Automated Docking with Selective Receptor Flexibility. J. Comput. Chem.30 (16), 27852791. 10.1002/jcc.21256

  • 33

    NanjanM. J.MohammedM.KumarB. R. P.ChandrasekarM. J. N. (2018). Thiazolidinediones as Antidiabetic Agents: A Critical Review. Bioorg. Chem.77, 548567. 10.1016/j.bioorg.2018.02.009

  • 34

    OzerG.KeyesT.QuirkS.HernandezR. (2014). Multiple Branched Adaptive Steered Molecular Dynamics. J. Chem. Phys.141 (6), 064101. 10.1063/1.4891807

  • 35

    OzerG.QuirkS.HernandezR. (2012a). Adaptive Steered Molecular Dynamics: Validation of the Selection Criterion and Benchmarking Energetics in Vacuum. J. Chem. Phys.136 (21), 215104. 10.1063/1.4725183

  • 36

    OzerG.QuirkS.HernandezR. (2012b). Thermodynamics of Decaalanine Stretching in Water Obtained by Adaptive Steered Molecular Dynamics Simulations. J. Chem. Theor. Comput.8, 48374844. 10.1021/ct300709u

  • 37

    OzerG.ValeevE. F.QuirkS.HernandezR. (2010). Adaptive Steered Molecular Dynamics of the Long-Distance Unfolding of Neuropeptide Y. J. Chem. Theor. Comput.6, 30263038. 10.1021/ct100320g

  • 38

    PatelS. (2015). Cerebrovascular Complications of Diabetes: Alpha Glucosidase Inhibitor as Potential Therapy. Horm. Metab. Res.48 (2), 8391. 10.1055/s-0035-1565181

  • 39

    PopovićD. M.StuchebrukhovA. A. (2004). Electrostatic Study of the Proton Pumping Mechanism in Bovine Heart Cytochrome C Oxidase. J. Am. Chem. Soc.126 (6), 18581871. 10.1021/ja038267w

  • 40

    Quezada-CalvilloR.SimL.AoZ.HamakerB. R.QuaroniA.BrayerG. D.et al (2008). Luminal Starch Substrate "brake" on Maltase-Glucoamylase Activity Is Located within the Glucoamylase Subunit. J. Nutr.138 (4), 685692. 10.1093/jn/138.4.685

  • 41

    RenL.QinX.CaoX.WangL.BaiF.BaiG.et al (2011). Structural Insight into Substrate Specificity of Human Intestinal Maltase-Glucoamylase. Protein & Cell2 (10), 827836. 10.1007/s13238-011-1105-3

  • 42

    RíosJ.FranciniF.SchinellaG. (2015). Natural Products for the Treatment of Type 2 Diabetes Mellitus. Planta Med.81 (12-13), 975994. 10.1055/s-0035-1546131

  • 43

    RosenbergR.BoughalebY.NitzanA.RatnerM. (1986). Effective Potentials from Langevin Dynamic Simulations of Framework Solid Electrolytes. Solid State Ionics18-19, 127135. 10.1016/0167-2738(86)90099-8

  • 44

    SantosC. M. M.FreitasM.FernandesE. (2018). A Comprehensive Review on Xanthone Derivatives as α-glucosidase Inhibitors. Eur. J. Med. Chem.157, 14601479. 10.1016/j.ejmech.2018.07.073

  • 45

    SatohT.ToshimoriT.YanG.YamaguchiT.KatoK. (2016). Structural Basis for Two-step Glucose Trimming by Glucosidase II Involved in ER Glycoprotein Quality Control. Scientific Rep.6, 20575. 10.1038/srep20575

  • 46

    SatoruS.HisakazuN.KitaroK.HajimeH. (2015). Review: Miglitol Has Potential as a Therapeutic Drug against Obesity. Nutr. Metab.12, 51. 10.1186/s12986-015-0048-8

  • 47

    SchäferG. (1983). Biguanides. A Review of History, Pharmacodynamics and Therapy. Diabete & metabolisme9 (2), 148163.

  • 48

    SimL.Quezada-CalvilloR.SterchiE. E.NicholsB. L.RoseD. 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), 782792. 10.1016/j.jmb.2007.10.069

  • 49

    SimL.WillemsmaC.MohanS.NaimH. Y.PintoB. M.RoseD. R. (2010). Structural Basis for Substrate Selectivity in Human Maltase-Glucoamylase and Sucrase-Isomaltase N-Terminal Domains. J. Biol. Chem.285 (23), 1776317770. 10.1074/jbc.M109.078980

  • 50

    SinghB.BulusuG.MitraA. (2015). Understanding the Thermostability and Activity of Bacillus subtilisLipase Mutants: Insights from Molecular Dynamics Simulations. J. Phys. Chem. B119 (2), 392409. 10.1021/jp5079554

  • 51

    StephenC.DavidM.LeiterA. L.SiewP. C.GiorgioS.MichelM. (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, 114. 10.1016/j.diabres.2018.05.028

  • 52

    SugitaniR.StuchebrukhovA. 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), 11401150. 10.1016/j.bbabio.2009.04.004

  • 53

    TashiroM.StuchebrukhovA. A. (2005). Thermodynamic Properties of Internal Water Molecules in the Hydrophobic Cavity Around the Catalytic Center of CytochromecOxidase. J. Phys. Chem. B.109 (2), 10151022. 10.1021/jp0462456

  • 54

    ThornberryN. A.GallwitzB. (2009). Mechanism of Action of Inhibitors of Dipeptidyl-Peptidase-4 (DPP-4). Best Pract. Res. Clin. Endocrinol. Metab.23 (4), 479486. 10.1016/j.beem.2009.03.004

  • 55

    TuyenD. T.YewG. Y.CuongN. T.HoangL. T.YenH. T.Hong ThaoP. T.et al (2021). Selection, Purification, and Evaluation of Acarbose−an α-glucosidase Inhibitor from Actinoplanes Sp. Chemosphere265, 129167. 10.1016/j.chemosphere.2020.129167

  • 56

    UsmanB.SharmaN.SatijaS.MehtaM.VyasM.KhatikG. L.et al (2019). Recent Developments in Alpha-Glucosidase Inhibitors for Management of Type-2 Diabetes: An Update. Curr. Pharm. Des.25 (23), 25102525. 10.2174/1381612825666190717104547

  • 57

    Vahedi-FaridiA.LichtA.BulutH.ScheffelF.KellerS.WehmeierU. 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), 709723. 10.1016/j.jmb.2010.01.054

  • 58

    Van BeersE. H.BüllerH. A.GrandR. J.EinerhandA. W. C.DekkerJ. (1995). Intestinal brush Border Glycohydrolases: Structure, Function, and Development. Crit. Rev. Biochem. Mol. Biol.30 (3), 197262. 10.3109/10409239509085143

  • 59

    WangJ.WolfR. M.CaldwellJ. W.KollmanP. A.CaseD. A. (2004). Development and Testing of a General Amber Force Field. J. Comput. Chem.25 (9), 11571174. 10.1002/jcc.20035

  • 60

    YuJ.ZhouY.TanakaI.YaoM. (2010). Roll: a New Algorithm for the Detection of Protein Pockets and Cavities with a Rolling Probe Sphere. Bioinformatics26 (1), 4652. 10.1093/bioinformatics/btp599

  • 61

    YuS.BojsenK.SvenssonB.MarcussenJ. (1999). alpha-1,4-glucan Lyases Producing 1,5-Anhydro-D-Fructose from Starch and Glycogen Have Sequence Similarity to Alpha-Glucosidases. Biochim. Biophys. Acta1433 (1-2), 115. 10.1016/s0167-4838(99)00152-1

  • 62

    ZhuJ.LiY.WangJ.YuZ.LiuY.TongY.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. 10.3389/fchem.2018.00437

  • 63

    ZhuJ.QiR.LiuY.ZhaoL.HanW. (2019). Mechanistic Insights into the Effect of Ligands on Structural Stability and Selectivity of Sulfotransferase 2A1 (SULT2A1). ACS Omega4 (26), 2202122034. 10.1021/acsomega.9b03136

Summary

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

Volume

9 - 2021

Edited by

Sonja Grubisic, University of Belgrade, Serbia

Reviewed by

Dragan M. Popovic, University of Belgrade, Serbia

Francesca Mocci, University of Cagliari, Italy

Updates

Copyright

*Correspondence: Wannan Li, ; Weiwei Han,

This article was submitted to Theoretical and Computational Chemistry, a section of the journal Frontiers in Chemistry

Disclaimer

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.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics