Abstract Heparan sulfate (HS) is a sulfated polysaccharide that plays a key role in morphogenesis, physiology and pathogenesis. The biosynthesis of HS takes place in the Golgi apparatus by a group of enzymes that polymerize, epimerize and sulfate the sugar chain. This biosynthetic process introduces varying degrees of sulfate substitution, which are tightly regulated and directly dictate binding specificity to different cytokines, morphogens and growth factors. Here, we report the use of molecular dynamics simulations to investigate the dynamics of substrate recognition of two glycosaminoglycan (GAG) sulfotransferases, N-deacetylase-N-sulfotransferase and 2-O-sulfotransferase to the HS chain during the biosynthetic process. We performed multiple simulations of the binding of the sulfotransferase domains to both the HS oligosaccharide substrate and sulfate donor, 3′-phosphoadenosine-5′-phosphosulfate. Analysis of extended simulations provide detailed and useful insights into the atomic interactions that are at work during oligosaccharide processing. The fast information matching method was used to detect the enzyme global dynamics and to predict the pairwise contact of residues responsible for GAG–enzyme binding and unbinding. The correlation between HS displacement and the location of the modified GAG chain were calculated, indicating a possible route for HS and heparin during sulfotransferase processing. Our data also show sulfotransferases contain a conserved interspaced positively charged amino acid residues that form a patch which controls the protein–GAG binding equilibrium. Together, our findings provide further understanding on the fine-tuned complex mechanism of GAG biosynthesis. Our findings can also be extrapolated to other systems for calculating rates of protein–GAG binding. glycosaminoglycan binding, Heparan sulfate, molecular dynamics simulation, sulfotransferases Introduction Heparan sulfate proteoglycans (HSPGs) are ubiquitously expressed at the cell surface and in the extracellular matrix of most living organisms, with convergence coinciding with the emergence of eumetazoa (Medeiros et al. 2000; Gesteira, Coulson-Thomas, Ogata, et al. 2011). HS is a polysaccharide that interacts with a wide variety of bioactive ligands, such as growth factors, morphogens (e.g., fibroblast growth factor, FGF) and their receptors (e.g., FGF receptor), extracellular matrix proteins (e.g., collagen), proteases, cytokines, chemokines and enzyme inhibitors (Bishop et al. 2007; Tecle et al. 2013; Coulson-Thomas et al. 2014; Xu and Esko 2014). The fine structure of HS/heparin (HEP) chains dictates their binding affinity to a plethora of biological ligands, thereby defining the specific roles HS/HEP plays in development and pathology and dictating their physiological function. More specifically, distinctive sulfation patterns along HS/HEP chains confer the specificity of HS/HEP to different ligands. The biosynthesis of HEP/HS glycosaminoglycans (GAGs) starts in the endoplasmic reticulum with the addition of a xylose residue to specific locations along a core protein. This protein is then transported to the Golgi apparatus where N-acetyl-glucosamine (GlcNAc) or glucuronic acid (GlcA) residues are sequentially added to the growing HS/HEP chain by polymerizing enzymes located in the Golgi, the exostosin (EXT) family of glycosyltransferases (Busse et al. 2007; Presto et al. 2008; Multhaupt and Couchman 2012) initiating HS chain elongation. The growing chain undergoes further substoichiometric and grouped modifications which introduce sulfate groups by means of specialized enzymes, namely, N-deacetylase/N-sulfotransferases (NDSTs) and 2-, 6- and 3-O-sulfotransferases (OSTs), and epimerize glucuronic residues to iduronic forms by C5-epimerase. The specificity and expression levels of these modifying enzymes dictate the fine structure of HS/HEP, producing characteristic domains which vary in size and number within each chain, dictating the biological functions of HS/HEP. NDST, the first modifying enzyme, removes the acetyl group from GlcNAc and adds a sulfate group, generating N-sulfoglucosamine (GlcNS, Figure 1A). The generation of GlcNS introduces N-sulfated (NS) domains along the HS chain. When these NS domains are flanked by unmodified N-acetylated domains (NA domains) NA/NS domains are generated. This organization specifies different protein-binding motifs; for example, it has been shown that FGFs specifically bind to NS domains (Kreuger et al. 2001), whereas interleukin-8 and vascular endothelial growth factor bind via interactions involving multiple spaced, modestly sulfated NS domains (Spillmann et al. 1998; Robinson et al. 2006). The precise mechanism that determines which locations along the HS chains should be modified by NDST, thereby regulating the presence of these NS domains, remains elusive. Sulfate transfer is dependent on a universal sulfate donor, 3′-phosphoadenosine-5′-phosphosulfate (PAPS) (Figure 1B), and despite its significant resemblance to phosphoryl-transfer reactions does not require divalent cations for activity (Skelton et al. 1991; Hooper et al. 1995; Lo-Guidice et al. 1995; Shailubhai et al. 1999). Fig. 1. View largeDownload slide GAG NST and 2OST-binding motifs. (A) and (C): Proteins depicted as a cartoon, and PAPS and oligosaccharides are shown as sticks. Reactive residues are labeled (Glu642 for NST and His142 for 2OST) and sugar acceptor atoms are indicated by an arrow. Panels (B) and (D): Sketch illustrating catalytic reactions at the active site. Coupled proton transfer takes place from the amine to Glu642 for NST (B) or His142 for 2OST (D) to allow sulfate transfer. GAG, glycosaminoglycan. Fig. 1. View largeDownload slide GAG NST and 2OST-binding motifs. (A) and (C): Proteins depicted as a cartoon, and PAPS and oligosaccharides are shown as sticks. Reactive residues are labeled (Glu642 for NST and His142 for 2OST) and sugar acceptor atoms are indicated by an arrow. Panels (B) and (D): Sketch illustrating catalytic reactions at the active site. Coupled proton transfer takes place from the amine to Glu642 for NST (B) or His142 for 2OST (D) to allow sulfate transfer. GAG, glycosaminoglycan. Even though great strides have been made in unveiling the biosynthetic pathway of HS/HEP, there is still a scarcity of data concerning the specific features of each enzyme which define the sulfation pattern, including the temporal control of these modifications and the potential formation of a Golgi biosynthesis complex named GAGosome (Esko and Selleck 2002). Interestingly, excess of PAPS has been shown to increase the incidence of NS domains demonstrating that PAPS is a limiting factor in the sulfation process (Carlsson et al. 2008) and thus may add a level of control to the complex modification process of HS and HEP chains (Bhattacharya et al. 2009) and other GAGs (Dick et al. 2008). Subsequent modifications to HS chains, namely, the epimerization of GlcA by C5-epimerase and sulfation of both N-glucosamine (GlcN) and hydroxyl groups by distinct sulfotransferases, are believed to take place in a cooperative sequential manner, with evidence showing enzyme interaction and the formation of an enzyme complex. These sulfate modifications are a prerequisite for the epimerization of HSPGs. The epimerization of C5 GlcA to iduronic acid not only increases HS chain flexibility but also determines the binding specificity for ligand recognition (Jia et al. 2009). The epimerization of HS occurs in sequence with the formation of 2- and 6-O-sulfated domains (Gama et al. 2006). OSTs also use the universal sulfate donor PAPS (Figure 1C and D), with key residues responsible for PAPS binding organized into two regions that are conserved among all sulfotransferases, named the 5′-phosphosulfate binding (5′PSB) and 3′-phosphate binding (3′PB) motifs (Kakuta et al. 1997, 1998). Up until now, emphasis has been generally placed on predicting the most likely binding processes for GAG-binding molecules that are predictive of high binding affinity (Babik et al. 2017; Samsonov et al. 2011). To date, the processing pathways of these enzymes have not been predicted, and the mean life span of the protein–GAG complexes has never been analyzed, despite being one of the most pertinent factors for sustained drug efficiency and safety aside from affinity. This property is intrinsically related to the time during which the GAG chains remain in the binding site, while the most distinctive feature of the unbinding process is the role played by the residues surrounding the open binding cleft. Many factors can influence the movement of oligosaccharides though the active site, including solvent interactions that promote unbinding by assisting in the breakage of shielded hydrogen bonds through the formation of water bridge interactions. This study aimed to extract the causality of correlated motions from GAG sulfotransferases molecular dynamics (MD) simulations, providing vital information on the conformational dynamics of enzyme-bound oligosaccharides during GAG sulfation assembly. Results Behavior of residues lining the sulfotransferase-binding cleft Using MD simulations of NST and 2OST complexed with PAPS and their respective oligosaccharides, we undertook a systematic investigation to elucidate the factors influencing uncoupling events. To date, the binding site for the oligosaccharide chain has been proposed within the large open cleft with a hydrophilic surface, a random coil (turn 5, residues 637–639) and an opposing alpha helix (alpha helix 6) forming the center of the cleft near the 5′-sulfate of the PAPS molecule for NST (Sueyoshi et al. 1998; Kakuta et al. 1999). The amphipathic random coil positions the putative charged active site residue side chains of Glu641 and Glu642 toward the center of the cleft, while the opposing region of the cleft contains side chains of residues required for binding, His716, Gln717 and His720. In unrestrained nanosecond MD simulations, subtle changes in oligosaccharide binding could be observed. An ample variation in the distance between amino acids involved in both catalysis and binding could be observed for both NST and 2OST simulations, termed coupling (approximation) and uncoupling (repulsion) events. Curiously, uncoupling events also occur in 2OST simulations. Despite structural differences between NST and 2OST, the 2OST catalytic domain also comprises a canonical sulfotransferase α/β domain that contains the PAPS-binding loop. The most noticeable difference is the occurrence of intermolecular contacts within the 2OST trimer, by extension of a C-terminal tail from each monomer across the active site of the neighboring monomer (Liu et al. 2014). A similar scaffold with a clamp-like orientation can be found in this structure, constituted by the discontinuous basic residues, K111, R184, R190 and K350, also forming the shape of a parallelogram with interior angles of ~100°. Residues Arg184 and Arg190 are disordered in the substrate-free structure but are found ordered in the bound ternary structure (Liu et al. 2014). Similar residues are conserved in all GAG sulfotransferases and also other biosynthetic proteins, such as C5-glucuronosyl transferase (Supplementary data, Figure S1). One interesting feature is the presence of a lysine (K350) from the adjacent monomer that helps to stabilize the bound oligosaccharide, from the reducing end, facilitating oligosaccharide binding to the clamp region that forms an antiparallel strand (N345–Y352) with the fifth strand (V104–T110) of the central beta sheet of the adjacent monomer (Liu et al. 2014; Xu et al. 2007). These studies reveal that Arg80, Lys350 and Arg190 of 2OST interact with the N-sulfo groups near the modification site, consistent with the dependence of 2OST on N-sulfation. Unlike NST, 2OST uses a histidine instead of a glutamate as the catalytic base. This residue is also present in cytosolic sulfotransferases. Analysis of trajectory clustering Trajectory clustering shows segments of generalized masked Delaunay (GMD) minima where conformational states represent opening regions formed of a long alpha helix (V182–L188—equivalent to alpha helix 6 in EST) and the fourth strand of their central beta sheet during the evolution of the simulation. These regions form the “buns” of the cleft (Figure 3) that sandwich the modified oligosaccharide and present a minimum of 1.1 during coupling simulations and 0.85 during uncoupling simulations. The structure of the minima of uncoupling simulations presents greater cleft conformational change, with again alpha helix 6 and the end of the fourth strand providing the main contributions to the transition events and protein activity, and, to a lesser extent, the loops formed by residues 280–292 (Figure 3A, C and E). Prediction of GAG uncoupling The structural transitions in NST and 2OST simulations were additionally characterized using a method based on higher order statistics using the TimeScapes software (Wriggers et al. 2009). This analysis is based on monitoring the time-dependent formation and breaking of side chain contacts in MD trajectories to pinpoint conformational rearrangements. We performed the analysis on each of the NST and 2OST coupling and uncoupling simulations, measuring GMD pair distribution-based system activity, a variable that is defined as the total number of contact changes over unit of time. From this analysis, we identified basins of minimal activity (GMD minima) in the trajectory. We then analyzed in detail the structures of each minimum identified in the GMD analysis paying attention to the conformation of residues surrounding NST and 2OST. TimeScapes analysis revealed that the clamp region “sandwiching” the oligosaccharides presents similar conformational propensities in both sulfotransferase simulations. Indeed, the GMD graph of uncoupling simulations showed that the events of high activity (Figures 2A and B and 3A and B for NST and 2OST, respectively) are associated with rearrangements in a region surrounding the 5′PSB motif, promoting the loss of recessed states in favor of a relaxed-like state. Additionally, in each GMD minima of both NST and 2OST simulations, the 5′PSB motif did not undergo large rearrangements, with low-activity events associated with uncoupling events and the maximum area of a triangular region formed by the PAPS sulfate and residues K111 and R188 for 2OST and H716 and T639 for NST. There are, however, often associated conformational changes in the loops composed of residues 261–285 for NST and 181–196 for 2OST (Figures 2 and 3C and D), which are associated with decreases in the area of bound conformations. TimeScapes analysis therefore complements the RMSF-based conclusions, confirming the high mobility of these loops. Fig. 2. View largeDownload slide A: Activity level during N-sulfotransferase uncoupling (A, C and E) and coupling (B, D and F) simulations. A contact metric based on higher order statistics (GMD) was employed to describe conformational changes in each NST simulation, measured as activity (refer to Materials and Methods). Basins of minimal activity (GMD minima, A and B) of the trajectories corresponding to maxima relaxation structures are reported as cartoons in C and D, with the reference structure starting at t = 0 (X-ray structure) and extended states. Structural states of turn 5 and opposing alpha helix 6 in GMD minima are represented by a cartoon colored according to the denoting timepoint in the GMD-activity plot. Structural states of turn 5 and opposing alpha helix 6 in GMD minima are represented by rectangular boxes colored according to their structural similarity with extended or recessed states measured by distance and area formed between T639 and H716 and the PAPS sulfate (putative binding cleft area E and F). Fig. 2. View largeDownload slide A: Activity level during N-sulfotransferase uncoupling (A, C and E) and coupling (B, D and F) simulations. A contact metric based on higher order statistics (GMD) was employed to describe conformational changes in each NST simulation, measured as activity (refer to Materials and Methods). Basins of minimal activity (GMD minima, A and B) of the trajectories corresponding to maxima relaxation structures are reported as cartoons in C and D, with the reference structure starting at t = 0 (X-ray structure) and extended states. Structural states of turn 5 and opposing alpha helix 6 in GMD minima are represented by a cartoon colored according to the denoting timepoint in the GMD-activity plot. Structural states of turn 5 and opposing alpha helix 6 in GMD minima are represented by rectangular boxes colored according to their structural similarity with extended or recessed states measured by distance and area formed between T639 and H716 and the PAPS sulfate (putative binding cleft area E and F). Fig. 3. View largeDownload slide A: Activity level during 2OST uncoupling (A, C and E) and coupling (B, D and F) simulations. The same contact metric higher order (GMD) employed for NST was used to describe conformational changes in each 2OST simulation. Basins of minimal activity (GMD minima, A and B) of the trajectories corresponding to maxima relaxation structures are reported as cartoons in C and D, with the reference starting structure at t = 0 (X-ray structure) and extended states. Structural states of turn 5 and opposing alpha helix 4 in GMD minima are represented by cartoons colored according to the denoting timepoint in the GMD-activity plot. These structural states in both GMD maxima and minima are represented by rectangular boxes colored according to their structural similarity with extended or recessed states measured by distance and area formed between R188 and K111 and the PAPS sulfate (putative binding cleft area E and F). Fig. 3. View largeDownload slide A: Activity level during 2OST uncoupling (A, C and E) and coupling (B, D and F) simulations. The same contact metric higher order (GMD) employed for NST was used to describe conformational changes in each 2OST simulation. Basins of minimal activity (GMD minima, A and B) of the trajectories corresponding to maxima relaxation structures are reported as cartoons in C and D, with the reference starting structure at t = 0 (X-ray structure) and extended states. Structural states of turn 5 and opposing alpha helix 4 in GMD minima are represented by cartoons colored according to the denoting timepoint in the GMD-activity plot. These structural states in both GMD maxima and minima are represented by rectangular boxes colored according to their structural similarity with extended or recessed states measured by distance and area formed between R188 and K111 and the PAPS sulfate (putative binding cleft area E and F). Dynamics of the GAG sulfotransferase interaction To further investigate how glycan binding may be altered by interaction with different amino acids at the GAG sulfotransferase interface, we analyzed bond formation distance from the sulfotransferase catalytic clefts. Enzyme residence times varied with a trend observed during uncoupling events, namely, increased bonding with the solvent toward the reducing end of the chain (Figure 4). An overall significant decrease in hydrogen bonding occurs in uncoupling simulations, which is accompanied by a slight increase in the bonding of reducing end monosaccharides, indicating imminent unbinding. NST simulations were performed with a previously modeled hexasaccharide (Gesteira et al. 2013) with the non-reducing end moiety (glucosamine) forming a hydrogen bond with R750 for 81.3% of the time in coupling simulations. This number falls to 32.9% for uncoupling simulations, while there is a parallel modest increase in hydrogen bonding between the reducing end GlcA and K831. 2OST displayed a similar pattern, with hydrogen bonding between the reducing GlcA and R354 being 90.4% during the course of coupling simulations and 41.7% during uncoupling simulations. A minor increase in the extent of the solvation shell can be noted in the vicinity of the active site monosaccharide closest to PAPS, possibly due to the overall forward motion of the oligosaccharide through the cleft. Evidence suggests that N-sulfation occurs from the non-reducing end toward the reducing end (Carlsson et al. 2008; Smeds et al. 2010; Sheng et al. 2011), while 2- and 3-O-sulfation occur in the same orientation to that of polymerization (reducing to non-reducing end) (Smeds et al. 2010; Prechoux et al. 2015). Fig. 4. View largeDownload slide Total number of monosaccharide hydrogen bonds during NST and 2OST simulations. Side chain oligosaccharide hydrogen bonds are shown as columns. Column bars represent average hydrogen bonds. Column color symbolizes the oligosaccharide chain, while the overlay line chart represents water solvation. The error bars represent the standard deviation of the average values calculated over three independent MD runs performed in this study. NRE, non-reducing end. Fig. 4. View largeDownload slide Total number of monosaccharide hydrogen bonds during NST and 2OST simulations. Side chain oligosaccharide hydrogen bonds are shown as columns. Column bars represent average hydrogen bonds. Column color symbolizes the oligosaccharide chain, while the overlay line chart represents water solvation. The error bars represent the standard deviation of the average values calculated over three independent MD runs performed in this study. NRE, non-reducing end. At the residue level, both NST– and 2OST–GAG interfaces show a few common hydrogen bond patterns. Hydrogen bond analysis for both coupling and uncoupling events demonstrates the variability of hydrogen bond formation with residues during both simulations. We verified that both NST and 2OST coupling and uncoupling simulations presented specific side chain hydrogen bond patterns which contributed toward the spacing and partial chain dislodgement of the oligosaccharide chain. Several hydrogen-bonded residue pairs were constant during the course of both coupling and uncoupling simulations, and this trend was not dependent on the length of neighboring contacts. GlcNA1 and GlcNA3 show hydrogen bond propensity in coupling simulations, while GlcNA3 hydrogen bonding decreased considerably and GlcA6 hydrogen bonding increased in uncoupling simulations (Fig. 4). 2OST-bound terminal GlcA also showed a more consistent hydrogen bonding pattern in coupling simulations, but a one-third drop was observed in uncoupling simulations. As observed in NST uncoupling simulations, the third residue from the active site showed increased hydrogen bonding in 2OST uncoupling simulations. MD simulations for NST and 2OST co-complexes showed an overall stable structure with minimal root-mean-square deviation (RMSD) from the average structures (only 3.2 and 3.4 Å, respectively—Supplementary data, Figure S2). Likewise, the two uncoupling co-complexes displayed deviations in a similar range (Supplementary data, Figure S2). Decomposed analyses of these differences in hydrogen bond formation were performed using the Visual Molecular Dynamics (VMD v.1.9.3) tool with a donor–acceptor distance cutoff of 3.5 Å and an angle cutoff of 60° (Humphrey et al. 1996). Figure 5 shows hydrogen bond consistency between active binding patch residues and ligands for both enzymes during the course of both coupling and uncoupling simulations. The results show that, overall, hydrogen bonds are more consistently maintained in coupling simulations than in uncoupling simulations, e.g., in the case of Lys676, Arg750 and Arg763, these occur in more than 75% of frames for NST coupling simulations and are clearly weaker in NST uncoupling simulations, suggesting atomic interaction disparities. The analysis revealed that Lys676 was the residue that most contributed toward hydrogen bonds in coupling simulations, whereas only a modest increase in oligosaccharide hydrogen bonding was observed for Lys833 and Arg835 during uncoupling simulations. An increase in bonding for residues at the top of turn 5, Lys833 and Arg835 (increase of 17% and 22%, respectively), was also observed. Fig. 5. View largeDownload slide Intermolecular hydrogen bond consistency across the MD simulation (final 10 ns) for NST and 2OST active site residues. The occurrence of oligosaccharide–protein hydrogen bonds between NST and 2OST residues (shown on the y-axis) and each oligosaccharide are shown for each fifth frame of each of the MD runs. (A) and (B) NST coupling and decoupling simulations. (C) and (D) 2OST coupling and decoupling simulations. Fig. 5. View largeDownload slide Intermolecular hydrogen bond consistency across the MD simulation (final 10 ns) for NST and 2OST active site residues. The occurrence of oligosaccharide–protein hydrogen bonds between NST and 2OST residues (shown on the y-axis) and each oligosaccharide are shown for each fifth frame of each of the MD runs. (A) and (B) NST coupling and decoupling simulations. (C) and (D) 2OST coupling and decoupling simulations. A similar analysis of the dynamics of 2OST coupling and uncoupling simulations shows that both Arg188 and Arg190 contribute almost equally to hydrogen bond occupancy, both losing more than 80% of contacts in uncoupling simulations. The remaining interface residues show <40% change in hydrogen bonding, and the percentage of long-ranged hydrogen bonds was greater in the coupling simulations in comparison to the uncoupling simulations. Another noteworthy feature was the long-lived hydrogen bonds that occurred in both coupling NST and 2OST simulations. During uncoupling simulations, however, small areas of consistent bonding could be observed, namely, Lys833 and Arg835 (Figure 5B, 47–62 ns) for NST and Lys111 (Figure 5D, 32–72 ns) for 2OST simulations. A quantitative analysis of the interplay of hydrogen bonds formed by the amino acids at the protein–GAG binding interface can explain the observed coupling and uncoupling conformations. We decomposed the hydrogen bonds into four groups: (i) antiparallel beta sheet and loop-oligosaccharide, (ii) antiparallel beta sheet and loop-opposing alpha helix, (iii) oligosaccharide-opposing alpha helix and (iv) oligosaccharide solvent (see Materials and methods for details). Figure 6A and B shows the hydrogen bonds and their average persistence as a chord diagram for NST groups (i), (ii) and (iii). During NST coupling simulations, numerous persistent hydrogen bonds were present between the beta sheet and the oligosaccharide, as this region defines the GAG-binding open cleft. Two slightly persistent hydrogen bonds were formed between the beta sheet (His716 and E642) and the opposing oligosaccharide chain (GlcA2 and GlcN3) (Figure 6A). Very few hydrogen bonds were found within the beta sheet and opposing alpha helix. For uncoupling simulations, a decrease in beta sheet-interacting residue hydrogen bonds was followed by an increase in opposing alpha helix hydrogen bonding. In addition to this shift in hydrogen bond formation, a decrease in hydrogen bond persistence was noted. This shift was also observed during 2OST coupling/uncoupling simulations. Hydrogen bond persistence was greater for His154 and Y94 during coupling simulations, while a shift toward hydrogen bonding between H191 and GlcNS4 was noted in uncoupling simulations. The average number of hydrogen bonds was also smaller compared to coupling simulations. In an analogous manner, a decrease in GAG hydrogen bonding to the antiparallel strand in the binding cleft occurred during 2OST uncoupling simulations. While most of the contacts were maintained in both simulations, the average number of hydrogen bonds formed between C-terminal residues and the antiparallel strand was inverted (Figure 6D–F). By the clamp opening, the GAG is released and forms interactions allowing the establishment of stable hydrogen bonds with C-terminal residues, as the decrease in hydrogen bonding is followed by a noticeable increase in hydrogen bonding with the C-terminal tail of the adjacent monomer (N345–Y352, Fig. 6D–F). Curiously, the residues surrounding monosaccharides in the active site, Tyr183, Lys191 and Tyr94, presented a significant decrease in hydrogen bonds in uncoupling simulations, while residues located toward the non-reducing end (GlcNS4 and GlcA5) were “pulled” by Lys350 and Tyr352. Fig. 6. View largeDownload slide Analysis of the hydrogen bond network. Chord diagram showing the hydrogen bonds formed by the amino acids at the protein–GAG binding interface in (A) NST coupling, (B) NST uncoupling, (D) 2OST coupling and (E) 2OST uncoupling. Hydrogen bond directionality (donor → acceptor) is represented by an arrow in the connection link. Links transparency indicates hydrogen bond persistence, with a lighter color being lower persistence and a grayscale color bar given as a visual guide. The average number of hydrogen bonds per frame is given between parentheses. Hydrogen bond analysis was performed on standard MD simulations. (C) and (F) show a cartoon representation of the region analyzed. Fig. 6. View largeDownload slide Analysis of the hydrogen bond network. Chord diagram showing the hydrogen bonds formed by the amino acids at the protein–GAG binding interface in (A) NST coupling, (B) NST uncoupling, (D) 2OST coupling and (E) 2OST uncoupling. Hydrogen bond directionality (donor → acceptor) is represented by an arrow in the connection link. Links transparency indicates hydrogen bond persistence, with a lighter color being lower persistence and a grayscale color bar given as a visual guide. The average number of hydrogen bonds per frame is given between parentheses. Hydrogen bond analysis was performed on standard MD simulations. (C) and (F) show a cartoon representation of the region analyzed. Binding free energies We calculated the theoretical binding free energy (ΔG) using the MM/PB(GB)SA method as implemented in Ambertools14 (Miller et al. 2012; Roe and Cheatham 2013). The binding energies calculated for NST coupling and uncoupling co-complexes were −92.1 ± 7.2 and −53.9 ± 3.8 kcal/mol, respectively. For 2OST, these binding energies were −121.1 ± 7.1 and −105.3 ± 6.9 kcal/mol. Here, it could also be seen that the overall similarity in binding energies for the coupling and uncoupling simulations did not result in identical interactions at an individual residue level. For example, Lys126 (−12.0 ± 3.7 kcal/mol) and Arg121 (−11.1 ± 1.4 kcal/mol) were identified as the prime affinity contributors for coupling events in NST and 2OST coupling simulations, respectively, whereas Gln135 (−10.2 ± 1.6 kcal/mol) and Arg121 (−9.0 ± 2.8 kcal/mol) were the dominant residues for NST and 2OST uncoupling simulations, respectively. In addition, differences in binding energy contributions from several other residues were evident between coupling and uncoupling events. Individual contributions (Eqs. (1)–(3)) to the free energy of binding of oligosaccharides to sulfotransferases are reported in Table I. Systematic errors have been estimated from the hysteresis between backward and forward alchemical transformations, with free energy contributions of all oligosaccharide positional and orientational restraints being smaller in the bound state, ~1.5 kcal/mol, in line with the soft harmonic potentials employed in the simulations. The hysteresis between the corresponding coupling and uncoupling simulations never exceeds at most a few 10th of kcal /mol (Supplementary data, Figures S5 and S6). The bulk contributions, ΔGbulkr,θ,ϕ,ψ and ΔGbulkΘ,Ψ,Φ, computed analytically, are about one order of magnitude larger than their bound-state complements. Free energy associated with coupling NST and 2OST to each respective oligosaccharide, ΔGsitealch, are larger than for 6OST (Table I). Backward and forward estimates of ΔGsitealch and ΔGsitbulkalch present a marginal hysteresis, resulting in very small errors (see Table I and Supplementary data, Figure S6). Convergence of all free energy perturbation (FEP) simulations was checked thanks to the ParseFep plugin54 of VMD (Supplementary data, Figure S6). Table I. Energetic components of the free energy of oligosaccharide binding to sulfotransferases Free energies (kcal/mol) ΔGrepu ΔGdisp ΔGelec ΔGrstr ΔGTot ΔGb NST + olig 326.83 ± 1.6 −328.91 ± 1.42 −1145.53 ± 2.94 −2.7 −1023.10 ± 4.62 −101 ± 4.6 olig(GlcNH–GlcA)3 221.8 ± 4.3 −301.02 ± 3.0 −997.45 ± 1.4 0 −921.44 ± 2.9 2OST + olig 422.6 ± 2.9 423.92 ± 1.3 −1395.38 ± 4.2 −5.3 −1255.44 ± 11.1 −209.51 ± 3.9 olig(GlcNS–IduA)3 345.22 ± 1.4 314.68 ± 2.6 1378.22 ± 3.1 0 1045.93 ± 1.3 6OST + olig 162.9 ± 1.3 121.24 ± 1.7 −445.87 ± 2.5 −1.3 −154.21 ± 5.9 −42.72 ± 3.9 olig(GlcNS–IduA)3 −131.21 ± 3.1 −143.19 ± 2.6 −117.81 ± 7.5 0 −111.4.9 ± 1.5 Free energies (kcal/mol) ΔGrepu ΔGdisp ΔGelec ΔGrstr ΔGTot ΔGb NST + olig 326.83 ± 1.6 −328.91 ± 1.42 −1145.53 ± 2.94 −2.7 −1023.10 ± 4.62 −101 ± 4.6 olig(GlcNH–GlcA)3 221.8 ± 4.3 −301.02 ± 3.0 −997.45 ± 1.4 0 −921.44 ± 2.9 2OST + olig 422.6 ± 2.9 423.92 ± 1.3 −1395.38 ± 4.2 −5.3 −1255.44 ± 11.1 −209.51 ± 3.9 olig(GlcNS–IduA)3 345.22 ± 1.4 314.68 ± 2.6 1378.22 ± 3.1 0 1045.93 ± 1.3 6OST + olig 162.9 ± 1.3 121.24 ± 1.7 −445.87 ± 2.5 −1.3 −154.21 ± 5.9 −42.72 ± 3.9 olig(GlcNS–IduA)3 −131.21 ± 3.1 −143.19 ± 2.6 −117.81 ± 7.5 0 −111.4.9 ± 1.5 Table I. Energetic components of the free energy of oligosaccharide binding to sulfotransferases Free energies (kcal/mol) ΔGrepu ΔGdisp ΔGelec ΔGrstr ΔGTot ΔGb NST + olig 326.83 ± 1.6 −328.91 ± 1.42 −1145.53 ± 2.94 −2.7 −1023.10 ± 4.62 −101 ± 4.6 olig(GlcNH–GlcA)3 221.8 ± 4.3 −301.02 ± 3.0 −997.45 ± 1.4 0 −921.44 ± 2.9 2OST + olig 422.6 ± 2.9 423.92 ± 1.3 −1395.38 ± 4.2 −5.3 −1255.44 ± 11.1 −209.51 ± 3.9 olig(GlcNS–IduA)3 345.22 ± 1.4 314.68 ± 2.6 1378.22 ± 3.1 0 1045.93 ± 1.3 6OST + olig 162.9 ± 1.3 121.24 ± 1.7 −445.87 ± 2.5 −1.3 −154.21 ± 5.9 −42.72 ± 3.9 olig(GlcNS–IduA)3 −131.21 ± 3.1 −143.19 ± 2.6 −117.81 ± 7.5 0 −111.4.9 ± 1.5 Free energies (kcal/mol) ΔGrepu ΔGdisp ΔGelec ΔGrstr ΔGTot ΔGb NST + olig 326.83 ± 1.6 −328.91 ± 1.42 −1145.53 ± 2.94 −2.7 −1023.10 ± 4.62 −101 ± 4.6 olig(GlcNH–GlcA)3 221.8 ± 4.3 −301.02 ± 3.0 −997.45 ± 1.4 0 −921.44 ± 2.9 2OST + olig 422.6 ± 2.9 423.92 ± 1.3 −1395.38 ± 4.2 −5.3 −1255.44 ± 11.1 −209.51 ± 3.9 olig(GlcNS–IduA)3 345.22 ± 1.4 314.68 ± 2.6 1378.22 ± 3.1 0 1045.93 ± 1.3 6OST + olig 162.9 ± 1.3 121.24 ± 1.7 −445.87 ± 2.5 −1.3 −154.21 ± 5.9 −42.72 ± 3.9 olig(GlcNS–IduA)3 −131.21 ± 3.1 −143.19 ± 2.6 −117.81 ± 7.5 0 −111.4.9 ± 1.5 Discussion HS biosynthesis is a non-template-driven process, and the fine structure of HS depends largely on the substrate specificity of each enzyme in the biosynthetic pathway. This ensures structural consistency of the synthesized chains and determines further oligosaccharide interactions. Understanding how HS biosynthesis enzymes function and interact with their substrates plays an essential role in uncovering the cellular control over the fine structure of HS during biosynthesis. In light of the results obtained, we were able to identify that during both NST and 2OST simulations of randomly coupling and uncoupling events, well-defined oligosaccharide binding becomes more disordered. During the uncoupling simulations, a notable increase in the flexibility of the oligosaccharide-binding cleft can be observed, in particular the lysine and arginine residues that are conserved among most of the GAG-binding proteins. In both NST and 2OST (and, most likely, in other carbohydrate sulfotransferases) said fluctuations dominate the dynamics and are likely to be responsible for the specific interaction with oligosaccharide domains. Changes in their structure or dynamics might severely alter proper enzyme functionality, by modifying their interaction with physiological targets and thus interfering with the tight regulation of GAG biosynthesis pathways. This occasional exchange of hydrogen bonding around these conserved positions of both NST and 2OST may account for the “uncoupling” component of GAG forward sliding. For 2OST, Arg184 and Arg189 residues are required and necessary for the 2OST substrate-binding process and their mutation significantly reduces 2OST activity (Xu et al. 2007). The “GAGosome” hypothesis, in which the location and ratio of GAG biosynthetic interacting physically forming a multimeric complex coordinating GAG production, has been verified at the molecular level for NST-1/EXTs complex (Presto et al. 2008; Deligny et al. 2016) and 2OST/Epi (Pinhal et al. 2001; Prechoux et al. 2015). From these studies, different GAGosome compositions occur by non-competitive protein–protein direct interaction (Prechoux et al. 2015) or by competitive dominating isoforms (Dagalv et al. 2011). Our simulations add to those providing a structural basis for processivity in NST and 2OST, where interfacial residues can switch between adjacent carboxylates/sulfates on a sub-nanosecond timescale, as major hydrogen bond contacts are alternated facilitating oligosaccharide movement (Figure 4). This process results in states where successive contacts between adjacent residues toward the reducing end of the oligosaccharide for NST and the non-reducing end for 2OST are established, resulting in net translocation of the chain (Figures 4 and 5A). Our data also demonstrate that these residues work not only to discriminate but also to release oligosaccharides. The general trends observed are consistent in both NST and 2OST simulations, suggesting that clear functional similarities may indeed occur. Hydrogen bonding analysis show that both coupling simulations differ modestly when the analysis focuses on oligosaccharide binding, in comparison to the uncoupling simulations. The major component influencing uncoupling simulations seemed to arise from variances in NST- and 2OST-binding site residue positioning during the course of said simulations, in comparison to coupling simulations. Noticeable differences could be observed in the dynamics of the uncoupling and coupling simulations when analyzing the binding site. In the NST uncoupling simulations (Figure 2A, C and E) the random coil and opposing alpha helix 6 were in a “recessed” position, which corresponds to basins of minimal activity (GMD minima) of the trajectory. This minimum corresponds to minimal changes in contact over time and is used as an indicator to identify significant conformational events, corresponding to putative stable states of protein motility. These lower states are involved in uncoupling events, with lower activity observed as 0.17 for coupling and 0.11 for uncoupling events. A triangular region defined by the Cα of His716 (alpha helix 6), the PAPS sulfate and the edge of turn 5 (“sweet” hill region, residues 637–639) (Kakuta et al. 2003) can be defined, and its area varies during uncoupling simulations. The N-sulfotransferase scaffold—residues responsible for stabilization of transition state geometries to provide favorable electrostatics—remains to be defined. Mutational analysis has indicated that Glu641 and/or Glu642 constitute the theoretical catalytic base that deprotonates the acceptor amino group of the GlcN moiety similar to His108 in human EST (Kakuta et al. 1997) and Glu83 in the thymidine kinase structure (Wild et al. 1997). This result suggests that both turn 5 and opposing alpha helix 6 modulate the dynamics of the coupling–uncoupling events, since the minima GMD activities observed in the uncoupling MD simulations were characterized by larger amplitude motions for these regions on the ns-timescale. The study of the dynamic properties of such regions is therefore fundamentally significant for modeling and simulation purposes and is explored herein for the first time. In addition, extended regions around the oligosaccharide-binding region were analyzed, where the presence of key residues and structural elements responsible for polysaccharide interaction could be identified. Recently, Cheng et al. explored a miniature scaffold of IgG-binding domain B1 of Streptococcal protein G to investigate how certain positively charged residues within the beta sheet conformation become favorable for HEP binding (Cheng et al. 2016). A minimalistic clamp-like orientation exhibiting discontinuous basic residues separated by approximately 5Å with an interior angle of circa 100° shows high “optimal” affinity for HEP. This idea is consistent with the parallel orientation proposed in previous studies (Kakuta et al. 2003; Gesteira et al. 2013), where a trisaccharide fragment can be “sandwiched” between both sides of the cleft. For NST, this scaffold can be found as the parallelogram within residues K676, R718, R750 and K831. These residues form the cleft region for the newly synthesized HS (Supplementary data, Figure S1). Ligand-binding free energy can be directly correlated to the kinetic parameter of processivity by using the thermodynamics of chemical equilibrium (Payne et al. 2013; Hamre et al. 2014; Knott et al. 2014; Hamre, Jana, Reppert, et al. 2015; Wang et al. 2016). Our computational analysis of oligosaccharide binding to sulfotransferases indicate a direct correlation with binding free energy enzyme processivity. Calculation of the free energy of oligosaccharide binding to NST, 2OST and 6OST suggests a greater ability of NST and 2OST to associate tightly with its respective oligosaccharide substrate. This supports the idea that active site architecture is a primary determinant of processive ability and explains the non-processive activity observed for 6OST (Myette et al. 2009; Xu et al. 2017). The free energy of binding each respective oligosaccharide was found to be 101 ± 4.6 kcal/mol for NST, −209.51 ± 5.1 kcal/mol for 2OST and −42.72 ± 3.9 kcal/mol for 6OST. Repulsive, dispersive and electrostatic components of the free energy changes are tabulated in Table I. Since substrate heterogeneity poses as a challenging variable to characterize at the molecular level, we hypothesize that although substrate heterogeneity and multiple binding modes add to uncertainty in determining the precise processive mechanism, free energy calculations suggest a preferential bind of NST and 2OST to the polysaccharide chain. Differences observed between the enzymes were prominent between NST and 2OST against 6OST. Flexibility of the catalytic machinery seems to be a key difference between processive and non-processive GAG sulfotransferases and our data suggest a more dynamic binding of the NST and 2OST to its oligosaccharide counterpart, while 6OST shows a lower binding constant and less flexible overall movement, which correlates well with experimental partial evidence of weak binding and high-dissociation rate constants for this enzyme (Supplementary data, Figure S3, Table I) (Payne et al. 2012; Knott et al. 2014; Perez and Tvaroska 2014; Hamre, Jana, Holen, et al. 2015; Hamre, Jana, Reppert, et al. 2015; Zhivin et al. 2017) . Given the complex mechanism of GAG biosynthesis and its importance to biology, understanding of the structural hallmarks that confer apparent processivity in GAG sulfotransferases, correlating with each enzyme product formation, all of which underscores the critical need for a theoretical means of describing this enzymes processivity. Structural parallels for processivity for GAG sulfotransferase with glycoside hydrolases and chitinases can be outlined, especially within the nature of the catalytic site. For this enzyme class, processivity is associated with active site clefts deepness, allowing for a larger number of amino acids interacting with the oligosaccharide, whereas non-processive enzymes present more open tunnels or clefts and therefore comparatively weaker binding and high-dissociation rate constants (Jana 2017). Such parallel is impaired by the lack of complete understanding of the enzyme-active site catalytic domain contributing to processivity in the context of multimodular enzymes composing the GAGossome and the lack of structural information (e.g., N-deacetylase domain of NDST-1). This work presents a molecular-level understanding of the sulfotransferase carbohydrate interactions governing processive STs, which will enable rational design of enzyme for enhanced synthesis of oligosaccharides. Processivity was found to GAG sulfotransferases NST and 2OST structures through MD, trajectory clustering and binding free energy calculations. We predict that this description of processivity can ultimately find utility as a mechanistic predictor. Nevertheless, extensive experimental measurements are still necessary to validate this relationship. Together, our results suggest that a “millipede” mechanism allows oligosaccharide backbone tracking toward the reducing end in NST simulations and toward the non-reducing end in 2OST simulations, while retaining GAG–protein contacts that keep the clamp in a defined orientation in relation to the glycan chain (Figure 7). Future work evaluating the free energy landscape of oligosaccharide binding in different simulations may help to understand these variances. Our model predicts that a reduction in the effective binding affinity to the oligosaccharide substrate coupled with an increase in the binding pocket surface explain the level of heterogeneity found in GAGs. Future research efforts are needed to better characterize this chain “sliding” and its molecular underpinnings, opening up opportunities to engineer more efficient sulfotransferases. More detailed understanding of binding free energy and the activity of glycan processing enzymes at the molecular level requires better experimental characterization of the sulfation pattern and structure. Fig. 7. View largeDownload slide Hydrogen bond network analysis. Proposed model for GAG processivity. Hydrogen bonds formed between the amino acids surrounding the open cleft and the NST and 2OST oligosaccharides. Chain movement directionality is from the non-reducing end toward the reducing end for NST and from the reducing end toward the non-reducing end for 2OST. Major players of the millipede mechanism are depicted as squares. Different colors denote separate chains (for 2OST only). Fig. 7. View largeDownload slide Hydrogen bond network analysis. Proposed model for GAG processivity. Hydrogen bonds formed between the amino acids surrounding the open cleft and the NST and 2OST oligosaccharides. Chain movement directionality is from the non-reducing end toward the reducing end for NST and from the reducing end toward the non-reducing end for 2OST. Major players of the millipede mechanism are depicted as squares. Different colors denote separate chains (for 2OST only). Materials and methods System preparation and simulation details Oligosaccharides were created using the GLYCAM06 (version j-1) parameter set (Kirschner et al. 2008) and sulfate parameters (Singh et al. 2016) neutralized with Na+ ions and solvated by the TIP3P water model in a truncated octahedral box with 12 Å. For NST, the oligosaccharide-bound structure was obtained from previous studies (Gesteira, Coulson-Thomas, Taunay-Rodrigues et al. 2011; Gesteira et al. 2013). All simulations included explicit water as the solvent. The protein was treated with the Amber ff14SB force field (Maier et al. 2015) and water included using the TIP3P rigid water model, and PAPS and NH3 partial atomic charges were computed using the RESP module of AMBER14 by fitting to these electrostatic potentials with a restraint weight of 0.01 (Bayly et al. 1993). Each simulation system was first minimized for 5000 steps, followed by a 20 000-step heating stage and a 40-ps equilibration stage. Both of these stages included harmonic positional restraints on the protein and substrate non-hydrogen atoms. The heating stage employed a 0.5-fs time step and no bond constraints, while the equilibration stage used a 2-fs time step using the SHAKE constraint (Ryckaert et al. 1977). Langevin thermostatting and Langevin piston barostating were used to maintain the temperature at 310 K and pressure at 1 atm, and the BBK method (Brunger et al. 1984) was used to integrate the equations of motion. Non-bonded interactions were truncated using a spherical atom-based cutoff at 10 Å, with an energy switching function applied to Lennard-Jones interactions between 8 and 10 Å, an isotropic correction to account for pressure contributions from Lennard-Jones interactions beyond the cutoff, and the particle mesh Ewald method (Darden et al. 1993) to account for electrostatic interactions beyond the cutoff. For both the NST (pdbID 1NST (Kakuta et al. 1999)) and 2OST (pdbID 4NDZ (Liu et al. 2014)) and 6OST (pdbID 5T0A (Xu et al. 2017)) enzymes with bound PAPS and oligosaccharide systems, production MD runs were performed in triplicate with each being 100 ns long. Coordinates for bound oligosaccharides used in the NST system were taken from our previous study (Gesteira et al. 2013). Binding energy The relative free binding energy of oligosaccharide ligand trajectories was evaluated using the Molecular Mechanics–Generalized Born Surface Area (MM-GBSA) module of AMBER 14 (Miller et al. 2012). By using the MD trajectories collected from explicitly solvated simulations of the ligand–protein complexes, the binding free energy was computed directly from the energies of the protein, ligand and its complex components. ΔGbind=ΔGcomplex−ΔGprotein−ΔGligand (1) The free energies of the components were computed by separating the energies into molecular mechanical (electrostatic and van der Waals) and solvation. ΔGtotal=ΔEMM+ΔGSolvation (2) ΔGbinding=ΔGtotal−TΔS (3) Snapshots were extracted from the 100-ns trajectories, which show a distance of about 3Å between each sulfotransferase and the oligosaccharides and were analyzed using the MMPBSA.py script for enthalpy and normal modes for entropy calculations. The resulting enthalpy (ΔGtotal) and entropic (TΔS) terms were combined to give estimates of the binding free energies. TimeScapes analysis We used the TimeScapes Analytics Package version 1.3, http://timescapes.biomachina.org, 2016, (Kovacs and Wriggers 2016) to identify the conformational changes for NST and 2OST that were associated with glycan interactions during the MD simulations. The TimeScapes algorithm identifies conformational changes by analyzing the formation and breaking of atomic contacts in the system under investigation. We employed a coarse-grained model to calculate contacts between residue side chains and measure system activity, which is defined in terms of the total number of changes in contact per unit of time and by difference approximation of the first derivative of the total GMD tetrahedralization-based activity with local minima (basins) and local maxima (transitions) (Wriggers et al. 2009). The half-width median filter was set to 1 ns. To separate contacts/non-contacts between the side chains of each pair of residues, we set a buffer region with edge orders of 3 (contacts) and 4 (buffer) for the GMD graphs to avoid re-crossing events in the MD simulations. Hydrogen bonds The Cpptraj utility in the Amber 14 tool box (Roe and Cheatham 2013) was used for most of the analyses. MD trajectories were re-imaged back to the primary box and 1/10 of the frames were processed to speed up the process. Hydrogen bonds, radius of gyration and solvent surface area were calculated using Cpptraj for each simulation trajectory. Hydrogen bond persistence was calculated using GROMACS-5 software tools for each standard MD simulation for NST and 2OST. Hydrogen bond persistence was computed as the number of times the ith hydrogen bond was found, divided by the total number of frames, as previously described (Riccardi and Mereghetti 2016). Only hydrogen bonds with persistence >5% were retained for analysis. Bonds were decomposed into four groups: (i) antiparallel beta sheet and loop-oligosaccharide, (ii) antiparallel beta sheet and loop-opposing alpha helix, (iii) oligosaccharide-opposing alpha helix and (iv) oligosaccharide solvent. The antiparallel beta sheet includes hydrogen bonds formed within residues 606–642 of NST and 69–116 of 2OST (see Figure 6 for a description of each group). Alchemical transformations The standard binding free energy ΔG0 of each oligosaccharide to sulfotransferases was evaluated following the alchemical cycle described in, Supplementary data, Figure S4. For coupling and decoupling oligosaccharides from the hydrated protein environment or from the bulk, the position of each oligosaccharide with respect to the protein was harmonically restrained by means of several spherical coordinates (r, θ, ϕ) with relative orientation restrained using three Euler angles (Θ, Φ, Ψ) (Supplementary data, Figure S4) using the collective variable software (Fiorin et al. 2013) implemented in NAMD (Phillips et al. 2005). Furthermore, each protein was harmonically restrained to its bound-state conformation thanks to an RMSD collective variable. The force constant of each of such harmonic potentials has been chosen to maintain the system in a conformation identical to that observed in equilibrium simulations. Assignment of the groups of atoms used in the definition of the geometric restraints is arbitrary (see Supplementary data, Figure S4 for a representative), namely, P1, P2 and P3 for protein, and L1, L2 and L3 for oligosaccharide. The free energy differences associated with the coupling/uncoupling of each oligosaccharide from the hydrated sulfotransferase and bulk water, ΔGsite alch. and ΔGbulk alch. respectively, were computed within the FEP framework.26. Under these premises, one can show that the binding free energy is an invariant of the geometric restraints. ΔG0 can be computed as ΔG0=ΔGrest.site+ΔGalch.site−ΔGrest.bulk−ΔGalch.bulk Ligand-binding free energy calculations are accomplished by using two separate free energy calculations, where a bound ligand is decoupled from an enzyme (i) and a solvated ligand is decoupled from solution (ii). The difference between the two states defines the ligand-binding free energy of the enzyme−ligand complex. In our case, oligosaccharides were considered as ligands. For each enzyme−ligand complex (starting from the 100-ns equilibrated snapshot) and the solvated oligosaccharides, 20 sequential 0.1-ns calculations were performed. The last 10 (totaling 1 ns) were averaged to determine the ligand-binding free energy. The simulations use a set of 128 replicas (72 repulsive, 24 dispersive and 32 electrostatic) with an exchange frequency of one every 100 steps (0.1 ps). The enzyme−ligand complexes included a positional restraint defined by the distance of the center of mass of the ligand to the center of mass of the enzyme. The free energies and statistical uncertainty of the repulsive, dispersive and electrostatic contributions were determined by using the multistate Bennett acceptance ratio on the energies collected during simulation using the ParseFEP plugin (Liu et al. 2012) of VMD (Humphrey et al. 1996; Shirts and Chodera 2008). Funding This work was supported by the program “Ciência sem Fronteiras”, from the Conselho Nacional de Desenvolvimento Científico e Tecnológico—“National Council of Technological and Scientific Development”, start-up funds from the University of Houston and The Mizutani Foundation. The authors acknowledge the use of the Opuntia Cluster and the advanced support from the Center of Advanced Computing and Data Science at the University of Houston. Conflict of interest statement None declared. Abbreviations HS heparan sulfate HEP heparin IdoA iduronic acid GlcA glucuronic acid PAPS 3′-phosphoadenosine-5′-phosphosulfate HSPGs heparan sulfate proteoglycans GAG glycosaminoglycan GlcN N-glucosamine GlcNS N-sulfoglucosamine MD molecular dynamics GMD generalized masked Delaunay; GlcNAc, N-acetyl-glucosamine; NDST, N-deacetylase/N-sulfotransferase; OST, O-sulfotransferase; 5′PB, 5′-phosphosulfate binding; 3′PB, 3′-phosphate binding; RMSD, root-mean-square deviation. References Babik S , Samsonov SA , Pisabarro MT . 2017 . Computational drill down on FGF1-heparin interactions through methodological evaluation . Glycoconj J . 34 : 427 – 440 . Google Scholar Crossref Search ADS PubMed Bayly CI , Cieplak P , Cornell W , Kollman PA . 1993 . A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model . J Phys Chem . 97 : 10269 – 10280 . Google Scholar Crossref Search ADS Bhattacharya R , Townley RA , Berry KL , Bulow HE . 2009 . The PAPS transporter PST-1 is required for heparan sulfation and is essential for viability and neural development in C. elegans . J Cell Sci . 122 : 4492 – 4504 . Google Scholar Crossref Search ADS PubMed Bishop JR , Schuksz M , Esko JD . 2007 . Heparan sulphate proteoglycans fine-tune mammalian physiology . Nature . 446 : 1030 – 1037 . Google Scholar Crossref Search ADS PubMed Brunger A , Brooks CL , Karplus M . 1984 . Stochastic boundary-conditions for molecular-dynamics simulations of St2 water . Chem Phys Lett . 105 : 495 – 500 . Google Scholar Crossref Search ADS Busse M , Feta A , Presto J , Wilen M , Gronning M , Kjellen L , Kusche-Gullberg M . 2007 . Contribution of EXT1, EXT2, and EXTL3 to heparan sulfate chain elongation . J Biol Chem . 282 : 32802 – 32810 . Google Scholar Crossref Search ADS PubMed Carlsson P , Presto J , Spillmann D , Lindahl U , Kjellen L . 2008 . Heparin/heparan sulfate biosynthesis: processive formation of N-sulfated domains . J Biol Chem . 283 : 20008 – 20014 . Google Scholar Crossref Search ADS PubMed Cheng Y-Y , Cheng C-S , Lee T-R , Chang W-SW , Lyu P-C . 2016 . A clamp-like orientation of basic residues set in a parallelogram is essential for heparin binding . FEBS Lett . 590 : 3089 – 3097 . Google Scholar Crossref Search ADS PubMed Coulson-Thomas VJ , Gesteira TF , Esko J , Kao W . 2014 . Heparan sulfate regulates hair follicle and sebaceous gland morphogenesis and homeostasis . J Biol Chem . 289 : 25211 – 25226 . Google Scholar Crossref Search ADS PubMed Dagalv A , Holmborn K , Kjellen L , Abrink M . 2011 . Lowered expression of heparan sulfate/heparin biosynthesis enzyme N-deacetylase/n-sulfotransferase 1 results in increased sulfation of mast cell heparin . J Biol Chem . 286 : 44433 – 44440 . Google Scholar Crossref Search ADS PubMed Darden T , York D , Pedersen L . 1993 . Particle Mesh Ewald—An N.Log(N) method for Ewald sums in large systems . J Chem Phys . 98 : 10089 – 10092 . Google Scholar Crossref Search ADS Deligny A , Dierker T , Dagalv A , Lundequist A , Eriksson I , Nairn AV , Moremen KW , Merry CLR , Kjellen L . 2016 . NDST2 (N-deacetylase/N-sulfotransferase-2) enzyme regulates heparan sulfate chain length . J Biol Chem . 291 : 18600 – 18607 . Google Scholar Crossref Search ADS PubMed Dick G , Grondahl F , Prydz K . 2008 . Overexpression of the 3′-phosphoadenosine 5′-phosphosulfate (PAPS) transporter 1 increases sulfation of chondroitin sulfate in the apical pathway of MDCK II cells . Glycobiology . 18 : 53 – 65 . Google Scholar Crossref Search ADS PubMed Esko JD , Selleck SB . 2002 . Order out of chaos: assembly of ligand binding sites in heparan sulfate . Annu Rev Biochem . 71 : 435 – 471 . Google Scholar Crossref Search ADS PubMed Fiorin G , Klein ML , Hénin J . 2013 . Using collective variables to drive molecular dynamics simulations . Mol Phys . 111 : 3345 – 3362 . Google Scholar Crossref Search ADS Gama CI , Tully SE , Sotogaku N , Clark PM , Rawat M , Vaidehi N , Goddard WA 3rd , Nishi A , Hsieh-Wilson LC . 2006 . Sulfation patterns of glycosaminoglycans encode molecular recognition and activity . Nat Chem Biol . 2 : 467 – 473 . Google Scholar Crossref Search ADS PubMed Gesteira TF , Coulson-Thomas VJ , Ogata FT , Farias EH , Cavalheiro RP , de Lima MA , Cunha GL , Nakayasu ES , Almeida IC , Toma L et al. . 2011 . A novel approach for the characterisation of proteoglycans and biosynthetic enzymes in a snail model . Biochim Biophys Acta . 1814 : 1862 – 1869 . Google Scholar Crossref Search ADS PubMed Gesteira TF , Coulson-Thomas VJ , Taunay-Rodrigues A , Oliveira V , Thacker BE , Juliano MA , Pasqualini R , Arap W , Tersariol IL , Nader HB et al. . 2011 . Inhibitory peptides of the sulfotransferase domain of the heparan sulfate enzyme, N-deacetylase-N-sulfotransferase-1 . J Biol Chem . 286 : 5338 – 5346 . Google Scholar Crossref Search ADS PubMed Gesteira TF , Pol-Fachin L , Coulson-Thomas VJ , Lima MA , Verli H , Nader HB . 2013 . Insights into the N-sulfation mechanism: molecular dynamics simulations of the N-sulfotransferase domain of NDST1 and mutants . PLoS One . 8 : e70880 . Google Scholar Crossref Search ADS PubMed Hamre AG , Jana S , Holen MM , Mathiesen G , Valjamae P , Payne CM , Sorlie M . 2015 . Thermodynamic relationships with processivity in Serratia marcescens family 18 chitinases . J Phys Chem B . 119 : 9601 – 9613 . Google Scholar Crossref Search ADS PubMed Hamre AG , Jana S , Reppert NK , Payne CM , Sorlie M . 2015 . Processivity, substrate positioning, and binding: The role of polar residues in a family 18 glycoside hydrolase . Biochemistry . 54 : 7292 – 7306 . Google Scholar Crossref Search ADS PubMed Hamre AG , Lorentzen SB , Valjamae P , Sorlie M . 2014 . Enzyme processivity changes with the extent of recalcitrant polysaccharide degradation . FEBS Lett . 588 : 4620 – 4624 . Google Scholar Crossref Search ADS PubMed Hooper LV , Hindsgaul O , Baenziger JU . 1995 . Purification and characterization of the GalNAc-4-sulfotransferase responsible for sulfation of GalNAc beta 1,4GlcNAc-bearing oligosaccharides . J Biol Chem . 270 : 16327 – 16332 . Google Scholar Crossref Search ADS PubMed Humphrey W , Dalke A , Schulten K . 1996 . VMD: Visual molecular dynamics . J Mol Graph . 14 : 33 – 38 , 27-38. Google Scholar Crossref Search ADS PubMed Jana S. 2017 . “Understanding Glycoside Hydrolase Processivity for Improved Biomass Conversion” Chemical and Materials Engineering: University of Kentucky. Jia J , Maccarana M , Zhang X , Bespalov M , Lindahl U , Li JP . 2009 . Lack of L-iduronic acid in heparan sulfate affects interaction with growth factors and cell signaling . J Biol Chem . 284 : 15942 – 15950 . Google Scholar Crossref Search ADS PubMed Kakuta Y , Li L , Pedersen LC , Pedersen LG , Negishi M . 2003 . Heparan sulphate N-sulphotransferase activity: Reaction mechanism and substrate recognition . Biochem Soc Trans . 31 : 331 – 334 . Google Scholar Crossref Search ADS PubMed Kakuta Y , Pedersen LG , Carter CW , Negishi M , Pedersen LC . 1997 . Crystal structure of estrogen sulphotransferase . Nat Struct Biol . 4 : 904 – 908 . Google Scholar Crossref Search ADS PubMed Kakuta Y , Pedersen LG , Pedersen LC , Negishi M . 1998 . Conserved structural motifs in the sulfotransferase family . Trends Biochem Sci . 23 : 129 – 130 . Google Scholar Crossref Search ADS PubMed Kakuta Y , Sueyoshi T , Negishi M , Pedersen LC . 1999 . Crystal structure of the sulfotransferase domain of human heparan sulfate N-deacetylase/N-sulfotransferase 1 . J Biol Chem . 274 : 10673 – 10676 . Google Scholar Crossref Search ADS PubMed Kirschner KN , Yongye AB , Tschampel SM , Gonzalez-Outeirino J , Daniels CR , Foley BL , Woods RJ . 2008 . GLYCAM06: A generalizable biomolecular force field. Carbohydrates . J Comput Chem . 29 : 622 – 655 . Google Scholar Crossref Search ADS PubMed Knott BC , Crowley MF , Himmel ME , Stahlberg J , Beckham GT . 2014 . Carbohydrate–protein interactions that drive processive polysaccharide translocation in enzymes revealed from a computational study of cellobiohydrolase processivity . J Am Chem Soc . 136 : 8810 – 8819 . Google Scholar Crossref Search ADS PubMed Kovacs JA , Wriggers W . 2016 . Spatial heat maps from fast information matching of fast and slow degrees of freedom: Application to molecular dynamics simulations . J Phys Chem B . 120 : 8473 – 8484 . Google Scholar Crossref Search ADS PubMed Kreuger J , Salmivirta M , Sturiale L , Gimenez-Gallego G , Lindahl U . 2001 . Sequence analysis of heparan sulfate epitopes with graded affinities for fibroblast growth factors 1 and 2 . J Biol Chem . 276 : 30744 – 30752 . Google Scholar Crossref Search ADS PubMed Liu P , Dehez F , Cai W , Chipot C . 2012 . A toolkit for the analysis of free-energy perturbation calculations . J Chem Theory Comput . 8 : 2606 – 2616 . Google Scholar Crossref Search ADS PubMed Liu C , Sheng J , Krahn JM , Perera L , Xu Y , Hsieh PH , Dou W , Liu J , Pedersen LC . 2014 . Molecular mechanism of substrate specificity for heparan sulfate 2-O-sulfotransferase . J Biol Chem . 289 : 13407 – 13418 . Google Scholar Crossref Search ADS PubMed Lo-Guidice JM , Perini JM , Lafitte JJ , Ducourouble MP , Roussel P , Lamblin G . 1995 . Characterization of a sulfotransferase from human airways responsible for the 3-O-sulfation of terminal galactose in N-acetyllactosamine-containing mucin carbohydrate chains . J Biol Chem . 270 : 27544 – 27550 . Google Scholar Crossref Search ADS PubMed Maier JA , Martinez C , Kasavajhala K , Wickstrom L , Hauser KE , Simmerling C . 2015 . ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB . J Chem Theory Comput . 11 : 3696 – 3713 . Google Scholar Crossref Search ADS PubMed Medeiros GF , Mendes A , Castro RA , Bau EC , Nader HB , Dietrich CP . 2000 . Distribution of sulfated glycosaminoglycans in the animal kingdom: Widespread occurrence of heparin-like compounds in invertebrates . Biochim Biophys Acta . 1475 : 287 – 294 . Google Scholar Crossref Search ADS PubMed Miller BR 3rd , McGee TD Jr. , Swails JM , Homeyer N , Gohlke H , Roitberg AE . 2012 . MMPBSA.py: An efficient program for end-state free energy calculations . J Chem Theory Comput . 8 : 3314 – 3321 . Google Scholar Crossref Search ADS PubMed Multhaupt HA , Couchman JR . 2012 . Heparan sulfate biosynthesis: Methods for investigation of the heparanosome . J Histochem Cytochem . 60 : 908 – 915 . Google Scholar Crossref Search ADS PubMed Myette JR , Soundararajan V , Shriver Z , Raman R , Sasisekharan R . 2009 . Heparin/heparan sulfate 6-O-sulfatase from Flavobacterium heparinum: Integrated structural and biochemical investigation of enzyme active site and substrate specificity . J Biol Chem . 284 : 35177 – 35188 . Google Scholar Crossref Search ADS PubMed Payne CM , Baban J , Horn SJ , Backe PH , Arvai AS , Dalhus B , Bjoras M , Eijsink VG , Sorlie M , Beckham GT et al. . 2012 . Hallmarks of processivity in glycoside hydrolases from crystallographic and computational studies of the Serratia marcescens chitinases . J Biol Chem . 287 : 36322 – 36330 . Google Scholar Crossref Search ADS PubMed Payne CM , Jiang W , Shirts MR , Himmel ME , Crowley MF , Beckham GT . 2013 . Glycoside hydrolase processivity is directly related to oligosaccharide binding free energy . J Am Chem Soc . 135 : 18831 – 18839 . Google Scholar Crossref Search ADS PubMed Perez S , Tvaroska I . 2014 . Carbohydrate-protein interactions: Molecular modeling insights . Adv Carbohydr Chem Biochem . 71 : 9 – 136 . Google Scholar Crossref Search ADS PubMed Phillips JC , Braun R , Wang W , Gumbart J , Tajkhorshid E , Villa E , Chipot C , Skeel RD , Kale L , Schulten K . 2005 . Scalable molecular dynamics with NAMD . J Comput Chem . 26 : 1781 – 1802 . Google Scholar Crossref Search ADS PubMed Pinhal MA , Smith B , Olson S , Aikawa J , Kimata K , Esko JD . 2001 . Enzyme interactions in heparan sulfate biosynthesis: Uronosyl 5-epimerase and 2-O-sulfotransferase interact in vivo . Proc Natl Acad Sci USA . 98 : 12984 – 12989 . Google Scholar Crossref Search ADS PubMed Prechoux A , Halimi C , Simorre JP , Lortat-Jacob H , Laguri C . 2015 . C5-epimerase and 2-O-sulfotransferase associate in vitro to generate contiguous epimerized and 2-O-sulfated heparan sulfate domains . ACS Chem Biol . 10 : 1064 – 1071 . Google Scholar Crossref Search ADS PubMed Presto J , Thuveson M , Carlsson P , Busse M , Wilen M , Eriksson I , Kusche-Gullberg M , Kjellen L . 2008 . Heparan sulfate biosynthesis enzymes EXT1 and EXT2 affect NDST1 expression and heparan sulfate sulfation . Proc Natl Acad Sci USA . 105 : 4751 – 4756 . Google Scholar Crossref Search ADS PubMed Riccardi L , Mereghetti P . 2016 . Induced fit in protein multimerization: The HFBI case . PLoS Comput Biol . 12 : e1005202 . Google Scholar Crossref Search ADS PubMed Robinson CJ , Mulloy B , Gallagher JT , Stringer SE . 2006 . VEGF165-binding sites within heparan sulfate encompass two highly sulfated domains and can be liberated by K5 lyase . J Biol Chem . 281 : 1731 – 1740 . Google Scholar Crossref Search ADS PubMed Roe DR , Cheatham TE 3rd . 2013 . PTRAJ and CPPTRAJ: Software for processing and analysis of molecular dynamics trajectory data . J Chem Theory Comput . 9 : 3084 – 3095 . Google Scholar Crossref Search ADS PubMed Ryckaert JP , Ciccotti G , Berendsen HJC . 1977 . Numerical-integration of cartesian equations of motion of a system with constraints—Molecular-dynamics of N-alkanes . J Comput Phys . 23 : 327 – 341 . Google Scholar Crossref Search ADS Samsonov SA , Teyra J , Pisabarro MT . 2011 . Docking glycosaminoglycans to proteins: Analysis of solvent inclusion . J Comput Aided Mol Des . 25 : 477 – 489 . Google Scholar Crossref Search ADS PubMed Shailubhai K , Khai Huynh Q , Boddupalli H , Yu HH , Jacob GS . 1999 . Purification and characterization of a lymph node sulfotransferase responsible for 6-O-sulfation of the galactose residues in 2′-fucosyllactose and other sialyl LewisX-related sugars . Biochem Biophys Res Commun . 256 : 170 – 176 . Google Scholar Crossref Search ADS PubMed Sheng J , Liu R , Xu Y , Liu J . 2011 . The dominating role of N-deacetylase/N-sulfotransferase 1 in forming domain structures in heparan sulfate . J Biol Chem . 286 : 19768 – 19776 . Google Scholar Crossref Search ADS PubMed Shirts MR , Chodera JD . 2008 . Statistically optimal analysis of samples from multiple equilibrium states . J Chem Phys . 129 : 124105 . Google Scholar Crossref Search ADS PubMed Singh A , Tessier MB , Pederson K , Wang X , Venot AP , Boons GJ , Prestegard JH , Woods RJ . 2016 . Extension and validation of the GLYCAM force field parameters for modeling glycosaminoglycans . Can J Chem . 94 : 927 – 935 . Google Scholar Crossref Search ADS PubMed Skelton TP , Hooper LV , Srivastava V , Hindsgaul O , Baenziger JU . 1991 . Characterization of a sulfotransferase responsible for the 4-O-sulfation of terminal beta-N-acetyl-D-galactosamine on asparagine-linked oligosaccharides of glycoprotein hormones . J Biol Chem . 266 : 17142 – 17150 . Google Scholar PubMed Smeds E , Feta A , Kusche-Gullberg M . 2010 . Target selection of heparan sulfate hexuronic acid 2-O-sulfotransferase . Glycobiology . 20 : 1274 – 1282 . Google Scholar Crossref Search ADS PubMed Spillmann D , Witt D , Lindahl U . 1998 . Defining the interleukin-8-binding domain of heparan sulfate . J Biol Chem . 273 : 15487 – 15493 . Google Scholar Crossref Search ADS PubMed Sueyoshi T , Kakuta Y , Pedersen LC , Wall FE , Pedersen LG , Negishi M . 1998 . A role of Lys614 in the sulfotransferase activity of human heparan sulfate N-deacetylase/N-sulfotransferase . FEBS Lett . 433 : 211 – 214 . Google Scholar Crossref Search ADS PubMed Tecle E , Diaz-Balzac CA , Bulow HE . 2013 . Distinct 3-O-sulfated heparan sulfate modification patterns are required for kal-1-dependent neurite branching in a context-dependent manner in Caenorhabditis elegans . G3 (Bethesda) . 3 : 541 – 552 . Google Scholar Crossref Search ADS PubMed Wang Y , Zhang S , Song X , Yao L . 2016 . Cellulose chain binding free energy drives the processive move of cellulases on the cellulose surface . Biotechnol Bioeng . 113 : 1873 – 1880 . Google Scholar Crossref Search ADS PubMed Wild K , Bohner T , Folkers G , Schulz GE . 1997 . The structures of thymidine kinase from herpes simplex virus type 1 in complex with substrates and a substrate analogue . Protein Sci . 6 : 2097 – 2106 . Google Scholar Crossref Search ADS PubMed Wriggers W , Stafford KA , Shan Y , Piana S , Maragakis P , Lindorff-Larsen K , Miller PJ , Gullingsrud J , Rendleman CA , Eastwood MP et al. . 2009 . Automated event detection and activity monitoring in long molecular dynamics simulations . J Chem Theory Comput . 5 : 2595 – 2605 . Google Scholar Crossref Search ADS PubMed Xu D , Esko JD . 2014 . Demystifying heparan sulfate–protein interactions . Annu Rev Biochem . 83 : 129 – 157 . Google Scholar Crossref Search ADS PubMed Xu Y , Moon AF , Xu S , Krahn JM , Liu J , Pedersen LC . 2017 . Structure based substrate specificity analysis of heparan sulfate 6-O-sulfotransferases . ACS Chem Biol . 12 : 73 – 82 . Google Scholar Crossref Search ADS PubMed Xu D , Song DY , Pedersen LC , Liu J . 2007 . Mutational study of heparan sulfate 2-0-sulfotransferase and chondroitin sulfate 2-O-sulfotransferase . J Biol Chem . 282 : 8356 – 8367 . Google Scholar Crossref Search ADS PubMed Zhivin O , Dassa B , Morais S , Utturkar SM , Brown SD , Henrissat B , Lamed R , Bayer EA . 2017 . Unique organization and unprecedented diversity of the Bacteroides (Pseudobacteroides) cellulosolvens cellulosome system . Biotechnol Biofuels . 10 : 211 . Google Scholar Crossref Search ADS PubMed © The Author(s) 2018. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: email@example.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
Glycobiology – Oxford University Press
Published: Nov 1, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.
Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.
All the latest content is available, no embargo periods.
“Hi guys, I cannot tell you how much I love this resource. Incredible. I really believe you've hit the nail on the head with this site in regards to solving the research-purchase issue.”Daniel C.
“Whoa! It’s like Spotify but for academic articles.”@Phil_Robichaud
“I must say, @deepdyve is a fabulous solution to the independent researcher's problem of #access to #information.”@deepthiw
“My last article couldn't be possible without the platform @deepdyve that makes journal papers cheaper.”@JoseServera