Simulation of the Structure of Sugar Chains of Glycoproteins

Deciphering the biological code for functions of oligosaccharides (glycans) represents the third frontier, after such characterization for proteins (oligopeptides) and DNA and RNA (oligonucleotides). At the level of glycan chemistry, i.e. composition, sequence, 3-D structure and biophysical properties, experimental "windows" are still relatively undeveloped. Although attachment of an oligosaccharide sidechain (glycosylation) to proteins is the normal state for most extracellular, and many intracellular, proteins, how such "decoration" affects the structure, dynamics, and functional properties of proteins is largely unknown. Indeed it is largely ignored, as nearly all high-resolution (NMR, x-ray crystallography) structural and other experimental molecular studies are performed on "undressed" protein, i.e. without attached glycans, for reasons of experimental tractability (avoid heterogeneous protein) or convenience (bacterial expression systems lack glycosylation machinery). For these reasons, computer simulation offers a particularly rich opportunity to complement the large amount of protein structural data, by studying the properties of proteins reclothed with their sugar chains. Simulation offers a second major advantage in being able to deal easily with glycan heterogeneity; glycans with variant chain length, sugar composition, branching pattern, and linkage type, as found experimentally by MS sequencing, can be simulated. Our interest started with prion protein (PrP), the unusual protein for which an aberrant conformational form (PrPSc) is associated with diseases such as mad-cow disease (BSE) and CJD in humans. Although the correlation of PrP glycosylation state with disease characteristics is well-known, NMR structures and other experimental molecular studies have all been done on unglycosylated protein, and, remarkably, interpreted without acknowledgement of their normal glycosylation. Such information is critical to understanding the potential of PrP to interact with other proteins and ligands, which underlies its normal, but as-yet-still-unknown, function, and its capacity to change into other conformational states with apparently neurotoxic or, unprecedented, infectious properties.

Principal Investigator

Dr Jill E. Gready
Computational Proteomics Group
John Curtin School of Medical Research


x04, d55

Facilities Used



Dr Johannes Zuegg
Computational Proteomics Group
John Curtin School of Medical Research
ANU; and
Alchemia, Brisbane
Mr Saburoh Sakai
RIKEN Genomics Science Center

RFCD Codes

250106, 250699, 270105, 270199

Significant Achievements, Anticipated Outcomes and Future Work

Molecular dynamics simulations have been undertaken to characterize the changes in the conformational flexibility and surface properties of PrP when sugar chains (N-glycans) are added to the two known glycosylation sites, and to define the properties of the sugar chains themselves. As PrP is also anchored, at its C-terminus, to the membrane by a GPI (glycophosphatidylinositol)-anchor, a GPI anchor was also added. Using NMR structures to develop starting coordinates, our initial work established suitable simulation conditions for unglycosylated PrP [1]; the model contained the structured C-terminal region of ~110 residues only as, unusually, the complete N-terminal region is shown to be unstructured by NMR. Defining appropriate simulation conditions was not straightforward due to the unusually open structure of PrP, which has a relatively low composition of formal secondary structure (only ~50%) and a very large number of charged residues on the protein surface; at effective pH 7, a net charge of -1 e.u. from 15 positively- and 14 negatively-charged residues. Simulations showed high sensitivity to the correct treatment of the electrostatic interactions, with stable trajectories achieved only by including all the long-range electrostatic interactions using the (computationally expensive) Particle Mesh Ewald (PME) method. These initial results on the importance of electrostatics and, particularly, stability of salt bridges had instructive implications for PrP structure and stability: PrP encounters a range of pHs physiologically during normal or abnormal recycling from the neutral pH membrane surface into endosomes or lysosomes (acidic pHs). This study was the first reported MD simulation for PrP.

Simulations were then undertaken on PrP models with attached sugars and GPI anchor, to probe their effects on the protein stability, electrostatics and dynamics and to clarify which regions of the protein might be occluded from solvent and approach of potential ligands or other proteins, either from the attached sugars themselves or by proximity to the membrane [2]. The N-glycans are highly branched, making up nearly 1/3 of the total mass. A large box of water containing >10,000 water molecules was required to solvate the total molecular system (~33,000 atoms in total), and trajectories were run for several nanoseconds (ns), i.e. these are computationally very demanding simulations. The results showed PrP is stabilised overall from addition of the glycans, but that the stabilization appears indirect, by reducing the mobility of the surrounding water molecules, and not from specific interactions such as H bonds or ion pairs. As the attached glycans have three negatively charged SiaLex groups each, the surface electrostatic properties of PrP are totally changed. A negative electrostatic field covers most of the surface, including the whole surface of two of the three helices, but, significantly, the unusual totally hydrophilic Helix-A is not affected. Modelling showed Helix-A could readily dimerise in anti-parallel fashion: this model was consistent with indirect evidence for how the experimentally observed dimer might form. Together with separate simulations of the GPI anchor in a membrane model, the results showed a highly flexible GPI anchor, which would maintain the protein at a distance between 9 and 13 Å from the membrane surface, with little influence on its structure or orientational freedom. The work was featured on the cover of Glycobiology, and the graphics also won the Todays Life Sciences award for best modelling picture in 2000.

As conversion between two different conformations of PrP, "normal" PrPC and "scrapie" PrPSc, is suggested to be part of the disease-causing mechanism, extended simulations (4 ns) were next conducted at elevated temperatures (effective 400 and 500 K, compared with previous 300 K) to investigate possible unfolding pathways [3]. Note that the simulation time is still much too short to simulate any transition state towards a PrPSc-like structure. Using Ca-distance matrices in combination with principal component analysis (PCA) we were able to identify regions of the protein which collapse first or, conversely, form fairly stable domains. Helix-C and Helix-A showed high intrinsic stability, even at 500 K, whereas Helix-B showed stability only in the glycosylated model, which might be attributed to its attached N-glycan. The core of the stable structure at 500 K for the glycosylated model, is formed mainly by Helix-B and Helix-C; contact maps and PCA clearly identifies this region to be the last to unfold. Although the simulations showed apparently different unfolding pathways for the glycosylated and glycosylated models, multiple trajectory simulations are necessary to clarify whether this is really the case; some such simulations have been undertaken in conjunction with the following study.

In this most recent study we have combined the methods and experience of this project with those of project x11, which is developing the linear-scaling semi-empirical QM method, MOZYME [4]. This is apt as electrostatics appear important for PrP structure and function, and MOZYME offers the unique opportunity to calculate molecular electrostatic potentials (MEPs) of the whole protein at a QM level. Also, as another study on a protein ion channel recently undertaken at the ANU [Bliznyuk et al. (2002) J Phys Chem B 105, 12673] had found artifactual results for EPs calculated from MM force fields compared with those calculated by MOZYME, we were interested to compare these two calculational methods applied to a different, and more taxing, protein problem. An extensive simulation series for models of wild-type and 13 mutants of human PrP, 11 of them disease-associated, and including some with glycosylation, was undertaken. The results did not show differences in MEPs which would seem likely to contribute to the disease association of the mutants. This negative result is consistent with a now-huge amount of experimental structural and biophysical data on these mutants which has similarly failed to show a link between the mutations and likely properties, such as reduced general or local stability or gross changes in surface or electrostatic properties. However, the study did show large differences between MEPs calculated by the several methods. Although the implications of these results require further analysis, in general they sound a warning against drawing quantitative conclusions about protein-ligand or protein-protein interactions from calculations where electrostatic (i.e. long-range) components dominate.

Future work includes a detailed study of the water mobility from the glycosylated and non-glycosylated PrP simulations in order to gain further insight into factors influencing both protein and oligosaccharide flexibility. We are also starting to consider the effect of N-glycans on the structure and properties of another extracellular domain.


Computational Techniques Used

The SANDER module in Amber 5, with the all-atom AMBER94 force field for the protein and an adapted GLYCAM93 force field for the glycans, was used for the simulations, using periodic boundary conditions and constant volume. The PME (Particle Mesh Ewald) method was used to calculate long-range electrostatic interactions. Both unglycosylated and glycosylated models were immersed in a rectangular box of pre-equilibrated TIP3P water molecules (~5,000 and ~10,000, respectively) giving total system sizes of ~15,000 and 33,000 atoms, respectively, with simulation times of ~1.5 to 4 ns for studies 1-3, and 400-500 ps for study 4. Simulation analysis was done for: trajectories with the program CARNAL in Amber 5; secondary structure with the DSSP program; hydrophobic interactions using the VORONOI method; PCA using the fast NIPLAS algorithm; cluster analysis using the Average Linkage method; surface area and surface properties using the program NSC; model manipulation and visual analysis using INSIGHTII; picture generation using MOLSCRIPT, RASTER3D, MOLMOL and POVRAY. MOZYME calculations were performed with MOPAC2000. For study 4, the MEP was calculated on a 65x65x65 grid, with spacing of 1.5 Å. For the AM1 method, a modified PMEP module in MOPAC2000 was used, while for AMBER charges the locally-written program GMIX was used. The modified PMEP module allows MEP calculation on a grid, as well as ability to add a contribution from the induced surface charges, for calculations on solvated molecules. Calculation of the MEP was also done with the DELPHI program, using default van der Waals radii and the same partial charges as for the Amber simulations, and also with partial charges obtained from the AM1 calculations in MOZYME.

The status of MOPAC on the VPP300, PC and APAC National Facility Compaq SC is given in the x11 report. Amber 5, particularly the SANDER module and PME option were highly vectorized (VPP300) or otherwise tuned for the VPP and PC. Some work has been done by National Facility staff to tune Amber (now) 6 for the Compaq SC, but further work is on hold pending release of a parallelized Amber 7 by its developers.


Publications, Awards and External Funding

1. J Zuegg, JE Gready. Molecular dynamics simulations of human prion protein: importance of electrostatic interactions for protein stability. Biochemistry 38, 1999, 13862-13876.

2. J Zuegg, JE Gready. Molecular dynamics simulations of human prion protein including the N-linked oligosaccharides and GPI anchor. Glycobiology 10, 2000, 959-974.

3. J Zuegg, JE Gready. Molecular dynamics simulations of human prion protein at high temperature. M/S in revision, 2002.

4. J Zuegg, AA Bliznyuk, JE Gready. Molecular electrostatic potential calculations on wild-type and mutants of human prion protein. M/S in preparation, 2002.

Study 4 was funded by a FECIT (Fujitsu European Centre for Information Technology) contract to ANUSF, as part of promotion activities for MOPAC2002 (demonstration of the capabilities and usefulness of MOZYME).