A Simple Method for Discovering Druggable, Specific Glycosaminoglycan-Protein Systems. Elucidation of Key Principles from Heparin/Heparan Sulfate-Binding Proteins.

Glycosaminoglycans (GAGs) affect human physiology and pathology by modulating more than 500 proteins. GAG-protein interactions are generally assumed to be ionic and nonspecific, but specific interactions do exist. Here, we present a simple method to identify the GAG-binding site (GBS) on proteins that in turn helps predict high specific GAG-protein systems. Contrary to contemporary thinking, we found that the electrostatic potential at basic arginine and lysine residues neither identifies the GBS consistently, nor its specificity. GBSs are better identified by considering the potential at neutral hydrogen bond donors such as asparagine or glutamine sidechains. Our studies also reveal that an unusual constellation of ionic and non-ionic residues in the binding site leads to specificity. Nature engineers the local environment of Asn45 of antithrombin, Gln255 of 3-O-sulfotransferase 3, Gln163 and Asn167 of 3-O-sulfotransferase 1 and Asn27 of basic fibroblast growth factor in the respective GBSs to induce specificity. Such residues are distinct from other uncharged residues on the same protein structure in possessing a significantly higher electrostatic potential, resultant from the local topology. In contrast, uncharged residues on nonspecific GBSs such as thrombin and serum albumin possess a diffuse spread of electrostatic potential. Our findings also contradict the paradigm that GAG-binding sites are simply a collection of contiguous Arg/Lys residues. Our work demonstrates the basis for discovering specifically interacting and druggable GAG-protein systems based on the structure of protein alone, without requiring access to any structure-function relationship data.

Introduction

Sulfated glycosaminoglycans (GAGs), such as heparan sulfate (HS), are nature’s most enigmatic biopolymers. Although made from a linear combination of simple saccharide rings, they display a staggering range of primary sequence diversity that surpasses the range possible for equivalent chains of other biopolymers. The occurrence of HS in planaria [1], the 2nd most primitive species of the animal kingdom, in a variably sulfated form closely matching that in humans suggests these biopolymers play critical roles in multiple fundamental biological processes. It is now recognized that GAGs bind to hundreds of human proteins with implications in physiological as well as pathological processes such as hemostasis and thrombosis, wound repair and inflammation, neuronal growth and amyloidogenesis, angiogenesis and cancer, defense against microbes and infection [2,3]. Yet, precious little is understood about how these unique, linear polysaccharides recognize and modulate their targets.

Considerable effort has been expended in understanding GAG–protein interactions at an atomic level. The earliest attempt to rationally deduce GAG binding sites (GBSs) on proteins was that of Cardin and Weintraub, who identified ‘XBBXBX’ and ‘XBBBXXBX’ as GAG-recognition domains, where B and X represent basic and hydropathic residues [4]. Later, these linear, α-helical or β-strand-like segments were extended to include other secondary structural elements [5,6]. However, such linear elements imply divergent evolution of GBSs, whereas the large structural diversity of GAG-binding proteins (GBPs) suggests exactly the opposite. In fact, most GBPs do not follow the simplistic Cardin and Weintraub rules. Sophisticated computational tools are being developed GBS identification [7-19]. While these approaches are successful within their set limits, a key question that remains unaddressed to date is the specificity of GAG-protein interactions [3,20,21], even though efforts have been made recently to do so [22,23]. For example, why are certain GAG–protein systems, e.g., heparin–antithrombin, highly specific [24,25], while others, e.g., heparin–thrombin [24,26,27], essentially nonspecific? More importantly, can ‘specific’ GAG–protein systems be more reliably predicted to help advance chemical biology and drug discovery?

Predicting the specificity of GAG–protein interactions is extremely challenging because of their reliance on long-range and non-directional Coulombic forces of attraction. Since multiple Arg/Lys and sulfates are involved in these interactions, many of which are redundant, most GAG–protein systems are traditionally forsaken as nonspecific. Yet, growing evidence suggests that GAGs display considerable non-ionic binding energy (10–40%) in recognizing their targets [24,28,29], which may arise from short-range and directional forces, such as hydrogen-bond(s), that induce ‘specific’ recognition. However, the exact origin of specificity among the multitude of interaction loci has been difficult to pinpoint.

The scientific community has historically focused on site directed mutagenesis of residues such as arginine, lysine and histidine, often present within the GBS () and which will likely possess a positive charge under physiological conditions. However, GBSs also possess polar residues such as asparagine and glutamine [29]. Histidine is a special case because it may or may not be charged under physiological conditions (pKa ~ 6.9) and its protonation state during GAG-binding isn’t always clear, but Asn and Gln are certainly never positively charged. Therefore, the reason for presence of Asn and Gln in GAG-binding sites remains unclear, even though it is firmly established that they interact with GAGs in many systems (). We asked a fundamental question: can these uncharged, polar residues have something to do with specificity?

We have uncovered that uncharged residues, such as Asn or Gln, help identify the GBS and also segregate specific GAG-protein systems from nonspecific ones. Our computational results in conjunction with structural and biochemical results show that an unusual constellation of ionic and non-ionic residues constituting the GAG-binding site is necessary for high specificity interaction. These principles can help parse proteins, including those that follow the Cardin–Weintraub rule [4] as well as those that do not, based on their specificity of interaction with GAGs. Based on our findings, we propose that Asn/Gln reduce desolvation penalties during the mostly electrostatic GAG-protein interactions, in addition to hydrogen bonding to sulfates and carboxylates on the GAG. Our results challenge the traditional paradigm that GAG-binding sites on proteins are located in contiguous segments, such as helices and/or turns and that specificity arises from Arg and Lys residues. Our work is expected to fundamentally change the landscape of discovery of ‘highly specific’ GAG-binding sites on proteins, which should greatly assist with identification of druggable GAG-protein systems for designing novel GAGs and GAG mimetics as drugs.

Results

Electrostatic potential alone at electropositive residues does not always identify a GAG binding site

Traditionally, a cluster of electropositive residues has been assumed to form a GBS, suggesting that the electrostatic potential (G ) at Arg/Lys residues should unequivocally identify GBSs. To quantitatively assess whether this expectation is correct, we calculated G on Arg/Lys residues of representative heparin-binding proteins including antithrombin, thrombin, FGF2, HS2ST1, HS3ST1 and HS3ST3A1. The G are represented on 2DSE plots, which we devised as new tools for easy quantitative visualization of energy at any locus, e.g., nitrogen donor atom of an Arg or Lys, on the protein surface in two-dimensions (see and Supplementary Methods). Briefly, the position of a locus (such as an atom) on a protein surface is projected from 3D space onto a 2D plane to obtain a scatterplot. The area of each data point in this scatterplot is scaled in size to reflect the G manifested at the atom by the rest of the structure. Therefore, the larger spots represent locations most conducive to binding negatively charged entities such as GAGs. Our expectation, as explained above, was that Arg/Lys of the GBS would demonstrate the highest (most negative) G .

Converse to our expectations, reveals that not all GBS Arg/Lys () carry a high G . This is in direct contrast to common empirical assumptions. For example, G for Arg46, Lys114 and Arg129 of antithrombin were significantly higher than other basic residues present in the heparin-binding site such as Arg47, Lys125, Arg132, Lys133, Arg136 and Lys275 (). In fact, Lys114 and Arg129 are known to be crucial for heparin binding [25,31,32], but Arg46 is not. Likewise, Arg93 stands out as the only residue with high G for thrombin, surpassing two residues (Arg97 and Arg101) known to play key role in heparin binding, () [26,33]. Arg120 of FGF2 is only one of several, similarly contributing residues of the GBS [34,35] but demonstrates significantly higher G (). Similar results were obtained in calculations for HS3ST3A1, HS3ST1 and HS2ST1, three enzymes of the heparan sulfate biosynthesis pathway. Although Arg166, Arg274 and Arg325 of HS3ST3A1 display high G (), the latter two are not known to be part of the GBS [36]. Likewise, Arg223 of HS3ST1 shows strong G () but is not known to be part of the GBS [37]. Finally, Arg184 of HS2ST1 demonstrates high G as well () but is not known to be important for binding GAGs [35,38]. This implies that not all Arg/Lys residues on a GAG-binding protein’s surface carry high electrostatic potential. More importantly, not all Arg/Lys possessing high electrostatic potential contribute to GAG binding. In fact, for some proteins, multiple loci of high G are observed in disparate locations, e.g., HS3ST3A1 () and HS3ST1 (), and not all are part of the GBS [36,37]. Thus, G at electropositive Arg and Lys residues alone does not always define GAG binding and specificity.

G at arginines and/or lysines does not identify the GBS on a protein.

G for Arg/Lys residues are mapped using 2DSE plots for multiple GAG-binding proteins including (a) antithrombin, (b) thrombin, (c) FGF2, (d) HS3ST3A1, (e) HS3ST1 and (f) HS2ST1. The maps reveal that GAG-binding site Arg/Lys residues may not always possess high G and not all Arg/Lys with high G on a protein are part of the GAG-binding site.

Desolvation is an important factor in GAG–protein binding

Considering that GES alone cannot define binding specificity, we set to identify other parameters that may contribute to specificity of GAG interactions. We hypothesized that to a first approximation the overall binding energy (ΔG) of a GAG–protein system is a composite of electrostatic (ΔGES) and desolvation (ΔGDS) free energy components. Additional terms, such as vibrational entropy and van der Waal’s energy [39,40], are likely to contribute but were not considered in this first approximation because it was our goal to identify GAG-binding sites from the structure of protein alone. Changes in vibrational entropy and van der Waal’s energy during binding will vary with nature of the cognate GAG binding partner, which precludes inclusion of these terms at this time. To test our hypothesis that ΔGES and ΔGDS are the major factors determining GAG-protein interactions, we studied the antithrombin–heparin pentasaccharide and thrombin–heparin systems, for which a large body of solution experimental data is available [41-45]. For both systems, ΔGES varies linearly (R2~0.6–0.7) with change in free energy of GAG binding (ΔΔGOBS) suggesting an important role for electrostatic forces to the interaction (; see Tables and ). Yet, while statistically significant (p<0.05), the correlation is moderate at best, which quantitatively confirms that electrostatics alone does not sufficiently address binding. To assess whether the release of solvent molecules is an important contributor in GAG-protein interactions, desolvation energy (ΔGDS) was calculated using the Poisson-Boltzmann Surface Area (PBSA) method [46] (see ). For both proteins, ΔGDS demonstrated an inverse correlation () with ΔΔGOBS (R2~0.6–0.7). Interestingly, electrostatic forces (ΔΔGES) and desolvation forces (ΔΔGDS) were found to be directly opposed to each other (R2 = 0.99, p<0.05) (). More importantly, the magnitude of desolvation effects was substantial, which implies that desolvation energy cannot be discarded and is an important driver of GAG-protein association.

Desolvation energy is critical for quantitative analysis of GAG–protein interaction.

Neither ΔGES (a and b) nor ΔGDS (c and d) alone explain the change in ΔGOBS for antithrombin (a and c) and thrombin (b and d) mutants studied to date. Any enthalpic gain due to electrostatics is opposed by desolvation (R2 = 0.99) in antithrombin (e) as well as in thrombin (f), suggesting that desolvation is critical for quantitative analysis of GAG-protein interactions. In all cases, the correlation was found to be significant at α = 0.05.

Since desolvation opposes electrostatic interactions, we hypothesized that residues capable of (a) producing non-ionic interactions such as hydrogen bonds (which also possess a partial covalent character to compensate for less favorable ionic interactions) with GAGs, and (b) lowering desolvation costs during such interactions, would play an important role in such interactions. Considering that GAGs possess several hydrogen bond acceptors, e.g., hydroxyls, carboxylates and sulfates, it is likely that asparagine and glutamine residues may serve as hydrogen bond donors.

GAG-protein sites are consistently identified by G at uncharged polar residues

It has been long known that GBSs do possess uncharged polar residues [29]. provides some typical examples. Their importance has mostly been construed from their presence in co-crystal structures and yet they have not been evaluated routinely by site-directed mutagenesis. However, whenever they have been, significant reduction in binding and activity are observed. For example, the Gln255Ala mutant of HS3ST3A1, a key enzyme that introduces a rare 3-sulfate group of the herpes simplex virus 1 recognition motif, demonstrates a >99% loss in activity [36]. Similarly, the Gln163Ala and Asn167Ala mutants of HS3ST1, the enzyme responsible for introduction of the 3-sulfate in unit F of DEFGH, show a ~60% drop in activity [37,47]. Likewise, Asn27, Gln123 and Gln134 of FGF2 have also been implicated in heparin binding, of which Asn27 makes a significant enthalpic contribution [34]. However, the reason for such observations has typically not been clear. Given the predominant roles played by Arg/Lys in GAG-protein binding, it seems unlikely that Asn/Gln alone would form motifs capable of binding GAGs. We predicted that a combined effect manifested by both, electropositive and uncharged residues, may provide us with clues about the roles of such residues in GAG-protein interactions.

Electrostatic interactions are effective at long range and therefore Arg/Lys may also affect the environment of other residues in a protein structure. For example, G at Asn and Gln residues may be significantly altered. demonstrates G at Asn and Gln residues (and other neutral hydrogen bonding donors) in the form of 2DSE plots. While G calculated at basic Arg and Lys residues failed to reliably identify the GBS (), the same calculated for neutral hydrogen bond donors identified the GBS consistently (). The highest G for any given protein structure is always present at the GBS, as observed for antithrombin, thrombin, FGF2, HS3ST3A1, HS3ST1, HS2ST1 and serum albumin (and ). The identification of the GBS was independent of specific or nonspecific nature of the protein. Quite clearly, the critical positioning of basic residues near uncharged hydrogen bond donors creates a unique environment which is preferentially bound by GAGs. Thus, G at neutral hydrogen bond donors can be called GGAG-binding. As will be demonstrated below, G can segregate specific GBSs from nonspecific ones.

2DSE plots for G at neutral hydrogen bond donors.

GAGs bind neutral hydrogen bond donors on the protein that possess significantly high GES. (a) Asn45 of antithrombin GAG-binding site possesses the highest GES within the structure. (b) In contrast, the nonspecific thrombin GAG-binding site demonstrates a diffused GES. Similarly, significantly high G are observed at (c) Asn27 of the FGF2 GBS; Asn27Ala mutation affects GAG-binding (ΔΔG~1.1 kcal/mol) almost as much as K125A (ΔΔG~1.7 kcal/mol), which had the largest effect, (d) Asn255 of the HS3ST3A1 GAG-binding site; the N255A mutant is inactive, and (e) Gln163 of HS3ST1; Gln163Ala mutant loses ~65% activity. (f) Diffused GES of HS2ST1 may represent its ability to bind low-sulfated GAGs. However, Asn91 and 112 of the HS2ST1 GAG-binding site possess a potential higher than His106, mutation of which is already known to affect GAG-binding.

Uncharged polar residues of specific GBSs possess significantly higher GGAG-binding than the rest of the protein

Antithrombin, FGF2 and sulfotransferases possess specific GBSs while thrombin and serum albumin are known to be nonspecific in nature [25,27,30-32,36,37,35,38,47-49]. All these proteins possess basic and uncharged polar residues in their binding site. We can correctly identify their GBSs by calculating GGAG-binding at uncharged hydrogen bond donors (vide supra), but can we segregate specific proteins from nonspecific ones? By first principles, a specific GBS will be an area on the protein surface that possesses significantly higher GGAG-binding compared to the rest of the protein. GAGs will therefore preferentially bind to that location rather than anywhere else.

For antithrombin, the location of Asn45 possesses a much larger GGAG-binding () than any other area on its surface, indicating that GAGs will preferentially bind to it. Replacement of Asn45 with Ala, Arg or Lys completely alters the unique electrostatic environment observed in the wt structure () suggesting Asn45 is a “hot spot” for binding the heparin pentasaccharide. In contrast, thrombin did not demonstrate any such unique neutral hydrogen-bond donor () despite the presence of many basic residues in its GBS. Likewise, bovine serum albumin demonstrated no outstanding neutral donor loci (). It is also known to interact with heparin nonspecifically [49].

An analysis of other GAG-interactors also confirmed our ability to identify such loci; the 2DSE plots for neutral hydrogen-bond donors of FGF2, HS3ST3A1 and HS3ST1 show focused G at highly localized regions () similar to AT (). The same observations are made at various dielectric constants and by employing the CHARMM charge model [50,51], (). Note that using the more accurate CHARMM charge model gave us a better representation of HS2ST1’s specificity (). Thus, GAGs prefer to engage loci containing a strategically positioned neutral hydrogen bond donors near positively charged residues. If the G at a few neutral hydrogen bond donors (G ) is significantly high in comparison with other regions of the protein surface, then it is a specific GBS. 2DSE plots are a convenient way of identifying the GBS and visualizing its specificity.

Further, a simple statistical distribution of GGAG-binding across any given protein can also be used to segregate specific GBSs from nonspecific ones, as we have done in . Here, GGAG-binding was visualized using violin plots. The specific proteins, antithrombin, FGF2, HS3ST1 and HS3ST3A1 possess at least one location where G is significantly higher (more negative) than the average over the entire surface, as delineated by the long tails in the plots. Conversely, the nonspecific proteins thrombin and serum albumin do not possess similar tails, which suggests G is rather uniform across most of their surfaces and GAGs will not prefer binding at one location more than any other. This explains why they are nonspecific.

Specific proteins demonstrate unique, non-uniform distributions of electrostatic potential across neutral hydrogen-bond donors.

Specific proteins such as antithrombin, FGF2, HS3ST1 and HS3ST3A1 demonstrate at least one location of electrostatic potential that deviates significantly from the mean. Nonspecific GAG-binding sites on proteins such as thrombin and serum albumin demonstrate a uniform, Gaussian distribution of the same, so no location is preferred significantly over another.

Pharmacophoric basis for targeting specific GBSs

What causes a biased GGAG-binding to exist in specific GAG-binding proteins? Close visualization of the antithrombin/FGF2-heparin interactions explains the origins of the unique environment surrounding Asn45 of antithrombin and Asn27 of FGF2 (). Both proteins possess positively charged basic residues very close to their respective neutral hydrogen bond donors, causing a large GGAG-binding at these locations. This is the unique feature responsible for GAG-binding specificity elicited in them. Similar structural features are also observed in the sulfotransferase enzymes (). Such features are a hallmark of specific GBSs.

The structural basis for existence of hot spots in GBSs.

Nature has designed specific GBSs by placing neutral hydrogen bond donors such as the ND2 and NE2 atoms of Asn and Gln respectively in close proximity to charged Arg or Lys residues, as seen in (a) antithrombin, (b) FGF2, (c) HS3ST3A1, (d) HS3ST1 and (e) HS2ST1. This close proximity maximizes the GES at these residues, thereby generating a specific GBS. Not all atoms are displayed, for the sake of visual clarity.

GAG sequences that target specific proteins can be identified by considering which mono- or disaccharides will bind these pharmacophoric features. The identification of GAG sequences that may bind with specificity to a given protein is an additional challenge in GAG-protein interactions. The structural diversity of GAGs requires docking of large GAG libraries onto protein surfaces, but can be used to reasonably represent GAG-protein interactions [8,10,13,15,16,18,40,52-54]. A single residue can be crucial for interactions in GAG-protein systems, e.g. antithrombin binds only heparin containing GlcNS3S [55] and FGF2 binds only IdoA2S-containing heparin [56]. Molecular modeling studies are also most successful at identifying binding partners for such systems where evidence of specific binding exists [18,52,54]. The Combinatorial Virtual Library Screening (CVLS) approach employs an “affinity filter” and a “specificity filter” to successfully identify GAG sequences that bind proteins [15,16,52]. Similar approaches can be refined to a great extent by focused searches targeting the pharmacophoric features identified by our method. First, our method can be used to identify hot spots, followed by docking only mono- or disaccharide libraries to identify preferred motifs that bind them. The CVLS algorithm can then be used to dock focused GAG oligosaccharide libraries enriched in those preferred mono- or disaccharides. Such a comprehensive approach would pave the path for design of chemical probes and drugs.

Discussion

Our computational results in conjunction with structural and biochemical studies indicate that (i) not all Arg/Lys carry equal and high electrostatic potential; (ii) not all Arg/Lys possessing high electrostatic potential are located in the GAG-binding site; and (iii) not all GAG-binding sites contain a neutral hydrogen-bond donor with significant electrostatic potential. On the contrary, high electrostatic potential is induced at neutral residues by nearby basic residues to generate GAG-binding sites with high specificity of recognition. In fact, we have found that considering the electrostatic potential that exists at neutral hydrogen bond donors not only depicts the correct GBS, but also reveals its specificity. Nature appears to engineer a few strategic neutral hydrogen-bond donors within the highly positively charged domain to induce specificity. At a fundamental level, this represents a major advance over the current paradigm of GAG–protein interactions, which are thought to originate only from Arg/Lys residues.

At a fundamental level, specificity of GAG–protein interaction can be of two types–‘biological’ specificity, which refers to how unique is the geometry of GAG-binding on the protein, and ‘chemical’ specificity, which refers to how unique is the sequence of a GAG that is recognized by the protein. Although most researchers assume that both biological and chemical specificities are identical for GAG-protein interaction, it is not so. For example, thrombin recognizes multiple GAG sequences, which can bind in different orientations on the protein [26] indicating that thrombin–GAG system is neither biologically specific nor chemically specific. In contrast, antithrombin recognizes a unique GAG sequence, which binds in a unique binding geometry indicating that antithrombin–GAG system is both biologically and chemically specific [24,25,28]. By the same token, proteins exhibiting intermediate forms of specificity, e.g., chemically nonspecific but biologically specific, and vice versa, may also exist. While concrete evidence of such systems is lacking, we predict that we will find such systems with increase in our ability to assess specificity of GBSs. The method we present here will certainly add to such abilities.

Operationally, this knowledge implies that high specificity GAG-binding sites on apo-protein structures will contain neutral hydrogen-bond donors with high (i.e., more negative) electrostatic potential compared to the rest of the protein. The advanced understanding can be used to develop a simple two-step algorithm for identification of proteins that can bind GAGs with reasonably high level of specificity ().

The two-step algorithm for identification of GBSs on proteins and elucidating their specificity.

The process involves preparation of protein; identification of neutral hydrogen bond donors in the structure; calculation of 2DSE plots for the protein; and evaluation of ‘hot spots’ for deduction of specificity of GAG–protein interaction.

The first step would involve calculation of electrostatic potential at every neutral hydrogen bond donor locus followed by statistical analysis of G gradients to identify ‘hot spots’ i.e., neutral loci among multiple Arg/Lys residues that presents high electrostatic potential. Hot spots can be identified easily using 2DSE plots. Such hot spots likely engineer high GAG-binding energy through ease of desolvation penalties and hydrogen bond formation, inducing unique GAG-binding orientation and specificity of interaction. Further biophysical investigation is required to confirm the importance of desolvation energy in eliciting specificity of binding, but our ability to elucidate specific and nonspecific GBSs is already a critical advancement because it allows identification of specific GAG-binding proteins. The presence of a specific GBS on a protein likely signifies a role for GAGs in regulating that protein. Mapping specific GBSs onto biological pathways will identify novel biological roles for GAGs. Furthermore, we now have the ability to rationally target biological pathways to design GAG-based chemical probes and drugs.

Our work is expected to fundamentally change the study of GAG-protein interactions because it minimizes reliance on approximations such as linear secondary structure sequences [4] or molecular modeling/docking on regions of high positive charge density [7-19]. Because the approach can work on most high resolution structures (1.5–3.0 Å) being reported today, it is likely to ease the discovery of ‘druggable’ GAG-binding proteins, which has been a major stumbling block. In fact, our method will be of major use in rationally guiding high-throughput docking of GAG sequences [15,16,52] for identifying novel therapeutically relevant GAGs and mimetics that can specifically bind to target protein(s) and also in computational analysis of the proteome to identify the GAG interactome.

Methods

2D-Surface Energy (2DSE) Plots

To the best of our knowledge, 2DSE plots have not been described in the literature. 2DSE plots enable quantitative visualization of 3D energy distribution on a protein surface in two dimensions. Here, the position of a locus on a protein surface is projected from 3D space onto a 2D plane to obtain a scatterplot, in which the area of the dot is scaled to reflect the G . The foundation of this novel graphical display is as follows. If ‘i’ is the ith atom at point (x,y,z) in 3D space and N the total number of atoms in the binding site, then total energy G is the sum of energies at each i th atom located at (x,y,z) coordinate. PDB files are a typical example where the atomic (x,y,z) coordinates are provided for proteins. Each atom can therefore be mapped onto the abscissa–ordinate plane by plotting (x +y ) against z. In this projection, the expression (x +y ) is the distance of an atom/residue from the origin on the x-y plane. Further, the size of each point in the plot can be scaled relative to each other to compare their contributions towards binding. Thus, the scatter plot helps visualize the energetics (e.g., G ) in 2D arising from loci in the 3D. Finally, each scatter point can be scaled using an exponential function to represent the change in energetics at individual loci in comparison to the standard, as defined by the well-established Eq 2.

Here, KD represents affinity (ΔG = RT ln K ) and can be substituted by a factor that can be correlated with affinity. G is one such factor known to influence GAG-protein interactions. Therefore, a direct comparison of affinities afforded by G at two different points becomes possible using Eq 2. Typical code for generating 2DSE plots is provided in Supplementary Methods.

GAG-Protein Interaction Energy and Desolvation Energy Calculation using Poisson-Boltzmann Surface Area (PBSA) Method

Electrostatic interaction energy (ΔGES) and desolvation energy (ΔGDS) involved in antithrombin-heparin or thrombin-heparin interactions were calculated using the Poisson-Boltzmann Surface Area (PBSA) method implemented in Openeye’s ZAP toolkit [46,57]. Bondi van der Waal’s radii [58], MMFF94 [59] charges and atom-centered Gaussians [46] were used for these calculations for which we assumed an inner dielectric constant of 2 [60] and an outer dielectric constant of 80. The pdb code 1tb6 [30] provided structural information for the calculations. ΔGES was calculated using the formula ΔGES = ΔGES,complex−(ΔGES,protein + ΔGES,GAG). Similarly, ΔGDS was calculated using ΔGDS = ΔGDS,complex−(ΔGDS,protein + ΔGDS,GAG).

Statistical Analysis

Statistical analysis was performed using SigmaPlot (sigmaplot.com), or R statistical environment (r-project.org), and the ggplot2 module (ggplot2.org). Statistical significance was calculated using a two-tailed ANOVA on data found to be normal using a Shapiro-Wilk test and demonstrating constant variance. In case of unequal variances, the non-parametric Spearman rank order correlation test was used.

Calculation of GES at Individual Residues

Multi-body G calculations were performed using the formula in implemented in the form of a Python script (Supplementary Methods). These calculations were performed at multiple dielectric constants, and also using the CHARMM charge model [50,51] to ensure that the results were similar.

2DSE plots for serum albumin.

(a) the G at Arg/Lys residues reflects the relatively hydrophobic nature of albumin (it is known to bind several hydrophobic ligands to reduce their bioavailability) and (b) the G at neutral hydrogen bond donors on serum albumin resembling the same map for thrombin ().

(TIF)

Click here for additional data file.

Mutation of Asn45 to Ala, Arg or Lys alters the electrostatic environment of the protein.

When compared to , these 2DSE plots clearly demonstrate that the “hot spot” at Asn45 cannot exist even on mutation to Arg/Lys.

(TIF)

Click here for additional data file.

The effect of various dielectric constants on potentials.

G on neutral H-bond donors were recalculated using dielectric constants (a) 2, (b) 3, (C) 4 and (d) 10. Clearly, the trend remains exactly the same.

(TIF)

Click here for additional data file.

Effect of using CHARMM charges on potentials.

CHARMM charges were used to calculate G on neutral hydrogen bond donors for all the cases. Clearly, while the G value may change, there is still a clear distinction between specific and non-specific proteins. However, HS2ST1 now seems far more specific than with MMFF94 charges, in line with expectations that the enzyme will possess a specific GBS.

(TIF)

Click here for additional data file.

An R script that creates 2DSE plots

(R)

Click here for additional data file.

The GAG-Binding Site Predictor written in python.

NOTE: This will require file 3 installed properly to function.

(PY)

Click here for additional data file.

A short molecule input/output module written in python.

(PY)

Click here for additional data file.

Table 1

Polar residues present in various GBSs.

These residues form direct interactions with GAGs, as evidenced by analysis of crystal structures.

Protein [PDB code]	Arg/Lys in GBS	Other polar residues in GBS
Antithrombin [1tb6, [30]]	Arg 46, 47, 129, 132, 136; Lys 114, 125, 133, 275	Asn 45
Thrombin [1tb6, [30]]	Arg 93, 101, 126, 165, 233; Lys 236, 240	Asn 184; Gln 256; His 87
Basic fibroblast growth factor (FGF2) [1fq9]	Arg 120Lys 26, 119, 125, 129, 135	Asn 27
3-O-sulfotransferase 3A1 (HS3ST3A1) [1t8u]	Arg 166, 190, 260, 370Lys 161, 162, 215, 259	Gln 255
3-O-sulfotransferase 1 (HS3ST1) [3uan]	Arg 67, 72, 197, 268, 276Lys 68, 123, 171, 173, 274	Asn 89, 167; Gln 163; His 92
2-O-sulfotransferase 1 (HS2ST1) [4ndz]	Arg 80, 184, 189, 190, 288Lys 111, 284, 289	Asn91, 108, 112; His 106, 140, 142

Table 2

Electrostatic interactions and desolvation energies for AT-heparin pentasaccharide complexes reported in the literature.

Mutant	ΔGES (kcal/mol)	ΔΔGES (kcal/mol)	ΔGDS (kcal/mol)	ΔΔGDS (kcal/mol)	Kd b (nM)	ΔΔGOBS (kcal/mol)
Wt	-199.8	0	184.6	0	6	0
K125Q	-125.1	74.7	116.6	-68.0	12	0.4
K136T	-200.0	-0.25	185.0	0.4	6	0
N135A	-198.9	0.8	183.7	-0.9	2	-0.6
N135A/R129Q a	-150.7	48.2	139.2	-45.4	1800	3.4
N135A/R129H a	-151.7	47.2	139.8	-43.9	820	2.9
N135A/K114A a	-85.4	113.5	80.2	-103.5	1800000 c	7.6
R132M	-156.7	43.0	145.2	-39.4	89	1.6
K133M	-159.4	40.4	147.0	-37.6	171	2.0

Calculated from a comparison with the N135A mutant.

Kd values were obtained from references [41–44]. ΔΔGOBS values were calculated using the formula ΔG = RTloge(Kd 1/Kd 2) where Kd 1 is for mutant and Kd 2 is for wt.

Calculated from data available in reference [43].

Table 3

Effect of electrostatic interactions on binding of various thrombin mutants with heparin.

Mutant	ΔGES (kcal/mol)	ΔΔGES (kcal/mol)	ΔGDS (kcal/mol)	ΔΔGDS (kcal/mol)	Kd b (nM)	ΔΔGOBS (kcal/mol)
Wt	-503.1	0	461.9	0	90	0
K169E	-275.9	73.8	392.4	-69.4	150	0.31
R175E	-401.6	101.5	366.8	-95.1	570	1.11
R233E	-337.3	165.8	311.5	-150.4	720	1.25
K236E	-275.9	227.2	255.5	-206.4	4800	2.38
K240E	-282.6	220.5	262.5	-199.4	1000	1.44
R233E-K240E a	-117.1	220.2	111.8	-199.7	16000	3.11

Calculated from comparison with the R233E mutant.

These values were reported in reference [45] as equilibrium constant for thrombin dissociating from the ternary antithrombin-thrombin-heparin complex. Since the mutations are at exosite II, where heparin binds (and not antithrombin), these quantify thrombin-heparin interactions. ΔΔG values were calculated using the formula ΔG = RTloge(Kd 1/Kd 2) where Kd 1 is for mutant and Kd 2 is for wt.