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.

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 (R_g), root-mean-square fluctuation (RMSF), and solvent-accessible surface area (SASA) values were calculated using VMD (Humphrey et al., 1996).

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: Δ G b i n d i n g , s o l v a t e d = Δ G c o m p l e x , s o l v a t e d − [ Δ G r e c e p t o r , s o l v a t e d + Δ G l i g a n d , s o l v a t e d ] . Then, the free energy change associated with each term of ΔG was calculated according to the following: Δ G s o l v a t e d = E g a s + Δ G s o l v a t i o n − T S s o l u t e , where ΔG _solvation 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: Δ G s o l v a t e d ≅ < E g a s > + < Δ G s o l v a t i o n > − T < S s o l u t e > = 1 N ∑ i = 1 N E i , g a s + 1 N ∑ i = 1 N Δ G i , s o l v a t i o n − T N ∑ i = 1 N S i , s o l u t e , where i is the index of a particular frame and N is the total number of frames analyzed. The gas-phase energies (E_gas) can be computed from the quantum mechanical (QM) calculations and used further as a part of the force field parameterization; therefore, the E_gas energies can be abstracted from the molecular mechanical (MM) energies and the corresponding force field (Miller et al., 2012).

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

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.

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

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 (β/α)₈ 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.

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 R_g 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 Å²; meanwhile, the other two systems were stabilized at 33,000 Å². 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.

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 (β/α)₈ 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.

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 (β/α)₈ barrel. Despiralization of these residues can enlarge the domain of the inserted loop 2, therefore, contributing to the architecture of the inhibitor binding site.

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

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

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 K_i value from Sim’s works (Sim et al., 2008).

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.

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

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.

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.