Insights on the structure and dynamics of glycosaminoglycans and their interactions with langering: NMR and computational studies

Tesis Doctoral: Insights on the Structure and Dynamics of Glycosaminoglycans and their Interactions with Langerin: NMR and Computational Studies Juan Carlos Muñoz García Sevilla, 2013 UNIVERSIDAD DE SEVILLA CONSEJO SUPERIOR DE INVESTIGACIONES CIENTÍFICAS Laboratorio de Glicosistemas Instituto de Investigaciones Químicas Insights on the Structure and Dynamics of Glycosaminoglycans and their Interactions with Langerin: NMR and Computational Studies por Juan Carlos Muñoz García Disertación presentada por la Universidad de Sevilla para aspirar al Título de Doctor en Química Juan Carlos Muñoz García Sevilla, 2013 Dr. PEDRO MANUEL NIETO MESA, Investigador Científico (CSIC), y Dr. JESÚS ANGULO ÁLVAREZ, Doctor Contratado “Ramón y Cajal” (CSIC) CERTIFICAN Que el trabajo titulado Insights on the Structure and Dynamics of Glycosaminoglycans and their Interactions with Langerin: NMR and Computational Studies se ha llevado a cabo bajo nuestra dirección y asesoramiento en los laboratorios del Instituto de Investigaciones Químicas (IIQ, CSIC-US) del Centro de Investigaciones Científicas Isla de la Cartuja (CICCartuja), Sevilla, constituyendo la memoria presentada por Juan Carlos Muñoz García (Licenciado en Química) para optar al grado de Doctor en Química. Sevilla, julio de 2013 Los directores Dr. Pedro M. Nieto Mesa Dr. Jesús Angulo Álvarez El trabajo de investigación recogido en la presente memoria se ha llevado a cabo fundamentalmente en el Insituto de Investigaciones Químicas (IIQ), que forma parte del Centro de Investigaciones Científicas Isla de la Cartuja, centro mixto del Consejo Superior de Investigaciones Científicas (CSIC) y la Universidad de Sevilla (US). El doctorando ha sido financiado mediante una beca predoctoral JAE-predoc del CSIC obtenida por concurrencia competitiva, durante el periodo 2009-2013. Dicha beca predoctoral también ha financiado dos estancias cortas del doctorando, la primera en los laboratorios del Prof. Robert J. Woods, en el Complex Carbohydrate Research Center (CCRC) de la Universidad de Georgia (Athens, Georgia, EEUU), y la segunda en el grupo liderado por la Dr. Anne Imberty, del Centre de Recherches sur les Macromolécules Végétales (CERMAV-CNRS), Grenoble (Francia). Algunos de los experimentos de RMN se han realizado en el Laboratorio de RMN de Barcelona (LRB), adscrito a los Centres Científics i Tecnològics (CCiT) de la Universidad de Barcelona, haciendo uso del espectrómetro de 800 MHz. Esto ha sido posible gracias a la financiación recibida por el antiguo Ministerio de Ciencia e Innovación. A Encarni A mis padres Index vii Index Nota al lector i Abbreviations iii Index vii Chapter 1. Introduction and objectives 1 1.1. Carbohydrates: complex macromolecules 1 1.1.1. Chemical structure 4 1.1.2. Synthetic carbohydrates 11 1.1.3. Biological importance 12 1.2. Glycosaminoglycans (GAGs) 14 1.2.1. Heparin and the singular role of L-iduronic acid 15 1.2.2. Heparin binding proteins: the acidic Fibroblast Growth Factor (FGF-1) case 22 1.2.3. Hyaluronic acid 24 1.3. C-type lectins 27 1.3.1. Structure of C-type lectin receptors (CLRs) 29 1.3.2. Structural features of glycan binding to CRDs 31 1.4. Langerin: a natural barrier to HIV-1 infection 33 1.4.1 Novel role of epithelial LCs and their transmembrane protein Langerin in 35 HIV-1 infection 1.4.2 Langerin 3D structure and known ligands 38 1.5 Objectives of the research 43 Chapter 2. Techniques and tools 47 2.1. Structure by NMR: the Nuclear Overhauser Effect (NOE) 49 2.1.1. Origin of the NOE 49 2.1.2. Rotating-frame NOE (ROE) 51 2.1.3. Measuring internuclear distances 53 2.2. Ligand-based NMR spectroscopy for binding studies 57 2.2.1. The equilibrium kinetics of binding: the fast exchange approximation 58 2.2.2. Transferred-NOESY experiment (tr-NOESY) 59 2.2.3. Saturation Transfer Difference spectroscopy (STD) 62 2.3. Molecular modelling 68 2.3.1. Modelling of glycans 69 2.3.2. Carbohydrates force fields 71 2.3.3. Molecular Dynamics (MD) simulations 76 2.3.4. Docking 79 Index viii Chapter 3. Structural studies of heparin-like oligosaccharides by NMR 89 and MD techniques 3.1 Library of sulphated trisaccharides 89 3.1.1 Background 89 3.1.2 Results and discussion 91 3.2 An inactive hexasaccharide sequence for the FGF-1 mitogenic activity 113 3.2.1 Background 113 3.2.2 Results and discussion 114 3.3 Methodology 127 3.3.1 Nuclear Magnetic Resonance 127 3.3.2 Modelling 129 Chapter 4. Structural features underlying Langerin interactions with GAGs 137 4.1 Calcium-dependent interactions 137 4.1.1. Sulphated GAGs: heparin-like trisaccharides 137 4.1.2 Non-sulphated GAGs: hyaluronic acid disaccharides 152 4.2 Non calciumdependent interactions 165 4.2.1 Results and discussion 165 4.3 Methodology 168 4.3.1 Nuclear Magnetic Resonance 168 4.3.2 Modelling 170 Capítulo 5. Conclusiones 178 References 182 Publications 197 Chapter 1 Introduction and objectives 1. Introduction and objectives 1 1.1 Carbohydrates: complex macromolecules Carbohydrates are the most abundant type of biomacromolecule existing in Nature, either alone or forming glycoconjugates with proteins (proteoglycans or glycoproteins) and lipids (glycolipids or lipopolysaccharides), and thus they are available in large quantities from natural sources. They function as energy storage (starch, glycogen), starting material in biosynthesis, structural constituents of plants (cellulose), and they are also the major components of shells, insects or crabs (chitin), bacterial cell walls (lipopolysaccharides), or virus capsids (HIV). Isolation, purification and chemical modification of carbohydrates have been areas of great interest and exploitation during the last decades. Most of the first studies on carbohydrates focused on plant polysaccharides (e.g. cellulose, starch, pectins) due to their wide range of applications. More recently, the role of carbohydrates in biological events was recognized[1], and glycobiology emerged as a new and challenging research area at the interface of biology and chemistry. In this context, carbohydrate-mediated recognition events are of key importance in biological phenomena, playing a pivotal role to the study of protein-carbohydrate interactions. Actually, the binding protein partners of carbohydrates encompass a wide variety of macromolecules involved in functions such as recognition, biosynthesis, modification, hydrolysis, and so on. A relevant part within the research field of glycobiology encloses the determination of the structures and functions of complex sugars (glycans). This has become a critical facet of postgenome science, proteomics in particular, since many proteins are post-translationally modified by glycosylation, and these modifications alter and regulate biological activities. Thus, glycans represent a major class of post-translational modifications that dramatically enhance the functional diversity of proteins (figure 1). During the last years, the scientific standpoint of glycobiologists has increasingly moved towards the concept of glycome, i.e., the complete set of glycan structures expressed by specific cells, tissues or organisms. This has in turn led to a need for analysis of larger numbers of glycan structures, which has accelerated the development of technologies with high-throughput potential[2]. The emerging omics domain of glycomics has dropped behind that of genomics and proteomics, mainly because of the inherent difficulties in analysing glycan structure and function[3]. The term glycomics is formed by the prefix glyco-, which means sweetness or sugar, followed by -omics (i.e., field of study) to be consistent with the naming convention established by genomics (which deals with genes) and proteomics (dealing with proteins). The definition of glycomics has evolved to cover a range of scientific disciplines that are applied to study the structure and function of carbohydrates (sugars) in biological systems. 1. Introduction and objectives 2 Figure 1. Glycome enhancement of the molecular and functional diversity of the proteome. Protein expression is based on a genetically encoded template, but post-translational modifications of proteins dramatically promote their functional diversity. The glycome represents the main class of posttranslational modifications, providing biological access to vast information space at minimum genetic cost. Source: Turnbull and Field 2007[2]. The nine common sugars found in mammalian cells (figure 2, top) can be combined in a myriad number of ways to form complex carbohydrate structures. The glycan collection (glycome) of a given cell or organism is thus many orders of magnitude more complex than the genome or the proteome. Thanks to the rapid development of enabling technologies such as high-throughput mass spectroscopy, glycan microarrays and carbohydrate chemistry, deciphering the complexity resulting from this diversity is increasingly possible (figure 2, bottom). In this regard, bioinformatics is a fundamental technology of growing importance in managing and integrating the diverse data sets from the different technologies (figure 2, bottom). The developing field of glycomics is earning its place alongside other established “omics” fields such as genomics and proteomics, as it was anticipated by experts in the field: “We envisage that the collective enterprise of glycomics over the next decade will begin the process of decoding the glycome, thereby yielding many new insights into its myriad functions and producing diverse advances in the biomedical arena”. Nature Chemical Biology (2007) 3; 74-77. “The knowledge gained from glycomics will be as important as a basis for the pharmaceutical industry as that discovered in the field of genomics and proteomics during the last 30 years”. Chem. Eur. J. (2005) 11; 3194 – 3206. 1. Introduction and objectives 3 Figure 2. (Top) The nine common sugar “letters” of mammalian glycomics. (Bottom) Combinations of cuttingedge technologies that aid in deciphering the glycocode, leading to new insights and biomedical applications. Currently, they are being exploited to endeavour large-scale analyses of the structure-function relationships of the glycome. Sources:Weiss and Lyer 2007[4] (top), Turnbull and Field 2007[2] (bottom). Many exciting applications of glycomics approaches have become evident, including diagnostics, new routes to glycotherapeutics and defined recombinant protein drugs. 1. Introduction and objectives 4 1.1.1 Chemical structure In contrast to other important biomolecules such as proteins and nucleic acids, carbohydrates can form branched structures by substitution of one or several hydroxyl groups, making them extremely complex and heterogeneous. The structure of an oligosaccharide is determined by the monosaccharide sequence, the glycosidic linkage sites, the stereochemistry of the glycosidic linkages ( or ) and the degree and type of substitution of hydroxyl groups (such as O-methylation or O-sulphation)[5]. Also, it is important to emphasize that water plays a central role in defining oligosaccharide conformation by influencing the geometry of the glycosidic linkages[6]. Carbohydrates, or less commonly named “hydrated carbons”, are compounds that often have the empirical formula Cn(H2O)n. A monosaccharide is an aldehyde or a ketone containing at least two additional hydroxyl groups. While two monosaccharides connected by a glycosidic bond are named a disaccharide, carbohydrates with three to ten units are typically called oligosaccharides, and larger structures are known as polysaccharides. Carbohydrates are chiral and optically active, and the majority of the naturally occurring ones present the D configuration[7]. Carbohydrates with five or six carbon atoms, pentoses and hexoses, respectively, can form intramolecular hemiacetals between the carbonyl group and the hydroxyl group on carbon 4 or 5. The resulting rings are called hexopyranoses (six-membered rings) or pentafuranoses (five-membered rings). The cyclic and the acyclic forms exist in equilibrium, with the hemiacetals (rings) being the most abundant forms. Upon hemiacetal formation, the former carbonyl carbon becomes a new stereocenter, called the anomeric center, with the hydroxyl group either equatorial or axial ( and  configuration in D-glucose, respectively). The equilibrium between the cyclic and acyclic forms allows the two anomeric forms to interconvert, a process called mutarotation[8]. Hexopyranose rings can exhibit different conformations. Figure 3 illustrates the conformational itinerary among the main canonical conformations (or puckers) of the hexopyranose ring. For the majority of carbohydrates, the energetically most favoured ring conformation is the chair conformation, which exists in two distinct forms, the 4C1 and the 1C4 chairs (C=chair), where the numbers refer to the atoms above (superscript) and below (subscript) a reference plane. The 4C1 conformation is generally favoured for D-sugars due to fewer non-bonding interactions between the ring substituents. For example, the 4C1 conformation is the only one observed by NMR spectroscopy for D-glucose[9]. Also, hexopyranoses can exist in boat (B), half chair (H) and skew-boat (S) conformations. For instance, L-iduronic acid presents three low-energy conformations in equilibrium: 1C4, 2SO and 4C1[10]. Glucopyranose exists in a 64:36 mixture between the  and  forms in aqueous solutions. Just considering steric hindrance, it would be expected for electronegative substituents to prefer the equatorial disposition; however, the opposite occurs. This higher than expected occurrence of the axial form is due to the endo-anomeric effect[11], which identifies the preference of the electronegative 1. Introduction and objectives 5 substituent at the anomeric carbon (C1 in aldoses) for an axial configuration (α-anomer), rather than for the equatorial orientation (β-anomer). This preference has its basis in the electronic structure of the O5–C1–O1 atomic sequence. The widely accepted justification for the anomeric effect is that it originates from nσ* hyperconjugation between the lone pair of electrons placed at the non-bonding orbital (n) of the ring oxygen atom (O5) and the antibonding σ* orbital of the adjacent C1–O1 bond , although there are other interpretations (figure 4 )[12]. For pyranoses, stabilizing hyperconjugation is maximized in the α-anomer[13]. On the other hand, the exo-anomeric effect manifests the preference of the Og-Cx glycosidic bond (Og: glycosidic oxygen; Cx: adjacent carbon to the right) to adopt a gauche orientation with respect to the C1-O5 bond (C1: anomeric carbon; O5: oxygen ring), thus it contains itself the the rotameric preference about the C1–Og bond (known as torsion ). The exo-anomeric effect arises again from hyperconjugation within the O5–C1–O sequence, but this time it is between the lone pair of electrons placed at the non-bonding orbital (np) of the Og oxygen and the antibonding σ* orbital of the O5–C1 bond (figure 4). Figure 3. Diagram of the pseudorotational itinerary of the pyranose ring according to Jeffrey and Yates[14], based on the Cremer-Pople ring-puckering coordinates θ and . The polar 4C1 (θ=0º) and 1C4 (θ=180º) chairs, together with the 12 equatorial puckers (θ=90º) are shown. The envelope (E) and half-chair (H) conformations are not pictured. 1. Introduction and objectives 6 Figure 4. Schematic representation of the stereoelectronic genesis of the endoand exo-anomeric effect. The diversity of carbohydrate structures results from the broad range of monomers (>100) of which they are composed and the different ways in which these monomers are joined (glycosidic bonds). Thus, even a small number of monosaccharide units can provide a large number of different oligosaccharides (also referred to as glycans), including branched structures, a unique feature among biomolecules. For example, the number of all possible linear and branched isomers of a hexasaccharide exceeds 1012[15]. Providing a structural basis for the multitude of biological roles played by carbohydrates, it is imperative to accurately determine their spatial (conformation) and dynamic properties in aqueous solution (the importance of dynamics in structural biology was highlighted by the prediction that ∼25% of mammalian proteins are fully disordered[16]. This goal promises to enable structure-based design of new medicines and materials, but its realization remains challenging due to experimental and computational difficulties of probing carbohydrate motions, which occur over a broad range of time scales. For instance, glycosidic linkages liberate on nanosecond time scales while pyranose ring conformational exchange (or puckering) and anomerization are microsecond[17] and millisecond time scale phenomena, respectively. The recognition mechanism of carbohydrates depends on their conformation, which is determined by 1) the sequence of the monosaccharides in the glycan, 2) the anomeric centres (i.e., α or β), 3) the linkage positions (i.e., 1–3, 1–4, 1-6), and 4) the chemical modifications to the core glycan (i.e., sulphation, phosphorylation, methylation, acetylation, etc.). The strength of this interaction is also determined by the carbohydrate conformation and orientation with respect to the binding site. Carbohydrates and their derivatives possess many hydroxyl groups and thus a large number of rotatable bonds. Due to the many hydroxyl groups, these compounds are usually highly water soluble, 1. Introduction and objectives 13 structure factors in extra-cellular matrices[29], and postor co-translational modifications of polypeptides[30]. Correct glycosylation patterns are essential for normal cell and organism function, and aberrant glycosylation is associated with numerous human diseases[31]. Polysaccharides make up a substantial part of bacterial cell-walls, the greatest part being peptidoglycans giving the membrane mechanical strength. Other polysaccharides, such as lipopolysaccharides, capsular polysaccharides and exopolysaccharides, are to a large extent covering the cell-walls of bacteria. Microbial surface polysaccharides play an important role in bacteria-host interactions. In Gram-negative bacteria the lipid bilayer contains lipopolysaccharides (LPS’s). LPS’s consist of three regions, the lipid-A, the core region and the O-antigen. When the lipid-A is released from the membrane it is toxic and harmful to mammals. The O-antigen consists of repeating units with 2-8 carbohydrate moieties, and is very diverse. The O-antigen is the part of the bacteria recognized by the immune system[5]. Bacteria can also produce capsular polysaccharides, which is an extracellular coat surrounding the bacteria associated with virulence[19]. An additional kind of polysaccharide produced by bacteria is exopolysaccharides (EPS). They are high molecular weight polymers generally composed of repeating units of D-glucose, D-mannose, D-galactose, L-fucose, L-rhamnose, and D-glucuronic acid[32]. Exopolysaccharides are used in a number of industrial products, e.g., as food additives and in medical applications[33]. Figure 12. Multiple roles of heparan sulphate (HS) in leukocyte entry into sites of inflammation. HS on activated endothelium binds L-selectin on leukocytes (granulocytes, monocytes and lymphocytes) during the rolling phase of leukocytes over the endothelium. Endothelial HS binds and presents chemokines to chemokine receptors on leukocytes, which leads to activation of leukocytes and movement of leukocytes towards the site of inflammation, while HS is also involved in the transport of chemokines across the endothelial cell barrier. 1. Introduction and objectives 14 1.2 Glycosaminoglycans (GAGs) Glycosaminoglycans (GAGs) are long unbranched polysaccharide chains consisting of repeating disaccharide units. By definition, one of the two sugars in this repeating unit is a D-hexosamine residue (i.e., D-glucosamine, D-GlcN, or D-galactosamine, D-GalN), giving glycosaminoglycans their name. The other monosaccharide is an uronic acid (D-glucuronic acid, D-GlcA, or L-iduronic acid, L-IdoA). Thus, GAGs are divided into different classes according to the nature of their repeating unit (table 1): hyaluronan (HA), keratan sulphate (KS), chondrotin sulphate (CS), heparan sulphate (HS), heparin (HEP) and dermatan sulphate (DS)[19]. HS and CS are synthesized in the Golgi apparatus, where the individual GAG chains are O-linked to a core protein, forming a large proteoglycan (PG)[34]. Keratan sulphate, on the other hand, can be either N-linked or O-linked to the core protein of the PG[35]. HA is not synthesized in the Golgi from the core protein but rather by an integral plasma membrane synthase, which secretes the nascent chain immediately[36]. The biosynthesis of GAGs is a complex non-template-driven process involving several enzymes that assemble the GAG polymer and then sulphate them at specific positions. The GAG attachment sites on the core protein of the PG have a consensus Ser-Gly/Ala-X-Gly motif. Importantly, this versatility of GAGs biosynthesis yields a large number of possible substitution patterns that make them into highly dense information carriers. The study of molecular recognition processes between carbohydrates and protein receptors has attracted considerable attention during the past years[37]. Among those processes, the interaction of glycosaminoglycans (GAG) and signalling proteins (figure 12), due to its biological relevance and large number of cases, has attracted the attention of many research groups and structural details of those interactions have been extensively studied[38]. There are many examples where the specificity of the interaction between the GAG and signalling proteins relies on the substitution pattern[39]. Among them it should be mentioned the interaction with the fibroblast growth factor (FGF) family[40], chemokines[41] and cytokines, antithrombin-III (AT-III)[42], lipases, apolipoproteins and ECM and plasma proteins[43]. 1. Introduction and objectives 15 Table 1. Types of GAGs and their disaccharide building blocks[44].* Category Disaccharide repeating unit Chemical substitutions Heparin (HEP) L-IdoA2X-α(1-4)-D-GlcNY,3X,6X-α(1-4) X: SO3-, Y: Ac or SO3Heparan sulphate (HS) D-GlcA-β(1-4)-D-GlcNY,3X,6X-α(1-4) X: SO3-, Y: Ac or SO3Chondroitin sulphate (CS) D-GlcA2X-β(1-3)-D-GalNAc,4X,6Xβ(1-4) X: SO3Dermatan sulphate (DS) L-IdoA2X-α(1-3)-D-GalNAc,4X,6Xβ(1-4) X: SO3Keratan sulphate (KS) D-Gal6X-β(1-4)-D-GlcNAc,6X-β(1-3) X: SO3Hyaluronic acid (HA) D-GlcA-β(1-3)-D-GlcNAc-β(1-4) None * Note that the abbreviations stand for: L-iduronic acid, L-IdoA; D-glucuronic acid, D-GlcA; D-glucosamine, D-GlcN; D-galactosamine, D-GalN; D-galactose, D-Gal. The acetyl (COCH3) and sulphate (OSO3−) groups are abbreviated using Ac and S, respectively. Also note that for HEP/HS and CS/DS pairs, only the major disaccharide repeating unit has been indicated, although the disaccharide sequence of its partner is also present. 1.2.1 Heparin and the singular role of L-iduronic acid Heparin is a linear polymer consisting of disaccharide repeating units of 1 4-linked hexopyranosyluronic acid and 2-amino-2-deoxyglucopyranose (glucosamine) residues[45]. The uronic acid residues typically consist of 90% L-idopyranosyluronic acid (L-iduronic acid, L-IdoA) and 10% D-glucopyranosyluronic acid (D-glucuronic acid, D-GlcA) 1 . Heparin has the highest negative charge density of any known biological macromolecule. This is the result of its high content of negatively charged sulphate (OSO3-) and carboxylate (COO-) groups[46]. Indeed, the average heparin disaccharide contains 2.7 sulphate groups. The most common structure occurring in heparin is the trisulphated disaccharide L-IdoA2S-(1-4)-D-GlcNS,6S shown in figure 13. However, other substitution patterns also participate, leading to the microheterogeneity of heparin. For example, the amino group of the glucosamine residue may be substituted with an acetyl or sulphate group or unsubstituted. Also, the 3and 6positions of the glucosamine residues can either be substituted with an O-sulphate group or unsubstituted. The uronic acid, which can either be L-iduronic or D-glucuronic acid, may also contain a 2-O-sulphate group. Glycosaminoglycan heparin has a molecular weight range of 5-40 kDa, with an average molecular weight of about 15 kDa and an average negative charge of approximately -75. The 1 The terms “hexopyranosyluronic” (idopyranosyluronic or glucopyranosyluronic) and “uronic” (iduronic or glucuronic) will be interchangeably used in the present Doctoral Thesis. The first term is, indeed, more restrictive as it exclusively refers to the  and  anomers in the ring form, whereas the second term allude to both the cyclic and acyclic sugars. Since GAGs are polymers formed by the repetition of cyclic sugars (rings), “hexopyranosyluronic” is indeed the most precise term. 1. Introduction and objectives 16 large structural variability of heparin, coming form its high heterogeneity and polydispersity, makes it an extremely challenging molecule to characterize. The structural complexity of heparin can be considered at several levels. At the proteoglycan (PG) level, different numbers of polysaccharide (or glycosaminoglycan) chains (possibly having different saccharide sequences) can be attached to the various serine residues present in heparin’s core protein. During their biosynthesis, heparin chains are attached to a unique core protein, serglycin, found only in mast cells and some hematopoietic cells. Tissue proteases act on this core protein to release peptidoglycan heparin, a small peptide to which a single long polysaccharide chain (100 kDa) is attached. This peptidoglycan is short-lived as it is immediately processed by a β-endoglucuronidase to a number of smaller (about 15 kDa) polysaccharide chains called glycosaminoglycan (GAG) heparin[47]. Most of the chemical and physical properties of heparin are related to GAG structure or sequence, conformation, chain flexibility, molecular weight and charge density. Heparan sulphate is structurally related to heparin but it is much less substituted with sulphate groups than heparin and has a more varied structure (or sequence). Like heparin, heparan sulphate is a repeating linear copolymer of a uronic acid 1 4 linked to a glucosamine residue[48]. Although D-glucuronic acid predominates in heparan sulphate, it can contain substantial amounts of L-iduronic acid. Heparan sulphate generally contains only about one sulphate group per disaccharide, but individual HS may have a higher content of this group. Also, heparan sulphate chains often contain domains of extended sequences having low or high sulfation[49]. While heparan sulphate contains all of the structural variations found in heparin (and vice versa), the frequency of occurrence of the minor sequence variants is greater than in heparin, making HS structure and sequence much more complex. HS chains are also polydisperse, but are generally longer than heparin chains, having average molecular weight of about 30 kDa ranging from 5 to 50 kDa[50]. Heparan sulphate is biosynthesized, as a proteoglycan, through the same pathway as heparin; however, unlike heparin, the HS GAG chain remains connected to its core protein. Heparan sulphate is ubiquitously distributed on cell surfaces and is also a common component of the extracellular matrix[49, 51]. Two types of core proteins, the syndecans (an integral membrane protein) and the glypicans (a GPI-anchored protein), commonly carry heparan sulphate GAG chains and correspond to the two major families of heparan sulphate PGs[51-52]. The HS chains on these heparan sulphate PGs bind a variety of proteins and mediate various physiologically important processes, including blood coagulation, cell adhesion, lipid metabolism and growth factor regulation[53]. Although structurally similar, heparin and heparan sulphate GAGs can often be distinguished through their different sensitivity towards a family of GAG-degrading microbial enzymes, the heparin lyases[54]. 1. Introduction and objectives 17 Figure 13. Major and minor disaccharide repeating units in heparin and heparan sulphate (X=H/SO3-, Y=Ac/SO3-/H). Adapted from Capila and Linhardt 2002[55]. Conformation Experimental structure determination methods such as X-ray crystallography[56], NMR spectroscopy[57], and fluorescence energy-transfer spectroscopy[58] have been applied in studies of carbohydrate conformation, either free or complexed with proteins. While NMR spectroscopy has been extensively used to characterize the dynamics of glycans in solution[59], interglycosidic linkage conformations are notoriously difficult to determine by NMR spectroscopy because of the paucity of Nuclear Overhauser Effects (NOEs)[60], the uncertainties in the Karplus-type equations employed to interpret scalar J-coupling constants[61], and the potential for the linkage to populate multiple rotamer states[62]. Moreover, NMR techniques employed to determine the structural properties of polysaccharides or protein–carbohydrate complexes are limited by molecular weight constraints. Alternatively, X-ray crystallography can be a powerful source of structural information. However, the presence of multiple glycoforms often prevents crystallization of glycoproteins, and the inherent flexibility of oligosaccharides is the presumed reason for the notable absence of X-ray structures for any but the smallest systems. Theoretical methods, such as Monte Carlo and molecular dynamics (MD) simulations, are employed increasingly to augment the experimental approaches in determining the conformational properties of carbohydrates, and biomolecules in general. The level of interest in applying classical simulations to oligosaccharides arises from experimental limitations and is demonstrated by the numerous force fields and parameter sets that have been derived for carbohydrates[63]. 1. Introduction and objectives 18 The conformation of heparin has been extensively investigated since its discovery[45a, 64]. Owing to its nature and topology, i.e., a rigid helix with a complete turn every four residues[65], discontinuous interactions with the same side of a protein surface are expected, grouping each three contiguous sulphate groups on opposite sides (figure 14)[66]. Additionally, the number and distribution of sulphate groups should play some role in the specificity of the interaction[67]. Another structural key aspect of the heparin or HS structure is that while it is very rigid from the backbone perspective (global conformation), at the same time it is quite flexible at the local level, i.e., when the conformational equilibrium of the iduronate ring is considered. Figure 14. Molecular conformation of heparin determined by NMR and molecular modelling from a dodecasaccharide representative of the regular region of heparin, with all the iduronate residues adopting either the 1C4 (top) or the 2SO (bottom) puckering[65]. Each three contiguous sulphate groups are marked within dashed-lined circles. L-iduronic acid, biosynthesized in the polymeric form through a single epimerization at the C5 position of D-glucuronic acid, confers unique properties to iduronate-containing biomolecules. In the manner of most of L-hexopyranoses, it could be expected L-IdoA residue to adopt a 1C4 chair conformation as its sole most stable conformer. However, more than one conformation is accessible[68] in solution and, furthermore, the equilibrium between them can be modulated[69]. This unique feature of the iduronate ring is related to its ability to adopt several conformations of comparable energies[66b, 70], which explains the particularly good ability of iduronate-containing GAGs to control the activity of proteins such as chemokines, growth factors or blood coagulation enzymes[71]. Since changes in the ring conformation alter both the dihedral angles between vicinal hydrogen atoms and the distances between them, one can employ NMR spectroscopy to track such ring puckers by monitoring the spin-spin vicinal coupling constants (3JHH) and the proton-proton NOEs. Thus, over the past decades, extensive solution NMR experiments[10, 72] as well as theoretical calculations[73] have been performed to better understand the conformational flexibility of the iduronate ring. As a result, the picture describing the L-IdoA ring puckering, initially thought to be depicted by the equilibrium 1. Introduction and objectives 19 between the 1C4 and 4C1 chair conformers, have been further completed when evidences about the 2S0 skew-boat pucker playing a critical role in the control of blood coagulation appeared[74]. This 2SO conformer (or pucker) can be easily identified by NMR because it produces an intra-ring exclusive NOE cross-peak corresponding to the close contact between H2 and H5 protons (figure 15), which are not at an NOE distance either in the 1C4 (figure 15) nor in the 4C1 chairs. Measured 3JHH couplings of iduronate as a monosaccharide[75], or as the non-reducing terminal of oligosaccharides[72, 75], indicated a mixture of both the 4C1 and 1C4 chairs with an additional contribution of the 2SO skew-boat. Specifically, an internal L-IdoA2S ring in heparin-like molecules shows a conformational equilibrium between the 1C4 chair and the 2SO skew-boat puckers, with a negligible or non-existent population of the 4C1 chair[75-76]. In this regard, it has been previously described that the 1C4 and 2SO conformers may interconvert with little changes to the geometry of the glycosidic linkages to adjacent residues in the polysaccharide chain[66b, 77] (C4-O4 and C1-O1 bonds present similar orientations in both forms), and thus anticipating the idea, later demonstrated, that the global conformation of the oligoor polysaccharide is independent of the iduronate conformational plasticity[77-78]. Concerning the iduronate flexibility, the existence of fast pseudorotational interconversion along the boats and skew-boats conformational space must also be considered (figure 3)[66b, 79]. The skew-boat 2SO occupancy in the L-IdoA and L-IdoA2S conformational equilibria is almost certainly biologically significant (inhibition of the coagulation cascade is thought to be initiated by antithrombin binding heparin with the iduronate residue in 2SO conformation, and synthetic heparins presenting 2SO-biased iduronate analogues are highly potent[80]. Recently, it has been reported the first complete exploration of the low-energy conformations for iduronate ring (L-IdoA and L-IdoA2S) by MD simulations (figure 16)[17]. This study predicted that iduronate residues undergo microsecond puckering equilibrium (1C4-4C1 conformations exchange of the iduronate ring on the microsecond time scale) and that this depends on substitution pattern (L-IdoA 2-O-sulfation stabilizes the 1C4 conformer) and epimerization (C5 epimerization leads to the 4C1 chair)[17]. These observations have been of fundamental importance as almost all historical carbohydrate simulations are sub-microsecond in duration. Furthermore, they have revealed how enzymatic chemical modifications (epimerization and sulphation) fine-tune the free energy landscape of the iduronate ring and thereby mediate protein selectivity. According to the authors, the theoretical free-energy landscape obtained for the iduronate ring (also for the D-GlcA residue; see figure 16), not amenable to experiment, provides a new route for the development of so-needed carbohydrate mimetic biomaterials and pharmaceuticals[17]. 1. Introduction and objectives 20 Figure 15. 3D representation of the 1C4-2SO conformational equilibrium of the iduronate ring. Note the variation of the H2-H5 and H1-H3 distances, outside the NOE range in the 1C4 chair conformer and within the NOE distance in the 2SO skew-boat pucker. Only the H4-H5 distance does not vary between both puckers. Figure 16. One-dimensional free energy (G) landscapes derived from equilibrium populations of each L-IdoA, L-IdoA2S and D-GlcA monosaccharide MD simulation. Source: Sattelle et al. 2010[17]. As hypothesis, it has been widely accepted that the driving force that determines the conformational equilibrium of the iduronate residue in heparin oligosaccharides is the electrostatic repulsion between anionic charges on neighbouring residues[69a]. Particularly for the regular region of heparin (iduronate in the L-IdoA2S form), due to the presence of three contiguous sulphate groups, NSO3-(GlcN)-2OSO3- (IdoA)-6OSO3-(GlcN), aligned on the same side of the helix (figure 14), it is tempting to speculate that the electrostatic “stress” existing on that part of the molecule might be at the origin of the singular conformational plasticity of the iduronate ring. Whereas the flexibility of the backbones of polysaccharides is usually associated only with rotation of the monosaccharide residues around the glycosidic bonds, the extra-flexibility induced by the presence of an equilibrium of two or more conformations of monosaccharide residues is a peculiar characteristic of iduronate-containing 1. Introduction and objectives 21 glycosaminoglycans which could contribute to their binding properties and biological “versatility”. These features contrast with the poor binding and biological properties of other glycosaminoglycans having approximately the same degree of sulphation and molecuIar weight, but the more rigid glucuronic acid as the major uronic acid. Different factors act as modulators of the conformational equilibrium of the iduronate ring. Thus, the type of counterion may shift the equilibrium towards one of the puckers by interactions in specific sites of the polysaccharide[75, 81]. This is the case of Ca2+ ions, which drive the equilibrium towards the 1C4 chair in heparin-like oligosaccharides[81b]. Furthermore, and more interesting in the framework of the present study, this equilibrium is highly sensitive to intramolecular factors such as the 2-O-sulfation of the iduronate residue and the sulfation pattern of the adjacent GlcN rings[68]. Thus, for internal L-IdoA or L-IdoA2S residues in heparin and heparan sulphate sequences, though only the 1C4 and 2SO puckers participate in the conformational equilibrium[75-76], this is displaced towards the 2SO conformation when IdoA2S is 4-O-substituted with a 3-O-sulphated GlcNS residue (54-69% according to 3JHH at room temperature[75]). On the other hand, the IdoA2S conformational equilibrium is shifted towards the 1C4 chair pucker as long as it is at the non-reducing terminal[75]. In addition, in the case of a terminal non-sulphated L-IdoA residue, the 4C1 form also contributes significantly to the equilibrium[75]. Regarding the global conformation of heparin, which is determined by the geometry of its glycosidic linkages (Φ and Ψ torsions), rigid and flexible behaviours have been observed for the GlcN-IdoA and IdoA-GlcN (figure 8) linkages, respectively[82]. In addition, the conformational space sampled by the Φ torsion is restricted by the conditions imposed by the exo-anomeric effect, and so it commonly presents a narrow distribution of values around the syn geometry (figures 8 and 17). On the other hand, the Ψ torsion, not limited by the anomeric effect, behaves rigidly (syn) within the GlcN-IdoA linkages, but may provide a significant flexibility to IdoA-GlcN linkages[82]. This additional flexibility of the Ψ torsion is characterized by the appearance of anti-ψ conformations (±180º) together with the more common syn-ψ disposition (figure 17). From the NMR viewpoint, while the syn-ψ geometry can be identified by the presence of the H1’-H4 and H1’-H6 NOE cross-peaks, the anti-ψ conformation produces the H1’-H3, H1’-H5 and H5’-H6 exclusive NOEs (figure 17). This allows to experimentally identify the existence of conformational flexibility around the IdoA-GlcN linkages in solution. 1. Introduction and objectives 22 Figure 17. Interglycosidic NOE distances in heparin-like fragments with the iduronate ring adopting both the 1C4 (left) and 2SO (right) conformations. Note the different of set NOEs observed for the minor and major anti-Ψ (bottom) and syn-Ψ (top) conformations, respectively, around the flexible IdoA-GlcN glycosidic linkage. Reducing-end GlcN, IdoA2S and the non-reducing rings are called a, b and c, respectively. 1.2.2 Heparin binding proteins: the acidic Fibroblast Growth Factor (FGF-1) case FGF-1 is a member of the Fibroblast Growth Factor family that interacts with heparin/heparan sulphate (HEP/HS) polysaccharides and the membrane receptors FGFRs, thus triggering a signal that leads to different cellular essential functions such as the regulation of embryonic development, homeostasis and regenerative disorders[83]. The formation of a FGF1-HEP/HS-FGFR2 ternary complex is the key step for the activation of the FGF signalling pathway. Dimerization of the receptors and subsequent autophosphorylation activates a mitogenic response through an enzymatic cascade[40, 84]. Previously, our group addressed the study of the factors that govern the activation of FGF-1 via heparin binding by measuring the induced mitogenic activities of synthetic oligosaccharides[67, 85]. As the helical structure of heparin drives the sulphate groups towards opposite sides of its molecular axis, the multimerization of FGF molecules might be, at first, favoured [86]. In fact, there is a crystallographic structure of the complex between heparin (hexasaccharide) and FGF-1 (PDB code 1AMX)[87], which corresponds to a FGF-1 dimer linked by a regular heparin chain. 1. Introduction and objectives 29 Table 2. Main features of some C-type lectins produced by DCs and LCs. C-type lectin Type Amino acids Production Ligand/s Function/s Key antibodies MMR (CD206) I 1456 DCs, LCs, Mo, M, DMECs Man, Fuc, sLex Antigen uptake[103] MG38, anti-human[104] DEC-205 (CD205) I 1722 DCs, LCs, actDCs, thymic ECs ? Antigen uptake[105] Dectin 1 II 247 DCs, LCs β-glucan[106] T-cell interaction[107] Dectin 2 II 209 DCs, LCs ? Antigen uptake[108] Langerin (CD207) II 328 LCs Man, Glu, Gal6S, Fuc, heparin* Formation of Birbeck granules[109], HIV-1 barrier[110] DCGM4, anti-human DC-SIGN (CD209) II 404 DCs HIV-1 (gp-120), SIV, mannan, ICAM-2, ICAM-3 T-cell interaction[111], HIV-1 pathology[112], migration[113], antigen uptake AZN-D1, anti-human Abbreviations: actDCs, activated dendritic cells; DMECS, dermal microvascular endothelial cells; Mo, monocytes; M, macrophages; sLex, sialyl Lewis X; ECs, endothelial cells. *To date, heparin is the first Langerin ligand reported to interact in a calcium-independent manner[114]. 1.3.1 Structure of C-type lectin receptors (CLRs) As we have introduced above, the common feature of all CLRs is that they possess at least one compact globular structure with a characteristic fold designated “C-type lectin-like fold” or “C-type lectin-like domain (CTLD)” that is unusual to any other known proteins. For the majority of CLRs that function as PRRs, the CTLDs bind sugars, usually in Ca2+-dependent manner, and therefore this domain is commonly called a “carbohydrate recognition domain” (CRD). 1. Introduction and objectives 30 All CRDs possess a characteristic “double-loop” fold. The whole domain can be regarded as a loop with two flanking α helices (α1 and α2) and two antiparallel β-sheets: Nand C-terminal β strands β1 and β5 constitute the basal β-sheet, and the top β-sheet is formed by strands β2, β3, and β4 (figure 22). The long loop region enters and exits the core domain at the same location, and is involved in Ca2+-dependent carbohydrate binding and, for some CRDs, in domain-swapping dimerization. Four highly conserved cysteine residues form two disulphide bridges at the bases of the loops: C1-C4 bridge links α1 and β5, and C2-C3 bridges β3 strand and a loop upstream the β5 strand (figure 22). The long loop region among different CRDs varies, and those that possess it are designated “canonical”, while those that lack it are called “compact”. The presence or absence of a short extension at N-terminus, a β1-hairpin, further subdivides CRDs to long or short forms, respectively. Two additional cysteine residues at the beginning of CRDs sequence are characteristic for the long form CRDs. The corresponding disulphide bridge C0-C0’ stabilizes the β-hairpin (figure 22). There may be up to four Ca2+-binding sites in the CRDs, and their occupancy depends on the sequence of a particular CRD. Sites 1, 2, and 3 are located within the long loop, while the fourth site (Ca-4) participates in the salt bridge formation between helix α2 and β1/β5 sheet[115] (figure 22). Figure 22. Cartoon representation of a common CRD structure in CLRs (DC-SIGN CRD; PDB code 1K9I). The long loop is shown in blue and the disulphide bridges in yellow sticks. The Ca2+-binding sites 1, 2 and 3 existing in this lectin are shown as cyan spheres, while the location of the fourth Ca2+-site (Ca-4, absent in this particular structure) has been drawn and appears as a cyan circle. Source: Doctoral Thesis of Ieva Sutkeviciute[116]. 1. Introduction and objectives 31 1.3.2 Structural features of glycan binding to CRDs Amino acid residues with carbonyl side chains involved in Ca2+ coordination in site 2 form two characteristic motifs in the CLR sequence that, together with the calcium ion itself, are directly involved in monosaccharide binding (figure 23A). The first group of residues, the “EPN motif, is contributed by the long loop region and contains two residues with carbonyl side chains separated by a proline in cis conformation (figure 23B). The carbonyl side chains provide two Ca2+-coordination bonds, form hydrogen bonds with the monosaccharide and determine binding specificity. The cisproline is highly conserved and maintains the backbone conformation that brings the adjacent carbonyl side chains into the positions required for Ca2+ coordination. The second group of residues, the “WND motif” (figure 23, B and C), is contributed by the β4 strand. Although only the asparagine and aspartate residues of this motif are involved in Ca2+-coordination, the tryptophan amino acid immediately preceding them is a highly conserved contributor to the hydrophobic core[117] and is a useful landmark for detecting the motif in a sequence. In the MBP-A structure shown in figure 23B, ASN205 and ASP206 provide three Ca2+-coordination bonds (two from the side chains, one from the backbone carbonyl of Asp) and also form hydrogen bonds with the sugar. Additionally, another carbonyl side chain is involved in site 2 formation, which belongs to the residue preceding the second conserved cysteine at the end of the long loop region (Glu193 in MBP-A; see figure 23B) and forms one coordination bond with the Ca2+ ion. The overall network of the hydrogen-bond donors and acceptors in binding site 2 (figure 23) determines the binding orientation of the carbohydrate and also which hydroxyls of the carbohydrate it can accept, i.e. the monosaccharide specificity. The EPN motif has a configuration that accommodates mannose-type monosaccharides (figure 23B), whereas the QPD motif determines specificity for galactose-like monosaccharides (figure 23C). In both of these motifs the cis configuration of the two carbonyl sidechains separated by proline is crucial for Ca2+-coordination and sugar binding. Besides the restrictions imposed by the H-bond network, other structural elements in the binding sites introduce selectivity to particular ligands within the mannose or galactose groups. As no other Ca2+-binding site except for site 2 is known to be involved in sugar binding, and as the site 2 motifs can be confidently detected in the sequence, it is common in the literature to associate the predicted Ca2+-dependent carbohydrate binding properties of an uncharacterized sequence with the presence of these motifs[118]. Although this is a useful simplification, it should be noted that the absence of the motifs associated with Ca2+-binding site 2 does not indicate that the CTLD is incapable of binding Ca2+, as there are two independent sites (1 and 4). Also, the presence of these motifs does not guarantee lectin activity for the CLR, as there are numerous examples of C-type lectins that contain the conserved motifs but are not known to bind monosaccharides. The other three Ca2+ binding sites (1, 3 and 4) play a structural stabilization role, as removal of Ca2+ increases susceptibility to proteolysis and changes physical properties of the domain. Ca2+ binding site 2 is also important for structural stability of the domain. It has been shown that pH-induced loss of 1. Introduction and objectives 32 Ca2+ causes the destabilization of the loops, which has an important physiological role for CRDs of endocytic receptors as internalization of ligand-bound receptor to acidic lysosomes and consequent Ca2+ loss leads to the release of the ligand for further processing, while receptor is recycled to the cell surface[115]. Figure 23. Ca2+-dependent monosaccharide binding by CLRs. (A) Schematic representation of a Ca2+-hexose-CLR complex. Two hydroxyl oxygens and the ring of the hexose are shown. The Ca2+ atom is shown as a large grey sphere, and oxygens as empty circles and ovals. Protein groups that act as hydrogen donors and acceptors are not shown. Black arrows show the direction of hydrogen bonds in “mannose-specific” CLRs, while red arrows indicate opposite directions in “galactose-specific” CLRs. (B) Mannose residue bound to MBP-A CRD (PDB code 2MSB). (C) GalNAc residue bound to MBP-A mutant CRD (PDB code 1BCJ). In (B) and (C), the coordination bonds are orange, and the hydrogen bonds where sugar hydroxyl acts as acceptor or donor are red or blue, respectively. The Ca2+ ions are shown as a cyan spheres. Adapted from the Doctoral Thesis of Ieva Sutkeviciute[116]. 1. Introduction and objectives 33 1.4 Langerin: a natural barrier to HIV-1 infection Langerin, a transmembrane type II C-type lectin with carbohydrate specificity, is almost exclusively expressed on epidermal Langerhans cells (LGs), i.e., a subset of dendritic cells (DCs) with very singular structural features. Due to its placement at the epithelium level, these cells were reported a few year ago to constitute the first recognition barrier to HIV-1 particles[110]. Interestingly, instead of the common antigen-presenting cell role of DCs, LCs internalize and subsequently degrade HIV-1 virions upon Langerin sequestration of the virus[110]. Thus, Langerin has been presented as as a natural barrier to HIV-1 transmission by LCs[110]. The host cell infection by HIV starts by binding of HIV envelope proteins “Env” (i.e, a trimmer of gp120 and gp41 heterodimers, where gp41 initially is hidden) to its primary receptors CD4+ T lymphocytes, a member of the immunoglobulin superfamily that enhances T-cell receptor (TCR)-mediated signalling (figure 24), being an absolute requirement for the infection. CD4+ T lymphocytes also express chemokine receptors CCR5 and CXCR4 that are exploited by HIV to enter the cells, hence also called HIV co-receptors (figure 24). Env binding event induces rearrangements in its gp120 subunit, which ultimately result in V3 loop repositioning and bridging sheet exposure that are essential for co-receptor engagement. The attachment of the virion can be relatively nonspecific, so that, for instance, HIV Env can interact with negatively charged cell-surface heparan sulphate proteoglycans[119]. More specific adhesion includes interactions between the envelope protein and α4β7 integrin[120] or CLRs such as DC-SIGN[112] and Langerin[110]. Either way of adhesion has been proposed to bring HIV envelope into close proximity with the host CD4 and a co-receptor, leading to the fusion of viral and target cell membranes[121] (figure 24). Subsequent binding to the co-receptors, CCR5 or CXCR4 depending on the virus strain R5 or R4, triggers the membrane fusion potential of Env, and usually is followed by the “surfing” of the virus particle to the site where productive membrane fusion may occur[122]. It is thought that HIV might usurp the host cell machinery to reach cell surface sites where membrane fusion can occur[123]. Besides, HIV may need to be endocytosed by the host cell for productive membrane fusion to occur[124]. Upon formation of Env-gp41 complex, the co-receptor undergoes conformational changes that expose its hydrophobic fusion peptide (figure 24), which then inserts into the host cell membrane and folds to form a six-helix bundle. The latter is the driving force that brings the opposing membranes into close proximity, resulting in the formation of a fusion pore[125] (figure 24). Once the virus enters the cell, it can start its replication and productive infection, which ultimately lead to the depletion of CD4+ T lymphocytes in the body, and thus the acquired immunodeficiency syndrome (AIDS). 1. Introduction and objectives 34 Figure 24. Schematic representation of HIV entry to a target cell. Source: Abbas et al. 2011[126]. Importantly, although the major target of HIV is CD4+ T cells, earlier studies have shown that DCs are crucial for HIV-1 infection enhancement and dissemination in mucosa, in the case of sexual HIV transmission, since the virus hijacks them to achieve the productive infection of the CD4+ T cells and the burst of the disease[127]. The stromal DCs that express a C-type lectin DC-SIGN (DC-Specific ICAM3 Grabbing Non-integrin) and reside in mucosa of vagina and ectocervix, have been repeatedly reported to be exploited by HIV-1 to enhance its infectivity of T cells. There exist several HIV infection mechanisms (figure 25) in which it is noticeable that DC-SIGN has a very important role in DC-mediated HIV transmission enhancement[112]. After virion binding to DC-SIGN, HIV can be endocytosed into DCs and, subsequently, the intact viral particles can be stored either in multivesicular bodies[112, 128], as integrated provirus (following productive infection of DCs) or as DC-SIGN-bound virions on the cell surface and protected from degradation[129] (figure 25). During infection there is an accumulation of intact viral particles on DC side while HIV receptors (CD4, CCR5) are presented on CD4+ T cell side. This situation greatly facilitates HIV-1 transfer from DCs to T cells. Indeed, it has been demonstrated that blocking DC-SIGN prevents HIV-1 binding and subsequent trans infection of CD4+ T cells[130]. While trans infection is responsible for the early stage of infection (24h after HIV exposure), there exists a different pathway of DC-SIGN-bound HIV transmission to T lymphocytes. This pathway is involved in long-term HIV transfer (72h after exposure) and it occurs as a cis infection of DCs by transfer of DC-SIGN-bound virus to canonical HIV entry receptors, CD4 and CCR5, which leads to productive infection of DCs and, in turn, to the presentation of increased viral load to the T cells[131]. 1. Introduction and objectives 35 Figure 25. Importance of DCs-T-cells interactions and DC-SIGN implication on HIV-1 transmission. Source: Hladik and McElrath 2008[132]. 1.4.1 Novel role of epithelial LCs and their transmembrane protein Langerin in HIV infection LCs constitute a subset of DCs, whose role is not completely clear, that are located in epithelium of mucosal tissues and epidermis[133] (figure 26A). Although initially they were assumed to function as antigen-presenting cells (APCs), like dermal or stromal DCs[134], the accumulating evidence of their fail to present antigens from various viruses and parasites to T cells activation, and the observation that LCs induce T regulatory cells, supported the hypothesis that these cells could have an inmunosuppresive tolerogenic role[135] . Notably, Langerhans cells are the only epidermal cells to constitutively express major histocompatibility complex class II molecules[136], CD1a molecules[133a], and Langerin[109] at their cell surface. LCs play a key role in the induction of immune responses against invading pathogens by capturing and processing foreign antigens and migrating to draining lymph nodes to present processed antigens to T cells[137]. A characteristic hallmark of LCs is the unique presence (in the cytoplasm) of a tennis racquetor rod-shaped membranous structures, known as Birbeck granules (BGs)[138] (figure 26, B and C), that are thought to be part of the endosomal 1. Introduction and objectives 36 recycling pathway (they are subdomains of the endosomal compartment)[109]. The protein Langerin has been shown to be the main component and responsible for the formation of these Birbeck granules[109]. Figure 26. (A) Langerhans cells (green) in the foreskin surrounded by nuclei of epithelial cells and other cell types (red). (B) Birbeck granules, site of HIV processing in LCs. (C) HIV-1 particles (black small circles) internalized into Birbeck granules of LCs (and mutant muLCs) upon capture by Langerin. Sources: (A, B) Schwartz 2007[139]. (C) de Witte et al. 2007[110]. As DC-SIGN, Langerin is a transmembrane type II C-type lectin, almost exclusively expressed in humans by epidermal LCs but also present on dermal CD103+ DCs and lymph node resident CD8+ DCs[109, 140]. It contains one calcium-dependent carbohydrate recognition domain with a short cytoplasmic tail with a proline rich motif[141], forms a trimer on the cell surface and, upon crosslinking with either a cell-bound or a soluble ligand, it induces the formation of BGs. As a C-type lectin, Langerin plays a role in pathogen recognition. However, only few pathogens have been demonstrated to interact with Langerin. Both HIV-1[110, 142] and Mycobacteria leprae[143] have been identified as pathogens that interact with human Langerin, whereas murine Langerin was shown to bind to Candida albicans[144]. LCs express receptors including CD4, CCR5 and Langerin[145]. It must be noted that, different from the other DCs (antigen-presenting cells), LCs play a significant preventive role in HIV infection process. In particular, it has been reported that epidermal LCs expressing Langerin efficiently bind HIV virions, which in turn are directed to Birbeck granules for degradation[110] (figure 27). Moreover, a recent study has shown that vaginal LCs may have low or no expression of Langerin, and thus they are 1. Introduction and objectives 37 susceptible to HIV infection[146]. These experimental evidences show an important protective function of Langerin in HIV invasion process. LCs reside in the epithelium (particularly in skin epidermis and mucosal epithelium) while DCs are placed at the sub-epithelium level (figure 27). Thus, LCs are the first cells encountered by HIV-1 particles, which are captured by the LCs membrane receptor Langerin, resulting in viral clearance and inhibition of HIV-1 transmission across the mucosal layer (figure 27). The disruption of the epithelial barrier through trauma or ulcerations enables HIV-1 to circumvent the LCs barrier and reach the DC-SIGN-expressing DCs, which act as antigen-presenting cells, thus mediating HIV-1 transmission to T cells (figure 27). On the other hand, the inhibition of Langerin function by drugs, co-infections or mutations, gives rise to the efficient transmission of HIV-1 to T cells but, in this case, through infection of LCs (figure 27). Hence, two C-type lectins such as DC-SIGN and Langerin present very different roles in HIV invasion despite the high homology of their CRDs and the overlap of ligand specificities[147]. While DC-SIGN promotes HIV transmission and infection, Langerin fights against HIV invasion. Based on the carbohydrate recognition specificity of Langerin CRD for mannose[148], fucose[148], N-acetyl-glucosamine[148] and 6-O-sulphate-galactose[148b] monosaccharides (see next section), it is likely that Langerin has a broader specificity for pathogens than it has been thought so far. Also, it has been very recently demonstrated that the trimeric extracellular domain of Langerin (LgECD) bind GAGs (specifically heparin) in a novel non-calcium dependent binding site[114] (see next section). Figure 27. Schematic representation of the implication of LCs on HIV-1 transmission. Note that only LCs are present at the first or more external pathogen-recognition level (epidermis), therefore being the first entity of the immune system to encounter the virus. Source: de Witte et al. 2007[110]. 1. Introduction and objectives 38 1.4.2 Langerin 3D structure and known ligands Langerin is formed by 328 amino acids (37,5 kDa) and has an overall molecular structure similar to DC-SIGN and other CLRs. It consists of a N-terminal short cytoplasmic domain, a unique transmembrane domain, and a large extracellular domain (ECD) subdivided into a neck domain and a C-terminal carbohydrate-recognition domain, CRD, (figure 28)[109]. The CRD of Langerin contains the EPN motif characteristic of mannose-type specificity (also presents the WND motif; see figure 29). On the other hand, the cytoplasmic domain of Langerin presents a proline-rich signalling motif (WPREPPP), which could function as a docking site for signal transduction proteins, and indeed, it was demonstrated to be important for Langerin intracellular targeting[149]. Like many other CLRs, Langerin exists as an oligomer (active form), forming trimers stabilized by a coiled-coil of helices in the neck region (figure 28). Trimer formation is essential for binding to oligosaccharide ligands because, as is typical for C-type CRDs, the CRD of Langerin has only low affinity for monosaccharides[148a, 150]. In addition, it should be noted that oligomerization of C-type lectins is also important for determining selectivity for particular oligosaccharide structures. For example, in serum mannose-binding protein, three CRDs in the trimeric unit are held in a fixed position via interactions between the CRDs and an helical neck region so that the binding sites are arranged to interact with arrays of sugars in polysaccharides of bacterial cell walls, but not with mammalian high mannose-type oligosaccharides[151]. Although both DC-SIGN and Langerin are C-type II lectins with similar sequences and structures of their CRDs, several important differences exist between these two CRDs (figure 30). The overall structure of Langerin CRD is maintained by two conserved disulphide bridges, in contrast to the four S-S bonds in DC-SIGN. Unlike DC-SIGN, Langerin CRD has only one Ca2+ ion at the calcium site 2 (conventional sugar binding site; see figure 22), and the lack of other Ca2+ ions in sites 1 and 3 (figure 22) might be at the origin for the high flexibility of the β2-β2´ loop comprising residues 258-262, which also leads to the formation of a large groove specific to Langerin structure. However, the most significant difference between these two lectins is present in their sugar-binding site topologies. A unique feature of Langerin CRD is the presence of two lysine residues (LYS299 and LYS313) within the sugar binding site that provides it the basic character that allows Langerin to accommodate sulphated sugars in the binding site. On the other hand, while PHE313 is an important side chain in DC-SIGN sugar binding-site as it forms stacking interactions with the sugar ring[152], PHE315 residue is not placed at the appropriate position in Langerin as to participate in such sugar interaction (figure 30, upper). The overall topology of Langerin CRD composes only a small binding site strongly constrained by LYS299, while DC-SIGN has a potential to adapt more extended oligosaccharides with a widely open binding site (figure 30, down), and thus less extensive secondary contacts take place with Langerin CRD than with DC-SIGN CRD. Also, the linkage between CRD and neck domains of DC-SIGN and Langerin are also different. Thus, while the flexibility in the former is retained, Langerin Chapter 2 Techniques and tools 2. Techniques and tools 47 As it has been outlined in Chapter 1, the physiological processes taking place in cells are the result of highly regulated intermolecular protein–protein and protein–ligand interactions. In particular, the binding of low molecular weight molecules (ligands) to macromolecules such as proteins plays a major role in the regulation of biological processes, e.g., signal transmission and cellular metabolism. The analysis of protein–ligand interactions is crucial, not only for understanding the regulation of biological functions, but also for designing novel bioactive molecules that modulate protein function or inhibit protein-ligand interactions[158]. With respect to carbohydrates, the elucidation of the 3D structures and dynamics properties of oligosaccharides and glycoconjugates, both in the free state and bound to proteins, is a prerequisite for a better understanding of the molecular basis of their associations and interactions, and the relationships between structures and functions, which are involved in the biochemistry of recognition processes and the subsequent rational design of carbohydrate-derived drugs. These have been claimed to be the main challenges in structural glycoscience[159] and many efforts in this direction still have to be done. A large variety of biophysical techniques have been developed to characterize protein–ligand complexes, e.g., surface plasmon resonance (SPR), isothermal titration calorimetry (ITC), fluorescence polarization assay (FP), fluorescence resonance energy transfer (FRET), enzyme-linked immunosorbent assay (ELISA), differential scanning fluorimetry (DSF), microscale thermophoresis (MST) and electrospray ionization mass spectrometry (ESI-MS) to cite the most popular methods[158d, 160]. SPR and MST techniques permit to obtain kinetic interaction parameters, whereas ITC measures the thermodynamic properties of binding in solution. In silico approaches have been also applied to search for ligands for a protein target (virtual screening) or to propose 3D models of protein–ligand complexes (docking calculations)[161], whereas X-ray crystallography and nuclear magnetic resonance (NMR) spectroscopy are both experimental techniques for resolving atomic structures[162]. Since the scientific interests of our group are focused on the structural studies of protein-carbohydrate complexes at atomic resolution, it is the X-ray crystallography and NMR spectroscopy the existing tools to our aim. However, one of the disadvantages of X-ray crystallography applied to carbohydrates is that the oligosaccharides, either in their free form or as part of glycoconjugates, are inherently difficult to crystallize, and structural data from X-ray studies are sparse[163]. Even when succeeding in crystal formation, part or the whole glycan is, in most cases, not observed in the high-resolution electron density map[57d] due to the intrinsic high flexibility of carbohydrates. Furthermore, the experimental assessment of carbohydrate recognition by X-ray crystallography is impeded by difficulties of co-crystallizing proteins and carbohydrates. To overcome this limitations, it has been tried, for instance, to stretch the polysaccharide into an oriented fibre[164]. Also, it has been employed electron diffraction to study very small crystals, or needles, that can be obtained from polysaccharides[165]. In any case, the amount of data collected to date is that small that building a model by molecular mechanics is necessary to resolve the 3D structure. 2. Techniques and tools 48 On the other hand, high-field NMR spectroscopy in solution state is one of the most important techniques for probing intermolecular interactions. NMR spectroscopy detects and reveals protein-ligand interactions with a large range of affinities, and it is widely used in pharmaceutical research to identify hits from compound library screening in drug discovery[166]. Protein–ligand complexes are analysed using the so-called protein-observed and ligand-observed NMR experiments in which the NMR parameters of the protein and the ligand, respectively, are compared in their free and bound states[166]. In particular, ligand-observed methods are not limited by the protein molecular size and therefore have great applicability for analysing protein–ligand interactions. The use of these NMR techniques has considerably expanded in recent years, both in chemical biology and in drug discovery. In protein-observed methods, the chemical shift perturbations of the protein resonances observed upon ligand addition are identified to localize the ligand binding site. This enables one to immediately distinguish specific from non-specific binding. The 3D structure of the protein-ligand complex can be resolved via heteronuclear experiments performed on isotopically labelled (13C, 15 N, 2H) protein samples. The structure resolution requires molecular dynamics calculations with experimental NMR restraints resulting from chemical shifts, scalar couplings, nuclear Overhauser effects (NOEs), paramagnetic interactions or residual dipolar couplings[162a, 167]. The major drawbacks are the experimental time and the need for a highly stable and soluble protein. In addition, these methods are limited in routine practice to proteins with low molecular masses (less than 30 kDa) to avoid great effort with regard to both labelling strategies and resonance assignment. NMR parameters such as transverse, longitudinal, and cross-relaxation rates strongly depend on the molecular rotational correlation time τc, which is directly related to the molecular weight. Ligand-based NMR experiments rely on the modification of such size-sensitive NMR parameters for the ligand in the presence of a protein receptor[166b, 168]. Considering a diffusion controlled protein-ligand binding of weak to moderate affinity (dissociation constant, KD, typically between 108 and 10-3 M), the association-dissociation process is fast within the chemical shift time-scale, so that the NMR parameters observed are a simple population-weighted average between the free and bound states. In contrast to protein-observed experiments, ligand-observation is more sensitive with larger receptors and do not require the use of isotopically labelled proteins. Ligand-based methods can be used for the detection of interactions and the measurement of protein–ligand affinities, and can also provide pertinent structural information on the protein–ligand complexes. 2. Techniques and tools 49 2.1 Structure by NMR: the Nuclear Overhauser Effect (NOE) 2.1.1 Origin of the NOE An accurate definition of the Nuclear Overhauser Effect (NOE) is the change in intensity of one resonance when the spin transitions of a dipolarly coupled nucleus are somehow perturbed from their equilibrium populations. This perturbation is achieved by either saturating a resonance, i.e., equalising the spin population differences across the corresponding transitions (then it is called steady-state NOE), or inverting it by reversing the population differences across the transitions (transient NOE). Thus, the magnitude of the NOE observed for spin I when spin S is perturbed ( { }) is expressed as the percentage of relative intensity change between the equilibrium intensity (I0) and that in the presence of the NOE (I), so that I{S}=I-I0 I0 100 Eq. 1 The intensity changes caused by NOE can be either positive (I>I0) or negative (I<I0) depending on the motional properties of the molecule and the signs of the magnetogyric ratios of the spins involved. To facilitate the understanding of the origin of the NOE, we will consider a system only formed by two homonuclear spin-1/2 nuclei of 1H (positive magnetogyric ratio), I and S, contained in a rigid molecule that tumbles isotropically in solution, i.e., it does not show any preferential axis about which to rotate. In this idealistic system both protons are not scalarly coupled (JIS=0) but they are enough close in space as to share dipolar coupling, this is, magnetic interaction through space between two spins such that both of them are able to sense the presence of the other dipolar-coupled partner. Therefore, upon selective saturation of 1H nucleus S, the spin populations of nucleus I will be also perturbed and the system will try to come back to the initial equilibrium situation. Although equilibrium recovery takes place by different relaxation mechanisms it is only the cross-relaxation pathways, characterized by the W0 (zero-quantum) and W2 (double-quantum) transition probabilities (or rates), those responsible for the NOE development (figure 1A). It is important to note that the W0 and W2 cross-relaxation pathways always compete with one another, with the dominant mechanism dictating the sign of the observed NOE and being dependent on the reorientational dynamics properties of the molecule (figure 1B). Thus, a small molecule performing a rapid tumbling in solution, which corresponds to a short correlation time τc, will generate fluctuating local magnetic fields of high frequency. On the other hand, a macromolecule, featuring slow tumbling rates (long τc) will give rise to low-frequency magnetic fields. These fluctuating local fields are the responsible for inducing cross-relaxation when their oscillation frequencies correspond to the W0 and/or W2 transitions. Furthermore, those of low-frequency are significantly more efficient than the high-frequency magnetic fields in activating the cross-relaxation pathways, due to the 2. Techniques and tools 50 different spectral density function (J()) 2 that describe the fast, intermediate and slow motions (figure 1B). Therefore, for slow tumbling the lower energy W0 process predominates such that a negative NOE (I<I0) is efficiently developed (figures 1B and 2). On the contrary, W2 (higher energy) is the dominant mechanism acting during cross-relaxation of fast tumbling molecules, so that a positive NOE (I>I0) is observed (with lower intensity compared to negative NOE; see figure 2). For and intermediate tumbling rate (usually medium-size molecules), the NOE will be either positive or negative depending on the dominant cross-relaxation process. Also, it has to be noted that there are two limits of motion in terms of the magnitude of the NOE developed. Thus, the rapid tumbling of a small molecule in a low viscosity solvent will favour the W2 process to a large extent, displaying high positive homonuclear NOEs. This is called extreme narrowing limit. In contrast, the slow molecular tumbling of large molecules in high viscosity solvents will stay in the spin-diffusion limit, characterised by an enormously favoured W0 mechanism, and thus, highly negative homonuclear NOEs. Apart from the cross-relaxation mechanisms, only responsible for the NOE growth, single quantum relaxation pathways (W1) activate to re-establish the equilibrium population differences of the nonsaturated nucleus (I in the simplified model) as soon as the NOE begins to develop, so acting against the NOE build-up. Thus, if W1 relaxation happens to be rather more efficient than W0 and W2 pathways together, the macroscopic magnetization will probably come back to the equilibrium before a measurable NOE is developed and this will not be observed. The NOE therefore results from the balance between distinct competing relaxation pathways, with its sign depending on the W2-W0 difference and its magnitude on the three Wo, W1 and W2 rates (see Eq. 2, derived from the so-called Solomon equation). { } [ - 0 0 ] [ ] Eq. 2 So, the ideal conditions for the NOE to be observed are inefficient W1 processes and efficient W0 or W2 transitions. 2 J() is the frequency distribution of the fluctuating magnetic fields associated with molecular motion. It may be viewed as the probability of finding a component of the motion at a given frequency  (in rad/s). 2. Techniques and tools 51 Figure 1. (A) The six possible transitions in a two-spin system. Note the two cross-relaxation pathways, W0 and W2. (B) Evolution of the spectral density function (J(ω)) as a function of the frequency of motion (ω) in logarithmic scale (ωt corresponds to the frequency of an hypothetical spin transition). Source: Claridge 2009[169]. Figure 2. Variation of the maximum theoretical homonuclear steady-state enhancement ( max), for NOE (black dashed line), TROE (blue bold line) and ROE (red bold line) experiments, in a two-spin system as a function of molecular tumbling rates in logarithmic scale (defined by the dimensionless parameter ω0τc, with ω0 being the spectrometer observation frequency and τc the rotational correlation time). The region of fast motion is the extreme narrowing limit and that of slow motion is the spin-diffusion limit. 2.1.2 Rotating-frame NOE (ROE) The greatest problem associated with NOE experiments is the zero-crossing region around ω0τc ≈ 1 where the conventional (laboratory-frame) NOE observed via steady-state or transient techniques becomes vanishingly small. This typically occurs for mid-sized molecules with masses of around 10002000 daltons, depending on solution conditions and spectrometer frequency. With the increasing interest in larger molecules in many areas of organic chemistry research coupled with the wider 2. Techniques and tools 52 availability of higher field instruments, this is likely to be a region visited ever more frequently by the research chemists’ molecules. Other than altering solution conditions (such as changing temperature) in an attempt to escape from this, the measurement of NOEs in the rotating-frame provides an alternative solution. In this case, the cross-relaxation rate between homonuclear spins is given by an expression that remains positive for all values of τc, and the undeniable benefit of ROEs is, quite simply, that they remain positive for all realistic molecular tumbling rates. For small molecules, the magnitude of the ROE matches that of the transient NOE, whilst for larger molecules it reaches a maximum for homonuclear spins of 68%, but under no circumstances does it become zero (figure 2). Similarly, the NOE and ROE growth rates are identical for small molecules but differ for very large ones. For very large molecules, the ROE therefore grows twice as fast as the NOE, and has opposite sign[170]. In essence, ROEs develop whilst magnetisation is held static in the transverse plane, rather than along the longitudinal axis (hence they are sometimes also referred to as transverse NOEs). To generate the required population disturbance of the source spins, the target resonance is subjected to a selective 180º pulse prior to the non-selective 90º pulse, such that it experiences a net 270º flip and is thus inverted relative to all others. Transverse magnetisation is then “frozen” in the rotating frame by the application of a continuous, low power spin-lock pulse to prevent evolution (in the rotating frame) of chemical shifts. The experiment is more frequently performed as the 2D experiment where it is usually termed ROESY (rotating-frame NOE spectroscopy). The situation during the spin-lock may be viewed as the transverse equivalent of events during the transient NOE mixing time (figure 3). The action of the spin-lock is to maintain the opposing disposition of magnetisation vectors, which would otherwise be lost through differential chemical shift evolution, and so allows the ROE to develop through cross-relaxation in the transverse plane. Spin relaxation here is characterised by a time constant called T1ρ, of very similar magnitude to T2. In utilising the spin-lock, one has effectively replaced the static B0 field of the conventional NOE with the far smaller rf B1 field, and it is this that changes the dynamics of the NOE. Whereas γB0 typically corresponds to frequencies of hundreds of megahertz, γB1 is typically only a few kilohertz, meaning γB1 << γB0 and hence ω1 (the rotating-frame frequencies) << ω0. The consequence of this is that ω1τc << 1 for all realistic values of τc, and all molecules behave as if they are within the extreme narrowing limit. Thus, ROEs are positive, any indirect effects have opposite sign to direct effects and tend to be weak, and saturation transfer can be distinguished by sign from ROEs, regardless of molecular size and dynamics. Against these obvious benefits are a number of experimental problems, principally TOCSY transfers (particularly contributing to strongly coupled systems such as glycans), also occurring during the spin-lock, and signal attenuation from off-resonance effects. An alternative ROESY sequence is the T-ROESY experiment, which is effective at suppressing TOCSY transfer. This is achieved by substituting the low-power continuous wave lock field used in ROESY sequence by a mixing sequence of 180ºx 180º-x or 180ºx 180º-x 360ºx 360º-x 180ºx 180º-x pulses. However, cross-relaxation rates in T-ROESY are an equal mixture of ROE and NOE components, and 2. Techniques and tools 53 moreover, the cross-relaxation rates measured are four times lower than those in conventional ROESY[171] (figure 2). Figure 3. (Up) Scheme of the pulse sequence for observing rotating-frame NOEs. The ROEs develops during the long spin-lock pulse that constitutes the mixing period τm. (Down) Situation of the magnetization during the spin-lock. The rotating-frame NOE experiment can be viewed as the transverse equivalent of the transient NOE experiment. Source: Claridge 2009[169] 2.1.3 Measuring internuclear distances As it has been commented in the previous section, the NOE is the result of the dipolar coupling of two spins which are close in the space to each other, so that this is a distance dependent effect. However, the exact NOE dependency with distance and which NOE experiment we must carried out to obtain accurate distance measurements are not straightforward issues. Thus, for instance, steady-state NOEs cannot readily be translated into internuclear separation because they result from a balance between the influences of all neighbouring spins (only relative distances are obtained from these NOEs). It has also been shown that for molecules that exhibit negative enhancements, steady-state measurements may fail to provide any reliable information of spatial proximity, and here one is forced to consider the kinetics of the NOE. Thus, saturation of the target resonance for periods that are far less than those needed to reach the steady-state would allow some NOE to appear, which is then sampled. Repeating the experiment with progressively incremented saturation periods allows the build-up to be mapped. Owing to the use of shortened saturation periods, the enhancements observed with this method are termed truncated driven NOEs or TOEs. Although once popular, this experimental approach is rather less used nowadays and as such shall be considered no further. The more common approach to obtaining kinetic data is to instantaneously perturb a spin system not by saturation but by inverting the target resonance/s (i.e. inverting the population differences across the corresponding transitions) and then allowing the NOE to develop in the absence of further external interference. In this case, the 2. Techniques and tools 54 NOE is seen initially to build for some time but ultimately fades away as spin relaxation restores the equilibrium condition; these enhancements are thus termed transient NOEs. The measurement of transient NOEs gained widespread popularity, initially in the biochemical community, in the form of the 2D NOESY experiment, which remains an extremely important structural tool in this area and increasingly in the analysis of smaller molecules. The 1D transient NOE experiment, also referred to as 1D NOESY, is also widely used within the chemical community as a gradient-selected sequence capable of providing high-quality NOE spectra. Transient experiments, whether 1D or 2D, are more commonly used qualitatively as ‘single-shot’ techniques, providing an overview of enhancements within a molecule rather than being employed to map the growth of the NOE. Unlike the steady-state enhancements, the transient enhancements are influenced by only a single internuclear separation (r) with r–6 dependence whilst the so-called initial rate approximation is valid. Under this approach, the two cross-relaxing spins initially behave as if they were an isolated spin pair and the growth of the NOE has a linear dependence on mixing time. As longer mixing periods are used, the relaxation of spin I begins to compete with cross-relaxation between I and S, so the build-up curve deviates from linearity and the NOE eventually decays to zero (figure 4). Thus, for the initial rate approximation to be valid, mixing times significantly shorter than the T1 relaxation time of spin I must be used. Only under these conditions is meaningful distance measurement possible. If, on the other hand, the goal is to qualitatively identify through-space correlations, as is more often the case in routine work, mixing periods comparable to T1 provide maximum enhancements. Since transient NOEs develop in the absence of an external radiofrequency field, they tend to give rise to weaker positive values (38% maximum) than the steady-state effects (50% maximum; see figure 2), so careful choice of timing is crucial to the success of transient experiments with fast-tumbling molecules. With respect to the slow-tumbling molecules, owing to the domination of efficient cross-relaxation both the transient and steady-state experiments give rise to a maximum negative NOE of -100%. The primary reason why transient NOE methodology is used, instead of the steady-state experiment, is the complications arising from the existence of spin diffusion for molecules within this regime (slow-tumbling). Thus, in this case there is also a demand for the use of short mixing periods. 2. Techniques and tools 61 To avoid it, one possibility is to use short mixing times. In addition, the tr-ROESY experiment has been proposed[176] to distinguish direct from indirect NOE cross peaks. Figure 6. (A) Shematic representation of a NOESY (left) and tr-NOESY (right) spectra. Cross-peaks are of the opposite sign to the diagonal peaks (positive NOEs) for a small molecule in the free state. Upon addition of the receptor, a sign change of cross-peaks takes place (same sign as the diagonal peaks; negative NOEs). (B) Nuclear Overhauser enhancements (NOEs) and tr-NOEs for α-L-Fuc-(16)-β-D-GlcNAc-OMe in the absence (filled symbols) and presence (open symbols) of Aleuria aurantia agglutinin, measured at 600 MHz as a function of the mixing time τm. Circles and diamonds refer to proton pairs H6proRGlcNAc-H6proSGlcNAc and H1Fuc-H6proSGlcNAc, respectively. Sources: Doctoral Thesis of Cinzia Guzzi[177] (A) and Claridge 2009[169] (B). The setup of transferred NOE experiments is identical to the setup of ‘‘normal’’ NOE experiments. The only difference is the preparation of the sample since the intensity of transferred NOEs strongly depends on the excess of ligand over protein. Applied to glycan systems, depending on the size of the carbohydrate ligand, three regimes may be distinguished: (a) The molecular weight of the carbohydrate ligand leads to correlation times ranging in the order of tens to hundreds of picoseconds, and therefore NOEs of the free ligand are positive. At 500 MHz, this is usually the case up to the size of trisaccharides. If charges are present as, for example, in sialic acid residues, the tumbling of the molecule is slower and one may observe negative NOEs already for a trisaccharide. (b) If the molecular weight is such that the correlation time approaches zero crossing conditions, no NOEs will be observable. For uncharged carbohydrates at 500MHz, this is usually the case for tetraand pentasaccharides. (c) Larger carbohydrates have correlation times of several nanoseconds and, therefore, display negative NOEs at frequencies of 500 MHz and higher. 2. Techniques and tools 62 In cases (a) and (b), the discrimination of transferred NOEs from free ligand NOEs is straightforward because at carbohydrate-to-protein ratios in which the equation 10 is fulfilled, the sign of the NOE changes upon binding from positive to negative. At the same time, the mixing time at which a maximum NOE is observed is reduced and in the range of 200 ms, as compared to 600-1000 ms for the free ligand. Because of this change in sign, the experiment has also been used to identify binding in mixtures of low molecular weight compounds[178]. In case (c), discrimination is less straightforward and usually requires the acquisition of NOESY experiments with different mixing times. 2.2.3 Saturation Transfer Difference spectroscopy (STD) The STD NMR experiment[179] is another spectroscopic technique to study the interactions, in solution, between a large molecule (receptor) and a medium-small sized molecule (ligand), and, alike tr-NOESY, it is based on the Nuclear Overhauser effect and the observation and analysis of the resonances of the ligand protons. The experiment is carried out by first registering a spectrum under conditions of thermal equilibrium with the irradiation frequency set at a value that is far from any ligand or protein signal (e.g. 40 ppm), i.e, the so-called off-resonance spectrum (figure 7, top), which is used as reference with signal intensities I0. A second experiment is then recorded, in which the protein is selectively saturated (on-resonance spectrum; see figure 7, middle), giving rise to ligand signals with Isat intensities. In general, the selective irradiation consists of a cascade of Gaussian-shaped pulses (low power) that saturate only a region of the spectrum that contains a few protein resonances (but not ligand signals), e.g., the aliphatic (from 0 to -1 ppm) or aromatic region (around 7 ppm), for a specific period of time (saturation time; typically from 0.5 to 5-6 seconds). The selective saturation is transferred to the whole protein via spin diffusion through the vast network of intra-molecular 1H-1H cross-relaxation pathways (intra-molecular NOE; see figure 7, middle), being a quite efficient processes due to the typical large molecular weight of the receptor. Also, saturation is transferred from the protein to the bound ligand via spin diffusion through inter-molecular NOEs. The dissociation of the ligand will then transfer this saturation into the bulk solution where it accumulates during the saturation time of the experiment, as a result of the much slower relaxation in the unbound that the bound state. In particular, as in fastexchanging protein-ligand systems the enthalpic relaxation (R1) of fast-tumbling molecules (small) in the free state is much slower than the kinetic off-rate constant of binding (koff >> R1), the accumulation of ligands molecules containing some of their resonances perturbed (NOE of large molecule) results in the macroscopic detection of transferred saturation on the ligand signals in the saturated STD NMR spectrum (Isat; figure 7, middle). 2. Techniques and tools 63 Figure 7. Scheme of the STD-NMR experiment showing the protein in surface representation and the non-exchangeable protons of the ligand as spheres. (Top) A 1D standard NMR experiments show only equilibrium intensities of the ligand in the free state (I0). (Middle) Upon selective saturation of some receptor signals, this is efficiently spread throughout the protein (yellow surface) by spin diffusion (intra-molecular NOEs). The fast exchange (transient binding) between the free and bound ligand states allows the transfer of magnetization (inter-molecular NOEs) from the receptor to the ligand protons in contact with the protein surface (salmon spheres, on-resonance spectrum). (Bottom) The difference spectrum (I0-Isat) only contains the ligand signals perturbed upon binding, whose intensities reflect the proximity of each proton to the protein surface. Furthermore, for those hydrogen atoms of the ligand establishing close contacts to the protein surface (4-5 Å) these Isat values will be lower than the I0 intensities, i.e, negative inter-molecular NOE, due to the transfer of the relaxation properties of the macromolecule to the small ligand in the bound state (figure 7, middle). By subtracting the off-resonance from the on-resonance spectrum (I0-Isat) the difference or STD spectrum is obtained, which will just contain the proton signals of the ligand in close contact to the protein surface (figure 7, bottom), and where any signal coming from non-binding 2. Techniques and tools 64 compounds is cancelled out. So, if a non-binder is present in solution its resonances will not appear in the STD spectrum. The signal intensities exclusively coming from the saturation transfer can be quantified (ISTD = I0-Isat), showing the proximity (high ISTD) or distance (low ISTD) of the proton to the receptor surface. Also, a blank experiment must be carried out to assure the absence of direct irradiation of the ligand. A sample containing the receptor at low concentration and a large molar excess of the ligand (1:50 up to 1:1000) is usually employed in STD NMR experiments. This precludes the perturbations of absolute STD intensities due to rebinding effects (i.e, a ligand already saturated experiences another association process, without previous full relaxation), which would impede to correctly determine the group epitope mapping (figure 8). The theory of NOE indicates that the magnetization transferred from receptor to ligand protons by intermolecular NOE depends on the inverse sixth power of their distances in the bound state. Thus, the shorter the protein-ligand proton-proton distance (bound state), the stronger the intensity of the corresponding STD signal. So, by normalizing all the measured STD intensities (I0-Isat/I0) against the most intense signal (which is arbitrarily assigned a value of 100%), the so-called “group epitope mapping” is obtained (expressed as percentages; see figure 8). This represents the fingerprint of protein-ligand contacts in the bound state, such that it illustrates which chemical moieties of the ligand are key for molecular recognition in the binding site. Figure 8. (Left) STD growth curves (absolute values) for the recognition of Man-α(1,2)-Man-α-O(CH2)NH2 by the anti-HIV-1 human antibody 2G12. (Right) Ligand group epitope mapping, or binding epitope for Man-α(1,2)-Man-α-O(CH2)NH2 compound. Note that proton H4A receives the highest saturation (thus the 100% is arbitrarily assigned) and is used as reference to calculated the STD percentages of the other protons. Residue of mannose-A makes the main contacts with the protein in the bound state. Sources: Doctoral Thesis of Pedro M. Enríquez-Navas[180], and Enríquez-Navas et al. 2011[175c]. Binding epitopes were commonly obtained for a given saturation time, assuming that the resulting ligand group epitope mapping did not depend on the chosen saturation time. However, significantly different R1 relaxation rates of the ligand protons can produce artefacts in the epitope definition[181]. In 2. Techniques and tools 65 particular, protons with slower R1 relaxation enable a more efficient accumulation of saturation in solution such that their relative STD intensities may be significantly overestimated at long saturation times, thus also overrating the proximity of those protons to the protein surface. Indeed, the structual information that the binding epitope provides can be affected by 1) differences in R1 relaxation rates of ligand protons, 2) the extent of saturation received in the first place, and 3) the kinetics of binding. Since these sources of distorsions are the consequences of differences in the ability to accumulate saturation in the free state, we can cancel them by deriving STD intensities close to zero saturation time, this is, when virtually no accumulation of saturated ligand molecules is taking place. This is usually referred to as initial slopes of the build-up curves (figure 9). In addition, under this approximation, possible artefacts comming from intra-molecular spin diffusion (bound state) can also be minimized. To calculate the initial slopes, Mayer and James proposed fitting the experimental build-up curves to the mono-exponential function[182]. ( ) ( ( )) Eq. 11 where STD(tsat) is the observed STD intensity, STDmax is the asymptotic maximum of the build up curve, tsat the saturation time, and ksat the rate constant related to the relaxation properties of a given proton that measures the speed of the STD build-up. ksat and STDmax are derived by least-squares fit, and the initial slope (STD0) of the curve is obtained as ( ) | Eq. 12 These STD0 values are then used to characterize the binding epitope independently of T1 and rebinding effects. Note that for ⁄ ( ) ( ) Eq. 13 Thus, if we plot the normalized STD factor versus the saturation time (tsat), ksat can be approximated as the inverse of the interpolated tsat value from a STD factor of 0.63 (figure 9). 2. Techniques and tools 66 Figure 9. Normalized STD build-up curves fitted to a mono-exponential function to obtain the initial slope (STD0) from the multiplication of STDmax and ksat. By interpolating a value of 0.63 for the STDmax factor, the ksat is obtained as the inverse of tsat. To get protein-ligand association curves from STD NMR experiments, Mayer and Meyer introduced the conversion from observed experimental intensities (I0-Isat/I0), which depend on the fraction of bound ligand LB f ][ ][ ][][ ][ ][ ][ PK P PLL PL L PL f DT LB     Eq. 14 , to STD amplification factors (STD-AF) by multiplying the observed STD by the molar excess of ligand over protein[183] [ ] [ ] Eq. 15 so that the STD amplification factor (STD-AF) depends on the fraction of bound protein (fPB). ][ ][ ][][ ][ ][ ][ LK L PLP PL P PL f DT PB     Eq. 16 Therefore, a plot of STD-AF values at increasing ligand concentrations will give rise to the protein-ligand binding isotherm, from which the dissociation constant KD can be derived. However, in eq. 16 the concentration of free ligand [L] is not known. This can be solved by employing a very low total concentration of protein, such that [ ] , and an excess of ligand. Under these conditions [ ] [ ] , so that 2. Techniques and tools 67 [ ] [ ] Eq. 17 and thus the dissociation constant KD can be easily determined by plotting the normalized STD-AF values at increasing total ligand concentrations, i.e, the binding isotherm (figure 10). This can be done graphically by interpolation as the KD coincides with the amount of free ligand [ ] that is necessary to reach a 50% of fPB ( 5.0 PB f ) [ ] [ ] [ ] Eq. 18 Therefore, the weaker the interaction, the higher the ligand concentration required to saturate half of the protein molecules. Figure 10. Example of STD binding isotherms, at different saturation times, of the WGA protein (46 µM) titrated with chitobiose. STD-AF values appear normalized against their corresponding plateau values. Note that different KD values (KD’, KD’’) are obtained at different saturation times. Adapted from Angulo et al. 2010[184]. The determination of KD from STD-NMR experiments is affected by different experimental parameters, as it has been thoroughly studied in our group[184]. In particular, the saturation time (tsat) employed (figures 10 and 11), the STD intensity of the signal (figure 11) and the fraction of bound ligand are key factors that must be chosen wisely for the accurate determination of dissociation constants. The apparent KD increases monotonically with tsat, thus understimating the protein-ligand affinity. Similarly, the higher the STD intensity of a resonance the more overstimated KD value will be obtained. Also, if the fraction of bound ligand is modified by increasing the receptor concentration, the apparent dissociation constant will be larger. 2. Techniques and tools 68 Figure 11. Effect of the experimental factors (saturation time and monitored proton) on the determination of the apparent binding constant for the BSA-L-tryptophan system on a sample of 20 μM of BSA. 2.3 Molecular modelling A force field consists of the combination of a mathematical formula and associated parameters that are used to describe the overall potential energy of a molecular system as a function of its atomic coordinates. Currently, it is widely accepted that force fields are critical to molecular simulation in many aspects of life sciences research. Understanding, analysing, and predicting 3D structural models of molecular systems, including their conformations, binding affinities and related properties, all depend on accurate atomic force fields. For this reason, there has been a great deal of effort devoted to the development and improvement of potential energy functions and their parameters, which are the two features that define a force field. Energy minimizations and MD simulations are often limited by inadequate description of the various force field parameters for the systems of interest. For instance, if a crystal structure is minimized without including penalty terms for the structure factor, then the deviations from the experimental structure are often much larger (0.5 to 1.5 Å) than the expected error[185], reducing confidence in molecular mechanics analyses. Thus, the key factor affecting the quality of MD simulations is the force field accuracy, which determines the goodness of the conformational sampling and dynamics performance in reproducing the experimental observables. Computational methods such as docking and molecular dynamics (MD) simulations (also homology modelling and computational mutagenesis) provide complementary tools, indispensable in many cases, to fully understand both X-ray and NMR data. Importantly in the context of the present thesis, 2. Techniques and tools 69 they perform particularly well (specially MD simulations) in characterizing the structure and dynamics of glycans and glycoconjugates[186] (see next section). The force fields commonly employed consist of a combination of bonded and non-bonded energy terms[187]. Given atomic positions and velocities, forces are calculated on the fly, as derivatives of the potential energy (V (r1, . . ., rN)), at specific time steps. The overall potential energy for a molecular system can be written in classical mechanics as ( ) ∑ ( ) ∑ ( ) ∑ ∑ [ ( )] ∑ ∑[ ] Eq. 19[63] where the parameters shown in red must be known from experiments or derived from QM calculations, and included in the force field. of Σbonds, Σangles and Σdihedral terms refer to the potential energy associated with bond-stretching, angle-bending, and proper (and improper) dihedral angle rotations, respectively, whereas ΣCoul and ΣLJ represent the pairwise electrostatic interaction and the Lennard-Jones (LJ) repulsion-dispersion potential energy terms, respectively. Thus, classical force fields are defined by both the functional form of the different terms contributing to the global potential energy and by the set of parameters that each term requires. The different energy terms are given by empirical formulae and/or harmonic functions penalizing deviations from ideal values, these being determined from high-resolution crystallographic or spectroscopic data and/or from calibration to QM calculations[187a, b, 187d, 188]. The accurate determination of these empirical parameters, such that when introduced in equation 19 lead to the correct potential energy landscape of the molecular system, is a crucial, meticulous and challenging task in force field development. Furthermore, due to the coupling between many of the force field terms (e.g. torsions and electrostatics), the parameterization process inevitably requires testing multiple sets of calculations for optimization. In this regard, the better the force field refinement protocol, the more probably the resultant parameters will be broadly applicable[189]. It is well known that hydrogen bonds formation is the driving force to many phenomena, including the generation and stabilization of secondary structures[190], protein folding and stability[191], molecular recognition[192], and drug binding and enzymatic reactions that involve transfer of protons[193]. Therefore, the detailed understanding of hydrogen bonds geometry and their incorporation into accurate potential functions is of fundamental importance, although many efforts in this direction still have to been done. 2. Techniques and tools 70 The resolution of X-ray crystallography data for proteins rarely covers beyond 1.0 Å. For this reason, studies of hydrogen bonds in proteins have been mostly limited to the coordinates of non-hydrogen atoms. Thus, the hydrogen bond selection criteria used in some of these studies are based on the distance between the potential donor and acceptor atoms[194]. However, the hydrogen bond geometry can be better understood in terms of the angle and distances involving the positions of hydrogen[195], even if the hydrogen positions are modelled only implicitly from their heavy atom neighbours. For these reasons, several strategies have been used to account for hydrogen bonding in crystallographic refinement and molecular simulations. The hydrogen bond potential is often implicitly parameterized as a combination of Lennard-Jones (L-J) and electrostatic terms. In force fields that use an explicit hydrogen bonding term, this is typically included as a distance-dependent function without any directional component. The functional form may be a L-J 6–12[196], a L-J 10–12[187a, b, 197], a L-J 6–9[196a, 198], or a Morse type potential[199]. On the other hand, CHARMM[187b] and MM3[200] force fields now include a cosine directional term. In these implementations, the hydrogen bond energy is minimized when the hydrogen bond N-H•••O is linear, with an angle of 180°. However, these do not reflect the non-linear directional preferences of hydrogen bonds at the acceptor molecule, conferred by their covalent component[194a, 195, 201]. In any case, it has to be noted that displacement of water molecules competing for hydrogen bonds is not accounted for in any force field. 2.3.1 Modelling of glycans As it has been commented in Chapter 1, carbohydrates encode an amount of potential information that is several orders of magnitude higher than in the case of any other biological macromolecule, and this complex encoding capacity arise from their enormously diverse (sequence), complex (ramifications) and flexible (local and/or global conformation) structures. This structural complexity makes the characterization of the 3D structure of oligosaccharides, their conjugates, and analogues, particularly challenging for traditional experimental methods. Thus, computational methods provide a basis for interpreting sparse experimental data and for independently predicting conformational and dynamic properties of glycans, which eventually can contribute unique insights into the relationship between oligosaccharide structure and biological function (figure 12). In the early years of biomolecular modelling, MD simulations were technically limited to small biological systems, e.g., small proteins[202], short DNA helices[203] and mono-[204] or disaccharides[205], and for short simulation times. Nevertheless, the growing evidences of the dynamical nature and fundamental role on biological functions of biomolecules acted as a catalyser for the development of computer modelling tools applied to structural biology[206]. Today, advances in computer technology and software algorithms enable us to sample the conformational space and dynamics of biomolecular systems for simulation times that go from hundreds of nanoseconds (e.g. most internal motions in glycans[6a, 207] to several microseconds (e.g. timeframe for the conformational equilibrium of the iduronate ring in GAGs[17]). However, it should be noted that very long timescales are not always necessary to obtain useful information. For instance, not too long MD simulations can be very 2. Techniques and tools 77 motion conserve energy and thus provide a suitable scheme for calculating a microcanonical ensemble. However, of more practical application is the canonical ensemble since it can readily be performed by coupling the molecular system to a constant-temperature bath, which rescales the atomic velocities according to the desired temperature. In a similar manner, constant pressure simulations can be performed. The common procedure for running stable MD simulations of complex molecular systems (inclusion of explicit solvent, counterions, etc) is shown in figure 14. Several algorithms have been developed for MD simulations that predict the time evolution of a system for a limited time. Thus, since physically observed properties are computed as the corresponding time averages on the individual microstates ensemble, for the results to be meaningful, the simulations must be sufficiently long so that the important motions are statistically well sampled. However, it must be considered that the longer the simulation the higher will also be the possibility of force field deviations to appear[6a]. Thus, the time scale of the phenomenon aimed to investigate should be considered in each particular case to decide the simulation time that interests us. Experimentally accessible spectroscopic and thermodynamic quantities can be computed, compared, and related to microscopic interactions. It should be noted that MD is severely limited by the available computer power. With currently available clusters and computing algorithms, it is feasible to perform a simulation with several thousand explicit atoms for a total time that goes from hundreds of nanoseconds to several microseconds (the use of GPUs is providing cutting-edge velocities). However, it may be possible that the carbohydrate molecules undergo dynamical events on longer timescales and/or that the accessible computational resources are not updated, so that the time scale of interest cannot be investigated with standard MD techniques. For those cases, another way is to use high temperature dynamics to allow the molecule to assume high-energy conformations. However, this approach has to be used with caution since it can force the molecules to adopt unrealistic conformations (artefacts). Figure 14. General scheme of the MD simulation protocol commonly followed. 2. Techniques and tools 78 Time-averaged restrained molecular dynamics simulations (tar-MD) As the inter-conversion between the L-IdoA2S conformers in GAGs is rapid on the NMR time scale, the observed NMR resonances reflect an average of both conformations. Similarly, NOES from both conformations are observed simultaneously. Thus, the use of a MD methodology involving “instantaneous” experimental constraints[237] would generate structures that simultaneously satisfy both sets of experimental data. However, there may be no single conformation in agreement with the whole set of NOEs, and even if one is generated, it may be highly strained and physically unrealistic (a so-called virtual conformer[238]). Under these considerations, it is not correct to treat the NOE data as affording a fixed distance boundary. Instead, NOE distance information should be used to enforce an average distance limit through time. This can be achieved by imparting particles with a memory of their history with respect to internuclear distances. At the same time, to truly model the physical nature of the NOE, it is necessary to account for the nonlinear dependence of the measured NOE intensity on the internuclear distance. Based on these grounds, the methodology time-averaged restrained molecular dynamics (tar-MD)[239] includes the presence of a penalty term, Epenalty, in the total potential energy equation. For each experimental restraint (distance or coupling constant), six keywords (r1 to r4 distances or couplings constants, and rk2 and rk3 force constants) are defined. These parameters delimit the shape of the restraining potential as follows (figure 15) ( ) ( )( ) ( ) ( ) ( ) ( )( ) Eq. 20 Figure 15. Representation of the penalty energy (Epenalty) applied in tar-MD simulations as a function of the value of the experimental parameter. 2. Techniques and tools 79 where R is the time-averaged value between atoms and r1 to r4 are the different limits defined below (r1 and r2) and beyond (r3 and r4) the experimental distance. The values of r1 to r4 have to be specified in Å or Hz, depending on the restraint type. rk2 and rk3 are defined in kcal·mol−1·Å−2 for distance restraints and in kcal·mol−1·rad−2 for dihedral restraints. The experimental key proton-pair distances, as well as coupling constants, are usually implemented as structural restraints with a 10% margin, using a flat well potential. Since the NOE arises from dipolar interactions between nuclei, the intensity of a NOE signal grows as r-6 (as long as the simulation time is higher than the correlation time for overall molecular tumbling; otherwise, it has been shown that r-3 averaging is necessary)[240]. Thus, for distance constraints the expression for R included in MD force fields is 〈 〉 (∑ ( ) ( ) ) Eq. 21 with being the characteristic time for the exponential decay or exponential decay constant, and t’ the total simulation time. The exponential decay constant “helps” the calculation to converge and it is commonly set to a value 10 times smaller than t’. This methodology has been extensively and successfully applied to the study of different molecular systems[241], iduronate containing carbohydrates included[78]. 2.3.4 Docking Molecular docking is a computational procedure that aims at predicting the preferred orientation and conformation of a ligand bound to its target protein. In order to perform computational protein-ligand docking calculations, the 3D structure of the receptor must be known. Each docking program operates slightly differently, but they share common features that enable them to (1) search for locations on the protein surface that lead to favourable interactions with the ligand, (2) sample the conformational space of the ligand, and (3) compute the interaction energy between the protein and ligand (score of “binding affinity” or scoring function). For instance, glide (grid-based ligand docking with energetics) uses a series of hierarchical filters to search for possible locations of the ligand in the active-site region of the receptor (figure 16). 2. Techniques and tools 80 Figure 16. Glide docking “funnel” showing the protocol followed to generate docked poses. Glide docking algorithm approximates a complete systematic search over ligand positions, orientations, and conformations in the receptor site, with increasingly demanding tests applied as the search space is reduced. Glide generates and docks many core conformations, but treats the rotamer groups sequentially, rather than combinatorially, which speeds up the calculation. The hierarchical protocol followed by glide, shown in figure 16, can be briefly described as follows: 1. Site-point search - Generate a 2-Å grid of site points in the active site. - Pre-compute histograms of distances between site point and receptor surface in grid setup. - Compare site point – receptor surface histograms with the ligand centre–ligand surface histogram. - Reject mismatched site points. 2. Dimensional tests and rough scoring - Diameter test: check steric clashes of atoms near ligand diameter for ~300 pre-specified orientations of the ligand diameter (figure 17). - Subset test: rotate about ligand diameter in 15º increments, and score atoms capable of establishing hydrogen bonds or ligand-metal interactions. - Greedy scoring: score all atom positions ±1 Å in x,y,z directions and use best score. - Refinement: move the whole ligand ±1 Å in x,y,z directions, re-score and reduce ~5000 poses to ~400 for energy minimization. 2. Techniques and tools 81 3. Energy minimization - Use pre-computed OPLS-AA electrostatic and van der Waals grids. - Anneal from soft-to-hard potential: smoothing reduces large initial energy/gradient terms from close contacts, permits freer movement. - Also optimize torsional angles when doing flexible docking. - Use Monte Carlo moves to explore nearby torsional minima for a small number of low-energy poses. 4. Final scoring. - Choose best pose(s) based on Emodel, which is a combination of the Coulomb-vdW energy, the GlideScore (enhanced version of ChemScore) and internal strain energy. - Final scoring based on GlideScore, consisting of ChemScore terms, the Coulomb-vdW energy and terms that penalize non-physical interactions Figure 17. Definition of ligand diameter and ligand centre parameters according to glide searching algorithm. Source: Thomas A. Halgren, Schrödinger. Protein interaction with the ligand relies on both the protein backbone fold and the orientation of the side chains in the binding site region. One of the most significant limitations in docking is that it is generally performed while keeping the protein surface rigid, which prevents the consideration of the effects of induced fit within the binding site. These difficulties are mostly due to the high number of degrees of freedom characterizing a protein–ligand system, which increases the computational cost of docking calculations. Thus, several approximations about the flexibility states may be introduced in molecular docking. The simplest approximation (rigid docking) considers only the three translational and three rotational degrees of freedom of the protein and those of the ligand, treating them as two distinct rigid bodies. However, the most widely used algorithms at present enable the ligand to fully explore its conformational degree of freedom in a rigid-body receptor[242]. 2. Techniques and tools 82 Induced Fit Docking (IFD) As we have introduced above, in standard docking studies ligands are docked into the binding site of a receptor where the latter is held rigid and the ligand is free to move. While this approximation present the advantage of reducing the computational cost, it however may give rise to misleading results, since in reality many proteins undergo side-chain or backbone movements, or both, upon ligand binding. These changes allow the receptor to alter its binding site so that it better adapts to the shape and binding mode of the ligand. This is often referred to as Induced Fit and is one of the most challenging features to model in structure-based drug design. Thus, a good Induced Fit Docking (IFD) protocol should both generate an accurate complex structure for a ligand known to be active but that cannot be docked in an existing (rigid) structure of the receptor and also rescue false negatives (poorly scored true binders) in virtual screening experiments, where instead of screening against a single conformation of the receptor, additional conformations obtained with the IFD protocol are used. Since glycans interactions present dissociation constants (KD) in the weak binding regime, i.e., poorly scored true binders from the point of view of docking scoring functions, the IFD method may result specially convenient to obtain an accurate description of their interactions. Furthermore, poor binders exhibit a significant dependence on the initial input conformation during docking due to the use of complex grid-based potentials and the practical limitations of thoroughly sampling the docked poses. Also, docking algorithms usually generates new ligand conformations through torsional variations only, so any differences in bond lengths and angles in the input ligand structures will persist through the docked poses, resulting in scoring and pose differences. For these reasons, docking several conformations of each ligand with variations in bond lengths and bond angles is a reasonable strategy to reduce input dependence. Grid generation Docking algorithms represent the shape and properties of the receptor on a grid by several different sets of fields that provide progressively more accurate scoring of the ligand poses. Glide allows to define the receptor structure by excluding any co-crystallized ligand that may be present, determine the position and the extent of the region for which receptor grids will be calculated (30Å-sided cube maximum) and set up constraints. In any docking job using these receptor grids, ligands are confined to the enclosing box (figure 18). Also, the ligand centre can be set during grid generation in glide. The ligand centre of a ligand is defined, in glide, as the midpoint of the longest line segment that can be constructed between any two atoms in the ligand. Furthermore, the ligand diameter midpoint box is the region in which the diameter midpoint of each docked ligand must remain. Each dimension of this box can be modified from its default value of 10 Å to the 6-14 Å range. When doing so, the enclosing box also changes its size to make the distance between faces of the enclosing box and the ligand diameter midpoint box alike. 2. Techniques and tools 83 A larger ligand diameter midpoint box can be useful to allow ligands to find unusual or asymmetric binding modes in the active site. Conversely, if the default ligand diameter midpoint box allows ligands to stray into regions you know to be unfruitful, you can confine their midpoints to a smaller box, eliminating some of the less useful poses and saving calculation time. Figure 18. 3D image showing the grid box (purple-lined cube) and the ligand diameter midpoint box (green-lined cube) in glide. Docking Algorithms The docking algorithms can be grouped into deterministic and stochastic approaches. While deterministic algorithms are reproducible, stochastic algorithms include random factors that do not allow the full reproducibility. The most widely used algorithms in docking simulations are described below. Incremental Construction Algorithms These algorithms consist of the division of a ligand into rigid fragments. One of the fragments is selected and placed in the protein binding site. The reconstruction of the ligand is then carried out in situ, adding the remaining ligand fragments. For example, DOCK[243] uses incremental construction algorithm to treat ligand flexibility. It generates points (sphere centres) that fill the binding site and try to capture the binding site shape properties for identifying favourable regions in which the ligand atoms may be located. The ligand is divided along each flexible bond to generate rigid segments. An 2. Techniques and tools 84 anchor fragment is then selected from all the rigid pieces and oriented in the active site by matching ligand atoms with sphere centres. After, fragments are added and all possible placements are scored on the basis of their interactions with the protein using the energetic scoring function. Then, the best anchor fragments are used for completing the construction of the ligand in the protein-binding site. Finally, the best scored poses of the complete ligand are selected. Genetic Algorithms Genetic algorithms are stochastic searching approaches that use techniques inspired by evolutionary biology to find reliable results. It mimics the process of evolution by manipulating a collection of data structures called chromosomes. AutoDock[244] uses this algorithm for obtaining reliable docking results. First, the protein is placed inside a cube with a predefined size, characterized by a defined number of points (grid points). In the second step, probes corresponding to the different atom types of the ligand are then moved through the cube and, in particular, at each point, protein–probe interaction energies are calculated and stored in affinity maps. Thirdly, a conformational search of the ligand is performed by applying the Lamarckian genetic algorithm[245]. At this stage, a minimization or local search is performed, and the new conformation is then considered as input for a new iteration of the genetic algorithm cycle. Hierarchical Algorithms It uses an exhaustive systematic search for discovering the most favoured ligand conformations in the protein active site, with a screening based on progressively restricted energetic cutoffs. A grid and a molecular surface containing information of the protein receptor properties are calculated before the algorithm search. Then a set of initial ligand conformations is produced and screens are performed over the whole phase space available to the ligand to locate promising ligand poses in the respective receptor fields. Afterwards, ligands are minimized in the field of the receptor using a standard molecular mechanics energy function[227, 233]. Finally, the lowest-energy poses are subjected to a Monte Carlo procedure that examines torsional minima and a composite scoring function is then used to select the correct docked poses. The algorithm used in glide[246] can be defined as a hierarchical algorithm. Scoring Functions Energy scoring functions are necessary to evaluate the free energy of binding ΔG (affinity) of ligand-receptor interactions. The Gibbs free energy equation describes the ligand-receptor free energy of binding as Eq. 22 where ΔG represents the energetic changes between the bound and unbound states of both ligand and receptor, ΔH is the enthalpy, T the temperature expressed in Kelvin, and ΔS the entropy of the system. 2. Techniques and tools 85 Furthermore, ΔG is related to the binding association and dissociation constants (Ka and Kd, respectively) as Eq. 23 (where R is the constant for ideal gases), which allows to obtain an estimate of binding affinity. Some sophisticated techniques for predicting binding free energies are currently too slow to be used in molecular docking of large sets of compounds[247]. Thus, fast scoring functions have been developed. Empirical scoring functions use a set of parameterized terms describing properties known to be important in protein–ligand binding to construct an equation for predicting binding affinities. Multilinear regression is used to optimize these terms using a set of known protein–ligand complexes. These terms usually describe polar–apolar interactions, loss of ligand flexibility (entropy), and desolvation effects. For instance, GlideScore 2.5 scoring function[246a] is a regression based empirical scoring function of the form ∑ ( ) ∑ ( ) ( ) ∑ ( ) ( ) ∑ ( ) ( ) ∑ ( ) Eq. 24 The first term describes the lipophilic and aromatic interactions, whereas the polar terms are included in the second, third and fourth terms (hydrogen bonds separated into differently weighted components that depend on the electrostatic properties of donor and acceptor atoms), and the Cmax-metal-ion term, which includes the anionic(ligand)-metal(receptor) interactions. The seventh term rewards instances in which a polar but non-hydrogen bonding atom is found in a hydrophobic region. Also, Coulomb and van der Waals interaction energies between the ligand and the receptor are evaluated as well as the solvation effect. GlideScore[246] has been optimized for docking accuracy, database enrichment and binding affinity prediction, and can be used as an empirical scoring function that approximates the ligand binding free energy. GlideScore should be used to rank poses of different ligands, for example in virtual screening. Other scoring functions different than GlideScore are also used in glide. Therefore, Emodel scoring function has a more significant weighting of the force field components (electrostatic and van der Waals energies), which makes it well-suited for comparing conformers but much less so for comparing chemically-distinct species. Glide uses Emodel to pick the "best" pose of a ligand (pose selection) and then ranks these best poses against one another with GlideScore. This means that if you save multiple poses per ligand, the apparent ranking of poses for a given ligand (by GlideScore) will not reflect the 2. Techniques and tools 86 actual ranking that glide used for pose selection. So, the value of Emodel scoring has to be considered to determine the highest ranked pose for a ligand. On the other hand, force-field-based scoring functions (e.g. those used by AutoDock or DOCK) are based on the non-bonded terms of the classical molecular mechanics force fields. In AutoDock[244, 248], the implemented scoring function presents five terms with coefficients empirically determined using linear regression analysis from a set of protein–ligand complexes with known binding constants. A 12-6 Lennard-Jones potential and a Coulomb term taken from the AMBER force field[196c] describe the van der Waals and electrostatic interactions, respectively. In addition, hydrogen bonding is described with a 12-10 Lennard-Jones term with Goodford directionality[249]. Also, a desolvation energy potential and an empirical measure of the unfavourable entropy of ligand binding due to the restriction of conformational degrees of freedom are included. Other docking approaches use knowledge-based scoring functions based on statistical observations of intermolecular close contacts in protein–ligand X-ray databases, which are used to derive potentials of mean force. This methodology assumes that the frequency of close intermolecular interactions between certain ligand and protein atoms contribute favourably to the binding affinity. In this approach, no fitting to experimental affinities is required and solvation and entropic terms are treated implicitly[250]. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 93 2D-NOESY experiments (278 K), from which cross-peaks corresponding to an anti-Ψ geometry around the IdoA-GlcN linkage could be identified (H1b-H3a, H1b-H5a and, in some cases, H5a-H6a+a’), together with larger H1b-H4a and H1b-H6a+a’ NOEs (with exceptions), the latter showing a major syn conformation (figure 3). The anti-Ψ exclusive H5b-H6a+a’ NOE was only observed, with a very weak intensity, for Tri1, Tri3, Tri5 and Tri6. However, it is noticeable that the growth of this NOE is surely affected by the loss of magnetization due to a very short longitudinal relaxation time (T1), as a consequence of the efficient (in terms of T1) relative reorientation of the H5b and H6a+a’ coupled protons (free rotation of methylene protons and IdoA2S conformational plasticity acting in the nanosecond time scale as relaxation mechanisms). Thus, although the anti-conformation is present, the H5b-H6a+a’ NOE intensity may not be observed (or very weakly) due to fast relaxation. This is why, to estimate the contribution of the anti-Ψ conformers (IdoA-GlcN linkages) it is preferable to consider the H1b-H3a and H1b-H5a NOEs (specially the former). A table containing the normalized H1b-H3a cross-relaxation rates for Tri1-Tri8 is included in the Appendix, showing that the GlcNS-IdoA2S-GlcNAc sequence (Tri8) contributes to the largest extent to the presence of anti- conformations around the IdoA2S-GlcN linkage. Figure 2. 1D-NOESY growth curves (at 278 K) corresponding to the H2b-H5b distance of trisaccharides Tri1 (6-OSO3in the reducing terminal) and Tri5 (6-OH in the reducing terminal). Note the significantly higher NOE intensity (shorter H2b-H5b distance) for the trisaccharide containing the 6-sulphate group (Tri1). 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 94 Table 1. Comparison between the experimental (exptl) and theoretical (tar-MD and free-MD) most relevant distances of the library of trisaccharides. The experimental values were derived from 1D-NOESY experiments at 278 K. The tar-MD derived results were calculated as <r-6>-1/6 over the 8000 frames of tar-MD simulation at 278 K. The free-MD calculated were calculated as <r-6>-1/6 and weighted on the populations of iduronate conformers at 278 K. The tar-MD and free-MD H1b-H6a distances represent the r-6 average over the H1b-H6aproR and H1b-H6aproS values. *Values derived from 2D-NOESY experiments at 278 K and 600 ms mixing time. The H1c-H2c NOE was used as reference (as in the case of 1D-NOESY derived distances). Comp Method Distance (Å) GlcN-IdoA2S linkage L-IdoA2S IdoA2S-GlcN linkage H1c-H3b H1c-H4b H1b-H3b H2b-H5b H1b-H3a* H1b-H4a H1b-H6a* Tri1 free-MD 2.3 2.8 3.1 2.6 2.7 2.4 3.5 tar-MD 2.4 2.6 3.3 2.9 4.2 2.3 2.9 exptl 2.6 2.7 3.0 3.0 3.3 2.5 2.7 Tri2 free-MD 2.3 2.9 3.1 2.6 2.8 2.4 3.4 tar-MD 2.6 2.5 3.2 2.8 4.2 2.3 2.9 exptl 2.6 2.6 3.0 2.9 3.0 2.4 2.6 Tri3 free-MD 2.3 2.8 3.1 2.6 3.9 2.3 3.3 tar-MD 2.5 2.5 3.3 3.0 4.3 2.4 2.8 exptl 2.6 2.7 3.0 2.9 3.4 2.6 2.7 Tri4 free-MD 2.3 2.8 3.1 2.6 3.8 2.3 3.2 tar-MD 2.6 2.5 3.3 2.9 4.1 2.3 3.1 exptl 2.7 2.7 3.0 2.9 3.2 2.5 2.7 Tri5 free-MD 2.3 2.8 3.4 2.9 3.8 2.3 3.5 tar-MD 2.3 2.6 3.6 3.3 4.1 2.3 2.9 exptl 2.6 2.6 - 3.2 2.8 2.6 2.5 Tri6 free-MD 2.3 2.6 3.3 2.8 2.9 2.4 3.5 tar-MD 2.4 2.5 3.6 3.1 4.2 2.3 2.9 exptl 2.6 2.6 3.2 3.2 3.7 2.6 2.6 free-MD 2.3 2.7 3.4 2.8 3.7 2.3 3.3 Tri7 tar-MD 2.4 2.5 3.3 3.0 4.3 2.3 2.9 exptl 2.6 2.6 3.2 3.2 3.1 2.5 2.6 free-MD 2.3 2.7 3.2 2.7 3.9 2.3 3.2 Tri8 tar-MD 2.5 2.5 3.4 2.9 4.2 2.3 3.1 exptl 2.6 2.6 3.1 3.0 2.6 2.6 - 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 95 Regarding the conformation of the iduronate ring, the analysis of the accurately measured coupling constants (table 2), the absence of the H5c-H5b NOE signal (exclusive of the iduronate 4C1 chair conformation), together with the observation of the 2SO-exclusive H2b-H5b and H1b-H3b contacts (NOE, ROE, T-ROE; see NOESY spectrum in figure 3), indicated the coexistence in solution of the 1C4 chair and the 2SO skew-boat conformers. This is a characteristic feature of the internal iduronate ring in heparins and heparin-derived oligosaccharides[69a]. From the NOE-derived distances obtained (table 1), we focused our analysis on those defining the local (distances H2-H5 and H1-H3 of the L-IdoA2S ring) and global (H1c-H3b and H1c-H4b distances of the GlcN-IdoA linkage, and H1b-H3a, H1b-H5a, H1b-H4a and H1b-H6a+a’ of the IdoA-GlcN glycosidic linkage) conformations, comparing them along the Tri1-Tri8 library. First, dealing with the local geometry (L-IdoA2S residue), the calculated values for the H4b-H5b distance (between 2.4 and 2.6 Å) are in agreement with the canonical 1C4 (2.5 Å) and 2SO (2.4 Å) conformations (PDB code 1hpn). Regarding the 2SO-exclusive H1b-H3b and H2b-H5b distances (table 1), the results showed shorter values in the Tri1-Tri4 ensemble ( 3.0 Å) compared to the Tri5-Tri8 one ( 3.2 Å). Since this distances, according to the canonical conformations, are much shorter in the 2SO conformer (2.9 and 2.4 Å for the H1b-H3b and H2b-H5b distances, respectively; 4.3 and 4.0 Å in the 1C4 pucker), this observation was in agreement with a higher population of 2SO pucker in the Tri1-Tri4 ensemble, thus when the 6-OSO3group is present in the reducing end GlcN residue (figure 4). Figure 3. Expansion of a NOESY experiment for Tri1 showing the signals corresponding to the anti and syn rearrangements of the Ido(B)-GlcN(A) glycosidic linkage. The most relevant exclusive NOE peaks of the anti-Ψ conformation are marked with a star symbol. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 96 Table 2. Comparison of the experimental (exptl) and theoretical (free-MD and tar-MD) proton-proton vicinal coupling constants (3JHH) for the L-IdoA2S ring in the Tri1-Tri8 compounds. The experimental values were measured at 278 K. The tar-MD derived results represent the average over the 8000 frames of tar-MD simulation at 278 K. The free-MD calculated values are weighted on the populations of iduronate puckers at 278 K. a Over 1 Hz of difference with respect to the corresponding experimental value. Comp Method H1-C1-C2-H2 H2-C2-C3-H3 H3-C3-C4-H4 H4-C4-C5-H5 free-MD 3.9 5.6 3.4 2.3 Tri1 tar-MD 3.1 5.1a 3.7 2.4 exptl 3.0 6.2 4.0 3.1 free-MD 4.1 5.9 3.5 2.3 Tri2 tar-MD 3.5 5.8 3.8 2.5 exptl 3.2 6.5 3.9 2.9 free-MD 4.0 5.9 3.6 2.4 Tri3 tar-MD 3.3 6.0 3.9 2.6 exptl - 6.3 3.7 3.3 free-MD 4.1 5.6 3.5 2.2 Tri4 tar-MD 2.9 5.2a 4.3 2.8 exptl 3.1 6.3 3.9 3.4 free-MD 3.0 4.1 2.9 1.9 Tri5 tar-MD 3.0 3.0a 3.1 1.9 exptl 2.3 4.8 - 2.7 free-MD 3.3 4.5 2.9 2.0 Tri6 tar-MD 2.2 3.7a 4.2 2.6 exptl 2.5 5.1 - 3.0 free-MD 3.0 4.3 2.9 2.0 Tri7 tar-MD 3.6a 3.8 3.6 2.4 exptl 2.1 4.8 4.0 2.7 free-MD 3.6 4.9 3.1 2.0 Tri8 tar-MD 2.8 4.9 3.7 2.4 exptl 2.8 5.3 3.5 3.1 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 97 Figure 4. Influence of 6-O-sulphation (reducing end GlcN) on the L-IdoA2S conformational equilibrium. The absence of the 6-O-sulphate group (Glc(A)) is indicated with a dashed line circle. About the global conformation, the distances defining the geometry around the GlcN-IdoA2S glycosidic linkage (H1c-H3b and H1c-H4b; see table 1) are pretty similar (from 2.6 Å to 2.7 Å) for all the trisaccharides, thus indicating a rather rigid syn conformation (no correlation with sulphation pattern). This observation is in agreement with the results obtained from MD simulations (see next section). With respect to the IdoA2S-GlcN glycosidic linkage, similar short H1b-H4a and H1b-H6a+a’ distances (2.4-2.6 and 2.5-2.7 Å, respectively), corresponding to syn-Ψ conformers, have been determined. Differently, the anti-Ψ exclusive H1b-H3a distance shows significant variations upon sulphation pattern, going from longer (3.7 Å for Tri6) to shorter (2.6 Å for Tri8) distances. Although these differences are not clearly correlatable to the substitution pattern, interestingly, for the less sulphated trisaccharide Tri8, both the syn-Ψ and anti-Ψ exclusive H1b-H4a and H1b-H3a distances, respectively, present an equally short value (2.6 Å, intense NOE peaks; see Tri8 NOESY spectrum in Appendix), thus indicating a very similar contribution of both conformers in solution. This result suggests that Tri8 sulphation pattern, i.e, just Nand 2-O sulphation at the non-reducing glucosamine and the internal iduronate ring, respectively, enhances or facilitates the presence of the anti-Ψ conformations around the IdoA-GlcN glycosidic linkage. Probably, this is due to the reduced electrostatic repulsion forces existing in Tri8 (less sulphated trisaccharide of the library). Modelling Previous 3JHH-based studies on the conformational equilibrium of the L-IdoA2S ring in heparin derivatives are in agreement with our observations based on NOESY experiments (see NMR section), i.e., provided that the iduronate rings are not present at the terminal positions, the 4C1 chair conformation does not participate in it [68-69] (or its contribution is too low as to be detected in NOESY 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 98 experiments). Thus, the conformational sampling of this ring is restricted to the 1C4 chair and the 2SO skew-boat conformations, with the latter being involved in the characteristic pseudorotational equilibrium of the hexopyranose ring (see Chapter 1, figure 3). Furthermore, the balance of the chair to skew-boat equilibrium in internal iduronate residues depends on both their 2-O-sulphation and the substitution pattern of adjacent glucosamine residues[262]. In this regard, our experimental results for Tri1-Tri8 suggest (see NMR section) that the population of 2SO conformer is strongly modulated by the presence or absence of 6-O-sulphation at the reducing terminal (at low temperature). Molecular dynamics is a powerful and widely employed technique to study the molecular conformation and dynamics, allowing a better interpretation of the experimental observables. Combined to NMR, it is especially useful in estimating distributions of conformers in equilibrium[78, 85c, 263]. For the case of L-IdoA2S, since the characteristic inter-conversion rate between the 1C4 and 2SO conformers occurs within the microsecond time scale[17], and therefore, far from being achieved by currently accessible simulation time, two alternative approaches can be used: either consider two different starting geometries for the iduronate residues and run two independent molecular dynamics simulations (unrestrained; free-MD), or use experimental observables as constraints in a single time-averaged restrained MD simulation (tar-MD). Based on these grounds, we have studied the conformation and dynamics of Tri1-Tri8 by both unrestrained and time-averaged restrained molecular dynamics simulations (free-MD and tar-MD, respectively), with explicit TIP3P[264] water molecules. Thus, we run a 20-nanosecond long free-MD simulation for each trisaccharide and for both L-IdoA2S conformations (1C4 and 2SO) as starting geometries, resulting in 16 independent MD simulations. Furthermore, an 8-nanosecond long tar-MD simulation with the iduronate ring adopting an initial 1C4 chair conformation was also accomplished for each trisaccharide, using the NOE-derived H2b-H5b distance (exclusive NOE of the iduronate 2SO pucker) as a sole constraint. It has to be noted that by the time we carried out the unrestrained MD simulations, there were not any force field for carbohydrates which included specific parameters and set of charges for sulphate and/or sulphamate groups. To overcome this technical limitation, the strategy followed was to combine the available parameters for sulphates and sulphamates (Altona´s[265]), which include a explicit hydrogen bond term with a Lennard-Jones 10-12 type potential[197, 266], with other sets of parameters for the carbohydrate moiety, water molecules and counterions developed under the same philosophy for consistency. Thus, the force fields Parm91[267] of Amber (for the water molecules and counterions) and Glycam93[268] (for the carbohydrate moiety) were used. On the other hand, the tar-MD simulations were carried soon after the parameters and partial charges for sulphate and sulphamate groups had been released in the framework of Glycam06[220] force field, the latter entailing a significant improvement compared to Glycam93[268] force field (first version of GLYCAM). Thus, Glycam06[220] together with Amber99SB[269] parameters and partial charges were employed (both of them lack the explicit Lennard-Jones 10-12 term for hydrogen bonding). 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 99 It is important to highlight that, to simplify, from now on we will use the term equatorial conformers (or puckers) to refer to the set of conformers of the equatorial region of the Cremer-Pople sphere of an hexopyranose ring (2SO, 2,5B, B3,O, 3S1, etc.; see Chapter 1). On the other hand, when we say pure 2SO we will mean just the 2SO conformation. Local conformation: plasticity of the L-IdoA2S ring Regarding the free-MD approach, to quantitatively determine the populations of 1C4 and 2SO puckers of the iduronate rings, we first monitored their four vicinal proton-proton dihedral angles (H1-C1-C2-H2, H2-C2-C3-H3, H3-C3-C4-H4 and H4-C4-C5-H5) along the 40000 frames of each 1C4 and 2SO trajectory obtained (when a conformational transition was observed, only the previous frames were considered; see Methodology). These values were turned into vicinal proton-proton coupling constants, 3JHH, by using the Haasnoot-Altona equation[270], which takes into account both the electronegativity and the orientation of the substituents on the H−C−C−H fragment, and averaged for each of the models (L-IdoA2S in 1C4 or in 2SO conformation; see Appendix). Next, to obtain the populations of conformers (1C4 and 2SO) of the L-IdoA ring in each trisaccharide, we carried out an iterative fit of the theoretical and experimental 3JHH according to 4 equations (one for each J-coupling; eq. 1) of the form     )()( 2 )1( 3 4 1 )1( 3exp )1( 3 2 4 1O MD mm S MD mm C mm SJfCJfJ O  Eq. 1 In this expression, f(1C4) and f(2SO) are the unknown molar fractions of each conformer,     MD mm J)1( 3 are the averages from MD and m is an index that runs from 1 to 4. Therefore, the experimental 3JHH coupling constants were considered as averages of the MD-derived ones for each conformer, weighted on the molar fraction of each. As the theoretical values were averages from MD simulations, they implicitly reflected the fluctuations around canonical conformations, which must be considered for this flexible hexopyranose ring[271], particularly to account for the pseudorotational conformational space in the case of the skew-boat conformer (2SO). In addition, the experimental measurements of 3JHH values at five different temperatures (278, 288, 298, 308 and 318 K) allowed us to monitor the population of conformers as a function of temperature (figure 5). 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 100 Figure 5. Temperature dependence of the populations of equatorial conformers of the central L-IdoA2S ring in the eight trisaccharides, classified by two chemical series, Tri1-Tri4 (empty red symbols, red lines), and Tri5-Tri8 (empty blue symbols, blue lines), and obtained by iterative fit of the NMR-derived and free-MD-calculated 3JHH. Also, the equatorial puckers populations predicted by tar-MD simulations at 278 K are shown (filled symbols). Note that the Tri1-Tri4 compounds showed significantly higher populations of equatorial conformers with both methods. From the tar-MD approach, the populations of 1C4 and equatorial conformers (only at 278 K) were determined by directly tracking the evolution of the Cremer-Pople puckering coordinates θ and  over time, with θ undergoing transitions from the south pole (180º, 1C4) to the equator (90º, equatorial puckers) and  fluctuating according to the pseudorotational equilibrium of the iduronate equatorial conformers ( is undefined in the poles). In principle, we thought that this approach was more accurate than the free-MD one as the latter may give rise to important deviations in the calculated values because it is subjected to 1) the experimental uncertainty of the coupling measurements, 2) the higher force field deviations of Glycam93 force field, and 3) the goodness of the least-squares fit.  The results obtained from both methodologies are shown in figure 5, classified by the two series of trisaccharides (Tri1-Tri4, in red; Tri5-Tri8, in blue) and method (tar-MD, filled symbols; free-MD, empty symbols). Thus, the data revealed clear differences in the distribution of populations of conformers of the L-IdoA2S ring between both series, with Tri5-Tri8 ensemble showing dependence with temperature (free-MD). Explicitly, conformational differences (figure 5) that are specific to the presence (Tri1-Tri4 series) or absence (Tri5-Tri8 series) of the 6-O-sulphate group on the reducing GlcN ring were observed. Particularly, the substitution of this bulky charged group with a neutral hydroxyl group (Tri5-Tri8 series) made the populations of equatorial puckers very sensitive to 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 101 temperature (figure 5). At low temperature (273 K), the 1C4 conformer was favoured in a larger extent in the Tri5-Tri8 ensemble (65–80 % from free-MD; 55-60 % from tar-MD), while the equatorial conformers were promoted in the Tri1-Tri4 series (47–43 % from free MD), with tar-MD simulations predicting their majority presence in solution at 273 K (53-58 %). As the temperature increased, the differences between both ensembles started to vanish (figure 5, free-MD). Focusing on the tar-MD data, they indicated that, as long as a 6OSO3group was present at the reducing end GlcN residue, the equatorial puckers mostly populated the conformational space sampled by the iduronate ring in heparin derivatives (figure 5, Tri1-Tri4), with the pure 2SO conformer being in all cases the predominant among the other puckers of the equator (table 3). Thus, whereas the presence of this functional group enhanced the 2SO population above 50%, its lack pushed the equilibrium towards a majority of the 1C4 chair pucker (figure 5, Tri5-Tri8). It is noticeable that, in all cases, the populations of equatorial conformers obtained were higher than the 3JHH derived ones (free-MD), particularly in the Tri5-Tri8 series (figure 5). Furthermore, when comparing the impact of 6-O-sulphation (Glc(A)) on the equatorial puckers populations of the L-IdoA2S ring (Tri1-Tri4 versus Tri5-Tri8 series), tar-MD simulations predicted smaller differences in pairs Tri1:Tri5 (12%) and Tri3:Tri7 (14%) compared to the results obtained from 3JHH fit (20 and 23 %, respectively; see Appendix). Interestingly, Tri8 tendency towards equatorial conformers compared to its partners Tri5, Tri6 and Tri7 was in agreement with free-MD data, i.e, Tri8 presented the highest population of equatorial puckers among the non-6-O-sulphated (at the reducing end GlcN residue) trisaccharides. In addition, and again in agreement with the 3JHH derived results, Tri1 presented the lowest population of equatorial conformers among the 6-O-sulfated trisaccharides (Tri1-Tri4). So, whereas the D-GlcNS,6S-L-IdoA2S-D-GlcNS,6S sequence (Tri1) showed the lowest tendency within its group (Tri1-Tri4) to populate the equatorial puckers, for the D-GlcNS-L-IdoA2S-D-GlcNAc sequence (Tri8) the opposite conformational behavior was observed. These correlations with the sulphation pattern, obtained from both tar-MD and unrestrained MD approaches, indicated that as long the reducing end GlcN residue was 6-O-sulphated, the simultaneous presence of Nand 6-O-sulphation in the reducing and non-reducing terminal, respectively, promoted in some extent the 1C4 chair conformer (figure 5, Tri1). On the other hand, when the reducing GlcN ring contained a 6-OH group (Tri5-Tri8), the population of equatorial puckers was enhanced provided that both the Nand 6positions at the reducing and non-reducing terminals, respectively, were not sulphated either. Thus, although the reducing end GlcN 6-O-sulphation clearly shifts the conformational equilibrium of the L-IdoA2S ring towards the 2SO skew-boat conformer, the other substituted positions are also playing a role in the modulation of this equilibrium. In this regard, it is significant the small difference observed between Tri1 and Tri8 (6%), within the same range as with other partners of their Tri1-Tri4 and Tri5-Tri8 ensembles, respectively. For these reasons and because of the discrepancies between tar-MD and free-MD data regarding specific N- (GlcN(A)) and 6-O-sulphation (GlcN(C)), we suggest that higher level of theory calculations are needed to unambiguously and accurately determine the “subtle” contributions of the N-sulphation and 6-O-sulphation of the reducing end GlcN and the non-reducing end GlcN residues, respectively, on the conformational equilibrium of the L-IdoA2S ring. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 102 From what we have discussed above, it seems reasonable to think that at the origin of the enhancement of the equatorial conformers upon 6-O-sulphation (Glc(A)) the change in the internal dynamics of the L-IdoA2S ring should play a role. Based on this hypothesis, the presence or absence of the 6-O-sulphate (Glc(A)) should be reflected on the reorientation properties of the torsions defining the L-IdoA2S puckering. Thus, we analyzed the internal auto-correlation functions (Cint(t)) for the vectors between vicinal protons of the L-IdoA2S residue (H1-H2, H2-H3, H3-H4 and H4-H5). We took fragments of the trajectories with the L-IdoA2S ring in 1C4 conformation and just before a conformational transition towards the equatorial puckers (we might expect to observe more evident variations on Cint(t) prior to a conformational change). Comparing Tri2 and Tri6, the results (Figure 6; see also Tri1:Tri5 pair in Appendix), indicated a significantly higher flexibility (lower value for the plateau) for the H2-H3 vector when the 6-O-sulphate group (Glc(A)) was present (Tri2), being the most flexible among the proton-proton vectors for Tri2. This correlates to the fact that the H2-C2-C3-H3 torsion participates in the 1C4  2SO conformational change to the largest extent. Therefore, 6-O-sulphation at the reducing-end GlcN residue seem to induce some strain on the H2-C2-C3-H3 torsion of the L-IdoA2S ring so that it promotes some additional flexibility on it (Figure 6; see also Appendix) that might be responsible of the higher tendency of Tri1-Tri4 to populate the equatorial conformers. Interestingly, the H4-H5 vector is, on the contrary, significantly more rigid when a 6-O-sulphate group (Glc(A)) is present, so that it seems to compensate the higher flexibility of the H2-H3 one, in agreement with the law of equipartition of energy states[272]. Figure 6. Comparison of the internal correlation function (Cint(t)) for the vicinal proton-proton vectors of the IdoA2S ring, in Tri2 (red line) and Tri6 (dark blue line) compounds. Only parts of the simulations with the iduronate ring in 1C4 chair conformation have been used. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 109 with an average occupancy of 57%, an average lifetime of 2.6 ps, an average distance of 2.4 Å and an average angle of 32º (table 4). Furthermore, it was not influenced by the conformation of the L-IdoA2S residue or the sulphation pattern. Also, another conserved hydrogen bond was observed between oxygen O5 of the non-reducing end GlcN and the hydrogen H3O of the L-IdoA2S residue. Nevertheless, this hydrogen bond presented, in average, low lifetimes and percentages of occupancy, below 1 ps and 20 %, respectively (table 4). Note that both the rdf and hydrogen bonds analysis have been based on the free-MD results. Tar-MD simulations were not consider for the analysis of solvation and hydrogen bonding properties due to the short lifetimes of the L-IdoA2S puckers (fast conformational transitions), which did not permit to reliably study the differential effect of iduronate puckering on them. 0 2 4 6 8 10 12 0,0 0,5 1,0 Radial distribution function Distance (Å) L-IdoA2S 1C4 L-IdoA2S 2SO Figure 8. Comparison of the solvation profile of the inter-glycosidic oxygen of the IdoA2S-GlcN linkage (O4c) for both conformations of the L-IdoA2S residue. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 110 Table 4. Set of averaged values for the percentage of occupancy, distance, angle and lifetime of the hydrogen bond O5(L-IdoA2S)-H3O(reducing-end GlcN) for the 16 trisaccharide models. Note: the hydrogen bond distance represents that between heteroatoms. The percentage of occupancy is defined as the number of frames, out of a hundred, of the trajectory for which a hydrogen bond exists under the chosen criteria (distance and angle cutoff of 3.0 Å and 120º, respectively, in this case). Compound L-IdoA2S conf. Occupancy (%) Distance (Å) Angle (deg) Lifetime (ps) Tri1 1C4 62.9 2.8 ± 0.2 31.2 ± 12.7 3.5 ± 6.8 2SO 24.6 2.9 ± 0.3 32.5 ± 14.0 2.1 ± 4.1 Tri2 1C4 64.1 2.8 ± 0.2 31.3 ± 12.9 3.2 ± 6.0 2SO 45.7 2.9 ± 0.2 32.7 ± 13.8 2.3 ± 4.4 Tri3 1C4 60.7 2.9 ± 0.3 32.3 ± 13.7 2.3 ± 4.8 2SO 64.5 2.9 ± 0.3 33.0 ± 13.8 2.2 ± 4.2 Tri4 1C4 67.6 2.8 ± 0.2 31.7 ± 12.9 3.2 ± 6.1 2SO 64.4 2.9 ± 0.2 33.0 ± 13.6 2.3 ± 4.4 Tri5 1C4 59.3 2.9 ± 0.2 32.2 ± 13.8 2.3 ± 4.5 2SO 60.1 2.9 ± 0.3 32.9 ± 14.0 2.1 ± 4.1 Tri6 1C4 68.1 2.9 ± 0.2 32.0 ± 13.1 2.9 ± 5.4 2SO 24.4 2.9 ± 0.3 33.1 ± 14.0 2.0 ± 3.8 Tri7 1C4 67.5 2.9 ± 0.2 31.1 ± 13.0 3.1 ± 6.1 2SO 55.3 2.9 ± 0.3 32.5 ± 13.9 2.2 ± 4.4 Tri8 1C4 68.4 2.8 ± 0.2 31.6 ± 12.8 3.2 ± 6.1 2SO 61.9 2.9 ± 0.2 32.9 ± 14.0 2.1 ± 4.1 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 111 Hydrodynamics For polysaccharides in general and, in particular, for glycosaminoglycans, the molecular motion and flexibility must be considered when discussing the molecular conformation, since the frequency, amplitude and geometry of such motions directly affect the measured NMR spectroscopic parameters[66b]. To formally describe a molecular motion, the use of a reorientational correlation function, or its corresponding density function, is required. These functions represent the time-dependent loss of orientational “memory”, for any given vector in the molecular framework, with respect to its initial orientation, as a result of the internal motions and the overall molecular reorientation (the latter typically occurring on the picosecond to nanosecond timescale). While they take a value of 1 (maximum probability) at time zero, over time they decay to a “plateau” which can take a minimum value of zero. If motion is spatially restricted, the correlation function will decay to a value higher than zero, so that values closer to one are associated to rigid motions and those closer to zero are related to flexible reorientations. This is called the parameter of order S2, which reflects the amplitude of motion (or freedom of reorientation) for the associated vector. Another parameter frequently used to deal with molecular motion is the correlation time, which is defined as the area under the correlation function. On the other hand, the spectral density function represents the frequency spectrum analogue of the reorientational correlation function, so that knowledge of one implies knowledge of the other. In our research field, it is usual to employ the model free approach of Lipari and Szabo[275] to analyse fast internal molecular motions. Within this approach, internal motions give rise to the exponential decay of the correlation functions to a “plateau” value less than one. The validity of this model relies on a much faster timescale for any internal motion than for the overall reorientation (τint<<τ0), so that the slower process of overall tumbling in solution subsequently causes a further loss in the orientational correlation. We determined the overall correlation time (τ0) for each trisaccharide from the analysis of the internal molecular motions considering the model free approach[275] (see Methodology for details). The results (Figure 9) showed a very significant increase of τ0 for the Tri1-Tri4 ensemble compared to the Tri5-Tri8 one. Thus, the removal of the 6-O-sulphate group on the reducing GlcN ring gave rise to the reduction of the global correlation time and, therefore, a faster reorientation in solution for the trisaccharides belonging to the Tri5-Tri8 ensemble. Interestingly, we observed that distinct sulphate groups exert a different influence on the overall correlation time. In particular, comparing 6-O-sulphation at the non-reducing and reducing terminal indicated that the former substitution augments τ0 to a smaller extent than the latter provided that the reducing end GlcN residue is N-sulphated (Tri1: Tri3 and Tri5:Tri7 pairs). On the contrary, when this is N-acetylated (Tri2:Tri4 and Tri6:Tri8 pairs) the overall correlation time is similarly enhanced. Regarding N-sulphation at the reducing terminal, its impact on τ0 was lower in all cases compared to the 6-O-sulphation at the same terminal (Tri1:Tri2, Tri3:Tri4, Tri5:Tri6 and Tri7:Tri8 pairs). 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 112 Figure 9. Population-weighted (at 278 K) overall correlation times (τ0) calculated from unrestrained-MD simulations for the differently substituted Tri1-Tri8 compounds. The exact τ0 values determined are shown on each point. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 113 3.2 An inactive hexasaccharide sequence for the FGF-1 mitogenic activity 3.2.1 Background FGF-1 is a member of the Fibroblast Factor family that, by forming a ternary assembly with heparin/heparan sulphate (HEP/HS) and the extracellular domain of membrane receptor FGFR2, triggers a signal that leads to different cellular essential functions (regulation of embryonic development, homeostasis and regenerative disorders). The key step for the activation of the FGF-1 signalling pathway is the formation of such a ternary complex (FGF1-HEP/HS –FGFR2) which gives rise to the dimerization of the receptors and, subsequently, the autophosphorylation that activates a mitogenic response through an enzymatic cascade. Previously in our group, in the context of a wider research programme about the factors that regulate the activation of FGF - FGFR signalling pathway by glycosaminoglycans, some hexaand octasaccharides containing the GlcN-IdoA repeating unit of the major sequence of heparin with diverse substitution patterns were prepared, tested and solved their structures[85]. Among them, we will focus our analysis in this chapter on the three hexasaccharides shown in Figure 10. While one of them represents the heparin regular region (Hexa1), the other two, Hexa2 and Hexa3, display a non-axially symmetric sulphate distribution. On one side, the size of these molecules was chosen as the minimal chain length to be expected to stimulate FGF-1-induced mitogenic activity that could be obtained with a reasonable synthetic effort. On the other hand, the sulphation pattern was varied to obtain distinct distributions of electrostatic potential and assuming that Hexa1-Hexa3 compounds would adopt a helix-like conformation in solution as found for heparin-like GAGs[66b]. This was later confirmed by NMR spectroscopy and MD simulations for Hexa1[85a] and Hexa2[276]. However, the conformational analysis of Hexa3 has not been carried out until now. Hexa2 was designed based on the X-ray structure published by DiGabriele et al.[87], which indicated that for heparin oligosaccharides to interact to FGF-1 the formation of a trans dimer was necessary (each heparin side interacting with one FGF-1 molecule). Thus, Hexa2, displaying sulphate groups on just one side of the helix, would bind FGF-1 in a much less extent (or would not interact) compared to Hexa1, so that the induced mitogenic activity would be greatly diminished. Surprisingly, Hexa2 activated FGF-1 as effectively as an octasaccharide of the heparin regular region[67]. Furthermore, it was later demonstrated that GAGs induced FGF-1 dimerization either in a cis or trans disposition with respect to the heparin chain is not an absolute requirement for biological activity[88]. On the other hand, Hexa3, which presents a pseudo-palindromic relationship with Hexa2, was synthesised[85d] with the aim to interact simultaneously with FGF-1, at sub-site a, and the receptor FGFR (see Chapter 1 and figure 17). Thus, the mitogenic activity would be maximized, as it was proposed by Pellegrini from the analysis of several crystallographic structures of FGF-1 and/or 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 114 FGFR-heparin oligosaccharides complexes[277]. However, when Hexa3 was subjected to the biological assay, it resulted to be inactive[85d]. Interestingly in contrast to the highly active Hexa2 sequence, this result demonstrated that the presence of the recognition sites in the oligosaccharide sequence is not enough to trigger the biological process, putting on evidence the complexity of FGF-1 activation. Furthermore, a previous study employing different synthetic oligosaccharides demonstrated that small variations in the sequence, size or sulphation pattern may dramatically impact their capacity to induce mitogenic activity[67]. Since the original work with Hexa3[85d] did not include a detailed structural analysis[85d], we decided to obtain its structure with the highest possible resolution to analyse the potential reasons for its unexpected lack of activity. To do so, the combination of NMR spectroscopy and molecular dynamics calculations was determinant. Figure 10. Structure (right) and schematic (left) representation of 3 different hexasaccharides presenting the heparin regular region (Hexa1), sulphate groups on just one side of the helix (Hexa2), and a pseudo-palindromic relation with the latter (Hexa3). The 3D structure of heparin helix (PDB code 1HPN[65]) has been considered in this 2D structure representation (right), showing the relative disposition of the sulphate groups. 3.2.2 Results and discussion Nuclear Magnetic Resonance It is well known that heparin oligomers longer than tetrasaccharides are anisotropic[278], exhibiting an hydrodynamic top rotor behaviour. Therefore, they rotate with different correlation times along the transversal (short) or the longitudinal (long) molecular axis[278a, 279]. Both Hexa1 and Hexa2 compounds exhibit such hydrodynamic anisotropic behaviour originated in the rigidity of their glycosidic bonds that it could reflect the electrostatic repulsion of the negatively charged groups (sulphates and carboxylates). In previous works, this behaviour was studied in terms of the differences between the perpendicular (τ) and parallel (τ||) correlation times derived from NOE experiments and 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 115 complete 13C relaxation analysis for Hexa1 and Hexa2[279]. About Hexa3, given its high chemical similarity with those, we should expect it to behave anisotropically as well. In an anisotropic molecule, the NOESY or ROESY based distances calculated using the Isolated Spin Pair Approximation (ISPA)[261] using a single distance (and therefore a single correlation time) as a reference are not accurate because they depend on the angle between the interprotonic vector and the molecular axis, which governs the correlation time of the vector. As an alternative method, it was proposed to use several reference distances selected to cover the range of possible orientations with respect to the rotation axis[85a-c]. However, in this case the accuracy needed to distinguish some of the characteristic features of heparin-derived oligosaccharides was lost. Thus, the off-resonance ROESY methodology [172] appeared as a more precise method to calculate distances from dipolar relaxation data without a pre-assigned model of motion (out from the isolated spin pair approximation (ISPA)[261]), allowing to simultaneously obtain the correlation time and distance for each pair of protons. This method relies on the determination of several off-resonance ROESY values by varying the tilted angle of the effective ROESY spin-lock field to achieve enough amount of independent data as to extract the correlation time for each vector and, from this, to calculate each interprotonic distance. For the reasons above mentioned, we employed off-resonance ROESY spectroscopy[172] as experimental approach to study the hydrodynamics properties of Hexa3. Therefore, we could extract the cross relaxation rates σNOESY and σROESY and the effective correlation times, τeff, by measuring a linear combination of NOE and ROE effects controlled by the spin lock offset. This way, we were able to accurately derive proton-proton distances for the anisotropic Hexa3 molecule independently of their relative orientation with respect to the molecular axis. Hexa3 spectra were assigned by standard procedures, i.e., identifying the spin systems of each ring by scalar coupling and interconnecting them via interglycosidic NOEs. On the other hand, arrays of several series of off-resonance ROESY experiments at several mixing times for different tilted angles (6, 10, 20 and 30 kHz) were recorded to determine the distances. All the growth curves for each proton and each tilted angle were linearly fitted. Then, calculating the growth rate of several series of off-resonance ROESY experiments corresponding to different spin lock offsets, the σNOESY, σROESY and τeff parameters were independently obtained for each proton pair (table 5; see Methodology, eq. 3 and 7), and then used to calculate the experimental interprotonic distances (see Methodology, eq. 8). The diffusional anisotropy has also been studied for Hexa3 in terms of the anisotropy factor, i.e., the quotient between the parallel (τ||) and perpendicular (τ) correlation times (τ/τ||). The most precise method to calculate these parameters (τ and τ||) is by COmplete Matrix Relaxation Analysis (CORMA)[280] based on 13C T1 and T2 and heteronuclear NOE measurements. However, this method is time consuming and, reasonably, it can be replaced by considering the relationship between the correlation times of two orthogonal interprotonic vectors of similar distances, with one of them being aligned with respect to the molecular axis. For Hexa3, we have chosen the vectors H1-H2 and H2-H4 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 116 of glucosamines C and E because 1) they are roughly orthogonal to each other and 2) vectors H1-H2 are nearly parallel to the molecular axis (Figure 11). Therefore, the effective correlation times, τeff, of H1-H2 and H2-H4 proton pairs (derived from off-resonance ROESY experiments; see table 5) approximately represent the τ and τ|| correlation times, respectively, associated to Hexa3. Thus, we calculated the anisotropy factor as the ⁄ ratio. The average anisotropy factor obtained for Hexa3 was 1.2 (τ/τ|| equal to 1.2 considering both C and E GlcN rings). The comparison (of the ⁄ ratio) to those previously reported for Hexa1 (τ/τ|| = 1.5 )[85a] and Hexa2 (τ/τ|| = 1.5 )[276], obtained by a similar methodology, clearly indicated that, although Hexa3 exhibits an anisotropic behaviour, it is slightly less anisotropic than Hexa1 and Hexa2 compounds. Figure 11. 3D structure of the average conformation (obtained over 500 ns of MD trajectory) of Hexa3, indicating with arrows the H1-H2 and H2-H4 vectors for ring C (D-GlcNS), roughly parallel and perpendicular to the molecular axis, respectively. NOESY experiments were also registered for Hexa3. The spectra showed a varied behaviour of the iduronate rings (B, D and F), which depended on their position and substitution along the hexasaccharide chain (Figure 12). For instance, the non-reducing end L-IdoA2S residue (ring F) did not give rise to the exclusive H2F-H5F NOE characteristic of the presence of the 2SO skew-boat conformation (Figure 12). On the other hand, two 2SO-exclusive H2-H5 NOE cross-peaks of medium-weak and weak intensity were identified, corresponding to the internal L-IdoA2S (B) and L-IdoA2OH (D) rings, respectively (Figure 12). Thus, we observed that, in Hexa3, the 2SO skew-boat conformer only participates in the conformational equilibrium of the non-terminal iduronate residues (B and D), although to a less extent (weaker NOE intensity) in that of the non-sulphated iduronate ring (D, L-IdoA2OH). Also, the NMR-derived (NOESY and off-resonance ROESY) H2-H5 and H1-H3 intra-ring distances of the iduronate ring (table 6), both in the NOE range for the 2SO conformer, showed a shorter distance for the L-IdoA2S ring B than the L-IdoA2OH ring (D), thus indicating the highest contribution of the 2SO conformer in the equilibrium of the former. This 2-sulphation effect on the conformation of the iduronate residue has been previously reported[10, 281]. In conclusion, the 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 117 NOESY experiments showed that the proportion of iduronate 2S0 conformer decreases along the Hexa3 chain from the reducing to the non-reducing end, being larger in B than in D ring, and absent or negligible in F residue. By combining NMR and MD-derived 3JHH couplings, these differences were translated into populations of 2SO pucker, which are discussed in the next section (table 7). Finally, regarding the geometry around the glycosidic linkages that define the global conformation in oligosaccharides, the presence of the pairs of intense H1’-H3&H1’-H4 and H1’-H4&H1’-H6 interglycosidic NOEs corresponding to the GlcN-IdoA and IdoA-GlcN linkages, respectively, indicated a major syn rearrangement (figure 13). Furthermore, the absence of the H5’-H6, H1’-H5 and H1’-H3 NOEs for the IdoA-GlcN linkages, exclusive of the anti-Ψ rearrangements, confirmed the unique presence of syn-Ψ conformations, and consequently, the rigidity of the Hexa3 backbone (figure 13). Table 5. σNOESY and σROESY cross relaxation rates and effective correlation times (τeff) of the intraand inter-residue proton pairs of Hexa3, calculated from off-resonance ROESY[172] measurements. σNOESY (s-1) σROESY (s-1) eff (ns/rad) H1B-H4A -0,09 0,61 0,50 H1B-H6A* -0,03 0,17 0,51 H1C-H3B -0,05 0,20 0,59 H1C-H4B -0,07 0,40 0,53 H1D-H4C -0,03 0,18 0,50 H1D-H6C -0,07 0,41 0,51 H1D-H6'C -0,05 0,16 0,75 H1F-H4E -0,08 0,69 0,46 H1F-H6E -0,13 0,83 0,50 H1F-H6'E -0,04 0,12 0,85 H1B-H3B -0,04 0,17 0,57 H4B-H5B -0,11 0,54 0,57 H1C-H2C -0,02 0,07 0,65 H2C-H4C -0,01 0,07 0,54 H1D-H2D -0,03 0,19 0,51 H1E-H2E -0,01 0,06 0,62 H2E-H4E -0,01 0,07 0,50 H4F-H5F -0,05 0,47 0,45 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 118 * Values assigned to each diastereotopic hydrogen (proR and proS) under the approximation that the most populated conformer around the GlcN(C) ω torsion is gg (in agreement to MD results, see below). Figure 12. Expansion of a NOESY experiment registered at 400 ms mixing time for Hexa3 showing the signals corresponding to the 2SO-exclusive NOE cross-peaks of the iduronate residues (B, D and F rings). The absence of the H2-H5 cross-peak for the non-reducing end L-IdoA2S residue (H2F-H5F) is indicated with a cross symbol, in red. Also note the very low intensity of the 2SO-exclusive H1-H3 NOE for ring F compared to those of the internal iduronate rings B and D. Figure 13. Expansion of a NOESY experiment registered at 400 ms mixing time for Hexa3 showing the crosspeaks corresponding to a syn rearrangement around the GlcN-IdoA (H1’-H3 and H1’-H4 NOEs, in green) and IdoA-GlcN glycosidic linkages (H1’-H4 and H1’-H6 NOEs, in bold black), respectively. The non-observed H5’-H6 NOEs, exclusive of the anti-Ψ rearrangements around the IdoA-GlcN linkages, are labelled in red. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 125 this hypothesis, the most representative conformation of Hexa3 obtained from MD was manually superimposed (backbone) to the NMR-resolved structure of Hexa2 bound to FGF-1 (PDB code 2ERM). All the different possible superimpositions were performed, starting with Hexa3 trying to fit its five sulphate groups located at the same side of the helix in both sub-sites of FGF-1 by inverting the directionality of the hexasaccharide chain (see Appendix). However, the resulting structures of inverted directionality showed unavoidable steric clashes, indicating that Hexa3 global conformation does not allow to maximize the interactions with FGF-1 (see Appendix), as Hexa2 does. Also, we superimposed the heavy atoms of the GlcNS-IdoA2S-GlcNAc,6S triad, contained in both molecules, in the two possible directionalities for Hexa3, this is, from the non-reducing to the reducing end (shown in figure 17) and its reversed mode (see Appendix), i.e., exchanging the positions of the GlcNAc,6S and GlcNS residues within the triad. Interestingly, this can be done, in principle, because the distances between the sulphate groups are similar. However, the reverse mode (GlcNAc,6S-IdoA2S-GlcNS superimposition) did not fit in sub-site a (steric hindrance) and, even more, the other part of the chain fell away from the binding site (see Appendix). Differently, in the non-reducing to reducing end orientation the GlcNAc,6S-IdoA2S-GlcNS triads of both Hexa2 and Hexa3 compounds presented the same geometry for their glycosidic linkages (because of the same directionality of the linkages), thus allowing the 3 sulphate groups of Hexa3 triad to be properly oriented to fully occupy sub-site a, while remaining the secondary sub-site (b) unoccupied. Therefore, the mode of interaction shown in figure 17 for Hexa3 complexed to FGF-1 is the only possible, thus confirming that this compound, whose chemical design aimed to simultaneously interact to FGFR and FGF-1 through its non-reducing and reducing terminal, respectively, was correctly devised. According to this model, the lack of binding site occupancy at sub-site b would attenuate the interaction between the GAG chain and FGF-1, explaining the lower IC50 value obtained for hexasaccharide 3. On the other hand, it was demonstrated that FGF-1 interaction with Hexa2 made more rigid the residues involved upon ligand binding (entropic cost) and more flexible those amino acids participating in the interactions with the FGF-1 receptor (entropy increase to compensate the entropic cost of ligand binding). Under these premises, we propose that the number of accessible FGF-1 conformations (or orientations of the side chains) able to interact to a FRFR, or number of “active microstates”, would be higher for the FGF-1 in the bound state (FGF-Hexa2), so that the enthalpy-entropy balance for the FGF-FGFR interaction would be favoured and, consequently, the biologically relevant FGF-Hexa2-FGFR ternary complex stabilized (high induced mitogenic activity). Following the same reasoning, for the interaction (weaker) of Hexa3 with just the sub-site a of FGF-1 (figure 17), we hypothesize that the number of accessible FGF-1 “active microstates” for FGFR recognition would be considerably lower, thus destabilizing FGF-FGFR interaction and, as a result, the formation of the FGF-Hexa3-FGFR ternary complex necessary for the signalling pathway to be triggered. We think this provide a reasonable hypothesis for the negligible induced mitogenic activity observed for hexasaccharide 3[85d]. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 126 Figure 17. 3D (left) and schematic (right) representation of the complexes of FGF-1 (ribbons) with Hexa2 (purple sticks; PDB code 2ERM) and Hexa3 (red sticks), the latter manually superimposed (through the sulphate groups) to Hexa2 at sub-site a. The two sub binding sites of FGF-1 (green) are labelled a and b. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 127 3.3 Methodology 3.3.1 Nuclear Magnetic Resonance Library of sulphated trisaccharides (subchapter 3.1) NMR experiments were performed on Bruker DRX 500 MHz spectrometer equipped with 5 mm inverse triple-resonance probe. NMR samples were prepared at pH* 7 in 500–600 or 200 mL in 5 mm or 3 mm tubes, at 2 mm and 6 mm, respectively, in 99.9% D2O and at several temperatures varying from 278 to 318 K. Sizes of acquisition matrices were 2 K_512 for COSY-dqf, gradient selected, experiments and 1 K_256 for TOCSY with mixing time of 80 ms. HSQC were recorded in gradient enhanced versions using echo-antiecho detection both with or without decoupling during acquisition. When it was required presaturation was applied by low power irradiation at water frequency. The preliminary results of the NOESY were unsatisfactory because the molecules were close to the zero crossing point and give very weak peaks. ROESY sequences were also applied but the strong coupling of the protons of the L-IdoA2S residue biased the results. No better results were obtained using T-ROESY. Therefore, we recorded all the NOESY experiments at 278 K to increase the correlation time in order to obtain negative NOE peaks. Additionally, the use of lower temperatures allowed us to exploit the increase of the population of the lowest energy conformations. The build-up experiments were acquired with 1D sequence selected with gradients spin echo (dpfgse)[286] and two spin echoes flanked by bipolar gradients during the mixing time (200, 300, 400, 500, 600, 700, 800, 1000, 1300, 1500 ms). An inactive hexasaccharide sequence for the FGF-1 mitogenic activity (subchapter 3.2) A NMR sample was prepared dissolving 2.5 mg of Hexa3 in 600 μL of D2O 100%, adjusting the pH* to 7.0. 2D experiments were recorded using z-gradients pulses when possible: gs-DQF-COSY,[287] TOCSY[288] using 80ms of spin lock, NOESY,[286, 289] ROESY,[288] HSQC using gradient selected with sensitivity enhanced versions[290], and coupled HSQC. In order to estimate the longitudinal and transversal cross relaxation rates, σNOESY and σROESY, compensated off-resonance ROESY[172] experiments were acquired using different mixing times (100, 200, 300, 400 and 500 ms) and several radiofrequency offsets (6, 10, 20 and 30 kHz) for the spin-lock. The spin locking field was 7.14 kHz and the experiments were carried out at 298 K. The off resonance pulse locks the spin along the effective field, which makes a  angle with the z axis. This angle is defined as follows: )/arctan( 1offset   Eq. 2 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 128 The pure longitudinal and the pure transverse cross relaxation rates (σNOESY and σROESY) were calculated for each proton pair from the dependence of the NOE cross peaks versus the  angle, according to the expression:  22 cos sen ROESYNOESYobs  Eq. 3 On the other hand, NOE and ROE cross relaxation rates are related to interproton distances and it can be expressed as[174a]:   )0()2(6 6JJrISNOESY    Eq. 4   )(3)0(2 6  JJrISROESY   Eq. 5 Transforming the cross relaxation rates in terms of interproton distances (rIS) and effective correlation times for each proton pair (τeff) requires an assumption on the behavior of the spectral density function J(nω), which involves that the motion of two interacting protons can be described by a single exponential. The spectral density can be then written as: 222 42 2 0 1 104 )( eff eff n nJ               Eq.6 where γ is the magnetogyric ratio of the nuclei, ħ is the reduced Planck constant, µ0 is the vacuum permeability and ω is the Larmor frequency (the later related with the magnetic field of the spectrometer). For each proton pair, correlation times and, therefore, proton-proton distances can be calculated from the ratio ROESYNOESY  / by solving the following equations[291]:     42 36493221 2 12    eff Eq. 7 6/1 22 42 2 0 1 3 2 104                         eff ROESY IS r       Eq. 8 We have obtained internuclear distances and correlation times disregarding equal mobility between different proton pairs. Thus, no model of motion was assumed a priori. This method reduces the intrinsic error resulting from the use of an internal reference (e.g. in the Isolated Spin Pair Approximation[261]) and is appropriate for anisotropic molecules. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 129 3.3.2 Modelling Library of sulphated trisaccharides (subchapter 3.1) A. Unrestrained MD (free-MD) Input preparation In all cases, the starting geometries were generated from the available data[65] deposited in the Protein Data Bank (PDB code 1HPN) and modified accordingly. The topologies were built with PREP-LINK-EDIT-PARM module of Amber 5.0, employing the residues and the set of partial charges published by Perez et al[236b] (the latter developed under the context of the set of parameters for carbohydrates PIM[292]) and the force fields parm91[266b] of Amber and glycam_93[268] together with the set of Altona parameters for sulphates[265]. Two independent starting geometries of each heparin-like trisaccharide structure were built, one with the IdoA2S residue in the chair 1C4 conformation and one with the IdoA2S in the 2SO skew boat geometry. Each of these models was immersed in a 41Å-sided cube of pre-equilibrated TIP3P water molecules. Molecular dynamics MD simulations were run on the Finis Terrae cluster belonging to the Centro de Supercomputación de Galicia (CESGA), Spain, taking advantage of the prioritized computing time we were awarded (ICTS-2010-ID119). To equilibrate the system we followed a protocol consisting of 10 steps. Firstly, only the water molecules were minimized, and then heated to 300 K. After, the water box together with the sodium ions were minimized and then followed by a short MD simulation (3 ps). At this point, the whole system is minimized by four consecutive steps imposing positional restraints on the solute, with a force constant decreasing step by step from 20 to 5 kcal/mol. Finally, an unrestrained minimization (100 steps) was carried out. The production dynamics simulations were accomplished at a constant temperature of 300 K (by applying the Berendsen coupling algorithm[293] for the temperature scaling) and constant pressure (1 bar). The Particle Mesh Ewald Method[267, 282b] (to introduce long-range electrostatic effects) and periodic boundary conditions were also turned on. The SHAKE algorithm for hydrogen atoms, which allows using a 2 fs time step, was also employed. Finally, a 9 Å cutoff was applied for the Lennard-Jones interactions. MD simulations have been performed with the sander module of Amber 6.0, with explicit treatment of the 10 12 hydrogen bond potential, in agreement with the parameters set for sulphates we use (Altona). A total of 16 MD simulations of 20 ns each were obtained. The trajectory coordinates were saved each 0.5 ps. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 130 The data processing of the 16 generated trajectories were done with the ptraj module of Amber 9.0, except for the Cremer-Pople puckering coordinates, which were calculated with the Carnal module of Amber 5.0. The final theoretical 3JHH values were obtained as averages for each of the models (L-IdoA in 1C4 or in 2S0) according to the expression     )()( 2 )1( 3 4 1 )1( 3exp )1( 3 2 4 1O MD mm S MD mm C mm SJfCJfJ O  ( Eq. 1) . Thus, to obtain the populations of conformers (1C4 and 2S0) of the L-IdoA2S ring in each trisaccharide, we performed an iterative fitting of theoretical and experimental J-coupling data (Eq. 1; the 4C1 conformation was disregarded as no experimental support was obtained for it, particularly the exclusive H5c-H5b NOE was not observed). In this equation   4 1C f and   O S f2 are the molar fractions of each conformer,     MD mm J)1( 3 are the averages from MD, and m is an index that runs from 1 to 4. Therefore, the experimental 3JHH coupling constants were considered as averages of the MD-derived 3JHH for each conformer weighted on the molar fraction of each one. As the theoretical values were averages from MD simulations, they implicitly reflected the fluctuations around canonical conformations, which must be considered for this flexible hexopyranose ring[271], particularly to account for the pseudorotational conformational space in the case of the skew boat conformer (2SO). In addition, the experimental measurements of 3JHH values at five different temperatures (5, 15, 25, 35 and 45 ºC), allowed us to monitor the population of conformers as a function of temperature. The overall correlation times (τ0) for the 16 trisaccharide models (8+8 with the iduronate ring adopting the 1C4 chair and 2SO skew-boat conformation, respectively) were calculated according to the model-free approach of Lipari and Szabo[275] from the auto-correlation functions for each proton-proton vector (between vicinal hydrogens), which were derived with the ptraj module of Amber 9.0. Since both the internal and overall motions act on the correlation function of each H-H vector, we first eliminated the translational and rotational components of the molecular tumbling at the global level to obtain the internal auto-correlation functions, for which only the internal motions contribute. This was done by RMS fit the backbone coordinates of each frame on the starting ones, prior to calculate the auto-correlation functions. Thus, the internal auto-correlation functions ( ( )) were fit to the Lipari and Szabo expression ( ) ( ) Eq. 9 This allowed us to obtain the parameters of order (S2) and the internal correlation times (τint). Since the correlation function describing the global motion (assuming an isotropic tumbling), C0(t), is given by the equation ( ) Eq. 10 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 131 the correlation function (C(t)) can be broken down by its different contributions (as long as both the global and internal motions are not correlated), so that ( ) ( ) ( ) Eq. 11 or ( ) ( ) Eq. 12 Therefore, the overall correlation times τ0 were determined by iterative fit to Equation 12, using the S2 and τint values derived from Equation 9. Furthermore, since they corresponded to independent values for the 1C4 and 2SO models, we calculated a representative τ0 value for each Tri1-Tri8 compound by weighting on the populations obtained for both puckers at 278 K (figure 8). B. Tar-MD Input preparation In all cases, the initial coordinates were based on the NMR-resolved dodecasaccharide structure of natural heparin[65] (PDB code 1HPN), following the protocol we have previously described[76]. The starting conformation of the L-IdoA2S unit of each trisaccharide was the 1C4 chair. The topology and coordinates files of every system were built with the tLEAP module of AMBER 11[294] package. All the trisaccharides were neutralized with sodium ions and then immersed in a TIP3P[264] water box, giving rise to systems of about 4000 atoms. The Glycam06g-1 parameters[220] were used to model the sugar moiety, including the sulphate and sulphamate moieties. For the water molecules and sodium ions, the Amber99SB parameters[269] were employed. Furthermore, the partial charges of GLYCAM06[220] were employed for the sugar moiety, adjusting the partial charge on the Oand Natoms bound to the SO3groups according to GLYCAM philosophy for charge development. For the O-isopropyl group, partial charges were derived from the molecular electrostatics potential (MEP) using the RESP method[295] with a constraint of 0.01, for consistency with the procedure employed in GLYCAM06[220] force field development. The HF/6-31G* level of theory was used for both the structure optimization and the MEP calculation. Detailing the procedure employed, a methyl-O-isopropyl and a D-Glc-OMe were built and charge constraints imposed as follows: the total charge of both molecular models is set to 0, the methyl group in the D-Glc-OMe must have a charge of +0.194 whereas the charge of the oxygen involved in the glycosidic linkage is set to -0.194, and both methyl groups are set to be equivalent and to be removed during the last step of charge derivation, and, thus, being both compounds merged to form a D-Glc-O-Isopropil. Additionally, the partial charges on aliphatic hydrogens and on the O-isopropyl group were constrained to 0 and -0.194, respectively, in agreement with GLYCAM philosophy. The standard error obtained was 0.005. It is noticeable that a similar protocol has been successfully applied for the development of parameters for different sugar derived compounds[296]. The quantum mechanical calculations and the RESP[295] procedure were carried out with ante-R.E.D 2.0 and R.E.D IV of the R.E.D web server[297]. 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 132 Molecular dynamics MD simulations run on the Finis Terrae cluster belonging to the Centro de Supercomputación de Galicia (CESGA), Spain, taking advantage of the prioritized computing time we were awarded (ICTS-2011-ID162). MD simulations were carried out with AMBER 11[294]. Prior to include the constraints, we performed an equilibration protocol consisting of an initial minimization of the water box (20000 steps), followed by a minimization of the whole system (10000 steps). Then, the TIP3P[264] water box was heated at constant volume until 278 K using a time constant for the heat bath coupling of 1 ps. The equilibration finished with 200 picoseconds of molecular dynamics simulation without restraints, at constant pressure (1bar) and turning on the Langevin temperature scaling with a collision frequency of 1 ps. Furthermore, non-bonded interactions were cut off at 8.0 Å and updated every 25 steps. Periodic Boundary Conditions and the Particle Mesh Ewald method[282] were turned on in every step of the equilibration protocol to evaluate the long-range electrostatic forces (the grid spacing was approximately 1 Å). The time-averaged restraints molecular dynamics were run with the same settings used in the last step of the equilibration protocol. A sole NOE-derived distance, that between H2 and H5 protons of the IdoA2S residue (table 2), was imposed as time-averaged constraint, applying a r-6 averaging. The equilibrium distance range was set to rexp-0.1Å ≤ rexp ≤ rexp+0.1Å. Trajectories were run at 278 K, with a decay constant of 800 ps and a time step of 1 fs. The force constants k2 and k3 used in each case go from 25 to 45 kcal·mol-1·A-2 (Supporting Information, table S5). The overall simulation length for the simulations was 8 ns. The coordinates were saved each picosecond, thus, obtaining MD trajectories of 8000 frames each. Convergence within the equilibrium distance range was obtained in all cases. The analysis of the tar-MD trajectories has been carried out with the ptraj module of AMBER 11[294], except for the Cremer-Pople coordinates, which were determined with an in-house script (see acknowledgements). The auto-correlation functions (Cint(t)) shown (Figure 9) were obtained following the procedure described in the previous section (A). In this regard, it has to be noted that the Cint(t) functions have been only obtained for the Tri1:Tri5 and Tri2:Tri6 pairs since the others (Tri3:Tri7 and Tri4:Tri8) showed too fast conformational transitions as to obtain the decay of the internal auto-correlation functions with the L-IdoA2S ring keeping the 1C4 conformation. An inactive hexasaccharide sequence for the FGF-1 mitogenic activity (subchapter 3.2) The initial coordinates were based on the NMR-resolved dodecasaccharide structure of natural heparin (PDB code 1HPN). The topology and coordinates files were built with the tLEAP module of AMBER 11[294] package. The system was neutralized with sodium ions and then immersed in a TIP3P[264] water box, giving rise to a molecular system of 4526 atoms. The GLYCAM_06h[298] 3. Structural studies of heparin-like oligosaccharides by NMR and MD techniques 133 parameters were used to model the sugar moiety, including the sulphate and sulfamate groups, and the AMBER ff12SB parameters[269] for the water molecules and calcium ions. The partial charges of GLYCAM06 were employed for the sugar moiety, adjusting the partial charge on the Oand Natoms bound to the SO3 groups according to GLYCAM philosophy for charge development. For the O-isopropyl group, partial charges were derived from the molecular electrostatics potential (MEP) using the RESP method[295] with a constraint of 0.01, for consistency with the procedure employed in GLYCAM06 force field[220] development. The HF/6-31G* level of theory was used for both the structure optimization and the MEP calculation. Detailing the procedure employed, a methyl-O-isopropyl and a D-Glc-OMe were built and charge constraints imposed as follows: the total charge of both molecular models is set to 0, the methyl group in the D-Glc-OMe must have a charge of +0.194 whereas the charge of the oxygen involved in the glycosidic linkage is set to -0.194, and both methyl groups are set to be equivalent and to be removed during the last step of charge derivation, and, thus, being both compounds merged to form a D-Glc-O-Isopropyl. Additionally, the partial charges on aliphatic hydrogens and on the O-isopropyl group were constrained to 0 and -0.194, respectively, in agreement with GLYCAM philosophy. The standard error and relative root mean square error were, respectively, 0.005 and 0.188. It is noticeable that a similar protocol has been successfully applied for the development of parameters for different sugar derived compounds[296, 299]. The quantum mechanical calculations and the RESP procedure were carried out with ante-R.E.D 2.0 and R.E.D IV of the R.E.D web server[297]. The molecular dynamics simulations have been run on a 4-node AMD Opteron Interlagos cluster (2,3 GHz, 16 cores per node). MD simulations were carried out with AMBER 12[300]. The equilibration protocol consisted of an initial minimization of the water box (20000 steps), followed by a minimization of the whole system (10000 steps); finally, the system was heated (40000 steps) at constant volume until 300 K using a time constant for the heat bath coupling of 1 ps. The production dynamics have been carried out at a constant temperature of 300 K, by applying the Langevin thermostat[301] with a collision frequency of 5 ps-1, and at constant pressure (1bar), applying Periodic Boundary Conditions (PBC) and using the Particle Mesh Ewald Method[282] (PME) to account for the long range electrostatic effect (the grid spacing was approximately 1 Å). The SHAKE algorithm[302] was also employed, thus, allowing a 2 fs time step, and non-bonded interactions were cutoff at 8.0 Å and updated every 25 steps. The equilibration protocol and production dynamics have been performed with the sander.MPI and pmemd.MPI modules of AMBER 12[300], respectively. One MD simulation of 500 ns have been performed, saving the trajectory coordinates each picosecond. The analysis of the MD trajectory has been done with the ptraj module of AMBER 12[300]. The Cremer-Pople coordinates were calculated with a script from R.J. Wood´s Group (see acknowledgements).