Molecular Basis of Chemokine CXCL5-Glycosaminoglycan Interactions*

Chemokines, a large family of highly versatile small soluble proteins, play crucial roles in defining innate and adaptive immune responses by regulating the trafficking of leukocytes, and also play a key role in various aspects of human physiology. Chemokines share the characteristic feature of reversibly existing as monomers and dimers, and their functional response is intimately coupled to interaction with glycosaminoglycans (GAGs). Currently, nothing is known regarding the structural basis or molecular mechanisms underlying CXCL5-GAG interactions. To address this missing knowledge, we characterized the interaction of a panel of heparin oligosaccharides to CXCL5 using solution NMR, isothermal titration calorimetry, and molecular dynamics simulations. NMR studies indicated that the dimer is the high-affinity GAG binding ligand and that lysine residues from the N-loop, 40s turn, β3 strand, and C-terminal helix mediate binding. Isothermal titration calorimetry indicated a stoichiometry of two oligosaccharides per CXCL5 dimer. NMR-based structural models reveal that these residues form a contiguous surface within a monomer and, interestingly, that the GAG-binding domain overlaps with the receptor-binding domain, indicating that a GAG-bound chemokine cannot activate the receptor. Molecular dynamics simulations indicate that the roles of the individual lysines are not equivalent and that helical lysines play a more prominent role in determining binding geometry and affinity. Further, binding interactions and GAG geometry in CXCL5 are novel and distinctly different compared with the related chemokines CXCL1 and CXCL8. We conclude that a finely tuned balance between the GAG-bound dimer and free soluble monomer regulates CXCL5-mediated receptor signaling and function.

Chemokines, a large family of signaling proteins, are highly versatile, play crucial roles in defining innate and adaptive immune responses by regulating trafficking of leukocytes, and also play a key role in various aspects of human physiology (1,2). Common to these diverse functions is the directed movement of various cell types to distal and remote locations. Cellular trafficking must be highly coordinated to elicit the required biological function, and its dysregulation could result in disease. Humans express ϳ50 chemokines that are broadly classified on the basis of their conserved cysteines into two major (CXC and CC) and two minor subfamilies (CX 3 C and C) (3,4). CXCL5/ENA78 belongs to a subset of CXC chemokines that function as agonists for the CXCR2 receptor ( Fig. 1) (5). CXCL5 is selectively expressed in alveolar epithelial type II cells, white adipose tissue macrophages, and cardiomyocytes and mediates diverse functions, from neutrophil trafficking to cancer cell migration and promoting obesity (6 -10). Elevated CXCL5 levels have been detected in various disease states, including rheumatoid arthritis, pancreatitis, inflammatory bowel disease, sunburn, and several cancers (11)(12)(13)(14).
Chemokine function is regulated via interactions with glycosaminoglycans (GAGs). 2 GAGs such as heparan sulfate (HS) are highly sulfated polysaccharides, are ubiquitously expressed by nearly all types of cells, and are anchored to the cell surface by covalent attachment to membrane proteins. GAGs can also exist as non-covalent macromolecular complexes in the extracellular matrix (15)(16)(17). Animal models and cellular studies have established that GAG interactions dictate the duration and steepness of chemokine gradients and that these gradients orchestrate cellular trafficking (18 -20). GAGs are acidic and contain carboxylate and sulfate moieties, whereas chemokines are basic or contain clusters of basic residues, thereby engineering electrostatic/H-bonding complementarity, resulting in binding.
The solution structures of human CXCL5 and its dimerization properties are known ( Fig. 1) (21). However, the structural basis and the molecular mechanisms by which CXCL5 binds GAGs remain unknown. To address this complete lack of knowledge, we characterized CXCL5-GAG interactions using NMR spectroscopy, isothermal titration calorimetry (ITC), and molecular dynamics (MD) simulations. Our studies show that the CXCL5 dimer is the high-affinity GAG ligand, that the GAG binding residues form a contiguous surface within a monomer, and that GAG-binding residues are also involved in receptor interactions, indicating that GAG-bound CXCL5 cannot acti-vate the receptor. Further, we observe that the GAG geometry in CXCL5 is novel and distinctly different than in the related chemokines CXCL1 and CXCL8 (22,23). MD simulations indicate that the roles of the individual lysines are not equivalent and that helical lysines play a more prominent role in determining binding geometry and affinity. We conclude that a delicate balance between the GAG-bound dimer and free soluble monomer regulates CXCL5-mediated receptor activation for orchestrated cellular trafficking to the target site.

Results
We characterized the structural basis of heparin interactions using solution NMR spectroscopy. We used three oligosaccharides of increasing length from a tetrasaccharide to a 14-mer (dp4, dp8, and dp14). For all oligosaccharides, a series of HSQC spectra was collected until essentially no change in chemical shifts was observed. We used heparin oligosaccharides as surrogates for heparan sulfate.
HS has a modular structure with sulfated sequences (defined as NS domain) separated by non-sulfated regions containing acetylated sequences (defined as NA domain). Both HS and heparin share a repeating disaccharide unit composed of glucosamine and a hexuronic acid. The structure of heparan sulfate is more diverse and complex because of differential N-deacetylation, N-sulfation, O-sulfation, and epimerization of glucosamine and hexuronic acid residues. Heparin, on the other hand, is more uniformly sulfated, has a lower molecular weight, can be considered a variant of HS, and has been shown to capture endogenous interactions (24).
CXCL5 exists as monomers and dimers (monomer-dimer equilibrium constant, ϳ1 M). To characterize the relative affinities of the monomer and dimer, we initially characterized heparin binding using an ϳ8 M protein sample where peaks corresponding to both the monomer and dimer could be observed in the HSQC spectrum (21). Upon adding heparin dp8 or dp14, peaks corresponding to the monomer disappeared, indicating that the dimer is the high-affinity GAG ligand (Fig. 2). Therefore, all subsequent titrations were carried out at ϳ100 M concentrations where the protein is essentially a dimer.
Binding of Heparin Oligosaccharides to the CXCL5 Dimer-For all oligosaccharides, we observed only one set of peaks during the titrations, although significant line broadening was observed for longer oligosaccharides at molar ratios of around 1:1 to 1:2, indicating that binding is in the fast-intermediate exchange regime on the NMR timescale. The CSP profiles are shown in Fig. 3. The pattern of CSP profiles for both dp8 and dp14 was similar. Chemical shift changes were observed for both solvent-exposed and buried residues, indicating that both direct and indirect interactions mediate the binding process. Perturbed residues are predominantly located in the C-terminal ␣ helix (C-helix), 40s turn, ␤ 3 strand, and N-loop (Fig. 4). In particular, the basic residues His 23 and Lys 25 from the N-loop, Lys 49 from the 40s turn, Lys 52 from the ␤ 3 strand, and Lys 64 , Lys 65 , Lys 69 , and Lys 76 from the C-helix show significant perturbation. In addition, chemical shifts of buried and acidic residues, Met 26 in the ␤ 1 strand and Val 66 and Asp 72 in the C-helix, were also perturbed. In the case of dp4, the same basic residues from the N-loop, 40s turn, ␤ 1 strand, and C-helix were perturbed; however, the extent of perturbation was lower, suggesting weaker interactions. Indeed, the apparent K d from NMR titrations was determined to be ϳ100 M compared with ϳ20 M for dp8 and dp14 (supplemental Fig. 1).
The CSP of basic residues can be attributed to direct interactions between the positively charged NH 3 ϩ group of lysines and sulfates/carboxylates of the GAG. However, describing the CSP of acidic and buried residues is less straightforward, as it could be due to direct or indirect interactions. To more stringently distinguish between direct and indirect interactions, we carried out 1 H, 15 N heteronuclear NOE experiments in the free and dp14-bound (Fig. 5) and dp8-bound (supplemental Fig. 2) forms. Heteronuclear NOEs are sensitive to motions in the picosecond-nanosecond timescale. Structured residues tend to have high NOE values (Ͼ0.8) and are minimally influenced by GAG binding. On the other hand, less structured or dynamic residues will have lower NOE values. These residues could become more structured on GAG binding, and, if so, their NOE values would increase, possibly reaching values comparable with those observed for structured residues. The data show higher NOE values for a number of residues in the bound form and, in particular, for the buried and acidic residues from the N-loop and C-terminal helix. These observations suggest that these residues are more structured in the bound form and that their CSP and altered NOEs must be due to indirect and not direct GAG interactions. A large change in NOE value for Lys 76 is striking considering that it also showed large CSP. However, it is relatively remote from other GAG binding residues (Fig. 4), and structural models also show that Lys 76 is not involved in binding, indicating that its CSP and NOE changes must be due to binding-induced structural changes (Fig. 5).
Stoichiometry of CXCL5-GAG Binding-We characterized the binding of heparin oligosaccharides to CXCL5 using iso-thermal titration calorimetry. ITC measurements, in addition to providing enthalpy, entropy, and binding affinity, also provide stoichiometry (25). The binding isotherms of dp8 and dp14 binding to CXCL5 are shown in Fig. 6. The data fit best to a single binding site model, yielding a stoichiometry of two GAG chains per CXCL5 dimer, and binding affinities (K d ) of 8.0 and 6.5 M for dp8 and dp14, respectively. These values were similar to those observed in the NMR experiments. Low heat release for dp4 indicates weak interactions, and so reliable K d or stoichiometry could not be determined.
Structural Basis of Receptor Interactions-Receptor activation involves interactions between the chemokine N-loop/␤ 3 strand and receptor N-domain residues (defined as site I) and between the ligand N-terminal and receptor extracellular/ transmembrane residues (site II) (26). Site I functions as a critical docking site, and we have shown previously that the structural basis of site I interactions can be studied outside of the context of the intact receptor by characterizing chemokine binding to N-domain peptides (27,28). Therefore, we characterized the binding of CXCL5 dimer to the CXCR2 N-domain peptide by NMR spectroscopy. CSP shows a prominent role for the N-loop and ␤ 3 strand residues (Fig. 7). Comparison of the GAG and receptor binding profiles indicates that residues that are perturbed on receptor binding are also perturbed on GAG binding, suggesting that GAG-bound chemokine is unlikely to bind the receptor. A schematic of the GAG-and receptor-binding residues and the extent of overlap is shown in Fig. 7B.
Structural Models of CXCL5-Heparin Oligosaccharide Complexes-The CXCL5 structure reveals that the basic residues identified from NMR experiments form a contiguous surface and a binding groove within a monomer. Moreover, an additional binding mode that spans the dimer interface across the antiparallel C-helices could be envisioned (Fig. 4). To gain conclusive insights into the structure of the complex, we carried out HADDOCK-based calculations that use CSP data as AIRs and shape complementarity and energetics to drive the docking process. The CXCL5 dimer and heparin dp14/dp8 were used as starting structures for the docking process. We observed similar results for both oligosaccharides and so confine our discussion to the CXCL5-dp14 models.
We performed three different HADDOCK runs to ensure that the input constraints did not bias specific structural models and that all possible binding geometries within a monomer and across the dimer interface were considered. In run 1, we modeled binding of one heparin with constraints given to both monomers of the dimer. In run 2, we modeled binding of two heparins with constraints given to both monomers of the dimer. In run 3, we modeled the binding of one heparin with constraints given to only one monomer of the dimer. Run 1 resulted in two major clusters. In one cluster, heparin bound to one monomer of the dimer with no evidence of interactions across the dimer interface (defined as model I) (Fig. 8). In the other cluster, heparin spanned the dimer interface and bound both monomers of the dimer (defined as model II) (Fig. 8). Run 2 resulted in only one major cluster with two heparins binding to two monomers of the dimer, similar to model I (Fig. 8). Run 3 also resulted in only one major cluster with heparin binding to a monomer of a dimer, similar to model I. NMR and ITC data are consistent with only model I. In model I, heparin engages all of the basic residues identified from NMR CSP and NOE measurements. However, in model II, interactions involve residues from the N-loop and C-helix of both monomers of the dimer, but interactions with the 40s turn and ␤ 3 strand residues are completely missing. Further, whereas model I allows binding of two heparins per dimer, model II allows binding of only one heparin per dimer that is inconsistent with the ITC data. HADDOCK results for dp14 alone revealed that the structures from model I could be further grouped into three subfamilies that are related through translation by a disaccharide unit (supplemental Fig. 3). In all three  15 N HSQC spectra showing the overlay of CXCL5 in the free (black) and dp14-bound (red) forms. Arrows indicate the direction of the peak movement. B-D, histograms of chemical shift changes in the CXCL5 dimer on binding heparin dp14 (B), dp8 (C), and dp4 (D). The basic residues Lys and His are shown in blue, and buried/acidic residues are shown in red. The horizontal lines at 0.19 ppm for dp14/dp8 and 0.10 ppm for dp4 represent the cutoff for a residue to be considered perturbed. The spectra were collected using a 100 M CXCL5 sample in 50 mM sodium phosphate (pH 6.0) at 25°C. families, basic residues from the N-loop, 40s turn, ␤ 3 strand, and C-helix from one monomer of dimeric unit are involved in binding. All of the binding poses had relatively similar energies. Further, there was no specific preference for the directionality of the heparin chain. Both orientations in each family of structures showed very similar energies, emphasizing that the structural plasticity of the binding surface is able to accommodate different sets of sulfate and carboxylate interactions. For both oligosaccharides, there was no evidence that C-helix Lys 76 is involved in binding, indicating that the CSP and NOE changes must be due to structure induction, as discussed earlier.
Analysis of the structures from HADDOCK docking showed that dp4 can bind in different orientations and locations and engage sulfates and carboxylates in different combinations. Major clusters correspond to a subset of C-helix and N-loop residues (model A) and a subset of N-loop, 40s turn, and ␤ 3 residues (model B) within a monomer of the dimer and a subset of C-helical residues from both monomers of the dimer (model C) (supplemental Fig. 4). Considering that a given cluster cannot engage all binding residues because of the dp4 length, multiple populations must exist in dynamic equilibrium with no preference for a particular model.
MD Simulation of CXCL5-Heparin Oligosaccharide Models-HADDOCK modeling studies for dp8 and dp14 oligosaccharides indicated that model I is the preferred mode of heparin binding. To understand residue-specific binding contributions, we carried out MD simulations of model I for dp8 heparin from HADDOCK docking. Simulations were carried out for 50 ns in NPT ensemble to explore the conformational behavior of negatively charged heparin over the positively charged domains of CXCL5 in the presence of explicit solvent. Throughout the simulation, lysine side chains were observed to be dynamic, and a given lysine residue could engage more than one sulfate/carboxylate group of the neighboring sugar residues. Simulations showed that there was considerable fluctuation of heparin and side chain-interacting residues in GAG binding epitopes on the protein surface (Fig. 9A). The heparin chain was also observed to undergo considerable translational motion within the observable span of 50 ns (supplemental Fig.  5).
To understand stability and specificity, residues participating in stable intermolecular H-bonds were calculated using a VMD H-bond plugin with a donor-acceptor distance   Of these residues, Lys 65 exhibited the strongest interaction, which persisted from initial to final state of simulation (supplemental Fig. 6).
GAG-protein co-crystal structures indicate that water molecules mediate binding interactions (29,30). These water-mediated interactions can further contribute to the stability of the complex (31-33). Our MD simulations also reveal the presence of such water-mediated interactions (supplemental Fig. 7), and a summary of these interactions is shown in Fig. 9C. Although Lys 65 and Lys 69 displayed consistent water-mediated interactions throughout the simulation (ϳ80 -120% occupancy), His 23 , Lys 25 , and Lys 64 displayed about half the occupancy (ϳ50 -60%). Other positively charged residues either did not or only minimally participated in these interactions.
We also calculated the relative single residue energy decomposition values to understand and also rank the contribution of the individual CXCL5 residues for heparin interactions. Fig. 9D shows that the direct H-bond-forming residues contributed relative free energy in the range of Ϫ6.1 to Ϫ1.2 kcal/mol. Rank order in terms of stability can be clustered into C-helical residues (Lys 64 , Lys 65 , and Lys 69 ) that provide the highest, N-loop residues (His 23 and Lys 25 ) that provide moderate, and 40s turn and ␤ 3 strand residues (Lys 49 and Lys 52 ) that provide the least contribution.

Discussion
Our studies indicate that the molecular basis of GAG heparin binding to CXCL5 is strikingly different compared with other closely related ELR chemokines (22,23,34). In CXCL5, lysine residues from the N-loop, 40s turn, ␤ 3 strand, and C-helix form a contiguous surface within a monomer of the dimer, with two GAGs binding two monomers with no interactions across the dimer interface. In stark contrast, heparin-binding residues in CXCL1 are present in two non-overlapping domains, with each domain spanning the dimer interface (22). On the other hand, GAG-binding residues in CXCL8 can be classified into core and peripheral residues that result in different GAGbinding geometries, with both interactions within a monomer and across the dimer interface (23). Residues that mediate GAG binding in CXCL5 also mediate GAG interactions in CXCL1. However, in CXCL1, additional interactions are observed from the N-terminal Arg 8 and ␤ 3 strand Lys 49 . Interestingly, N-terminal Arg is absolutely conserved among all ELR chemokines and, in fact, is part of the ELR motif. The residue corresponding to ␤ 3 Lys 49 is present in some chemokines and is actually a Glu in others, including in CXCL5 (Fig. 1C). These observations suggest that a few basic residues and their locations in the context of quaternary structures dictate very different GAG binding geometries.
Our ITC data indicate a stoichiometry of two oligosaccharides per CXCL5 dimer. We propose a "clamp" model for the CXCL5 dimer, similar to what was proposed earlier for CXCL1, in which the chemokine dimer is sandwiched between two GAGs (Fig. 10). Heparan sulfate has a modular structure consisting of sulfated regions (NS) interspersed with non-sulfated acetylated regions (NA). Therefore, it is possible that two NS regions on the same chain bind two monomers of dimer (Fig. 10A). Binding of two NS regions from a single HS has also been previously proposed for CXCL1 and CXCL8. However, the binding geometries are distinctly different and, in the case of CXCL1, differ by an ϳ90°rotation about the 2-fold axis (Fig.  10B). In the case of CXCL8, biochemical studies have proposed a horseshoe model where GAG lies in the same plane with NS domains binding parallel to the helix that is linked by the NA domain (35). Proteoglycans have multiple GAG chains, and so it is also possible that CXCL5 interacts with NS domains from two adjacent GAG chains because of spatial proximity on the cell surface (36).
Molecular dynamics studies show that rapid, transitory, and plastic conformational fluctuations mediate CXCL5-heparin binding and that C-helical residues actively engage the GAG chain. This plasticity is further enhanced by water-mediated interactions with surface-exposed residues. Single-residue decomposition studies suggest that C-helical residues, compared with the basic residues from other domains, play a prominent role in defining the binding affinity and geometry. These observations need to be confirmed with mutational studies. Interestingly, structure-function studies for CXCL1 show no such preference, with mutations in all domains resulting in similarly reduced GAG interactions.
Previous studies have shown that CXCL5 function is modulated by proteolytic cleavage (37). For instance, metalloproteases from Pseudomonas aeruginosa inactivate CXCL5 by cleaving the C-terminal helix at residues between Pro 61 and Phe 62 , which results in impaired neutrophil trafficking. Considering that these residues are buried in the GAG-bound form, CXCL5 would be much less susceptible to proteolysis, indicating that GAG interactions also indirectly play a role in evading microbial infection. Similar observations of proteolytic protection have been reported for chemokines CXCL1 and CXCL12 in the GAG-bound form (38,39).
Animal model studies have shown that the ability of chemokines to reversibly exist as monomers and dimers regulates neutrophil recruitment (24,40,41). During active cellular trafficking, CXCL5 concentration can vary by orders of magnitude, and so it may exist as a monomer, dimer, or both at different times and locations. Our results suggest that CXCL5 predominantly exists as a dimer in the GAG-bound form and also that receptor-binding residues are occluded in the GAG-bound For receptor interactions, residues that showed significant CSP and are solvent-exposed and N-terminal Glu 10 , Leu 11 , and Arg 12 (of the ELR motif) are considered involved in binding.
form. Interestingly, it has also been shown that GAG-bound CXCL1, CXCL8, and CXCL12 dimers are also occluded from receptor binding (22,23,39). Collectively, these observations indicate that, although the GAG geometries are quite different, GAG-bound chemokines are not accessible for receptor binding, indicating that only the free monomer activates the receptor. Differences in binding geometry could be intricately intertwined with optimal GAG-binding affinity, receptor interactions, and/or distribution between free and GAG-bound forms of the monomer and dimer. Importantly, the latter dictates GAG-bound haptotactic gradients and soluble chemotactic gradients. We propose that GAG interactions are highly regulated for each chemokine, which, in turn, influences the makeup and strength of the gradients, for orchestrating directed trafficking of different cell types to various locations in normal physiological functions and that dysregulation in these processes results in disease.

Materials and Methods
NMR Experiments-NMR spectra were recorded at 25°C using Bruker Avance III 600-and 800-MHz spectrometers equipped with QCI and TCI cryoprobes, respectively. Titration of heparin-derived oligosaccharides (dp4, dp8, or dp14) to 15 N-CXCL5 was carried out in 50 mM sodium phosphate (pH 6.0). Expression and purification of 15 N-CXCL5 were carried out as described previously (21). Heparin oligosaccharides were purchased from Iduron. According to the manufacturer, the oligosaccharides were purified using high-resolution gel filtration chromatography, the main disaccharide unit is IdoA,2S-GlcNS,6S (ϳ75%), and they show some variation in sulfation pattern and contain uronic acid at the non-reducing end and a C4-C5 double bond as a result of the heparinase endolytic action.
The protein concentration was ϳ100 M for all NMR experiments. Aliquots of 10 mM dp4, dp8, and dp14 prepared in the same buffer were added to the protein sample, and a series of 1 H, 15 N HSQC spectra was collected until essentially no changes in chemical shifts were observed. For all titrations, the final protein:ligand molar ratio was 1:6. The observed chemical shift perturbation (⌬␦ obs ) was calculated as a weighted average chemical shift change of 1 H (⌬␦ H ) and 15 N(⌬␦ N ). The apparent dissociation constant (K d ) was determined by fitting the binding-induced chemical shift changes (⌬␦ obs ) as described previously (42). In the case of receptor titrations, a series of 1 H, 15 N HSQC spectra was collected on adding aliquots of the CXCR2 N-domain (1 mM) to the ϳ70 M CXCL5 sample. The final  CXCL5:CXCR2 molar ratio was 1:4. To determine relative GAG affinities of the CXCL5 monomer and dimer, dp14 was titrated to ϳ8 M CXCL5 in 50 mM sodium phosphate (pH 7.5) at 40°C. At this concentration, both monomer (9%) and dimer (91%) peaks were observed. 1 H, 15 N NOE Experiment-Steady-state 15 N heteronuclear NOEs were measured using a gradient-selected, sensitivity-enhanced pulse sequence (43). The "NOE on" experiment used a 4-s proton presaturation period with a 2-s relaxation delay. In the "NOE off" experiment, presaturation was applied 4 MHz off-resonance with the same recycle delay. The heteronuclear NOE values were calculated as a ratio of peak intensities with and without proton saturation.
ITC-Theinteractionsofdp4,dp8,anddp14heparinoligosaccharides with CXCL5 were studied at 25°C using a Microcal VP titration calorimeter. The protein and GAG solutions were centrifuged and degassed under a vacuum before use. Titrations were performed by injecting 30 ϫ 6 l aliquots of 0.5 mM heparin oligosaccharides into 20 M CXCL5 solution in 50 mM sodium phosphate (pH 6.0). All titrations were carried out at least twice, and the raw data were corrected using buffer and protein controls and analyzed using Origin software supplied by Microcal.
Docking of the CXCL5-GAG Complexes-Molecular docking of heparin to the CXCL5 dimer was carried out using high ambiguity driven biomolecular docking (HADDOCK) (44,45). The docking was carried out using NMR chemical shift perturbations (CSPs) as ambiguous interaction restraints (AIRs) between the interacting partners to drive the docking process. Docking simulations were carried out using the CXCL5 dimer NMR structure (PDB code 2MGS) and oligosaccharide structures generated from the NMR structure of heparin 12-mer (PDB code 1HPN) as the starting structures (21,46). CSP data from NMR titration experiments were used as the AIRs to drive the docking process. Residues with a chemical shift perturbation greater than the average plus 1 S.D. and an accessible surface area Ͼ 30% were selected as "active" residues. Residues that were surface neighbors of active residues were selected as "passive" residues, and both active and passive residues were used as AIRs. Optimized parameters for liquid simulation (OPLS) were used for the non-bonded interaction. The topology and parameter files for heparin oligosaccharides were generated using the PRODRG server (47). A total of 3000 complex structures were generated during the initial rigid body docking. The best 1000 structures in terms of intermolecular energies were further subjected to semiflexible simulated annealing. During this step, semiflexibility was introduced for the side chains of the binding interface residues. This step was followed by explicit solvent refinement, where the oligosaccharide and the protein interface residues were allowed to be flexible. The pair-wise "ligand interface RMSD matrix" over all structures was calculated, and final structures were clustered using a cutoff value of 7.5 Å. The clusters were sorted using RMSD and "HADDOCK score" (weighted sum of a combination of energy terms).
MD Simulations-CXCL5-heparin dp8 complexes obtained from HADDOCK runs, which were consistent with both NMR and ITC data, were used as starting geometries for unrestrained MD simulations for 50 ns. We settled upon 50 ns based on the simulation trajectory that indicated sufficiency of information regarding the binding interactions. Analysis of data collected for the final half of the 50-ns run showed that all heavy atoms of CXCL5, including those belonging to C-terminal lysine residues that interact with the GAG sequence, reached a stable state (RMSD Ͻ 1.5 Å). This stability of the CXCL5 structure enabled quantitative analysis of GAG motion on the chemokine surface. The GAG sequence was found to exhibit translational motion on the positively charged surfaces. These motions arose from interaction of sulfate groups with more than one lysine/ histidine partners so that the GAG sequence would return back to the original set of interactions, suggesting reestablishment of interactions. This implied that all key events were being sampled within the time frame of 50 ns.
Each docked complex was methodically prepared using the Leap module of AMBER 14 using AMBER-ff12SB force field and GLYCAM_06J-1 force field (48,49) parameters for CXCL5 and heparin dp8, respectively. The total charge of the complex was then brought to zero by adding an adequate number of counterions. The charge-neutralized complex was centered in a three-point water box (TIP3P) with a minimum distance of 12 Å between the wall and the nearest atom in a complex, followed by generation of initial coordinates and parameters for the starting solvated complex (50). Initially, the system was relaxed to achieve a minimum energy state. Energy minimization was then carried out in two steps with the non-bonded cutoff of 10 Å. In the first step, the solute atoms, including the counterions, were restrained by a harmonic potential with a force constant of 100 kcal/(mol⅐Å 2 ). The water molecules were relaxed using 500 cycles of steepest descent and 2000 cycles of conjugate gradient methods. In the second step, the whole system was relaxed using conjugate gradient minimization of 2500 cycles without any restraints. Equilibration and the simulation process were validated using the physical observables of the system, which confirmed that the system obeyed NPT ensemble. The MD simulations were performed over 50 ns, with snapshots captured every 1 ps. Throughout the simulation, the sugar puckering for IdoA2S and GlcNS6S was maintained at 2 S O and 4 C 1 conformations, respectively.
The binding free energy of the CXCL5-heparin complex was computed using the post-processing MM/PB(GB)SA method from the reordered trajectories (51). MM/PB(GB)SA employed single-residue energy decomposition to estimate energy contribution of the individual CXCL5 residues responsible in the bound state. Both energy calculations were performed with default parameter settings by employing the Python version of MM/PB(GB)SA module from Amber Tools13 (52).