Computational Studies on the Potency and Selectivity of PUGNAc Derivatives Against GH3, GH20, and GH84 β-N-acetyl-D-hexosaminidases

β-N-acetyl-D-hexosaminidases have attracted significant attention due to their crucial role in diverse physiological functions including antibacterial synergists, pathogen defense, virus infection, lysosomal storage, and protein glycosylation. In particular, the GH3 β-N-acetyl-D-hexosaminidase of V. cholerae (VcNagZ), human GH20 β-N-acetyl-D-hexosaminidase B (HsHexB), and human GH84 β-N-acetyl-D-hexosaminidase (hOGA) are three important representative glycosidases. These have been found to be implicated in β-lactam resistance (VcNagZ), lysosomal storage disorders (HsHexB) and Alzheimer's disease (hOGA). Considering the profound effects of these three enzymes, many small molecule inhibitors with good potency and selectivity have been reported to regulate the corresponding physiological functions. In this paper, the best-known inhibitors PUGNAc and two of its derivatives (N-valeryl-PUGNAc and EtBuPUG) were selected as model compounds and docked into the active pockets of VcNagZ, HsHexB, and hOGA, respectively. Subsequently, molecular dynamics simulations of the nine systems were performed to systematically compare their binding modes from active pocket architecture and individual interactions. Furthermore, the binding free energy and free energy decomposition are calculated using the MM/GBSA methods to predict the binding affinities of enzyme-inhibitor systems and to quantitatively analyze the contribution of each residue. The results show that PUGNAc is deeply-buried in the active pockets of all three enzymes, which indicates its potency (but not selectivity) against VcNagZ, HsHexB, and hOGA. However, EtBuPUG, bearing branched 2-isobutamido, adopted strained conformations and was only located in the active pocket of VcNagZ. It has completely moved out of the pocket of HsHexB and lacks interactions with HsHexB. This indicates why the selectivity of EtBuPUG to VcNagZ/HsHexB is the largest, reaching 968-fold. In addition, the contributions of the catalytic residue Asp253 (VcNagZ), Asp254 (VcNagZ), Asp175 (hOGA), and Asp354 (HsHexB) are important to distinguish the activity and selectivity of these inhibitors. The results of this study provide a helpful structural guideline to promote the development of novel and selective inhibitors against specific β-N-acetyl-D-hexosaminidases.

These studies on the crystal structures of GH3, GH20, and GH84 β-N-acetyl-D-hexosaminidases in complexes with PUGNAc and its derivatives partially clarified the binding mechanisms of related inhibitors with different GH β-N-acetyl-D-hexosaminidases. Nevertheless, the crystal structure of the complexes of GH20 and GH84 β-N-acetyl-D-hexosaminidases bound to N-valeryl-PUGNAc and EtBuPUG has not been reported to date. Furthermore, few in-depth studies have been published about the dynamic changes of GH3, GH20, and GH84  a Results determined previously . b Results determined previously (Macauley et al., 2005). c Results determined previously (Stubbs et al., 2006). β-N-acetyl-D-hexosaminidases bound to PUGNAc, N-valeryl-PUGNAc, and EtBuPUG. Such dynamic mechanisms could be closer to the models of inhibitor binding with target enzymes in vivo, thus it may provide interesting to investigate the details of the protein-ligand interaction (Anderson et al., 2005;Zhang et al., 2005;Falck et al., 2008;Vergara-Jaque et al., 2012;Soumana et al., 2016). Given the profound physiological functions of GH3 VcNagZ, GH20 HsHexB, and GH84 hOGA on the regulation of biological processes, in this paper, PUGNAc, N-valeryl-PUGNAc, and EtBuPUG (which possess different inhibitory potency against these three enzymes) were selected as model compounds. Then, possible inhibitory modes of selected inhibitors against the specific β-N-acetyl-D-hexosaminidases were investigated using molecular docking, molecular dynamics (MD) simulations. Finally, comparing these established binding patterns, clarifies the reasons for both the potency and selectivity of PUGNAc and its derivatives against GH3, GH20, and GH84 β-N-acetyl-D-hexosaminidases, which could further guide the discovery of highly effective and specific inhibitors.

Molecule Preparation and Optimization
PUGNAc, N-valeryl-PUGNAc, and EtBuPUG molecules were constructed using the Sketch Molecule module of SYBYL 7.3 software package (Horsch et al., 1991;Stubbs et al., 2006Stubbs et al., , 2007Tripos Associates, 2006). To achieve corresponding low energy conformations, these compounds were then fully optimized at the B3LYP/6-31G(d) basis set, using density functional theory (DFT) in Gaussian 16 (Frisch et al., 2016). The optimized conformation was used for the subsequent docking study.

Molecular Docking
Molecular docking can effectively predict and characterize predominant binding sites of a ligand located in the protein.
In this study, the X-ray crystal complex structures of GH3 VcNagZ (PDB code: 5UTQ) (Vadlamani et al., 2017), GH20 HsHexB (PDB code: 1NP0) (Mark et al., 2003), and GH84 hOGA (PDB code: 5M7T) (Roth et al., 2017) were obtained from the Protein Data Bank (http://www.rcsb.org/pdb/home/ home.do). These were selected as the starting structure for docking, using the Surflex-Dock program of the SYBYL 7.3 software package (Tripos Associates, 2006). Surflex-Dock uses an empirically derived scoring function that is based on the binding affinities of protein-ligand complexes and on their respective X-ray structures. Several preliminary steps were required before docking, such as removing water molecules and adding hydrogens. In the following, the ligand mode was applied to generate the protomol, a suitable putative ligand pose, which was then used for the docking of the small molecules PUGNAc, N-valeryl-PUGNAc, and EtBuPUG. During the docking studies, ring flexibility was considered, and default settings were chosen for the remaining parameters. The binding affinities of the protein-ligand complexes were obtained with the Surflex scoring, which considers hydrophobic, polar, repulsive, entropic entropic, and solvated effects (Welch et al., 1996).

MD Simulations
The enzymes of VcNagZ, HsHexB, and hOGA and three compounds of PUGNAc, N-valeryl-PUGNAc, and EtBuPUG, constituted nine systems, which were used for the MD calculations in the Amber14 program (Case et al., 2014). For each of the ligand-protein system simulations, the force field parameters for proteins and inhibitors were represented by AMBER03 (Duan et al., 2003) and GAFF (Wang et al., 2004) force field, respectively. Then, each complex structure was solvated in a truncated octahedral box in which TIP3P water molecules extended 10 Å from the complex. Counterions (Cl − or Na + ) were added to achieve electrostatic neutrality. Two successive energy minimizations were realized using the Sander module prior to MD simulations. First, the hydrogen atoms and water molecules were minimized with the first 2,500 steps of steepest descent and the last 2,500 steps of the conjugate gradient algorithm. Then, 2,500 cycles of the steepest-descent and 2,500 cycles of the conjugated gradient algorithm were used to minimize all atoms of the systems without restriction. Following minimization, these systems were gradually heated from 0 to 300 K in the NVT ensemble during 100 ps with a restrain force constant of 5 kcal/mol/Å 2 , which was applied to the backbone atom. Next, the 100 ps MD simulations with a 5 kcal/mol/Å 2 restraint on the backbone atoms were carried out at a pressure of 1 atm and 300 K to equilibrate the density using Berendsen's barostat. These systems were then equilibrated at 300 K and pressure of 1 atm without any constraints. Finally, MD simulations of 20 ns were performed at a constant temperature of 300 K and a pressure of 1 atm employing the PMEMD module in Amber14 (Case et al., 2014). The SHAKE method (Ryckaert et al., 1977) was used to constrain hydrogen atoms and the time step was set to 2.0 fs. The particle mesh Ewald (PME) algorithm (Darden et al., 1993) was applied to treat the long-range electrostatic interactions with default cutoff of 8.0 Å under periodic boundary conditions.

MM/GBSA Calculations
The coordinate trajectories were saved every 10 ps and stable conformations, that were generated from the last 2 ns of simulations, were applied to binding free energy calculations and analyses using MM/GBSA methods implemented in Amber14 (Case et al., 2014). In MM/GBSA, the binding free energy ( Gbind) was calculated as follows: (1) where G complex , G receptor , and G ligand represent the free energy of complex, protein, and ligand, respectively. G bind was evaluated via the entropy term (T S), the solvation energy term ( G solvation ), and the gasphase interaction energy ( E MM ) between inhibitors and receptors. E MM represents the gasphase interaction energy, which includes internal ( E internal ), electrostatic ( E electrostatic ), and van der Waals energies ( E VDW ). G solvation accounts for the change of solvation energy that combines polar and nonpolar components of the desolvation free energy. T S is the change of conformational entropy in response to ligand binding. The non-polar contribution was estimated via the solvent accessible surface area (SASA) using the LCPO method implemented in Amber (Weiser et al., 1999). The polar component of the desolvation energy was determined by the GB (igb = 2) of Onufriev et al. (2004). The solute and outer dielectric constants were set to 80 and 1, respectively. All energy components were performed for 2000 snapshots and extracted from the last 2 ns MD trajectories. Then, the energy decomposition was used to calculate the contribution of each residue of the systems to the total binding free energies using the MM/GBSA method in Amber14 (Case et al., 2014).

RESULTS AND DISCUSSION
Molecular Docking of PUGNAc, N-Valeryl-PUGNAc, and EtBuPUG Into VcNagZ, HsHexB, and hOGA To investigate the basis for both the potency and selectivity of PUGNAc derivatives toward VcNagZ, HsHexB, and hOGA, the binding modes of PUGNAc, N-valeryl-PUGNAc, and EtBuPUG into these three enzymes were observed via molecular docking. Firstly, density functional theory (DFT) calculations at the B3LYP/6-31G (d) basis set were performed to optimize the PUGNAc derivatives using the Gaussian 16 program. The optimized molecules are shown in Figure 2A. Compared to PUGNAc and N-valeryl-PUGNAc, the oxygen on the C-6 in the pyranose ring of EtBuPUG was found to have become perpendicular to the paper plane. This phenomenon may be because the 2-isobutamido group of EtBuPUG is relatively larger and changes the conformation of its pyranose ring. Furthermore, a blue positive area was observed near the pyranose ring of these compounds, as well as a slightly yellow negative area close to benzene, which complemented the charge in the active pocket of the enzymes ( Figure S1). Thus, good electrostatic interaction between small molecules and enzymes may be beneficial for protein-ligand binding. Subsequently, these three optimized inhibitors were docked into the active sites of VcNagZ, HsHexB, and hOGA, respectively. Among docking score results, total_score is expressed as the binding constant Log(Kd), where a high value of the total score represents a good binding ability. Crash indicates the degree of inappropriate penetration by the ligand into the protein and of interpenetration (selfclash) between ligand atoms (Pham and Jain, 2006). Crash scores close to 0 are favorable and beneficial for protein-ligand binding. Polar indicates the contribution of polar interactions to the total score. Similarity indicates the Surflex-Sim similarity between the scoring pose and the ligand (Jain, 2000).
As shown in Table 2, the total scores of compounds were consistent with the ranking of their experimental activities, and the similarity scores all exceeded 0.6, which indicated that the results of docking were trustworthy. Furthermore, the polar scores of these systems of PUGNAc combined with VcNagZ, HsHexB, and hOGA appeared to be very high, which indicates that PUGNAc have strong polar interactions with the three β-Nacetyl-D-hexosaminidases. The crash absolute value of EtBuPUG was largest, which may be because its bulky 2-acetamido substituent group leads to a relatively larger steric hindrance.
The detailed binding modes of PUGNAc, N-valeryl-PUGNAc, and EtBuPUG toward VcNagZ, HsHexB and hOGA are demonstrated in Figures 2B-D. The pyranose ring of these three inhibitors was found to be bound to the−1 subsite (active pocket) of all the three enzymes and the benzene ring extends out from the active pocket. Figure 2B shows that the inhibitors are located  well at the VcNagZ active site pocket and formed electrostatic, hydrophobic and hydrogen bond interactions with the residues Asp65, Arg140, Lys170, His171, Asp253, Asp254, and Met257, as previously reported (Vadlamani et al., 2017). Figure 2C shows the inhibitors bound to the active site of HsHexB and their interactions with residues Arg211, Glu355, Tyr450, Asp452, Tyr456, and Glu491. Figure 2D shows the binding modes between the inhibitors and hOGA, and small molecules interact with residues Gly67, Tyr69, Lys98, Asp174, Asp175, Phe223, Asn280, Asp285, Tyr286, Asn313, and Trp645. The above results indicated that the obtained docking sites are reliable and conducive to the development of molecular dynamics.

MD Simulations
In an effort to further investigate the appropriate interaction mode of PUGNAc, N-valeryl-PUGNAc, and EtBuPUG binding to VcNagZ, HsHexB, and hOGA, 20 ns MD simulations were conducted to further consider the flexibility of the whole protein and observe the dynamic binding mechanism of inhibitorenzyme systems. As illustrated in Figure 3a, the dynamic convergences were all achieved at approximately 16 ns and the RMSD values of the backbone atoms were ultimately maintained at around 1.6-3.1 Å for the nine simulated systems. These data indicated that the complexes underwent appropriate conformational change (Figure 3a and Figures S2,S3).  (d), respectively, at 20 ns of MD simulations. The ligands 1, 2, and 3 are shown as sticks with pink, cyan and yellow-orange color, respectively. The residues that interact with the compounds are shown as lines with green. Figure S4 represent the superimposition of the conformations of the three inhibitors bound into VcNagZ, HsHexB, and hOGA, respectively, at 20 ns of MD simulations. Compared to N-valeryl-PUGNAc and EtBuPUG, PUGNAc is deeply-buried in the active pockets of these three enzymes. The conformation of N-valeryl-PUGNAc is similar to that of PUGNAc but it is buried shallower than PUGNAc. This may result in the activity of N-valeryl-PUGNAc being lower than that of PUGNAc (reduced activity by 7-, 6,111-, and 870-fold against VcNagZ, HsHexB, and hOGA than against PUGNAc). EtBuPUG bearing branched 2-isobutamido further made its volume larger and unstable to binding into the active pocket, which resulted in a reduced activity by about 10-fold against VcNagZ, HsHexB, and hOGA compared to N-valeryl-PUGNAc. Figure S5 shows the two-dimensional (2D) representation of the binding modes of these compounds in the active pockets of VcNagZ, HsHexB, and hOGA, respectively. The results indicate that PUGNAc forms strong hydrogen bond interactions with the catalytic residues Asp253 (VcNagZ), Asp254 (VcNagZ), Asp354 (HsHexB), and Asp175 (hOGA). Moreover, the Nphenylcarbamate of PUGNAc also forms π-π stacking and hydrophobic interaction with the surrounding residues Phe36 (VcNagZ), Tyr450 (HsHexB), and Phe223 (hOGA). However, N-valeryl-PUGNAc and EtBuPUG hardly interacted with the corresponding catalytic residues of these enzymes. Furthermore, the number of hydrogen bonds between PUGNAc and VcNagZ, HsHexB, and hOGA exceeded that between N-valeryl-PUGNAc, EtBuPUG, and these three enzymes. These results provide the basis for the different potencies of PUGNAc, N-valeryl-PUGNAc, and EtBuPUG against the three enzymes.

Figures 3b-d and
Relative to the inhibitory activity of the PUGNAc derivatives, selectivity is a more attractive issue for exploration. PUGNAc binds well to VcNagZ, HsHexB, and hOGA; therefore, it displays a lack of selectivity. Upon increasing the size from N-acetyl to the N-valeryl group, N-valeryl-PUGNAc, possesses 667-and 121fold selectivity for VcNagZ over HsHexB and hOGA, respectively. Further increasing the size of the substituent, EtBuPUG bearing a larger 2-isobutamido, shows 968-and 109-fold selectivity for VcNagZ over HsHexB and hOGA. The reason for these significant selectivities was shown by the MD simulations results. Firstly, the active site architectures of VcNagZ, HsHexB, and hOGA were found to be sufficiently different (Figure 3). The volumes of their active pockets were calculated using the software CASTp (http://sts.bioe.uic.edu/castp/index.html). The order of the sizes of their active pockets is VcNagZ (31.9 Å) > hOGA (25.8 Å) > HsHexB (18.2 Å). The VcNagZ has a large pocket with remarkable flexibility that can accommodate various inhibitors. HsHexB has a shallow substrate-binding pocket. hOGA has a medium pocket that more willingly than HsHexB accommodates rather voluminous substituents attached to the location of the substrate 2-acetamido group (Henrissat and Davies, 1997;Stubbs et al., 2007;Liu et al., 2012). Therefore, N-valeryl-PUGNAc may remain stable in the active pocket of VcNagZ and has different degrees of deviations in the active sites of HsHexB and hOGA. EtBuPUG is located well within the active pocket of VcNagZ and completely moved outside of HsHexB; however, it could still enter into the active pocket of hOGA in a special mode ( Figure S4). These results reveal show the causes of selectivity from the structural characteristics of the enzymes. Furthermore, interactions between inhibitors and binding residues of these nine systems showed apparent differences ( Figure S5). N-valeryl-PUGNAc binds well into the pocket of VcNagZ by forming eight hydrogen bonds with Ala180, Arg73, and Asn285, but moves far from the pocket of HsHexB, which leads to interactions of four hydrogen bonds with Asp354, Asp452, and Glu491. The interactions between N-valeryl-PUGNAc and hOGA were also found to be weaker than those of N-valeryl-PUGNAc and VcNagZ. EtBuPUG remained firmly anchored in the open pocket of VcNagZ by forming five favorable hydrogen bond interactions with Asp65, Arg140, and His171. However, EtBuPUG cannot enter into the pocket of HsHexB, which results in activity loss with HsHexB, and only forms two hydrogen bond interactions with Arg211 and Glu491. These results may explain why EtBuPUG exhibits an 968-fold selectivity against VcNagZ over HsHexB.
Moreover, the differences in the catalytic mechanisms of VcNagZ, HsHexB, and hOGA also led to the selectivity of these PUGNAc derivatives. VcNagZ was first reported by Vocadlo et al. using a two-step double displacement mechanism (catalytic residues: Asp253 and Asp254) (Vocadlo and Withers, 2005). However, HsHexB (catalytic residues: Asp354 and Asp355) and hOGA (catalytic residues: Asp174 and Asp175) were reported to adopt a substrate-assisted mechanism, in which the 2-acetamido group of the non-reducing end sugar acts as a nucleophile, which results in the formation of a bicyclic oxazolinium intermediate (Terwisscha van Scheltinga et al., 1994;Tews et al., 1996;Mark et al., 2001). These differences may enable PUGNAc, N-valeryl-PUGNAc, and EtBuPUG to be selective to HsHexB/VcNagZ or hOGA/VcNagZ, especially the specificity of EtBuPUG to VcNagZ than HsHexB. This may also explain why the selectivity of these PUGNAc derivatives against hOGA/HsHexB is not large, and does not exceed 10-fold.

Binding Free Energy Calculation
To further explore the cause of the inhibitory potency and selectivity of PUGNAc derivatives against VcNagZ, HsHexB, and hOGA, binding free energy analyses of the nine systems were performed by using MM/GBSA calculation methods. The results show that the absolute value of the binding free energy is consistent with the trend of inhibitory activities of compounds to the corresponding enzyme. As shown in Table 3, the value of the binding free energy ( G TOT ) of PUGNAc is more negative than that of N-valeryl-PUGNAc and EtBuPUG, which means that PUGNAc has better inhibitory activities against the three tested enzymes than N-valeryl-PUGNAc and EtBuPUG. The binding free energy value of the HsHexB-EtBuPUG system is the maximum, indicating the worst inhibitory activity of EtBuPUG against HsHexB. Besides, experimental inhibition constants K i (in µM) of PUGNAc derivatives measured against VcNagZ, HsHexB and hOGA (Macauley et al., 2005;Stubbs et al., 2006Stubbs et al., , 2007 were converted into log unit (pK i = -LogK i ). We found the experimental binding affinities (pK i ) is in good agreement with the binding free energies calculated using MM/GBSA method results, and a linear correlation was obtained that yielded a good correlation coefficie (R 2 = 0.895) ( Figure S6).
Furthermore, four individual energy components ( E vdw , E ele , G SA , and G GB ) were carefully analyzed to better understand the energy terms with the greater impact on binding affinity. The van der Waals ( E vdw ) and electrostatic ( E ele ) contributions are highly favorable for the inhibitor binding to enzymes, and the electrostatic interactions exceed the van der Waals interactions. The data exhibited here are consistent with the previously reported electrostatic environment in which the active pocket of glycosidase is beneficial for the elimination of barriers (Passos et al., 2011). Moreover, the order of van der Waals interactions relative to the specific enzyme was EtBuPUG > N-valeryl-PUGNAc > PUGNAc, which matches the size sequence of the compound. The larger substituent at the 2position of these compounds maybe closer to the hydrophobic region of the protein. The nonpolar desolvation energies ( G SA ) among these nine systems are not much different and range between −4.33 and −6.76 kcal mol −1 . However, the polar contributions of desolvation ( G GB ) are very unfavorable for the binding of PUGNAc derivatives and these enzymes, which could offset their electrostatic interaction. Thus, the nonpolar interactions ( E vdw + G SA ) in the system are more favorable for the binding of ligands to receptors.

Decomposition Analysis of the Binding Free Energies
To obtain a detailed quantitative analysis of the contribution of each residue to VcNagZ, HsHexB, and hOGA with PUGNAc, N-valeryl-PUGNAc, and EtBuPUG interactions, MM/GBSA free energy decomposition analyses were performed to decompose the total binding free energies into residues. Residues near 4 Å of the ligands are used to calculate the decomposition free energy, and the key residue contributions of VcNagZ, HsHexB, and hOGA against the PUGNAc, N-valeryl-PUGNAc, and EtBuPUG, respectively, are shown in Figure 4.
In the VcNagZ-inhibitor systems, the differences in interactions are mainly reflected in the six residues, Arg140, Lys170, Ala180, Asp253, Asp254, and Asn285 ( Figure 4A). Among these, the interactions of PUGNAc with the key residues Lys170, Asp253, and Asp254 are stronger than those of N-valeryl-PUGNAc and EtBuPUG. In particular, PUGNAc formed the strongest interaction with catalytic residue Asp253 and Asp254 (−4.57 and −5.61 kcal mol −1 ), while N-valeryl-PUGNAc only bound to Asp254 with moderate interaction (−3.29 kcal mol −1 ) and EtBuPUG only bound to Asp253 with lesser interaction (−2.39 kcal mol −1 ). These results could explain why the inhibitory activity of PUGNAc against VcNagZ was better than that of N-valeryl-PUGNAc and EtBuPUG. In addition, the movement of N-valeryl-PUGNAc in the VcNagZ pocket moved it closer to Ala180 and Asn285, thus forming more favorable interactions with them (−2.68 and −4.31 kcal mol −1 ). Similarly, the movement of EtBuPUG in the active pocket brings its pyranose closer to Arg140, resulting in greater interaction (−4.0 kcal mol −1 ).
These results may explain why the order of activity against hOGA was PUGNAc > N-valeryl-PUGNAc > EtBuPUG. Furthermore, the apparent differences between these three inhibitors and hOGA were mainly identified on the catalytic residue Asp175. Only PUGNAc possesses a strong interaction with the value of−2.62 kcal mol −1 . N-valeryl-PUGNAc and EtBuPUG could not interact with Asp175, resulting in lower inhibitory activity.
The decomposition analysis of the binding free energies explained why the inhibitory activities of N-valeryl-PUGNAc and EtBuPUG against VcNagZ was better than those of hOGA and HsHexB. The reason is that N-valeryl-PUGNAc and EtBuPUG cannot form an important binding interaction with the key residues Asp175 (hOGA) and Asp354 (HsHexB) in the active pockets. The contribution of the catalytic residues Asp253 and Asp254 (VcNagZ), Asp354 (HsHexB), and Asp175 (hOGA) in the enzyme-inhibitor systems were important for distinguishing the activity of these inhibitors.

CONCLUSION
A computational approach that combines molecular docking, and MD simulations was employed to provide insights into the activity and selectivity of PUGNAc derivatives against GH3, GH20, and GH84 β-N-acetyl-D-hexosaminidases (VcNagZ, HsHexB, and hOGA). Firstly, the molecular docking combined with 20 ns MD simulations were performed to investigate the differences of these three enzyme active pocket architectures as well as the possible inhibitory patterns of the nine inhibitorenzyme systems. The results show that the best-known inhibitor PUGNAc is located well at the active pockets of all three enzymes; therefore, it lacks selectivity. When the size of 2-acetyl of the sugar ring of PUGNAc is enlarged, the result is selectivity despite a degree of decrease in activity. Specifically, N-valeryl-PUGNAc, bearing a larger line 2-valeryl group, could be accommodated in the pockets of VcNagZ and hOGA. In particular, EtBuPUG, with the largest branched 2-isobutamido group, could only bind well to the VcNagZ pocket but was completely removed from the active pocket of HsHexB, thus lacking interactions with HsHexB. These distinctive differences in binding modes may lead to EtBuPUG exhibiting 968-fold selectivity to VcNagZ over HsHexB.
Analyses of MM/GBSA free energy calculations and MM/GBSA free energy decomposition demonstrate that the predicted energy of the PUGNAc-enzymes complexes is more favorable than that of the N-valeryl-PUGNAc-enzyme and the EtBuPUG-enzyme systems. Moreover, the van der Waals and electrostatic contributions are favorable to the combination of PUGNAc derivatives with the three enzymes. Furthermore, the catalytic residues Asp253 (VcNagZ), Asp254 (VcNagZ), Asp354 (HsHexB), and Asp175 (hOGA) are the most important residues for distinguishing the potential of the related inhibitors.
This is the first systematic study of the possible binding mechanisms for the dynamic interaction of GH3, GH20, and GH84 β-N-acetyl-D-hexosaminidases with inhibitors in solution. It is more in line with the interaction between drugs and target enzymes in the body. These computational study results indicate that positively charged functional groups (to bind to residue Asp) as well as an aromatic group at suitable positions of inhibitors would be helpful to improve the inhibitory potency against β-N-acetyl-D-hexosaminidases. On the other hand, the selectivity of inhibitors is mainly related to the pocket size of the enzymes. The bulky and branched group at the position of the 2-acetamido of the substrates would be beneficial for the selectivity against GH3 β-N-acetyl-D-hexosaminidases over GH20 and GH84 β-Nacetyl-D-hexosaminidases. Therefore, selectivity to different GH β-N-acetyl-D-hexosaminidases could be achieved by regulating the size of the group at this position. In summary, this work provides considerable structurally guided improvements for the development of novel, potent, and specific β-N-acetyl-Dhexosaminidases inhibitors.

AUTHOR CONTRIBUTIONS
The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.