Computer Simulations of Membrane–Sugar Interactions Jon Kapla

©Jon Kapla, Stockholm University 2016 ISBN 978-91-7649-363-2 Printed in Sweden by Holmbergs, Malmö 2016 Distributor: Department of Materials and Environmental Chemistry, Stockholm University The cover art was created with: VMD (http://www.ks.uiuc.edu/Research/vmd/), Inkscape (https://inkscape.org)

Abstract

Carbohydrate molecules are essential parts of living cells. They are used as energy storage and signal substances, and they can be found incorporated in the cell membranes as attachments to glycoproteins and glycolipids, but also as free molecules. In this thesis the effect of carbohydrate molecules on phospholipid model membranes have been investigated by the means of Molecular Dynamics (MD) computer simulations. The most abundant glycolipid in nature is the non-bilayer forming monogalactosyldiacylglycerol (MGDG). It is known to be important for the membrane stacking typical for the thylakoid membranes in plants, and has also been found essential for processes related to photosynthesis. In Paper I, MD simulations were used to characterize structural and dynamical changes in a lipid bilayer when MGDG is present. The simulations were validated by direct comparisons between dipolar couplings calculated from the MD trajectories, and those determined from NMR experiments on similar systems. We could show that most structural changes of the bilayer were a consequence of lipid packing and the molecular shape of MGDG. In certain plants and organisms, the enrichment of small sugars such as sucrose and trehalose close to the membrane interfaces, are known to be one of the strategies to survive freezing and dehydration. The cryoprotecting abilities of these sugar molecules are long known, but the mechanisms at the molecular level are still debated. In Papers II–IV, the interactions of trehalose with a lipid bilayer were investigated. Calculations of structural and dynamical properties, together with free energy calculations, were used to characterize the effect of trehalose on bilayer properties. We could show that the binding of trehalose to the lipid bilayer follows a simple two state binding model, in agreement with recent experimental investigations, and confirm some of the proposed hypotheses for membrane–sugar interactions. The simulations were validated by dipolar couplings from our NMR investigations of TRH in a dilute liquid crystal (bicelles). Furthermore, the assumption about molecular structure being equal in the ordered and isotropic phases was tested and verified. This assumption is central for the interpretation of experimentally determined dipolar couplings in weakly ordered systems. In addition, a coarse grain model was used to tackle some of the problems with slow dynamics that were encountered for trehalose in interaction with the bilayer. It was found that further developments of the interaction models are needed to properly describe the membrane–sugar interactions. Lastly, from investigations of trehalose curvature sensing, we concluded that it preferably interacts in bilayer regions with high negative curvature.

iii

Till Jannie och Margit

List of Papers

The following papers, referred to in the text by their Roman numerals, are included in this thesis. I. Molecular Dynamics Simulations of Membranes Composed of Glycolipids and Phospholipids Jon Kapla, Baltzar Stevensson, Martin Dahlberg and Arnold Maliniak, J. Phys. Chem. B, 116, 244-252 (2012). Specific contributions: Performed the simulations and all analysis. Produced all figures and contributed to the writing. II. Molecular Dynamics Simulations of Membrane–Sugar Interactions Jon Kapla, Jacob Wohlert, Baltzar Stevensson, Olof Engström, Göran Widmalm and Arnold Maliniak, J. Phys. Chem. B, 117, 6667-6673 (2013). Specific contributions: Planned and performed the simulations and most analysis. Produced all figures and contributed to the writing. III. Molecular dynamics simulations and NMR spectroscopy studies of trehalose–lipid bilayer systems Jon Kapla, Olof Engström, Baltzar Stevensson, Jacob Wohlert, Göran Widmalm and Arnold Maliniak, Phys. Chem. Chem. Phys., 17, 2243822447 (2015). Specific contributions: Planned and performed simulations and most analysis. Produced most figures and contributed to the writing. IV. Coarse Grained Molecular Dynamics Simulations of Membrane– Sugar Interactions Jon Kapla, Baltzar Stevensson and Arnold Maliniak, Submitted for publication. Specific contributions: Planned and performed the simulations and all analysis. Produced all figures and had a lead role in the writing.

Reprints were made with permission from the publishers. vii

Contents

Abstract

iii

List of Papers

vii

Abbreviations

xi

1

2

Introduction 1.1 Carbohydrates . . . . . . . . . . . . . . . . . . . . 1.2 Lipids . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Phospholipids . . . . . . . . . . . . . . . . 1.2.2 Glycolipids . . . . . . . . . . . . . . . . . 1.3 Lipid bilayers as model membranes . . . . . . . . 1.3.1 Membrane models in computer simulations

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

13 15 16 20 21 22 24

Computer simulations of lipid bilayers and carbohydrates 2.1 Computational chemistry . . . . . . . . . . . . . . . . 2.2 Molecular Dynamics simulations . . . . . . . . . . . . 2.3 Force fields for lipids and carbohydrates . . . . . . . . 2.3.1 Lipid force fields . . . . . . . . . . . . . . . . 2.3.2 Carbohydrate force fields . . . . . . . . . . . . 2.3.3 Choosing force field . . . . . . . . . . . . . . 2.3.4 Level of detail . . . . . . . . . . . . . . . . . 2.4 Software and simulation parameters . . . . . . . . . . 2.4.1 Periodic boundary conditions . . . . . . . . . 2.4.2 The isothermal–isobaric (NPT) ensemble . . . 2.5 Analysis of simulation results . . . . . . . . . . . . . . 2.5.1 Area per lipid . . . . . . . . . . . . . . . . . . 2.5.2 Orientational order . . . . . . . . . . . . . . . 2.5.3 Mechanical properties . . . . . . . . . . . . . 2.5.4 Distributions . . . . . . . . . . . . . . . . . . 2.5.5 Bilayer thickness . . . . . . . . . . . . . . . . 2.5.6 Diffusion . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

27 27 30 31 31 32 32 33 33 34 34 36 37 37 39 41 42 42

ix

. . . . . .

2.6

3

4

5

2.5.7 Time correlation functions . Potential of mean force calculations 2.6.1 Umbrella sampling . . . . . 2.6.2 Metadynamics . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

Membranes composed of glycolipids and phospholipids (I) 3.1 Force field construction . . . . . . . . . . . . . . . . . . 3.2 Simulated systems . . . . . . . . . . . . . . . . . . . . 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Structural properties . . . . . . . . . . . . . . . 3.3.2 Orientational order and NMR dipolar couplings . 3.3.3 Lateral diffusion of lipids . . . . . . . . . . . . .

. . . .

. . . .

. . . .

43 43 44 44

. . . . . .

. . . . . .

. . . . . .

47 48 50 51 51 53 55

Investigations of membrane–sugar interactions (II–IV) 4.1 Modelling of the trehalose interactions (Paper II) . . . . . 4.1.1 Choice of force field . . . . . . . . . . . . . . . . 4.1.2 Simulated systems . . . . . . . . . . . . . . . . . 4.1.3 Results . . . . . . . . . . . . . . . . . . . . . . . 4.2 Validating the results (Paper III) . . . . . . . . . . . . . . 4.2.1 DMPC dipolar couplings . . . . . . . . . . . . . . 4.2.2 TRH dipolar couplings . . . . . . . . . . . . . . . 4.3 Reaching longer timescales and larger systems (Paper IV) . 4.3.1 TRH interactions . . . . . . . . . . . . . . . . . . 4.3.2 Structural properties of the bilayer . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

57 59 59 60 61 64 64 65 68 69 71

Conclusions 5.1 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75 77

Sammanfattning

lxxix

Acknowledgements

lxxxi

References

lxxxiii

Abbreviations

AFM

Atomic Force Microscopy

CMC

Critical Micelle Concentration

CV

Collective Variable

DFT

Density Functional Theory

DGDG

Digalactosyldiacylglycerol

DHPC

1,2-dihexanoyl-sn-glycero-3-phosphocholine

DLPC

1,2-dilauroyl-sn-glycero-3-phosphocholine

DMPC

1,2-dimyristoyl-sn-glycero-3-phosphocholine

DPPC

1,2-dipalmitoyl-sn-glycero-3-phosphocholine

FRAP

Fluorescence Recovery After Photobleaching

GAL

β -D-galactose

GLU

Glucose

GUV

Giant Unilamellar Vesicle

LUV

Large Unilamellar Vesicle

MC

Metropolis Monte Carlo

MD

Molecular Dynamics

MGDG

Monogalactosyldiacylglycerol

MLV

Multi Lamellar Vesicle

MM

Molecular Mechanics

MSD

Mean Square Displacement

MTD

Metadynamics

NMR

Nuclear Magnetic Resonance

NOE

Nuclear Overhauser Effect

xi

NPT

Isothermal–isobaric (Gibbs) ensemble

NVE

Microcanonical ensemble

PW

Polarized water model in the Martini force field

QM

Quantum Mechanics

RDF

Radial Distribution Function

SANS

Small Angle Neutron Scattering

SPC

Simple Point Charge water model

SPT

Single Particle Tracking

SUC

Sucrose

SUV

Small Unilamellar Vesicle

TCF

Time Correlation Function

TRH

α,α-trehalose

US

Umbrella Sampling

W

Standard water model in the Martini force field

1. Introduction

From the complex membranes and organelles all the way down to the single molecules, atoms and electrons, the building blocks of the cell are essential for all living organisms. The complex structures and processes in the cells make room for many fields of research seeking to unravel the underlying mechanisms. All disciplines involved in trying to understand how it all ties together are collectively adding to the big picture, independently of the length and timescales of the performed research. The biological top down view of organisms and their function, the medical and pharmaceutical approach to cells and the physicochemical investigations of molecules, atoms and electrons, are all needed to get the pieces of the puzzle to fit. The interconnected processes in the cells of living organisms are very sensitive to defects. Many deceases are effects of faults in the cell functions, and to find effective treatments it is important to understand the whole chain of events. Altzheimer’s decease is an important example. The affected brain suffers a massive loss of neurons and disruptions in synaptic signals resulting in damaged memory functions and eventually death. The cause of the decease is not fully understood, but it is commonly accepted that the effects are related to an overproduction of the amyloid β -peptide. This peptide strongly interacts with copper ions and the resulting complexes are thought to be a cause of an increased oxidative stress on the brain cells. Thus, understanding the structural properties of amyloid β and the redox reactions related to its complexes is crucial for finding treatments for the decease [1; 2]. Similarly, research of the mechanisms of de- and rehydration of plant cells has proven useful for developing new methods for storing biological materials [3], and molecular knowledge about the processes involved in photosynthesis is used to develop artificial ways to turn captured CO2 into useful biosynthesis precursors [4], or to split water into hydrogen and oxygen for fuel purposes [5]. Most of the processes related to the living cells are in one way or another depending on the functions of the highly complex membranes that protect the interior of the cells from the world outside. Deep knowledge about the functions and structure of the cell membranes unlocks pathways to other functions of the cell. It is there molecules are transported into and from the cells, and it is also there the cells can initiate communication by the release and reception of signal substances. Lipids and proteins constitute the largest part of the cell 13

membrane structure and function. However, many processes in the membranes are also highly dependent on the presence of various carbohydrate molecules that can be found in and around the membrane structures. In my work, I have used molecular dynamics (MD) computer simulations to investigate carbohydrate molecules involved in cell membranes and their surroundings. Computer simulations of model membranes consisting of mixtures of phospho- and glycolipids (Paper I) give access to the structure and dynamics of the membrane at the molecular level. Together with direct comparisons with data from nuclear magnetic resonance (NMR) experiments of similar systems it is possible to investigate how the orientational order as well as molecular conformations and dynamics change with the lipid composition. Having access to good experimental datasets is crucial for interpreting results from computer simulations. The close collaboration with a research group at the Department of Organic Chemistry has resulted in several successful structural NMR investigations of carbohydrates in dilute liquid crystalline media (lipid bicelles) [6–10]. In the NMR measurements, lipid bilayer structures were used for alignment of the solutes in the magnetic field, but it is generally not known to what extent the interactions between lipids and solutes affect the outcome of the measurements. To address this uncertainty, we characterize the interactions between the disaccharide trehalose and a phospholipid bilayer in Paper II and III. From simulations at different sugar concentrations, we can interpret how the bilayer properties change with the amount of sugar, and by calculating the free energy of the interaction a simple model for the trehalose–bilayer interaction is suggested. The results are validated by direct comparison with dipolar couplings retrieved from NMR experiments. The last paper (Paper IV) is a continuation of the work with trehalose where we try to reach longer timescales and larger system sizes in our simulations by using a well known coarse grained force field to describe the molecular interactions with less details. We do this as an attempt to overcome some of the statistical problems that we encountered in the fully atomistic simulations. The disadvantage of a coarser model is that the detailed molecular information is sacrificed for the gain in simulation time and dynamics. The following sections of this introduction will give an overview of carbohydrates and lipid molecules in general, and the molecules I have been focusing this work on in particular. In Chapter 2, the key simulation and analysis methods used in my work are presented. In Chapters 3 and 4, important aspects of the published results (Papers I-IV) are summarized and discussed, followed by general conclusions and future aspects of the work in Chapter 5.

14

1.1

Carbohydrates

In their simplest forms, carbohydrates are used as energy storage and signal substances in living organisms and they can protect plants and animals from freezing and dehydration. Carbohydrate molecules are found in and around the cell membranes as free molecules, but also incorporated in the membrane structure in the form of glycolipids and glycoproteins. In addition, small carbohydrate molecules can be connected together to form complex polysaccharides such as cellulose, which is the major structural component of plant cell walls. Cellulose is also the main component of products made of wood, paper and fabric. Another important polysaccharide in biology, starch, is used by plants for storage of energy, similarly to the way animals turn excess carbohydrates into glycogen for storage in the liver. It is also used in industry in a large variety of products including food and pharmaceuticals. Carbohydrates have the general empirical formula Cn(H2O)n. The simplest carbohydrates are polyhydroxy aldehydes or ketones, but the general classification also includes derivatives of these with similar properties. The smallest carbohydrates are the three carbon long L- or D-glyceraldehyde, and dihydroxyacetone. These are called trioses, and with an increment of the number of carbons, the successive carbohydates are named tetroses, pentoses, hexoses and so on. In this thesis I will focus on the group of carbohydrates commonly called aldoses which originates from the aldehydes. Carbohydrates within this group are also called sugars or saccharides [11]. The L- and D-forms of the aldehydes (represented as L- and D-glucose in Fig. 1.1a and 1.1b) relate to the stereochemistry of the molecules and how they rotate plane polarized light. Simply put, the two forms, often referred to as enantiomers, are mirror images of each other with the same chemical composition. Interestingly, only the D-forms are abundant in nature. The variation within the group of aldoses does not solely come from the enantiomers and the possible stereochemical combinations of the hydroxyl groups. Many important carbohydrate molecules are aldehyde derivatives, and/or combinations of carbohydrate molecules. Sugar alcohols such as glycerol and sorbitol are produced by reduction of the carbonyl group in D-glyceraldehyde and D-glucose respectively. Amino sugars can be obtained by substituting a hydroxyl group with an amino group, and many carbohydrates are found in nature as phosphate esters. Altogether, a vast amount of aldose derivatives are possible to synthesize, but only a limited amount of these are actually found in nature [11]. Since the aldose carbohydrates contain both alcohol and aldehyde groups, the aldopentoses and aldohexoses have the ability to form intramolecular hemiacetals. The resulting ring forms are called furanoses and pyranoses, and most pentoses and hexoses take this form in aqueous solution. As an ex15

1 CHO

CHO

a HO

C

H

H

C

OH

H HO

HO

C

H

H

HO

C

H

H

ʟ-glucose

HO

OH

C

H

C

OH

C

OH

ᴅ-glucose

OH

6

C

6 CH2OH

CH2OH

c

4

b

d

OH

6

O HO

4 HO HO α-ᴅ-glucose

1

OH

O

4 HO

OH HO 1

β-ᴅ-glucose

Figure 1.1: Aldehydes exhibit different stereochemical isomers, here represented by different forms of glucose. The L- and D-enantiomers of glucose are shown in a) and b) respectively. In solution, the aldoses transform into rings. In c) and d) the D-glucose anomers α-D-glucose and β -D-glucose are shown.

ample, D-glucose (Fig 1.1b) forms α-D-glucopyranose (Fig 1.1c) and β -Dglucopyranose (Fig 1.1d), where the α and β forms differ in the orientation of the hydroxyl group on the C1 carbon. This kind of isomers are called anomers, and the anomer with the hemiacetal hydroxyl group on the same side as the exocyclic CH2OH is designated β . Similarly, if the hydroxyl group is on the opposite side it is designated α. The furanose and pyranose rings can further react with other rings (or functional groups) with their alcohol and hemiacetal hydroxyl groups to form di- tri- and oligosaccharides (more than three) [11]. Two carbohydrates are particularly important for this work: i) β -D-galactose (GAL, Fig 1.2a), a diastereomer to D-glucose which is commonly found in the headgroups of plant glycolipids, and ii) α,α-trehalose (TRH, Fig 1.2b) which is a disaccharide of D-glucose in which the hemiacetals of two α-Dglucopyranose rings have reacted to form a (1←→1) glycosidic linkage. It is known for its ability to protect plants from freezing and dehydration.

1.2

Lipids

Lipid molecules are important parts of living organisms. Oils and fats are used to store excess energy in plants and animals. Phospholipids, glycolipids and 16

a

OH HO OH

4 HO

b 6'

HO

4'

4

OH

6

6

OH

O OH 2

HO

OH

2

OH O

OH

O O

2' HO

Figure 1.2: The two carbohydrate molecules that has been involved in the work leading up to this thesis: a) β -D-galactose (GAL) and b) α,α-trehalose (TRH).

other types of amphiphilic molecules (i.e. molecules that have both hydrophobic and hydrophilic groups) are typically found in the cell membranes where they serve as the structural base [12; 13]. Even though lipids are the main components that make up the cell membrane structure, the function of the cell depends on a vast number of other molecules such as proteins, sterols and carbohydrates that are also embedded in the membrane. This creates a highly heterogenic and complex environment (Fig 1.3) which the experimental and computational techniques available today struggle to fully describe. In the present work we have focused on lipid bilayers consisting of glycerophospholipids and glyceroglycolipids. Therefore the focus of the introduction to lipid molecules will give a brief overview of lipids and their general behavior, with a focus on these two types. For a more in depth view of lipids and their functions in biological membranes, I refer to e.g. Biochemistry of Lipids, Lipoproteins and Membranes [14].

Figure 1.3: A schematic picture of the highly heterogenic cell membrane with the phospholipid bilayer as a structural base and incorporated biomolecules such as proteins and carbohydrates. Image taken from Wikimedia Commons. Released to the public domain by ’LadyofHats’ Mariana Ruiz 2007.

17

Lipid molecules come in many forms and sizes. When you speak about lipids in a scientific context, the first type that comes to mind is probably the phospholipid (or phosphocylglycerol). Due to its abundance in the cell membranes of mammals and bacteria [15] it has become the most studied lipid type, but phatty acids, triacylglycerols (fats), sphingolipids and glycolipids are also parts of this class of molecules. Common to all lipids is that their solubility in water is limited, but they are readily soluble in organic solvents [12]. Most lipids are amphiphilic (Fig 1.4), and this property is the reason why we find these molecules in bilayer formations in the cell membrane. In solution, above the critical micelle concentration (CMC), amphiphilic lipids group together and form micellar or lyotropic liquid crystalline phases (Fig 1.5). The driving force behind this spontaneous phenomenon is called the hydrophobic effect and it is mediated through the change of entropy in the whole system. Since the highly polar water molecules cannot create hydrogen bonds with large parts of the lipids, they are forced into states of lower entropy when they create solvation cages around the lipid molecules. At CMC, water molecules are released from the states of low entropy by the formation of lipid aggregates, in which the lipids face their polar groups towards water. The decrease in entropy of the lipids, and the enthalpy increase due to lipid head group repulsion in the aggregates, result in a positive contribution to the system’s free energy. However, due to the large increase in entropy of the free water molecules, the net outcome of this transition is a total decrease. Typical phases that are formed in lipid aggregation include micelles (L1 ), hexagonal (HI ) and lamellar (Lα ). In some cases, such as low water content, the lipids turn their heads in the opposite direction and form so called reversed phases (type-II [16]) such as reversed micelles or reversed hexagonal phases (HII ) [13; 17].

Non-polar hydrocarbon tails

Polar head group

Figure 1.4: A schematic picture of an amphiphilic lipid molecule. The polar headgroup is hydrophilic and the non-polar hydrocarbon tails are hydrophobic.

The abiltity of the lipids to form the liquid crystalline phases in solution is also affected by the molecular shape. Cylindrical or rod like molecules usually form the lamellar phases, while cone shaped molecules are more likely to take part in higher curvature structures such as the spherical micellar or hexagonal phases. A common measure for the lipid shape is the packing parameter, P = v(al)−1 [18], where v is the volume of the hydrocarbon tails, l is the tail length, and a is the equilibrium area per molecule at the aggregate surface. The unity 18

Lipid concentration

Lamellar

Hexagonal

Micelle

Figure 1.5: Three common phases of polar lipids in water. From the bottom with increasing lipid concentration: Micellar phase, Hexagonal phase and Lamellar (bilayer) phase. Intermediate and reversed phases may exist depending on temperature and the shape of the lipid molecules.

packing parameter, P = 1, indicates a cylindrical molecule. When P < 0.5, the molecule is shaped like a cone or a wedge, and is more likely to be found in the micellar or normal hexagonal phases. If P > 1, the molecule is also cone or wedge shaped, but has the tendency to form reversed phases due to the area in the hydrocarbon region being larger than the lipid head group [17; 19]. One example of a structure where the mixture of lipids with different shape is of great importance is the bicelle, which is a two-component bilayered micelle. The bicelle is a mixture of bilayer forming lipids and shorter lipids that typically forms micelles in solution. The resulting shape of the aggregate can be described as bilayer discs, with the shorter lipids enriched at the edges where the curvature is high (Fig. 1.6). Bicelles tend to align in a magnetic field and they are commonly used as alignment medium in NMR spectroscopy. This particular use of bicelles is described in more detail in Section 2.5.2. The liquid crystalline phases that the lipids can form in solution also exhibit internal phase transitions depending on the temperature and hydration level. For lipid bilayers, a transition from the fluid phase, Lα , to the gel phase, 19

B0

Figure 1.6: A schematic picture of bicelles aligning in a magnetic field. The cutout illustrates the composition of the bicelle, with longer bilayer forming lipids in the middle and shorter micelle forming lipids at the edges.

Lβ , takes place at the specific gel transition temperature, Tm , below which the lipid tails become highly ordered, often in a tilted orientation. The behavior of the bilayer in the gel phase differs substantially from the behavior in the fluid phase. As an example, the bilayer thickness fitted from X-ray scattering experiments of a fully hydrated phospholipid membrane with the saturated 1,2-dimyristoyl-sn-glycero-3-phosphocholine lipid was found to be more than 10% larger in the gel phase than in the fluid phase [20].

1.2.1

Phospholipids

The most abundant lipid molecules in a living cell membrane are the phospholipids. The typical phospholipid consists of a glycerol (or sphingosine) backbone onto which one phosphate group and two fatty acids are attached. The glycerol together with the phosphate group constitutes the hydrophilic part of the molecule, while the fatty acids, typically 12 to 18 carbons long, are hydrophobic. The acidic nature of the phosphate group makes it likely to react with other molecules such as etanolamine, choline or carbohydrate molecules. The phospholipids are often distinguished by the head group composition, and the length and saturation level (i.e. the number of double bonds) of the acyl chains. All work in this thesis has involved a phosphatidylcholine called 1,2-dimyri20

stoyl-sn-glycero-3-phosphocholine (DMPC, Fig 1.7). This particular lipid has fully saturated acyl chains, which is rather atypical for the lipids found in the cells, where double bonds in at least one of the lipid tails are more common. It is, however, one of the most studied lipids, mainly due to practical reasons since the gel transition temperature for this lipid is low (24 ◦C) and bilayers formed by DMPC remain fluid at ambient conditions. The packing parameter of a DMPC molecule is close to unity, and thus it tends to form flat bilayers in aqueous solution. O

Fatty acid tails

Phosphate O Glycerol Choline O O P O + N – O H O Polar head group O

Figure 1.7: The phospholipid DMPC. The polar, zwitterionic head group consists of the glycerol backbone with a phosphatidyl choline attached. The two fully saturated fatty acid tails are 14 carbons in length.

1.2.2

Glycolipids

In the cell membranes of plants and animals, carbohydrates are often found chemically attached to other molecules to form glycoproteins and glycolipids. The function of these molecules in the membranes vary, and are often related to cell signalling and cell identity. In the human body, glycolipids are found to be parts of the communication between cells, which is often constructed as the interaction between glycolipids or glycoproteins on one cell surface, with specific receptor molecules such as lecitins on the other [21]. Glycolipids are also major components in the thylakoid membrane in higher plants, where they act as a structural base. A glycolipid that we have studied in bilayers, in mixtures together with DMPC, is monogalactosyldiacylglycerol (MGDG, Fig 1.8). This molecule represents over 50% of the dry weight in the thylakoid membrane in plant chloroplasts, which makes it the most abundant polar lipid in nature [22]. The facts that MGDG molecules can be found close to the core of photosystem I [23] and that plants with defective MGDG synthase exhibits reduced production of chlorophyll [24], are both strong indications that this particular glycolipid takes active part in processes related to photosynthesis [25]. Additionally, MGDG has been assigned promising anti-inflammatory and antiproliferative* properties in animal tissue [26]. * Inhibition

of cell growth.

21

The headgroup of MGDG consists of a glycerol backbone that has a β -Dgalactose monosaccharide attached to it through a glycosidic bond. The unsaturated tails are 16-18 carbons long, with a total of 6 double bonds; three in each tail. The small galactose headgroup and the highly unsaturated tails result in a relatively high value of the packing parameter, and MGDG therefore tends to form a reversed hexagonal phase in water. O

Glycerol

HO O

O O Unsaturated tails

H

HO O

OH

OH

Galactose

O

Figure 1.8: The glycolipid MGDG. The headgroup consists of a glycerol backbone onto which a GAL is attached through a glycosidic linkage. The two unsaturated hydrocarbon tails are 16-18 carbons long with three cis-doublebonds each.

1.3

Lipid bilayers as model membranes

The diversity of molecules and the number of processes going on simultaneously in the cell and its surroundings, make direct studies of cell membranes difficult with available experimental and computer simulation techniques. Therefore the need for simpler systems, representative of the cell membrane, has emerged to assist the study of single proteins, lipids and mixtures thereof. Today there exist three main membrane models that are used throughout the experimental biosciences: Lipid monolayers, lipid vesicles and supported bilayers. In computer simulations, the sizes of the model membranes used in experiments are usually too large to be studied directly due to computer resource limitations, and as a result the most common models used in simulations are limited to smaller bilayer patches. The simplest of the three experimental membrane models, the lipid monolayer (Fig. 1.9a), is also the oldest. It first caught scientific attention in the late 18th century, when scientists noticed that the spread of oil films on water led to a reduction of ripples on the water surface. One of the most famous experiments from this time is reported in a series of communications between Franklin and coworkers in 1774 [27]. Inspired by tales from fishermen, Franklin observed that a teaspoon of oil poured onto the surface of a pond, spread out to an area corresponding to several square yards with a surface "as smooth as a looking glass". These observations, probably the oldest reported on oil films on water, eventually led to a new era of surface chem22

a

b

Figure 1.9: Schematic pictures of two model membranes used in experiments: a) a lipid monolayer at the air/water interface, and b) a lipid vesicle.

istry, beginning in the late 19th century with the work of Pockels, who used a surface film balance technique to observe impurities in the surface tension of water [28]. Langmuir further developed the surface balance and performed the first systematic studies of monolayers on water. His work led to a more scientific understanding of the films at the molecular level, and one of the most important findings was the preferential orientation of the fatty acid molecules in the films. He was awarded the Nobel prize in 1932 for his work on surface chemistry. Subsequently, together with his coworkers Blodgett and Schaefer, he showed that it was possible to transfer successive monolayer films of any thickness to a solid substrate, so called Langmuir–Blodgett films [29]. The use of Langmuir monolayers and Langmuir–Blodgett films are today standard procedures to characterize lipid molecules. Typically, changes of surface pressure and surface pressure area isotherms between lipid systems with different compositions are monitored with a Langmuir trough to gather information about lipid packing and lateral phase separations. It is possible to vary parameters such as temperature and lipid composition to mimic the cell environment, which makes these methods suitable to investigate the behavior of other molecules in membrane like environments [30–32]. Analysis techniques, such as Brewster microscopy, can be utilized directly on the monolayer/water interface for analysis of the lipid topology [33]. Other techniques, especially Atomic Force Microscopy (AFM) requires a transfer of the monolayer to substrates using the Langmuir–Blodgett approach [34; 35]. Vesicles, or liposomes (Fig. 1.9b), are spherical assemblies of lipids containing one or more bilayers that enclose a small volume of water. One large advantage of using vesicles as model membranes is that in contrast to monolayers they consist of lipid bilayers, which can closer mimic the complex environments of cell membranes. Vesicles can be prepared as multilamellar (MLV) with sizes 0.5-10 µm, small unilamellar (SUV) with sizes below 50 nm, large unilamellar (LUV) with sizes 100-500 nm and giant unilamellar (GUV) with sizes 5-100 µm [32]. Analysis of the vesicles can be obtained by several methods. The lamel23

larity may be detected by 31 P NMR targeting the phospholipid headgroups, or X-ray scattering methods. Size distributions can be obtained by light scattering methods, microscopy, tunneling electron microscopy and so on. To investigate the dynamics and properties of the lipid bilayer and incorporated molecules, it is common to use spectrophotometry or fluorescent spectroscopy techniques as well as dynamic light scattering and NMR techniques [36]. Lately, many studies of vesicles and liposomes utilize lasers to pull out narrow membrane nanotubes from the surface of GUVs. Investigations of tubes with differing diameters may yield information about mechanical and dynamical properties of the lipid membrane [37; 38]. To be able to use certain analysis techniques, or to accomplish stable membrane samples it is common to use so called supported bilayers. Langmuir monolayers, vesicles or micelles are transferred, in a Langmuir-Blodgett like fashion, to a solid substrate with a soft surface (typically mica), on which the resulting supported bilayer is flat. The advantages of using supported bilayers is that it enables some control over the distribution of lipids in the two layers [39; 40], and possiblities to use techniques such as AFM [41; 42] to investigate the properties of the bilayer. Apart from the three most commonly used model membrane systems mentioned above, there also exist other types of lipid assemblies that can be more suitable for certain experimental techniques. Examples of such systems are micelles and bicelles that are found particularly suitable for investigations of bilayer–peptide interactions in solution NMR [43]. In addition, similar studies of bicelles with intermixed glycolipids, show that these assemblies are suitable models for NMR studies of the structure and dynamics of lipid bilayer systems, and a possible route for investigations of glycolipid–protein interactions [44].

1.3.1

Membrane models in computer simulations

In this thesis, computer simulations are used to investigate properties and dynamics of lipid bilayers and molecules in their proximity. There are important differences between the experimental models and the models typically used in simulations. First and foremost, the size of computer generated membranes are substantially smaller than the experimental counterparts. The limiting factor is mainly the available computer power, and even though we have seen an extensive growth in computer development the last 30 years, the number of lipid molecules in typical simulated membranes are still well below 5 000 lipids for classical simulations of atomistic resolution, and much smaller for quantum chemical calculations where the electrons have to be taken into account. The bilayer systems in this thesis were simulated on supercomputers with thousands of processor cores. A typical simulation that used 192 cores to simulate 24

128 lipids and 10 000 water molecules produced around 50 ns simulation data per day. As a comparison, a unilamellar liposome or vesicle with a diameter of 100 nm contains around 100 000 lipid molecules† . Furthermore, the limits of system sizes also limits the shape range for simulated membranes. Most classical membrane simulations the last 15 years have used rectangular membrane patches of sizes up to a few hundred lipids per leaflet (Fig. 1.10). There are, of course, cases where state of the art computers have been used to simulate rectangular membranes of hundreds of thousands of particles [45], but for most applications the size limits, and the effects of them, is something that have to be considered when designing and analyzing simulations.

Figure 1.10: Typical bilayer patch with 128 DMPC molecules (64 in each leaflet) used in simulations. The two images show the patch from different angles. Created with VMD [46].

It is also important to recognize differences in timescales within which the experimental techniques and simulations operate. A good example is NMR, where measured quantities such as dipolar interactions are measured in Hz, and thus averaged over several seconds. Typical atomistic simulations lie within the nanosecond timescale (the simulations in this thesis typically spans 100500 ns), and slower processes may drastically affect the accuracy of calculated averages.

† If

we assume a lipid headgroup area of 0.6 nm2 and a bilayer thickness of 4 nm.

25

26

2. Computer simulations of lipid bilayers and carbohydrates

The modern era of computer simulations began around the end of World War II. The development of computers for military applications for code breaking and weapon development during the 1940’s paved the way for scientists such as Metropolis who 1953, together with his coworkers, published the well known algorithms for importance sampling using the Monte Carlo methods for statistical evaluation of mathematical functions [47]. In 1956, the first Molecular Dynamics simulations of hard spheres were performed by Alder and Wainwright [48], closely followed by simulations of real particles and fluids [49; 50]. Moore’s famous prediction in 1965 that the number of components in the computer microchips will double every 12 months [51], with the modification to every two years in 1975 [52], has so far proved highly accurate. The fast development of powerful computers has enabled extensive use of computer models to complement various experimental techniques, not least in biochemistry and biophysics. The complexity of systems such as the cell membrane, together with the fact that no single method is perfect for all applications, has forced the community to develop combinations of different methods at different length and time scales. These so called multiscale approaches typically combine quantum chemistry based calculations (ab initio) at the atomic level, classical dynamic or statistical methods (such as Molecular Dynamics and Monte Carlo) at the molecular level with coarser models approaching the macroscopic (Fig. 2.1). A milestone for computational chemistry was reached 2013, when Martin Karplus, Michael Levitt and Arieh Warshel were awarded the Nobel Prize in Chemistry for their work on multiscale modeling.

2.1

Computational chemistry

The computational methods used for simulations of chemical systems, can generally be divided into three groups: Quantum Mechanical Methods (QM), Metropolis Monte Carlo Methods (MC) and Molecular Dynamics Methods (MD). The last two are based on Molecular Mechanics (MM), and the atomic 27

Figure 2.1: The time and lengthscales in Multiscale Modelling.

interactions are described by sets of classical interaction potentials often renowned as force fields. Quantum mechanical approaches, also referred to as ab initio (from first principles), are suitable to investigate systems and processes where the electronic behavior of the atoms and molecules is of major importance for the accuracy of the calculations. This can include the calculations of chemical bonds, reactions and partial charges, as well as the optimization of chemical structures to their lowest energy configuration. These approaches solve the electronic Schrödinger equation numerically, and utilize the Born–Oppenheimer approximation assuming that the atomic nuclei are stationary in comparison to the electrons. The electronic orbitals are estimated by sets of basis functions which are optimized for lowest energy, in agreement with the variational principle of quantum mechanics. Commonly used methods are the Hartree–Fock based methods where molecular orbitals are approximated as linear combinations of atomic orbitals, and Density Functional Theory (DFT) methods where, instead of using electronic orbitals, the electrons are treated in a mean field approach based on the Hohnenberg–Kohn theorems relating the electron density to the wave function. The latter method is less computationally expensive but still highly accurate, and this is a reason for its popularity. A rule of thumb regarding the choice of simulation technique for a particular chemical system is that the more precise a technique is, the more computer 28

power is needed to perform the calculations. The Hartree–Fock based methods, as an example, have to take all electronic wavefunctions into account, and are therefore computationally expensive (standard Hartree–Fock scales as O(N 4 ), where N is the number of basis functions [53]). In practice it is usually not feasible to use the QM methods on systems above a few thousand atoms in size, and they are often used on small parts of larger systems such as binding sites of large proteins [54], or to calculate partial charges in lipids during force field development [55]. In MC and MD, electrons are usually* not taken into account, and bonds are often constrained to reduce the degrees of freedom. The atoms in these models behave classically, and interact with harmonic, coulombic or Lennard–Jones type pair potentials. The performance is therefore much higher, and allows for large systems and extensive sampling. Metropolis Monte Carlo, which most often is what people refer to with just Monte Carlo, is a statistical method that is used to randomly sample the potential energy landscape of a system according to the Boltzmann distribution of states [47]. If we imagine a system of a few particles, each particle is moved stepwise in an iterative fashion. The move is accepted if the potential energy of the system is decreased, otherwise the move is evaluated by drawing a random number between 0 and 1, comparing it to the Boltzmann factor, e−∆U/kB T , and accepting the move if the random number is smaller or equal. The procedure is repeated for all particles iteratively until the system reaches equilibrium. At equilibrium, properties of interest can be calculated as the ensemble averages. The MC method may seem like a simple and robust choice for most systems, but it is in fact rather unpractical when it comes to large condensed systems, e.g. lipid bilayers where molecular dynamics such as rotational and translational diffusion are important for the interpretation of the simulations. Due to its purely statistical approach, time is not present in MC simulations and dynamical properties cannot be obtained, even if there exist algorithms that take dynamics or kinetics into account in the form of known transfer rates [56]. A more straightforward method to achieve simulations with an available time axis, is to use MD simulations. This is what is most commonly used for large biochemical simulations, with no exception for this thesis, and it allows for systems up to millions of particles in size with modern scientific computer clusters and highly optimized software. The following sections of this chapter will introduce MD as a tool for simulations of lipid bilayer systems and how structural and dynamical properties of the bilayer and its associated molecules can be calculated from the resulting time trajectories of atomic positions and validated by existing experimental data.

* Except

in hybrid QM/MM methods.

29

2.2

Molecular Dynamics simulations

Molecular Dynamics Simulations are methods in which the classic Newtonian equations of motion are solved numerically to move the atoms through space and time. For detailed descriptions of the procedures involved, I recommend the book Understanding Molecular Simulations by Daan Frenkel and Berend Smit [57]. Here, only the most important factors for successful simulation of bilayer systems in combination with carbohydrates will be addressed. The main idea in MD is to let the atoms interact through sets of predefined potential functions, force fields, that are parametrized to certain properties based on experimental data or QM calculations. The forces acting on the atoms at any given time can then be calculated and used to solve the equations of motion utilizing computational algorithms referred to as integrators, of which the most well known are the Velocity-Verlet [58] and Leap-Frog [59]. Most force fields are additive, meaning that the sum of separate contributions from intramolecular (bonds, angles, torsions) and intermolecular (van der Waals interactions, electrostatics) interactions constitutes the expression for the potential energy of the whole system (eq 2.1). Utot = Ubond +Uangle +Utorsion +UL−J +Uelec =

∑ kb (b − b0 )2 + ∑ kθ (θ − θ0 )2 + ∑ 

+ ∑ 4εi j L−J

σi j ri j

12

 −

kφ (1 + cos(nφ − δ )

torsion

angle

bond

σi j ri j

6 !

(2.1)

qi q j elec 4πε0 ri j

+∑

where equilibrium parameters related to bond lengths (b0 ), angles (θ0 ), torsion angles (multiplicity n and phase shift δ ) are provided by the force field for any set of atom types, together with the corresponding force constants for the harmonic potentials (kb , kθ , kφ ). Non-bonded Lennard–Jones and electrostatic parameters are usually given per atom type as the self interaction strength (εii ), atom size (σii ) and charges (qi ). Cross interactions are treated with combination rules, of which the most common are the Lorenz–Berthelot rules: σii + σ j j 2 √ εi j = εii ε j j

σi j =

(2.2) (2.3)

Equation 2.1 is representative of many existing force fields, even though additional terms to stabilize molecular structure (often referred to as improper torsions) or more advanced non-bonded terms are common. The intramolecular interactions are described by harmonic potentials, which forces the user 30

to decide the molecular topology, such as chemical bonds and protonation state of acids, in advance. The development of accurate force fields for biochemical simulations is tedious work, and most force field development takes place within the framework of existing well known force fields, each that has its own methodology for parameter fitting and charge calculations. It is also worth noting that every biomolecular force field has been optimized for a certain water model. There exist many different parameter sets that describe water molecules [60], and the most commonly used, TiP3P [61], SPC [62] and SPC/E [63] are based on point charges. As a user, it is important to recognize which water model is suitable for a certain force field implementation, and keep the properties of the specific water model in mind when analyzing the results. To assure proper sampling and reliable ensemble averages, the system simulated with MD needs to be ergodic. The meaning of ergodicity is that all thermally available points in phase space (all possible configurations) have to be accessed. The ergodicity principle also states that, if a system is ergodic, the ensemble averages can be exchanged for time averages as the time approaches infinity. According to the Boltzmann distribution of states, the probability of accessing high energy states is low and the frequency by which the system accesses such states is subsequently low. This imposes a problem in simulations where states with high energy barriers are present, since there is no guarantee that the system visits these states within the timescale of the simulation. If the ergodicity criterion is not met, the resulting time averages are not representative of the expected thermodynamic ensemble.

2.3

Force fields for lipids and carbohydrates

The choice of force field is the most important factor that influences the results from MD simulations. There is no single force field that is general enough to describe all systems equally well, and most force fields are specialized within a certain area of research. One thing to note is that, as the methodology in parameter fitting and validation differs between force fields, mixing of parameters from different force fields is generally not recommended unless it can be assured that the levels of interaction between particles in the mix are preserved.

2.3.1

Lipid force fields

For lipids there exist a number of force fields at atomic resolution (atomistic) to choose from. Historically, lipid parameters have mainly been developed to support simulations of proteins in membrane environments, but with an increasing knowledge about the role of lipid molecules in cell membranes, more 31

and more effort has been invested in producing more accurate lipid parameters. Early attempts to simulate lipid bilayers often included a fixed area or surface tension, and there was an ongoing debate in the late 90’s about whether or not an applied surface tension is correct or not for simulations of small bilayer patches [64–66]. In 1997, Berger et al. [67] reported their parameters for lipids, compatible with the Gromos force field [68], which showed results in excellent agreement with experimentally determined properties such as the lateral bilayer area and acyl chain order parameters without using any preset nonzero surface tension. Today, the general view is that a flat bilayer exhibits zero or a very small surface tension, and that bilayer simulations should be performed accordingly. Most modern parameter sets are now stable at zero surface tension, and there exist several accurate sets of lipid parameters compatible with large well known force fields such as Amber [69] including Stockholm Lipids (SLipids) [70], OPLS-AA [71], Charmm [72–74] and Gromos [75–77].

2.3.2

Carbohydrate force fields

Carbohydrate force fields are difficult to parametrize correctly. Most carbohydrate molecules exhibit a broad conformational range that is hard to generalize into a transferable parameter set of interaction potentials. Most existing carbohydrate force fields are targeting molecules in aqueous solution and for simple monosaccharides the largest challenge is to balance the intramolecular interactions so that the distribution of anomeric ring conformations is correct. With an increasing number of interconnected carbohydrates (oligosaccharides, polysaccharides) and various exocyclic attachments, the problems related to correct conformational distributions increase. The first attempts to develop general carbohydrate force fields treated only small saccharides, but today there exist a number of parameter sets, namely within Glycam [78], Charmm [79–82], OPLS-AA [83; 84] and Gromos [85; 86], that have broader ranges of carbohydrate molecules to choose from [87].

2.3.3

Choosing force field

When first choosing a force field for a particular system of interest, it is reasonable to begin with something that has already been tested for similar systems. Since mixing of parameters from different development methodologies is discouraged, different molecules within the system should have parameters from the same source. As an example, if a certain carbohydrate molecule only exists within Glycam, any other molecules in the system should also be parametrized within Glycam, or in any force field that is compatible with the Glycam methodology such as Amber. If a molecule is not available in a certain 32

force field implementation, it should be parametrized using the proper methodology. Only if it can be properly motivated, a mix of force field parameters is eligible.

2.3.4

Level of detail

Another consideration when choosing a force field, is the level of detail. Basically, there exist three levels to choose from when it comes to classical MD; atomistic (or all atom), united atom and coarse grained. In an atomistic force field, all atoms are represented explicitly. To speed up the sampling and simulation time, a common modification to the atomistic case is to unite non-polar hydrogens with the corresponding carbon atoms. The resulting united atom force field produces a smoother potential landscape that allows for faster reorientation of hydrocarbon chains, thus increasing the sampling rate of possible conformations. The reduction in the number of atoms also reduces the computational demands. One of the most used united atom force fields is Gromos [68; 88]. Generally, the information available from simulations in a united atom lipid force field is the same as from atomistic simulations, except for properties where the explicit hydrogen positions are important (e.g. when calculating NMR order parameters). In such cases, the hydrogen positions have to be calculated from the positions of the neighboring carbon atoms. To increase the sampling rate, coarse graining procedures are often used to reduce the degrees of freedom in the system. Typically, this is accomplished by grouping atoms together into larger spherical interaction sites, or beads. The Martini force field for simulations of biomolecules [89] is an attempt to create a general, transferable set of parameters for coarse grained simulations of proteins, lipids and sugars. Based on the Gromos scheme of interaction potentials, groups of four atoms are coarse grained into single beads. The resulting force field contains a set of predefined bead types corresponding to different functional groups in various biomolecules. Thus, by choosing wisely between the predefined bead types, new molecules can be added to the force field without extra parametrization. Despite the high level of simplification, the Martini force field has been successfully deployed for many different applications [90].

2.4

Software and simulation parameters

To run the simulations, there exist a number of available software packages for biochemical simulations. The most used force fields, e.g. Amber, Charmm and Gromos, come bundled with their own proprietary simulation software, and there also exist some standalone software packages which can utilize several existing force fields, e.g MDynaMix [91], NAMD [92] and GROMACS [93; 94]. 33

All simulations in papers I-IV were performed with Gromacs. The developers of Gromacs claim it is the most optimized and actively developed of all simulation programs for simulations of biochemical molecules such as proteins and lipids. It is a collaboration project that started in Berendsen’s research group at the University of Groningen, Netherlands, and is now directed from research groups at Science for Life Laboratory in Sweden, with contributors from all over the world. The Gromacs code can run simulations in many of the existing force fields, and it is released as open source [95], which make it available for everyone with the possibility to change the code as needed. The force field developers and/or software developers usually provides recommended input parameters for their own or several simulation software, and best practice is to follow these recommendations as close as possible. Electrostatic interactions are most often treated explicitly up to a real space cutoff distance. At longer distances the Particle Mesh Ewald summation (PME) [96] or Poisson–Boltzmann reaction field [97] methods can be employed. The nonbonded Lennard–Jones interactions are calculated explicitly up to a certain distance after which they are sharply cut off or more smoothly shifted to zero, often with an extra van der Waals correction to account for neglected long range interactions. Normally, parameters related to the thermodynamic ensemble and performance tuning are left for the user to decide upon.

2.4.1

Periodic boundary conditions

The limits of system size in biomolecular simulations give rise to edge effects. The molecules at the edge of the simulation box do not exhibit bulk properties, and the number of edge molecules is substantial in a moderately sized simulation box (∼ 104 particles). To enable the calculations of bulk properties, a common method is to use so called periodic boundary conditions. In this approach, the box is replicated infinitely in all directions with the result that the molecules at the edge of the simulation box are interacting with molecules in the neighboring box (Fig. 2.2). In this way the edge effects can be held at a minimum. However, it is important to make sure that the simulation box is large enough, not to allow large molecules to directly interact with their mirror images in the neighboring boxes, when the distance approaches the real space cutoffs. This is normally referred to as the minimum image convention.

2.4.2

The isothermal–isobaric (NPT) ensemble

Integration of the Newtonian equations of motion, using a time independent system Hamiltonian (the expression for the total energy of the system), samples system states in the microcanonical ensemble (NVE, where N, V and E 34

Figure 2.2: A 2-D schematic of periodic boundary conditions. When the molecule with the arrow moves to the left from the box in the center, it comes back in again from the other side. This applies for all box replicas in all dimensions.

denote constant number of particles, volume and energy, respectively) by default. The most useful ensemble for reliable simulations of lipid bilayers is the isothermal–isobaric (Gibbs) ensemble, commonly denoted NPT, where the temperature and pressure are only allowed to fluctuate around predefined average values. Allowing for a fluctuating volume of the simulation box, this ensemble enables simulations of bilayers without a predefined area per lipid. To use this ensemble, the temperature and pressure need to be kept constant (or fluctuating around a constant value) during the course of the simulation. The temperature can be controlled using so called thermostats, that couple the system to an external heat bath. The Berendsen thermostat [98] accomplished a constant temperature by scaling the atomic velocities directly by a weak coupling algorithm, in relation to the reference temperature. Due to neglections of the local stochastic contributions to the kinetic energy, this method suppresses the temperature fluctuations with the result that the use of this thermostat does not sample a system in a proper thermodynamical ensemble [99]. To assure sampling in the proper ensembles, more advanced ways of dealing with the temperature coupling have been developed. Commonly used are the Nosé–Hoover thermostat [100; 101] which acts through a modification of 35

the system Hamiltonian, and the more recent velocity rescaling thermostat developed by Bussi et al. [102]. The latter is a modification of the Berendsen thermostat, adding a stochastic term to the velocity scaling to assure proper temperature fluctuations. The parameters related to constant temperature include the type of thermostat, the expected average temperature and an empirical coupling strength parameter, which normally differs between the implementations and give rise to a certain amount of ambiguity. In addition, parts of the system can be coupled independently. This is often used when the degrees of freedom of parts of the system (proteins, lipid bilayers etc.) differ from e.g. the solvent and it assures a uniform temperature across the whole simulation box. The pressure is controlled in a similar fashion adjusting the system volume with so called barostats. The straightforward way is to use direct scaling of the simulation box with a weak coupling algorithm to approach the preset average pressure. An implementation of this can be found in the Berendsen barostat [98] and since this method, similarly to the corresponding thermostat, suppresses volume fluctuations the system is not sampled in a proper thermodynamic ensemble. Several more advanced barostats exist, of which the most commonly used is the Parrinello–Rahman barostat [103]. This implementation uses a similar approach as the Nosé–Hoover thermostat, resulting in proper fluctuations of the box volume and thus a proper thermodynamic ensemble of states. The settings related to pressure includes the type of barostat (usually Berendsen or Parrinello–Rahman), the expected average pressure, the thermal compressibility of the system and a coupling strength parameter. Additionally, lipid bilayers are usually simulated with a semiisotropic pressure coupling, meaning that the box axes spanning the lateral plane of the bilayer are scaled simultaneously to keep the size ratio of the bilayer constant and prevent elongation of the system in one direction (which may break the minimum image convention in the other direction).

2.5

Analysis of simulation results

The results from an MD simulation are printed to file as time trajectories of atomic positions and velocities together with the potential and kinetic energies of the system. At each time frame the size of the simulation box is also available for analysis. To calculate time averages of wanted properties, care has to be taken to assure that the system is as close to equilibrium as possible. Usually this is accomplished by following the time evolution of a property and observe the convergence. If it converges to a fluctuation about a constant value the analysis is performed from the time of convergence to the end of the trajectory. If the property never reaches convergence, the simulation time is 36

probably too short for accurate evaluation of the trajectories.

2.5.1

Area per lipid

The area per lipid molecule (AL ) can be determined from simulations as well as experiments. It is a quantity that can give useful information about structural changes in bilayers, and it is extensively used to investigate if a lipid bilayer simulation has converged towards equilibrium. Since experimental values for many different lipid types are readily available, these are also used to validate lipid force field parameters during parametrization as well as lipid bilayers in simulation studies. The AL can be determined experimentally by several different methods including X-ray scattering techniques [104; 105], NMR [106– 108] and gravimetric methods [109]. From simulations of one component bilayers, AL can be calculated as the area of the simulation box in the plane of the bilayer, A0 , divided by the number of lipids in a single leaflet, NL /2: AL =

2 hA0 i NL

(2.4)

This is also often referred to as the projected area per lipid. The approach works for smaller bilayer patches of a few hundred lipid molecules but as the system size increases, bilayer undulations become more pronounced and the projected area underestimates the true bilayer area. In these situations the area fluctuations stemming from undulations can be corrected for by using additional terms with a dependence on the bending modulus of the bilayer and the temperature [110]. In mixed bilayers with two or more lipid types or other incorporated molecules such as sterols or proteins, other approaches are needed to determine AL . The most straightforward method is to use a grid in the plane of the bilayer and assign each lipid the nearest grid points. The area per lipid can then be determined from the fraction of grid points of a given molecule type, multiplied by A0 and divided by the number of molecules of the same type. This approach works rather well for smaller bilayer patches, but the number of grid points needed [111] to get a good estimate of AL increases rapidly with system size. A computationally more effective method is to use Voronoi tesselation [112] to assign a polygon to each lipid molecule in the bilayer.

2.5.2

Orientational order

Local orientational order in soft matter in general, and lipid bilayers in particular, can be expressed in terms of the order parameter S, which is the second 37

order Legendre polynomial P2 (cos φ ):  S = hP2 (cos φ )i =

 1 2 (3 cos φ − 1) 2

(2.5)

where the angular brackets denotes the average over all molecular orientations and φ the angle between a vector in the molecule and the director (the average orientational direction of the molecules, which in bilayer simulations corresponds to the bilayer normal). The order parameter measures the degree of order of a molecular vector, and the highest value, 1, corresponds to perfect alignment along the director. A fully disordered system exhibits order parameters close to 0, but it can also take this value in a special case when the vector is tilted at an angle of 54.7°, the so called magic angle. Negative values of the order parameter corresponds to alignment perpendicular to the director, and the minimum value is −0.5. The most frequently used experimental method to determine local orientational order parameters is the deuterium (2 H or D) NMR spectroscopy. In this technique, the quadrupolar interaction for the deuterium spins is monitored, and the order parameters for the CD-bonds, SCD , can be extracted. These order parameters, in particular for lipid tails, are often used in bilayer simulations to validate the interaction parameters. An alternative technique, that provides completely analogous information, are measurements of NMR dipole-dipole interactions (dipolar couplings) for CH-bonds which, in turn, provide SCH order parameters. Note that, in deuterium NMR both experimental procedures and the resulting spectra are relatively simple, whereas determination of dipolar couplings is a technically demanding procedure. On the other hand, isotopic labelling (required for deuterium experiments) can be both complicated and expensive. Simulations of deuterated lipid bilayer systems are rare, which results in the comparison between calculated SCH and experimental SCD , being the most common. The NMR dipolar couplings are mediated through space, and dependent on the anisotropic orientational distributions of the spin-spin vectors. These interactions depend on the inverse cubic distance between the interacting spins and provide important information about the orientational order and molecular structure. The latter through the strong dependence on intramolecular distances. Since the magnitude of dipolar couplings is proportional to the molecular order, all couplings average to zero in isotropic solutions. In principle, there are two techniques to measure dipolar couplings: i) solid state NMR and ii) solution (high resolution, HR) NMR. The former is used for investigations of systems where the orientational order is naturally present: Lipid bilayers possess a preferential orientation of the lipid molecules and therefore the dipolar couplings can readily be measured. For measurements of dipolar couplings in solution (using HR NMR) an orientational order has to be induced by aligning 38

the molecules using a liquid crystalline solvent (or polymer gel). A frequently used alignment medium consists of bicelles (bilayer discs, see Fig. 1.6), which are formed by a mixture of lipids such as DMPC, and a shorter lipid, e.g. 1,2dihexanoyl-sn-glycero-3-phosphocholine (DHPC), that can cover edges of the bilayer. Bicelles create the orientational order necessary for measurement of dipolar couplings, and simplify the spectra by aligning in the presence of the magnetic field of the NMR spectrometer. From simulations the dipolar couplings between atoms i and j can be calculated using µ0 γi γ j h¯ Di j = − 2 3 S (2.6) 8π hri j i where µ0 is the magnetic constant, γi and γ j the magnetogyric ratios of the interacting spins, ri j the distance between them and S is the order parameter.

2.5.3

Mechanical properties

The elastic properties of lipid bilayers are important quantities in the characterization of biological membranes. The response of a membrane to mechanical stress is strongly dependent on the compressibility and bending moduli. Elastic properties, and bending rigidity in particular, can be measured experimentally by several techniques; monitoring the thermal undulation spectrum by light microscopy [113; 114], analysis of the area change of vesicles under micropipette pressurization [115] (also known as micromechanical techniques [116]) or measuring the force needed to pull thin tubes, tethers, out of the membrane [117], to mension a few. From simulations the bending modulus is typically determined utilizing the Helfrich–Canham theory [118; 119], in which the membrane is modeled as a thin structureless sheet. The bilayer shape is monitored through height variations from a reference plane, and is typically analyzed by spectral analysis in Fourier space in terms of the reciprocal space vector [120]. In this thesis a molecular approach [121–123] for calculating the bending rigidity has been used. The model relates the bending modulus, kC , of the membrane to the splay modulus, which is derived from the potential of mean force for the splay deformation   P(α) PMF(α) = −kB T ln (2.7) sin α where P(α) is the the probability density for the angle, α, between vectors defining the long molecular axis of two different lipids. The lipid splay modij ulus, χ12 , corresponds to the free energy cost of splaying (Fig. 2.3) one lipid molecule (i) with respect to another ( j). 39

α

Figure 2.3: The splay of two lipid molecules is defined as the splay angle α between the length axes of the molecules.

According to Helfrich’s continuum theory the free energy of splay is related to the molecular tilt: χ12 PMF(α) ∝ (∇θ )2 (2.8) 2 in the region of small tilt angles, θ , which corresponds to the tilt of lipid molecules with respect to the bilayer normal. The key point behind the model is then that, if the splay distribution is calculated only for nearest neighbours with small tilt angles (θ < 10◦ for at least one of the lipids in a pair), the splay modulus can be determined from the best quadratic fit of the splay distribution in an interval of small splay angles, α. The bending modulus, kC , is calculated from ϕi j 1 1 = (2.9) ∑ ij kC ϕtotal i, j χ12 where ϕi j and ϕtotal is the number and total number of nearest-neighbor pairs, respectively. In the case of a one component bilayer, eq 2.9 simplifies to kC = χ12 . The area compressibility modulus, KA , for a lipid bilayer reflects the resistance of the bilayer to compression in the membrane plane. It can be obtained from simulations by the fluctuations of the area per lipid [124], σA2L KA =

2 hAL i kB T σA2L NL

(2.10)

Since the compressibility modulus is calculated through fluctuations in projected area per lipid, it contains contributions from undulations that typically 40

are corrected for in experimental reports. In analogy with AL , these contributions can be corrected for in the calculations, by additional terms to eq. 2.10 that are dependent on the bending modulus of the bilayer (kC ) and the temperature [110].

2.5.4

Distributions

A convenient method to examine the behavior of a certain property is to calculate the distribution of the property over all molecules and time. The common way to visualize the results is to distribute the calculated observations into a binned histogram (often called distribution) k

n = ∑ mi

(2.11)

i=1

where n is the total number of observations, k the number of bins and mi the number of observations that falls into a specific bin i. The histogram is most often normalized with the total number of observations and is in that case a discrete approximation of the probability density function, P(x). Typical uses for histograms in MD simulations are to monitor angles, torsions, positions and distances between atoms and molecules, and they are powerful tools for visualization of datasets for paired comparisons. Normalized histograms relate to the free energy as PMF(x) = −kB T ln(P(x))

(2.12)

where PMF(x) is the potential of mean force and P(x) is the probability density. Since the PMF is based on the distribution of states that the system has visited during the simulation, it only reflects the accessible regions of the free energy landscape and most likely omits states with low probability. In Section 2.6 two more elaborate ways of PMF calculation are presented, in which the system is forced into high energy states that may not have been accessed during normal simulation due to high energy barriers. A special form of distribution is the Radial Distribution Function (RDF), which can be used to investigate local atomic and molecular structure. It is calculated as the distribution g(r) of molecular pairs at the distance r around a chosen molecule or atom type. For a single molecule or atom, g(r) takes the form: n(r) g(r) = (2.13) ρ4πr2 ∆r where n(r) is the number density in a small slice at distance r, ρ is the mean particle density and ∆r is the thickness of the slices in which the pairs are 41

calculated. Since the function is normalized with the mean particle density so that it converges to 1 at large values of r, it can be interpreted as the local deviation from the uncorrelated bulk. The RDF can be used to directly compare lipid bilayer simulations with X-ray scattering data, since g(r) is related to the static structure factor S(Q) that can be experimentally determined [125; 126].

2.5.5

Bilayer thickness

The thickness of a lipid bilayer is typically measured experimentally by X-ray scattering, from which the electron density profile is retrieved. The thickness of a phospholipid bilayer is measured as the peak to peak distance corresponding to the distance between the phosphate groups [20; 127]. From MD simulations the bilayer thickness can similarly be retrieved either from the calculated electron density across the bilayer, or from the peak to peak distance in the atom distribution along the bilayer normal.

2.5.6

Diffusion

The diffusion of molecules in membrane systems can be measured by several different methods. One extensively used method is fluorescence recovery after photobleaching (FRAP), by which fluorescent labeled molecules in the bilayer are photobleached by a laser beam, and the recovery of fluorescence, depending on the influx of molecules from the surrounding areas, is subsequently monitored [128]. Other methods to determine diffusion in lipid bilayers include single particle tracking (SPT) [129; 130] and NMR techniques [131; 132]. The molecular diffusion in MD simulations can be determined from the mean square displacement (MSD) following the Einstein relation

 [r(0) − r(t)]2 = lim 2dDt t→∞

(2.14)

 where [r(t) − r(0)]2 is the MSD, d is the dimensionality of the system and D is the self diffusion coefficient. Thus, the self diffusion coefficient can be estimated by a linear fit to the MSD calculated from simulations. Since there is no a priori knowledge about sufficient simulation length, care has to be taken to assure that the simulations are long enough. Usually the estimated error of the diffusion coefficient is calculated through block averages of the linear fit, but in a recent publication [133] Pranami and Lamm stress the importance of repeated simulations related to the fact that different starting configurations may give substantially different values of the diffusion coefficient. In bilayer 42

systems, molecules attached to the bilayer (e.g. lipids and proteins) exhibit two dimensional (lateral) diffusion, and d = 2 in eq. 2.14.

2.5.7

Time correlation functions

Calculation of the time correlation function (TCF) related to an average property of interest is a powerful method to monitor convergence of the simulation, and to estimate the precision of the calculated average [110]. For the property hAi, the TCF is defined as: CA (t) = hA(0)A(t)i

(2.15)

where t is the time. The TCF often decays exponentially, in which case the correlation time τ can be obtained through a fit to one or several exponentials (C(t) = exp(t/τ1 )+exp(t/τ2 )+...). The correlation time, τ, gives information about how long time it takes for the property A to be uncorrelated its value at the starting time, and is valuable to estimate the simulation length needed to get statistically well determined time averages. A special case of the TCF is the second order Rotational Correlation Function, which give information of the time correlation of a molecular vector: C2 (t) = hP2 {cos[φ (t)]}P2 {cos[φ (0)]}i

(2.16)

where P2 = 21 (3x2 − 1) is the second order Legendre polynomial. At long times the events become uncorrelated and the function converges to the order parameter squared, S2 S2 = C2 (∞) = hP2 {cos[φ (∞)]}i hP2 {cos[φ (0)]}i

2.6

(2.17)

Potential of mean force calculations

The Potential of Mean Force (PMF) corresponds to the free energy profile along a reaction coordinate, and calculations of PMFs can give important information about chemical processes. Many processes in biology are defined by changes in free energy along a reaction pathway, and can thus be described by a PMF through the same pathway. Areas where the free energy landscape is of great importance include protein–ligand association [134], drug discovery [135], and interactions of various compunds with bilayers and membranes [136; 137]. As previously mentioned (Section 2.5.4), determining the PMF of a process through the probability density distributions (eq. 2.12) can, at best, only be a rough estimation. It is more rule than exception that there exist states in 43

a simulated molecular system that are not accessed during the length of the simulation due to high energy barriers (and thus low probability). To overcome these limitations, a number of methods have emerged to bias the system such that it is forced into states of low probability. Two of these methods have gained a lot of interest lately, and are now considered standard tools when performing molecular simulations: Umbrella sampling (US) and Metadynamics (MTD).

2.6.1

Umbrella sampling

In Umbrella Sampling (US) [138], the system is forced to remain in certain configurations, rN , through a biasing potential acting in the reaction coordinate (that e.g. can be a distance or a torsion angle). The bias is most often expressed as a harmonic potential u(rN ) = ku (rN − rN0 )2

(2.18)

where ku is the force constant, and rN is the reaction coordinate. The modification of the potential energy of the system can be expressed as U ∗ (rN ) = U(rN ) + u(rN )

(2.19)

A US calculation is often based on a set of simulations that represents windows, rN0 , along the reaction coordinate. To ensure sampling along the whole coordinate, the windows are chosen so that the harmonic bias potentials in each window overlap with the nearby windows, and the PMF can be evaluated from the bias potential using the weighted histrogram method [139]. The biased simulations in US are not sampled according to the Boltzmann distribution of possible states but, since the biasing potential is known, Boltzmann averages can be obtained by scaling the averages obtained from the biased simulation, with the corresponding Boltzmann factor related to the bias potential N hA(rN )eβ u(r ) i∗ hAi = (2.20) N heβ u(r ) i∗ where the subscript ∗ denotes that the average is determined according to the biased probability and β = 1/kB T .

2.6.2

Metadynamics

Metadynamics (MTD) [140] is a simulation method that, in analogy with US, is used to calculate the free energy of a system along one or several reaction coordinates, also renowned as collective variables (CVs). During a MTD simulation the system’s position in the space described by the CVs is monitored, 44

and as the system evolves through points in this space the system Hamiltonian is biased by added Gaussian functions. The Gaussian functions eventually fills the potential energy wells along the CVs one by one, forcing the system over higher barriers. When all the energy wells are filled up by Gaussian functions, the system exhibits free diffusion along the CVs, and the free energy can be evaluated from the sum of the biasing potentials. In normal MTD the added Gaussian functions are static, and it is difficult to determine when to stop the simulation, which can result in overfilling of the free energy landscape. There exist more elaborate methods to deal with convergence of the MTD simulations, and one of the most used is called Well Tempered MTD [141]. In this method, the Gaussian functions have an adaptive height, and the system is converging towards additions of Gaussians with zero height, thus making it easier to know when the system is converged. The Well Tempered MTD simulations used in this thesis were performed using the PLUMED2 [142] software, which can act as a plugin package for several simulation programs, including Gromacs.

45

46

3. Membranes composed of glycolipids and phospholipids (I)

Glycolipids are important molecules in biology. Apart from acting as structural components in various types of cell membranes [22], they have been found to play important roles in the communication between procaryotic cells [21] and in the processes of photosynthesis in higher plants [23; 25; 143]. Before 2005, simulation studies of glycolipids incorporated into membranes were scarce, with only a few published articles on the matter. The lack of generalized glycolipid parameters and the computer power needed to run simulations with explicit lipid molecules and solvent made it difficult to run any large scale simulations at that time. Most simulation studies during the 1990’s focused on structure determination based on data from NMR experiments [144–146], with MD simulations of single molecules on the picosecond timescale up to a few nanoseconds. One of the first attempts to model glycolipids in a membrane environment was conducted 1992 by Ram et al. [147]. In this study, glycolipid moieties were simulated in a dielectric representation of the membrane using the AMBER force field and structural preferences depending on the dielectric environment were noticed. In 2005, Róg et al. successfully parametrized and ran 25 ns simulations of fully hydrated bilayers of 128 1,2-di-O-palmitoyl-3-O-D-glucosyl-sn-glycerol (DP-GLUC) [148]. They used parameters within the OPLS force field, and it was probably the first simulation study of glycolipids in a larger membrane patch ever performed. Subsequently, in 2007, they published a more elaborate test of different type of glycolipid head groups and their effect on membrane properties [149]. At the same time the Glycam06 carbohydrate force field was extended with parameters for lipids, including glycolipids, to be run in the AMBER suit of simulation software [150]. The first investigation I was engaged in as a new PhD student in early 2010, was an ongoing simulation project with the aim of looking at structural and dynamical changes in bilayers with mixtures of DMPC (Fig. 1.7) and MGDG (Fig. 1.8). The project was based on an earlier in-house experimental NMR study of similar mixtures of the same lipid types, performed by Castro et al. [151], who recognized a continuous phase transition in the 31 P NMR spectra with increased amount of MGDG in the bilayers. The observed spectra 47

showed typical line shapes [152] corresponding to a phase transition from a lamellar Lα phase formed by DMPC to a reversed hexagonal, HII , formed by MGDG. The main idea behind the simulation study was to complement the NMR experiments with structural and dynamical information from MD simulations of similar systems.

3.1

Force field construction

The interaction model used for the simulations of MGDG (Paper I), was a mixture of parameters from different force fields. The GAL part of the head group was modeled using OPLS-AA parameters for carbohydrates [83], whereas the rest of the lipid was modeled with parameters from Berger et al. [67]. Apart from the lack of publicly available glycolipid parameters at the time (2008), the motivation for the mixing of parameters in this particular case, comes from the fact that the Berger lipids utilize atom types from OPLS internally. Exceptions are two hydrocarbon parameters that are parametrized from pentadecane to improve the lipid properties, and charges from Chiu et al. [153]. In the Berger parameter set, the non-bonded interactions with water and other Berger lipids are treated as OPLS/Berger, and interactions with all other molecules as Gromos. Thus, the combination of OPLS and Berger should not impose a large methodological problem, as long as no other Gromos molecules are present. Some fine tuning, based on comparisons of similar atom types in both force fields, of the molecular charges close to the glycosidic linkage had to be made to ensure a neutral molecule. The glycosidic linkage was treated with OPLSAA parameters on the GAL side of the glycosidic oxygen, and with Berger parameters on the other, including dihedrals originating therefrom. One problem that arises when combining force fields in this way is the treatment of non-bonded so called 1-4 interactions* . This type of interaction is typically scaled by a constant factor from the full interaction strength, usually implemented to balance the intramolecular interactions with the intermolecular to ensure stable molecule structures. The scaling factor may differ between force fields, which is indeed the case between Gromos and OPLS; which use scaling factors of 1.0 and 0.5, respectively. Berger et al. solved this problem by using Ryckaert-Belleman [154] dihedral parameters that includes the 1-4 interactions, together with explicit lists of 1-4 pair interactions. Mixing of 1-4 scaling factors has generally not been implemented in any available simulation software until recently (AMBER 11 [155]), but a trick referred to as the half-ε double pairlist method [156], utilizing additive interaction potentials in Gromacs to accomplish mixed scaling, can be used in some rare cases. Thus, to * The

48

interaction between atoms that are separated by three covalent bonds.

accomplish the mixed scaling this method was applied to non-bonded interactions originating from the glycosidic oxygen with atoms on the Berger side of the glycosidic linkage. The 1-4 interactions in GAL were explicitly scaled according to the OPLS methodology together with interactions between atoms in GAL with atoms described by Berger parameters further down in the molecule. Figure 3.1 shows the MGDG molecule with indications of the regions in which the two parameter sets are active . C4-C7

C18

O

g1

O O C18

C4-C7

g2

g3

H

HO O

HO

OH G4

O G6

G1

OH

O

Berger

OPLS

Figure 3.1: A schematic picture of the 18:3-18:3-MGDG molecule. The regions of the molecule treated with OPLS-AA and Berger parameters are indicated in the figure. The intermediate region requiring special treatments is shaded in grey.

A second problem that arises when mixing force fields is the choice of water model. The Berger force field was parametrized towards the simple point charge (SPC) [62] water model, while the OPLS force field is optimized for TIP3P [61]. To decide which water model to choose for the combination of the two, both SPC and TIP3P were used in 100 ns simulations after which bilayer properties such as the AL and orientational distributions in lipid molecules were calculated and compared. To illustrate the rather small difference, a comparison of some properties is presented in table 3.1. Both water models produce averages of structural properties within 2%, which is the reason for not showing results from both models in Paper I. Note that the SPC simulations were extended after this comparison, which is the reason for differences between the values in Table 3.1 and values reported in the paper. Table 3.1: Comparison of properties from simulations with different water models. The averages are calculated from systems corresponding to MG-45 (45% MGDG in a bilayer with DMPC) in Paper I. Here AL is the area per lipid determined by a grid approach (Section 2.5.1) and θP-N is the tilt angle of the vector between P and N atoms in DMPC with respect to the bilayer normal.

Water model SPC TiP3P Diff./%

ADMPC /nm2 L 0.62 ± 0.04 0.61 ± 0.04 1.6

AMGDG /nm2 L 0.61 ± 0.03 0.61 ± 0.05 0

θP-N /deg 70.0 ± 7.0 69.0 ± 9.0 1.4

49

3.2

Simulated systems

The simulated systems were tailored to be as close as possible to some of the measurements in the previous NMR study [151]. In the NMR experiments, the measurements were performed on 50:50-mixtures of 18:3-18:3-MGDG (with 18 carbons and three double bonds in each tail) and 18:3-16:3-MGDG. The 18:3,18:3-MGDG is the most abundant variant in nature [157]. Therefore, since a ternary mixture of lipids would be unpractical in the first MD simulation study of MGDG, only the 18:3,18:3 variant was used in the simulations. Three main bilayer systems were simulated; the reference system MG0 with 128 DMPC molecules (64 in each leaflet), MG-20 with 36 MGDG and 144 DMPC (20% MGDG) and MG-45 with 80 MGDG and 96 DMPC molecules (45% MGDG). The bilayers were hydrated with 25 water molecules per lipid, which could be considered close to fully hydrated if compared to the number of 29 water molecules per lipid reported by Högberg et al. [158] to be the limit for full hydration of a DMPC bilayer. The level of hydration directly affects the spacing between lipid molecules. A decrease in hydration leads to a decrease of AL as water molecules are removed from the bilayer interface. Following the decrease in lateral area is an increased order, which eventually will increase the gel transfer temperature Tm , but severe effects on the bilayer properties of phospholipid bilayers are not observed until the level of hydration decreases below 12 water per lipid, the limit of the lipid hydration shell [159]. In Figure 3.2, a snapshot of the simulation box corresponding to MG-45 is displayed.

Figure 3.2: A snapshot from the simulation of MG-45. DMPC is colored grey to highlight MGDG. The image was created with the VMD software [46].

50

3.3 3.3.1

Results Structural properties

The AL for DMPC and MGDG were calculated according to the grid based method (Section 2.5.1). The polygons in the resulting grid map of the lateral plane of the bilayer (Fig. 3.3a) was shown to correspond to individual lipid molecules if the number of grid points exceeds 200 × 200 [111] for a 128 lipid bilayer. Thus, a 250 × 250 grid was used in the area calculations. 250

a gy

0 250

b gy

0

0

gx

250

Figure 3.3: Plot of the 250 × 250 grid used for the calculations of area per lipid. The plot is a snapshot of the lower leaflet, taken from a single time frame in the MG-20 simulation. The lipid molecules are mapped onto grid points (gx , gy ) along the x- and y-axes of the simulation box. The polygons correspond to a) single lipids of any type and b) MGDG (the black background corresponds to DMPC).

The rather broad distribution of areas (Figure 2, Paper I) is a reflection of the surface roughness, since the projection to the grid only takes the xy-position of the center of mass of the lipid into account. Thus, neighbouring lipids situated at different z-positions may have overlapping areas, resulting in smaller values of AL for the lipid closest to the bilayer center. The average AL for DMPC (MG-0) was in reasonable agreement with experimental values [107; 160; 161] that are reported in the span 0.60-0.67 nm2 depending on the level of hydration and the temperature. The result is also consistent with simulation studies 51

O O P – O O O

180

O

90 θ/°

O

0

O

0

θ

H

P(θ)

Bilayer normal

0.8

N+

at similar hydration levels [158; 162]. The average area calculated for MGDG is more difficult to directly compare to experiments due to the lack of experimental reports of MGDG in similar environments, but it is reasonable to assume that the area increases with concentration as the lamellar phase undergo lateral stress [16] due to the hexagonal nature of MGDG. An experimental study [163] of MGDG monolayers report an area per lipid of 0.82 nm2 which is substantially larger than the values calculated in the simulations. Accordingly, upon increased amount of MGDG in the simulations, the AL increase by 2 nm2 for MGDG. The AL for DMPC, however, decrease by 1 nm2 . The decrease in area for DMPC was ascribed to a closer packing of the molecules in the bilayer. This is also indicated by the observation that the orientation of the vector between the P and N atoms in the DMPC headgroup becomes more parallell to the bilayer normal (Fig. 3.4).

Figure 3.4: Angle of the P-N vector in DMPC in simulations with 0% (green), 20% (red) and 45% (black) MGDG.

The distribution of molecules in the bilayer can be visually observed in the grid used for the area calculations. In Figure 3.3b, the distribution of MGDG is accentuated by plotting the molecular type instead of individual lipids related to the grid. Using visual inspections of simulation snapshots throughout the simulations, together with the RDF calculations (Fig. 5, Paper I) that show only weak preferential lateral ordering of the lipids, we concluded that no clear lateral phase separations of the lipids are taking place within the timescale of the simulation. A recent investigation [164] of the role of MGDG in the thylakoid membranes of higher plants indicate coexistence of the Lα and HII phases in the bilayer stacks typical for these membranes. The ratio of the two phases was found to be regulated by the level of hydration (low hydration increases the amount of HII ) and the ratio of MGDG and the Lα forming digalactosyldiacylglycerol (DGDG) lipids. The systems studied in Paper I correspond to MGDG concentrations which show lamellar characteristics in the NMR investigation by Castro et al. [151], and thus no HII phase should be observed in the simulations. In addition, NMR investigations [44] of fast-tumbling DMPC/DHPC52

Electron density / e nm-3

bicelles intermixed with MGDG, did not reveal any phase separations at any MGDG concentration (≤ 30%). The bilayer thickness was estimated as the distance between the peaks in the z-position distributions of the phosphorous atoms in DMPC and the center of mass of the GAL-groups in MGDG. The results 3.33, 3.56/3.20 and 3.69/3.48 nm are in reasonable agreement with experiments [20; 105; 160] and simulations [158] of phospholipid bilayers. Calculating the thickness as the distance between the peaks in the total electron density profiles (Fig. 3.5) results in similar values; 3.1, 3.2 and 3.3 nm for MG-0, MG-20 and MG-45, respectively. Worth noting is that MG-0, compared to MG-20 and MG-45, shows a higher electron density in the parts of the density distribution that correspond to bulk water. The cause of this difference is unknown, but it could be related to osmotic effects as a consequence of changed surface structure in the bilayer when MGDG is added. 450 400 350 300 MG-0 MG-20 MG-45

250 200

-4

-3

-2

-1

0 1 rz / nm

2

3

4

Figure 3.5: The electron density profile along the bilayer normal for MG-0 (green), MG-20 (black) and MG-45 (red).

3.3.2

Orientational order and NMR dipolar couplings

The C-H order parameter, SCH , is often calculated in simulations to investigate membrane properties and to validate the correctness of force field parameters. If the SCH profile is in good agreement with experimental results, other characteristics of the lipid bilayer such as the AL is likely to be correct [67]. In the NMR experiments, the sign of the dipolar couplings cannot be readily determined, whereas the sign can be determined directly from the simulations. If the conformational transitions of the molecules is fast enough in the simulations to give well determined dipolar couplings (i.e. the errorbars do not overlap both the positive and negative values of the experimental couplings), the correct sign can be assigned to the experimental couplings using the information from the simulations. 53

For DMPC (MG-0), the calculated order parameters show typical characteristics for a phospholipid bilayer, and the corresponding NMR dipolar couplings show a good agreement with the couplings reported in the NMR experiments [151]. Calculated couplings in the headgroup of DMPC are also in good agreement with experimental couplings in this region. This indicates that the DMPC molecules are well described in the simulations, and the sign of the experimental couplings can readily be assigned from the results in the simulations. The tail order parameters in DMPC increase slightly when MGDG is included in the bilayer. Order parameters in the acyl chains of MGDG are in good agreement with the experiments. The double bonds show typical characteristics with SCH on C10 close to zero (Fig. 3.6a), related to an orientation of the CH-vectors in this position close to the magic angle, in good agreement with both the NMR experiments by Castro et al. [151] and earlier investigations of unsaturated lipids [165]. Dipolar couplings in the GAL hadgroup of MGDG (Fig. 3.6b) are not well determined from the simulations. Large errorbars in this molecular segment make unambiguous assignments of the correct signs for the experimental couplings impossible. We concluded that this result is related to slow reorientation of the torsional angles around the glycosidic linkage, which also is demonstrated in the rotational correlation functions calculated for vectors in the GAL moiety in MGDG and the PC-hedgroup in DMPC. The slow dynamics may be an effect of an incorrect force field description of the glycosidic linkage and since other glycolipid force fields now are available, an appropriate route to investigate this problem would be to rerun the simulations in another well known force field such as Glycam [150].

4 D / kHz

-SCH

0.2 0.1

a

0

-8

(b) g3 (a) g3 g2 (S) G6 ) (R G6

G5

8 10 12 14 16 18 Acyl chain carbon

G4 G3

6

b

-12 G2

4

-4

G1

2

0

Figure 3.6: Tail and headgroup order in MGDG (MG-45) represented as a) acyl chain order parameters and b) calculated dipolar couplings (black circles) in the GAL headgroup. The red squares and blue triangles are the positive and negative experimental [151] values, respectively.

54

3.3.3

Lateral diffusion of lipids

The lateral diffusion coefficients, Dlat , of the lipids in the three systems were calculated by linear fits to the MSD of the lipids (see Section 2.5.6 for details). The diffusion of DMPC in MG-0 was calculated to 5 µm2 s−1 , which is in agreement with diffusion coefficients calculated from simulations of DMPC bilayers with the same force field [158]. Experimentally determined diffusion coefficients in DMPC bilayers, which typically ranges between 10 - 20, µm2 s−1 [132; 166], were however around 2 to 4 times larger than the calculated Dlat . Diffusion coefficients calculated from simulations are highly dependent on system size [167] and the lipid parameters used, and it is a quantity that is not fitted during force field development. Discrepancies between experimental and calculated diffusion coefficients are common, but still an important concern if the dynamics of a system is of interest, which was also pointed out by one of the referees during peer review of Paper I. In MG-20 and MG-45 both DMPC and MGDG are subjects do a decrease in Dlat with increasing amount of MGDG. In MG-45 the diffusion is roughly decreased by a factor of 2. Pulse field gradient (PFG) NMR measurements of lipid lateral diffusion in a DMPC bilayer [132] and the glycolipid rich Acholeplasma laidlawii membrane (consisting of 70-80% glycolipids and 20-30% phospholipids) [168] resulted in Dlat =16 µm2 s−1 and Dlat =5 µm2 s−1 resepectively. By assuming an ideal mixture of lipids in the bilayers (eq 3.1) and utilizing the results from PFG-NMR, we estimated Dlat for MGDG to be 2.2 µm2 s−1 , which is close to the value calculated in MG-45: Dtot = DDMPC xDMPC + DMGDG xMGDG

(3.1)

where xDMPC and xMGDG are the fractions of the lipids in the bilayer. The approximation is crude, and the fact that DMPC follows the decrease in Dlat by roughly the same amount as MGDG indicates that the assumption of an ideal mixture is far from true. However, the small size of the simulation box may accentuate collective effects of the lipids that are not accounted for in the MSD calculations. In a simulation of a one component bilayer consisting of 72 dipalmitoylphosphatidylcholine (DPPC) a three times larger Dlat than in a bilayer with 288 lipids was observed [167]. This phenomenon was ascribed to a finite size effect of correlated lipid motions. In addition a 25% reduction of the calculated lateral diffusion was observed when removing the center of mass (COM) motion of the individual leaflets for the simulation of 72 lipids. In the 288 lipid simulation the reduction was only 6%. Based on earlier reports of unrealistic COM motion of the individual leaflets in small size molecular simulations [169; 170], the simulations in Paper I utilized COM motion removal for the individual bilayer leaflets during the simulations. Whether or not this 55

was a good choice remains to be tested, and simulations of larger membrane patches would be necessary to give a more reasonable picture of the lateral dynamics of the lipids.

56

4. Investigations of membrane– sugar interactions (II–IV)

Trehalose (Fig. 1.2b) is subject for a lot of controversy in the scientific litterature. The fact that small sugars in general, and TRH in particular, show important abilities to protect organisms from cellular damage due to anhydrobiosis and freezing has been recognized for a long time [3; 171–174]. The processes involved at the molecular level are however still actively debated [175]. In their attempt to reconcile the opposing views of the molecular interactions involved at the membrane interface, Andersen et al. [176] recently summarized the field into two main interaction hypotheses: 1. A direct interaction model (Fig. 4.1) in which the TRH molecules interact directly with the lipid head groups at the membrane surface. The model is further divided into three plausible interaction scenarios: a) Preservation of lipid–lipid headgroup distances by intercalation of TRH into the surface and substitution of water molecules, b) strong interactions with both lipids and water, keeping the hydration shells around the lipids intact or c) vitrification (glass formation) of sugar molecules on top of the membrane surface. 2. An exclusion model (Fig. 4.2) in which the sugar molecules are residing outside the hydration region of the lipids keeping the original hydration level of the lipids intact as the water amount decrease during the anhydrobiotic stress. By combining results from Small Angle Neutron Scattering (SANS) experiments and dialysis of bilayer systems consisting of DMPC and small sugar molecules known to stabilize membranes namely TRH, Glucose (GLU) and Sucrose (SUC), Andersen et al. [176] conclude that the hypotheses mentioned above are not mutually exclusive. They base their conclusions on the concept of preferential binding, or affinity, which in practice means that at low concentrations the sugar molecules interact directly with the membrane surface (direct interaction model). At a certain concentration threshold, the surface becomes saturated with sugar molecules, after which all added sugar molecules will be excluded from the surface (exclusion model). In paper II, we performed 57

ge

a

ter Wa

an ch ex

Water entrapment

b

Vit ri

fic ati

on

c

Figure 4.1: Illustration of the direct interaction model: a) Water exchange and intercalation of TRH between lipid head groups, b) Water entrapment: TRH interacts strongly with both water and lipids and c) Vitrification of TRH on top of the bilayer.

Exclusion

Figure 4.2: Illustration of the exclusion model: The TRH molecules are excluded from the hydration region of the bilayer.

simulations of a DMPC bilayer at several TRH concentrations with results supporting their conclusions. The project started as a collaboration with a research group at the Department of Organic Chemistry at Stockholm University. To characterize sugar molecules they utilize NMR spectroscopy in solution to extract conformational and structural information from measured J- and residual dipolar couplings. The J-couplings are mediated through the intramolecular bonds, and are readily available from NMR measurements in isotropic solution. The magnitude of the J-couplings contain information about atom connectivity, bond lengths and conformations. The latter is usually obtained through the so called Karplus relations [177; 178] that relates J-couplings to dihedral angles in the molecule. Typically, the information from J-couplings is complemented by intramolecular distances obtained from measurements of the Nuclear Overhauser Effect (NOE) [179], which is a distance dependent cross relaxation effect of excited spins mediated through space between neighbouring atoms in the molecule. The residual dipolar couplings, on the other hand, is mediated through space at long range distances and depend on molecular orientation through the order 58

parameter S (section 2.5.2) which is effectively averaged to zero in isotropic media. Thus, to be able to measure the dipolar interactions, the samples must possess some orientational order which is commonly accomplished by using alignment media in the form of gel forming molecules or lipid bicelles (Fig. 1.6) at low concentrations. The analysis of the structural data originating from aligned samples relies on the important assumption that the conformational space of the measured molecules is unaffected by the medium used for alignment. Lipid bicelles were used in the NMR measurements of dipolar couplings in TRH, and the original aim of the MD simulation project (Papers II-IV) was to validate the experimental method by structural investigations of TRH in the presence of a phospholipid bilayer. In Paper III, residual dipolar couplings from the NMR experiments are fitted to order parameters calculated from the MD trajectories. The fits are used together with conformational calculations from MD simulations in both isotropic solution and with the DMPC bilayer to validate both the sugar force field and the assumption of unaffected conformational space of TRH in the prescence of bicelles.

4.1

Modelling of the trehalose interactions (Paper II)

The bicelles used as alignment medium in the NMR study of TRH consisted of DMPC and the shorter DHPC in a 3:1 ratio. Since the largest part of the bicelles is thought to be lamellar (DMPC) the natural choice for the simulations was to approximate the bicelles as a pure DMPC bilayer. At the time, the most recent experimental study of TRH interacting with a DMPC bilayer was the work by Andersen et al. [176]. Having reliable experimental data at hand is of great importance to validate the results from any simulation study, and thus we chose to simulate several TRH concentrations within the same concentration range used in the experimental study to enable a direct comparison with their results.

4.1.1

Choice of force field

The simulations were planned to be executed in Gromacs, and the first force field candidate was Gromos. A promising set of updated parameters for carbohydrates, compatible with the Gromos 53A6 [180] parameters, were released in 2011 by Hansen and Hünenberger [85], and lipid parameters were available from Kukol [76] and Poger [77]. Early simulation tests of a DMPC bilayer indicated problems with the lipid dynamics in both the Kukol and Poger parameters, with calculated lateral diffusion coefficients close to one order of magnitude lower than typical experimental coefficients (>10 µm2 s−1 ). Based on these early tests, we decided to instead use the Stockholm Lipids 59

(SLipids) [70] force field, another recently developed lipid force field showing promising results for properties related to both structure and dynamics. Since SLipids claims to be compatible with Amber parameters, a reasonable choice of parameters for TRH was to use the well established (and Amber compatible) Glycam [78] force field. The Glycam force field does not use any special scaling of 1-4 interactions, while Amber uses a scaling of 0.8333, so a similar trick [156] as used for MGDG was used for the TRH parameters to accomplish a mixed scaling of 1-4 interactions in Gromacs. Potential of mean force calculations (Section 2.6) of a TRH molecule approaching a DMPC bilayer, revealed similar interaction strengths for both Hansen/Kukol and Glycam/Slipids parameters, which implies that simulations of TRH/DMPC systems with these two combinations are probable to yield similar information, given that the solubility of TRH is similar.

4.1.2

Simulated systems

To investigate how the properties of a DMPC bilayer changes with the TRH content, a total of ten main systems were constructed. The bilayer consisting of 128 lipids was simulated with TRH concentrations ranging from 0 to 20% (w/w). The simulations were labelled TRHN, where N corresponds to the number of TRH molecules in a particular simulation and ranges between 0 and 260. Even at low concentrations clustering of TRH molecules is substantial, a tendency which previously have been reported from NMR studies [181] in isotropic solution. At higher concentrations the clustering of TRH, together with a strong attraction of TRH towards the bilayer surface, enabled the formation of aggregates connecting the upper and lower bilayer leaflets over the periodic boundaries. To overcome this effect, the TRH molecules had to be added to the boxes in a titration like manner, i.e. adding TRH to the simulation boxes in smaller amounts with equilibration in between. A problem with this approach is that TRH molecules interacting directly with the bilayer may be trapped beneath clustering TRH molecules on top of the bilayer surface, which can lead to problems with ergodicity if the rearrangement of molecules within the clusters is slow compared to the timescale of the simulation (∼ 300 ns). Such problems have recently been addressed [182] in bilayer systems with cholesterol, in which a substantial decrease in diffusivity of the system leads to a switch of timescale for the dynamics from nanoseconds to several microseconds or even longer. A typical simulation box (TRH80) is displayed in Fig. 4.3. Note that all TRH molecules are situated on the bilayer surface, either interacting directly with the lipid headgroups or with other TRH molecules. The amount of TRH residing in the water phase is very small throughout the simulation. 60

Figure 4.3: Snapshot of the TRH80 simulation (Paper II). The TRH molecules (green) interact strongly with both the DMPC bilayer and other TRH molecules resulting in clusters on the bilayer surface. For clarity, the water molecules have been removed from the picture. The blue rectangle illustrates the size of the simulation box. The image was created with the VMD software [46].

4.1.3

Results

To examine the effect of TRH on the structure and dynamics of the DMPC bilayer, the AL , bilayer thickness (dP−P ), lateral diffusion coefficient, Dlat and the area compressibility modulus, KA were calculated for all ten simulations. A clear trend of increasing AL and decreasing dP−P with increased amount of TRH could be observed. This is consistent with the results from SANS experiments [176] of similar DMPC/TRH systems, as well as electron spin echo experiments of systems with 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) and TRH at low temperature [183]. The trends are in line with the direct interaction hypothesis, and indicate that TRH molecules are intercalating the lipid headgroup region in accordance with the water exchange model. 61

rc

a

200

b

NB

PMF/kBT

5

0

-5

0

1

2 3 rz /nm

4

100

5

0

0

100

N

200

300

Figure 4.4: In a) the PMF used for the parameters in the interaction model (eq. 4.1) is displayed on top of the bilayer to indicate the binding region. The cutoff distance, rc , is indicated in the figure with a dashed vertical line. In b) the calculated amount of bound TRH together with the model is displayed as a function of total number of TRH molecules. The parameter a0 is estimated as 0.50 nm2 and the cutoff rc is 2.8 nm.

At the same time, a drastic slowdown in dynamics can be observed as a tenfold decrease of Dlat at high TRH concentrations. Following the decrease in lateral diffusion is a drastic increase in the area compressibility, KA ; an observation that come with a substantial amount of uncertainty since the slow dynamics increase the timescale needed to obtain reliable time averages of the area fluctuations. To further scrutinize the interactions of TRH with the bilayer, the PMF for a TRH molecule approaching the bilayer was calculated by the means of umbrella sampling (Fig 4.4a). The PMF clearly shows an attractive part of the interaction potential in the headgroup region of the bilayer. Utlilizing the PMF, together with assumptions about the size and volume of the TRH molecule in the bulk and on the bilayer surface, a simple two state binding model was suggested: " #   NB NB −vB /vF N = nF 1 − − 1− (4.1) nF nB where N is the total number of TRH molecules, NB is the number of bound molecules, nF the total number of free sites, nB the total number of binding sites, and vB and vF are the free volumes of the bound and free TRH states, respectively. The parameters nB and vB were estimated using the average area of the bilayer, the integral of the PMF in the attractive region, and an area estimation of the bound TRH molecule. The parameters, nF and vF were determined from the volume of the free water phase and the volume of a TRH molecule in solution which were derived from the solubility of TRH [184]. Probability density distributions of TRH along the axis perpendicular to the bilayer plane were used to calculate the amount of bound TRH molecules, NB . 62

Since the boundary between bound and free molecules (here free molecules can be interpreted as not bound directly to the lipids) is diffuse, cutoffs corresponding to the attractive part of the PMF were used to estimate the binding region (Fig. 4.4a). The calculated NB values (Fig. 4.4b) are closely following the interaction model (eq. 4.1), which is in agreement with dialysis experiments [176], a recent investigation of membrane–sugar interactions by Pulsed Electron Paramagnetic Resonance of spin labels [185] and a characterization of sugar interactions with a DPPC monolayer [186], which indicates specific lipid–sugar interactions, but a relatively weak thermodynamical driving force at room temperature. However, calculated partition coefficients P = [TRH]bound /[TRH] f ree = zb /z f

(4.2)

where zb and z f corresponds to the number of available bound and free states respectively, are ten times larger in the simulations than in experiments [176]. This difference can be ascribed to a one order of magnitude higher bilayer concentration in the simulations compared to experiments [176], which leads to a ∼ 10 times higher affinity of TRH towards the bilayer (Figure 4, Paper II). Note that the model (eq. 4.1) does not take TRH–TRH interactions into account, and since we observe a substantial clustering of TRH molecules on the bilayer surface this may have to be taken into consideration in future simulation studies of TRH–membrane interactions. The trends in AL and dP−P together with the observations of bilayer saturation at ∼ 1 TRH per lipid molecule are indications of direct TRH–bilayer interactions. Above the saturation point, added TRH molecules are excluded from direct interaction with the bilayer. Thus we can conclude, that both main interaction hypotheses suggested by Andersen et al. [176] can explain the membrane–sugar interactions in different concentration regions, and that they are not mutually exclusive. The intercalation of TRH between lipid molecules support the water exchange model of the direct interaction hypothesis, and the clustering of TRH molecules on the bilayer surface that is followed by a reduction in dynamics is a clear indication of glass formation (vitrification). Figure 4.5 shows TCFs of TRH reorientation at different distances, rz , from the bilayer center of mass. It is clear that the change in dynamics is drastic even for a single TRH molecule approaching the bilayer, and a strong TRH–TRH interaction at the interface is thus likely to result in highly immobilized aggregates. Vitrification of TRH is usually observed at reduced hydration [187], and it is still unclear if the high clustering of TRH in the simulations are physically correct or an artifact stemming from unbalanced interaction potentials in the force field. 63

1 rz = 1.5

C2(t)

rz = 2.0 0 rz = 2.5 0

0n m

nm

nm

rz > 3 nm

0 0

1

2

3

4

5

t / ns

Figure 4.5: Rotational time correlation functions calculated from umbrella sampling trajectories at different distances from the bilayer center of mass. The blue curve at rz = 2.0 nm corresponds to the minimum of the PMF in Fig. 4.4a.

4.2

Validating the results (Paper III)

In Paper III we return to the original aim of the simulation study: to validate the MD simulations and the assumption that the conformational space of TRH is unaffected by the dissolved bicelles. The combination of MD simulations and NMR experiments provide a powerful tool to investigate systems at the molecular level. Cases where both the magnitude and relative sign of the dipolar couplings can be retrieved from NMR are rare. In contrast, both the sign and magnitude of all possible dipolar couplings can readily be calculated from simulations and, if the force field is accurate, simulations can successfully be used to complement experimental data with such information. Thus, to be able to validate our simulations and possibly give some insight to the technique of using bicelles to accomplish anisotropy, a set of dipolar couplings was measured by NMR and subsequently compared to the simulation results from Paper II.

4.2.1

DMPC dipolar couplings

To validate the lipid force field (SLipids [70]), a set of dipolar couplings previously measured in DMPC by the solid state NMR [151] was compared to calculated couplings (eq. 2.6) from the TRH0 simulations. A good fit to the experiments was obtained for all but two couplings in the glycerol backbone of the lipid headgroup (Fig. 4.6). In these two sites the magnitudes of the dipolar couplings are too small compared to the experiments, and it has been ascribed to a shortcoming in the parametrization of the lipid force field [188]. The problem originates from inaccurate parameters for the torsions in the glycerol backbone, which result in a too fast reorientation around the g2 and g3 carbons. Private communication with the force field developers indicate that this 64

6

d / kHz

4 2 0 -2 -4 -6 C2g1B C2g1A g1 g2 g3 α β γ

Figure 4.6: Some of the dipolar couplings calculated in DMPC. The negative and positive values of the experimental values [151] are represented as black and blue circles respectively. The calculated couplings in g2 and g3 deviate substantially from the experiments due to shortcomings in the force field.

problem is about to be resolved.

4.2.2

TRH dipolar couplings

Calculated dipolar couplings in TRH from the TRH20 simulation are three orders of magnitude larger than couplings measured in our NMR experiments. Basically all sugar molecules are attached to the bilayer throughout the simulations, and to achieve order parameters close to experiments (in the order of 10−4 compared to 10−1 in the simulations) the contribution to the time average from bound TRH molecules would have to be decreased from 1/1 to around 1/1000. The discrepancy between simulations and experiments can only to some extent be explained by a difference in lipid concentration, since the water content in the NMR experiments is only around two times higher than the simulations. In Figure 4.7a, the calculated dipolar couplings are linearly fitted to the experimental dataset. The slope of the linear fit was subsequently used to scale the calculated couplings for a direct comparison with the experimental counterparts. In Figure 4.7b, the experimental couplings are plotted together with the scaled calculated couplings and it is clear that a linear scaling can describe the difference between them in most cases. The H6S-H6R coupling is one of few exceptions. This particular coupling is not converged, and it is thus excluded from the fit as an outlier. Overall, it seems that the one bond C-H couplings are all poorly averaged during the simulation, which is manifested by large errors of up to ∼ 3 kHz. However, the H6S-H6R coupling is the only coupling that affects the linear fit drastically. For examination of the structural information contained in the dipolar coup65

8

a

dsim / kHz

4 y(x)=1.35x

- 0.20

0

-4 -8

d / Hz

2

-1

-0.5

0

0.5 dexp / Hz

1

1.5

2

b

0 -2 C6-H4 C5-H6R C5-H4 C4-H6R C4-H3 C4-H2 C3-H4 C3-H2 C3-H1 C2-H4 C2-H1 C1-H3 C1-H2 C1-H1’ C6-H6S C4-H4 C2-H2 C1-H1 H6S-H6R H6R-H5 H4-H5 H3-H4 H2-H3 H1-H2

Figure 4.7: Dipolar couplings in TRH. a) Linear fit of the experimental dipolar couplings, dexp and the calculated dipolar couplings, dsim from the TRH20 simulation. The dipolar coupling marked in red, corresponding to H6S-H6R, is excluded from the fit. b) The experimental dataset (red circles), backcalculated couplings from a fit (eq. 4.3) to simulations (black crosses), and dipolar couplings directly calculated from the simulation TRH20 and scaled by the slope of the linear fit in a) (blue circles). The errorbars of the calculated couplings are left out for clarity.

lings, eq. 2.6 can be transformed into a relationship between a coordinate system fixed in the molecule (Fig. 4.8) and the vector connecting the interacting spins, through the elements (order parameters) of the Saupe matrix, S [189]:    µ0 γi γ j h¯ h  2 z 2 x 2 y di j = − × S 3 cos θ − 1 + (S − S ) cos θ − cos θ zz xx yy i j i j i j 16π 2 ri3j i + 4Sxy cos θixj cos θiyj + 4Sxz cos θixj cos θizj + 4Syz cos θiyj cos θizj (4.3) where θiaj (a = x, y, z) are angles between the spin–spin vector and the molecular axes (Fig. 4.8). The elements of S, given by Sab (a, b = x, y, z), are called the order parameters. Equation 4.3 assumes that all possible molecular conformations are accessed, since ri j , Sab , and θiaj become dependent on the molecular structure if the dipolar interaction, di j , occurs between spins in different 66

molecular fragments. The TRH molecule exhibit a C2 -symmetry axis through the glycosidic oxygen (corresponding to the z-axis in Fig. 4.8), which results in vanishing Sxz and Syz . Thus only the three first order parameters in eq. 4.3 are needed to determine di j .

Figure 4.8: The molecular coordinate system in TRH used to calculate the directional cosines in eq. 4.3.

Since the calculated dipolar couplings in TRH are subject to very large errors, they can not be used directly to validate the correctness of the sugar parameters used in the simulations of TRH. The directional cosines in eq. 4.3 were obtained directly from the simulations, and the order parameters Sab were fitted to the experimental data set by a least square fitting procedure. The fitted order parameters were then used to back-calculate the dipolar couplings for direct comparison with the experiments. To estimate the error of the fit, sets of 10 000 randomly normal distributed datasets within the experimental errors (±3σexp ) were generated and fitted to the same number of normal distributed directional cosines (±3σsim ). The errors for the back-calculated couplings were then estimated from the standard deviation of back-calculated sets of dipolar couplings from the the 10 000 fits. The back-calculated couplings, together with the corresponding errors, are included in Fig. 4.7b. The agreement with the experimental dataset is very good, which indicates that the sugar force field is correctly describing the conformational space of TRH, and that the DMPC bilayer is not affecting it. This is also confirmed by calculations of torsion angles within TRH from both the simulation TRH20 and isotropic simulations of TRH (Fig. 6, Paper III). The torsions are clearly unaffected by the TRH–DMPC interactions. 67

4.3

Reaching longer timescales and larger systems (Paper IV)

The slow dynamics encountered in simulations of TRH–DMPC at high TRH concentrations reduce the accuracy of the calculated bilayer properties. As an example, in Paper II the calculated compressibility modulus of the DMPC bilayer increases more than tenfold, going from the reference simulation (TRH0) to the simulation with highest amount of sugar (TRH260). It is difficult to estimate how much of the increase is actually coming from changes in bilayer mechanics, and how much is coming from insufficient sampling of the area fluctuations. In addition, problems with ergodicity may also have an effect on the calculated time averages. Coarse grain simulations can be a straightforward way to boost the dynamics and reduce the degrees of freedom of a specific system. In addition, the decrease of the computational effort that follows from the reduced number of atoms, allows for simulations of several times larger systems. In Paper IV we try to tackle the problems with slow dynamics by changing the force field from the atomistic Glycam/SLipids to the coarse grained Martini force field parameters for lipids and sugars [89; 190]. At the same time the system sizes are increased by a factor of nine, corresponding to an increase from 128 to 1152 lipid molecules. The lipid type was changed from DMPC to the shorter 1,2-dilauroyl-sn-glycero-3-phosphocholine (DLPC) as a result of the coarse graining procedure in Martini, in which four heavy atoms are grouped together into one large bead. The length of the DMPC tails (14 carbons in each chain) are unfortunately not evenly dividable by four. Coarse grained representations of the molecules are displayed in Fig. 4.9a.

Figure 4.9: a) Coarse grain representations of DLPC (left) and TRH (right). b) A snapshot of a CG simulation with DLPC, TRH and polarized Martini water (PW). The TRH molecules are colored green and the water is purple. The image in b) was created with the VMD software [46].

68

4.3.1

TRH interactions

Simulations with both the standard Martini water model (W) and the polarized water model (PW), labelled MA- and MP- respectively, revealed that large clusters of TRH, attached to the bilayer surface (Fig. 4.9b), were formed within 50 − 100 ns of simulation time. Potential of mean force calculations of a TRH molecule approaching the DLPC bilayer showed large differences from the atomistic PMF (Fig. 4.10a), here also illustrated in table 4.1 by a comparison of the integral Z rc

zB =

0

e−W (rz) /kB T drz

(4.4)

where W (rz ) is the potential of mean force, and rc is the cutoff to determine the border of the bilayer region of the PMF. The integral zB is related to the partition coefficient (eq. 4.2), and implemented as a parameter related to vB in the binding model (eq. 4.1) suggested in Paper II. It is clear that the Martini parameters underestimate the TRH–bilayer interactions in comparison to the atomistic force field (including the Gromos parameters, for which the corresponding zB = 30 nm). Table 4.1: Comparison of the integral zB (eq. 4.4) from PMF calculations for different force fields and water models.

Force field Glycam/Slipids Martini Martini

Water model TiP3P W PW

2 0

-2

-4 1.0

1400 1200 1000 800 600 400 200 0

NB

PMF / kBT

4

zB /nm 43 2 3

a 1.5

2.0 2.5 rz / nm

3.0

3.5

b 0

500 1000 1500 2000 2500 N

Figure 4.10: a) Potential of mean force profiles calculated using atomistic (AAPMF, red), Martini W (MA-PMF, blue), Martini PW (MP-PMF, green), and Martini (P3-PMF, black). b) Binding isotherm calculated according to eq. 4.1 with the same colors as in a). The black symbols correspond to the number of bound TRH molecules in the CG simulations. Red triangles indicate the scaled numbers of bound TRH molecules derived from the atomistic simulations with 20, 80 and 170 TRH molecules corresponding to 180, 720, and 1530 in the CG systems.

69

Free energy / kBT

Private communication with two of the Martini developers, S.-J. Marrink and C. López, resulted in suggestions for two force field changes in TRH, that affect the TRH–DLPC and TRH-TRH interactions: 1) A change from the polar bead type P4 to P3 to improve the interactions with the bilayer and 2) a change from the bead type P2 to the ring type SP2 which decreases the TRH-TRH attraction and possibly can reduce the clustering. Modification 1 resulted in a PMF almost identical to the atomistic PMF (Fig. 4.10a) and the corresponding binding isotherm (Fig. 4.10b). Additional simulations with each of the suggested force field modifications (labelled P3- and P2-), and a combination of the two (labelled PP-), were subsequently carried out with medium and high TRH concentration (720 and 1440 TRH molecules). An atomistic system (AA0) with the same size as the CG simulations was simulated as a reference, and two systems (labelled 00-TRHL/H where L and H corresponds to 720 and 1440 TRH molecules respectively) were simulated with all TRH–TRH interactions turned off. In addition to the larger simulations, three Well Tempered MTD simulations (Section 2.6.2) were carried out for atomistic (AA-MET), normal Martini (MP-MET) and modification 2 (P2-MET) to calculate the free energy of interaction between two TRH molecules in solution (Fig. 4.11). 4 2 0

-2

-4 -6 0.0

0.5

1.0 1.5 rij / nm

2.0

2.5

Figure 4.11: Free energy of interaction between two TRH molecules simulated with atomistic (AA-MET, red), Martini PW (MP-MET, green) and Martini modification 2 to decrease TRH–TRH attraction (P2-MET, black).

Since the first modification, resulted in a PMF (P3-PMF) in close resemblance to the atomistic PMF (AA-PMF), a similar behavior in the larger simulation as in the corresponding atomistic simulations was expected. Figure 4.10b gives a different picture, where the number of bound TRH molecules (calculated from distributions of TRH molecules along the axis parallel to the bilayer normal) all follows the binding isotherm corresponding to the unmodified Martini parameters independent of the force field modifications. From these results we concluded that the PMF is not a valid descriptor of the TRH–bilayer interactions in the CG simulations. In our atomistic study (Paper II) we saw that the binding of TRH to the DMPC bilayer is following a simple two state binding isotherm, which is consistent with recent experimental findings. In Paper IV 70

we suggest that the discrepancy from this behavior in the CG simulations is related to the balance between TRH–DLPC, TRH–TRH and TRH–Water interactions in the system, and that further tuning of the CG parameters is needed for a good description of the TRH–bilayer interactions. Only looking at the free energy of interacting TRH molecules in Fig. 4.11 is not sufficient to describe the clustering behavior in the larger simulations. The attractive well in the interaction free energy profile corresponding to the atomistic parameters is two times deeper than in the corresponding profile for normal Martini. One plausible explanation to the differences between atomistic and CG simulations is that the highly reduced free energy contribution from water entropy, related to hydrogen bonding, is prohibiting the interactions of TRH to be mediated in a realistic way when the solvent is described by the Martini parameters. A fact that further strengthens this explanation is the shape change of the TRH–DLPC PMF when the water model is changed from the spherical water molecules in the standard W model, to the more realistic structure of the PW model. In fact, in his review of the Martini force field [90], Marrink stresses that entropy driven phenomena may not be accurately described with the existing Martini water models.

4.3.2

Structural properties of the bilayer

To characterize the influence of TRH on the lipid bilayer in the CG simulations, the area per lipid and elastic properties of the bilayer were examined. In Fig. 4.12a, areas per lipid (Section 2.5.1) for the different simulations are displayed. The decrease in AL between the atomistic reference simulation (AA-0) and the Martini reference simulation (MP-0) is a reflection of the parameter change between the simulations. Experiments show that AL in pure DLPC bilayers are around 2% larger than the corresponding area for DMPC [191]. The calculated value for DLPC agrees well with the experimental value at T = 303 K, 0.61 nm2 (the simulations are performed at T = 310 K), which could be an indication of problems with the transferability of the Martini parameters between different temperatures [90]. When TRH is added, the AL decrease, which is opposite to the behavior observed in Paper II and recent experiments [176; 183]. The reason behind the decrease is not yet fully understood, and undulation correction factors for the areas (Section 2.5.1) did not change this picture. A clear trend of AL between the different parameter sets could not be observed. The compressibility and bending moduli, KA and kc , were also calculated (Section 2.5.3) for all parameter sets. The increase in KA between AA-0 and MP-0 (Fig. 4.12b) is a direct consequence of the change of lipid type, since DLPC are two carbons shorter per tail. The value of KA in AA-0, 235 ± 71

AL / Å 2

62.5

a

61.5 60.5

KA / mNm-1

59.5 340

270

b 200

H R -T 00 H R -T PP H R -T P2 H

H

TR

TR

P3

P-

M

-0

0 P-

M

A A

Figure 4.12: a) Area per lipid, AL , and b) Area compressibility modulus, KA , calculated from simulations with different force field parameter sets. Black diamonds correspond to simulations without TRH, red circles to low TRH concentration (720 TRH) and blue squares to high TRH concentration (1440 TRH).

34 mNm−1 , is in excellent agreement with the experimental value 234 ± 23 mNm−1 [115]. The MP-0 value 274 ± 5 mNm−1 is low compared to the experimental value for DLPC, 302 mNm−1 [20; 192]. However, a slightly lower value is expected due to the temperature difference between the experiments (303 K) and the simulation (310 K) [193]. The errors were estimated by division of the trajectories into 3-10 parts, each for which KA was determined (Fig. 4.13). The reported errors are the mean standard deviations over all trajectory divisions. The AA-0 simulation is the only simulation that has not converged with respect to the fluctuations in lateral area, which can be observed as an increase in KA when the number of trajectory slices are increased. Overestimation of the compressibilty modulus has previously been reported as an effect of too short simulation time [110], and AA-0 is the shortest of the main simulations (100 ns). No trend in KA related to parameter sets could be observed. Addition of TRH, however, results in an overall increase, which is consistent with our atomistic simulations in Paper II, but ihe increase at high TRH content is far 72

KA / mNm-1

380

a

310

240 KA / mNm-1

420

b

350

280

3

4

5 6 7 8 9 Number of slices

10

Figure 4.13: Area compressibility modulus, KA , calculated from the trajectories divided into slices: a) AA-0 (red hexagons), MP-0 (green stars) b) PP-TRHL (720 TRH, grey hexagons), PP-TRHH (1440 TRH, blue stars).

from as dramatic as in the atomistic case. The calculated kc values are in good agreement with experimental values for both AA-0 and MP-0. The values for CG is rather constant over the simulations, and no change can be observed when TRH is added, which is consistent with experiments [105] and previous simulations [110]. The calculated bending moduli were subsequently used to calculate undulation corrections [110] for AL and KA (Sections 2.5.1 and 2.5.3). The correction factors were exactly 2% for all area calculations and much less for the area compressibility moduli. Looking at snapshots of the CG simulations (such as Fig. 4.9b), we noticed bilayer bending around the TRH clusters. In the last section of Paper IV, we investigated curvature sensing of a single TRH molecule by simulating three degrees of bilayer buckling (flat, medium and high curvature). We concluded that the TRH molecule prefers to interact in regions with high negative curvature (Fig. 4.14). Even though the CG parameters are questionable for describing TRH–bilayer interactions, this result indicates that the approximation of a bicelle as a flat bilayer patch is reasonable.

73

z

x Low

P(x,z)

High

Figure 4.14: Probability density distribution for one single TRH molecule at coordinates (x,z), simulated with a buckled DLPC bilayer (x-axis compressed by a factor of 0.8). The TRH molecule preferably interacts in regions with high negative curvature.

74

5. Conclusions

This thesis has focused on the influence of certain carbohydrate molecules on the structure, order and dynamics of lipid bilayers. Molecular Dynamics (MD) computer simulations have been used to simulate lipid–carbohydrate systems which previously, have been studied by NMR experiments. Molecular Dynamics simulations is a powerful tool to complement NMR experiments with detailed molecular information. In addition, MD can be used to calculate the structural properties which cannot be determined experimentally. In Paper I, mixtures of the non-bilayer forming glyceroglycolipid MGDG and the bilayer forming phospholipid DMPC were investigated. Three bilayer compositions were simulated: 0%, 20% and 45% MGDG, which corresponds to three of the samples previously investigated by NMR [151]. Calculated C–H dipolar couplings in both lipids were in good agreement with the experimental datasets, for which the correct signs could be determined from the simulations. The exception was the GAL group of MGDG, where large errors in the calculated couplings, made an unambiguous sign determination impossible. The cause of the error was ascribed to a slow reorientation around the glycosidic linkage, which may be an effect of inaccurate force field parameters. The intermixture of MGDG into a DMPC bilayer resulted in several changes of the bilayer structure, most of which could be ascribed to an increased lateral stress of the bilayer caused by the tendency of MGDG to form reversed hexagonal phases. Investigations of the lateral distribution of lipid molecules reveal no separation into separate phases or domains, which is consistent with the NMR study in which the corresponding samples lie within the concentration range where the bilayer show lamellar characteristics. In Paper II, the influence of the disaccharide trehalose (TRH) on a DMPC bilayer was examined, with the original aim of validating the use of lipid bicelles as alignment medium in NMR experiments. It was found that TRH molecules bind directly to the bilayer at low to medium concentrations, which give a direct effect on bilayer properties such as area per lipid and thickness. It is clear from the Potential of Mean Force (PMF) calculations that TRH has a region of strong interaction at the bilayer surface. The calculated PMF was used to construct a simple two state binding model that was found to correctly describe the concentration dependent binding of TRH to the bilayer. The model shows that TRH molecules are interacting directly with the bilayer up to 75

a concentration threshold at which the bilayer becomes saturated. After saturation, added TRH molecules are excluded from the surface, seen as a plateau in the isotherm. This behavior has also been observed in experiments [176; 185], and strongly indicates that the hypotheses describing the anhydrobiotic properties of TRH, related mechanisms of direct interaction and sugar exclusion, are not mutually exclusive and can be valid for the same system within different concentration regions. In the TRH–DMPC simulations a substantial aggregation of TRH molecules at the bilayer/water interface was observed. This behavior was ascribed to be the cause of the observed reduction in lipid lateral diffusion and increase in the area compressibility modulus. Such glassy states of TRH have previously been suggested as one hypothesis to explain the ability of TRH to protect membranes during anhydrobiotic stress. In Paper III, it is found that the slow dynamics of the aggregates results in a three orders of magnitude larger orientational order of TRH molecules in the simulations, compared to our NMR experiments. Despite the difference in calculated dipolar couplings (kHz) to the experimental couplings (Hz), fits of the experimental dataset to the simulations could be performed by utilizing the relation between the dipolar couplings, the Saupe matrix of order parameters and internal molecular coordinates. The fits are in good agreement with the experiments, and it was concluded that the Glycam force field is correctly describing the internal molecular structure of TRH and that the conformational space of the simulated TRH molecules is unaffected by the presence of the bilayer. In addition, calculated torsional angles in TRH molecules, simulated in both isotropic solution and in the presence of the DMPC bilayer, were also unaffected by the bilayer interactions, indicating that the assumptions used for the analysis of dipolar couplings in NMR are correct. The slow dynamics of the TRH molecules that are attached to the bilayer imposes a methodological problem related to the ergodicity criterion in Molecular Dynamics simulations. If all of the configurational space of the system is not visited during the timescale of the simulations, the results are not representative of the intended ensemble. In Paper IV, we try to overcome this problem by simulating the TRH–bilayer systems with a coarse grain force field, Martini, to be able to reach a longer timescale and larger systems. Bilayers consisting of DLPC were simulated with TRH at two different concentrations. It turned out that the TRH–bilayer interactions in Martini (both with the standard and polarized water models) give substantially different results than the corresponding atomistic simulations from Paper II. The TRH molecules were clustering to an even greater extent, and calculated PMFs show a much weaker attraction towards the bilayer. The result is that large clusters of TRH are attached to the bilayer, and few TRH molecules are actually interacting directly with the 76

bilayer as opposed to the atomistic simulations. To tune the TRH–bilayer interactions towards the results of the atomistic simulations, some modifications of the TRH parameters were made. We found that even if the PMF corresponding to one of the modifications was close to the atomistic PMF in appearance, the overall behavior did not change when a larger system with many TRH molecules was simulated with the same modification. The conclusion from these results is that the balance of interactions in the coarse grain system is different from the atomistic, and that the PMF is not a good descriptor of the TRH–bilayer interactions. Since the entropic contributions to the total free energy in the Martini water models are limited [90], it is natural to think that further developments of these models are needed to fully describe processes that to a great extent are entropy driven. As a closing of Paper IV, we investigated the curvature sensing of TRH. From simulations of a single TRH molecule with buckled bilayers, we conclude that the TRH molecule preferably interacts in regions with high negative curvature.

5.1

Outlook

The simulation projects within this thesis inject some knowledge into areas of research that are very active. Even if there exist methodological problems in the simulations, as there are in any modelling, the results gained can help design new simulations and better parameters to describe the various dynamical and structural phenomena with more accuracy. Since the glycolipids in the thylokaid membranes in higher plants are closely connected to the processes of photosynthesis [25; 143], a deeper knowledge about these membranes and their phase behavior may contribute to the overall knowledge about photosynthesis as a whole. The glycolipid studied here, MGDG, is one of the main components of this membrane structure, and the tendency of MGDG to form the reversed hexagonal phase seems to be important for the membrane functionality. There are indications of lamellar and inverted hexagonal phase coexistence in the thylokaid membranes [164; 194], regulated by changes in the ratio of MGDG and the other main glycolipid in the Thylokaid membrane, the bilayer forming digalactosyldiacylglycerol (DGDG). Thus, it would be interesting to extend the scope of the MGDG project in this thesis, to perform large scale simulations of MGDG/DGDG mixtures to investigate mechanisms behind phase regulation and membrane properties in such mixtures. This would of course require that the glycolipid parameters are validated, possibly by a change of force field, to be certain that all aspects of the molecules are described correctly. Such investigations could in the long 77

run also provide useful knowledge to investigations of related processes in the cells where phase regulation is of importance for membrane morphology. The mechanisms behind the cryoprotecting properties of the disaccharide TRH on biological membranes are still an area of scientific debate. Even though recent experimental investigations show that TRH affects the lipid bilayers through direct interaction [176; 185], other investigations indicate the opposite [195; 196]. Our simulations show that the coexistence of mechanisms in different concentration regions, first proposed by Andersen et al. [176] and later supported by Konov et al. [185], is probable. However, the uncertainty of the TRH clustering and the slow dynamics that followed it, increase the ambiguity in our conclusions and there are aspects that need to be considered to increase the accuracy of the simulations. Firstly, the sugar/lipid interaction parameters need to be further investigated. As previously pointed out [181], it is difficult to know whether or not the substantial clustering we see in the simulations is physically correct, or merely an artifact of unbalanced interaction parameters. The discrepancy between the atomistic and coarse grain simulations increases the uncertainty. A comparison between different existing force fields and parameter variations would be a reasonable starting point. Secondly, experimental studies of the interaction strength between sugar molecules and lipids are needed to validate simulation results. Free energy calculations show a substantial attraction between TRH molecules and the lipid bilayer, but it is difficult to know if the results are reasonable due to the lack of experimental data to compare with. Possible routes for experimental determination of sugar affinity could be to utilize the Quartz Crystal Microbalance, AFM and other surface analysis techniques to quantify the binding of TRH to a supported bilayer, and monitor changes in structure, mechanics and viscoelastic properties. If we can increase the accuracy of the simulations by force field parameter tuning, extensive comparisons with existing experimental data and conducting simulations and experiments in collaboration, there is a good chance that we can further increase the overall knowledge about membrane–sugar interactions.

78

Sammanfattning

I den här avhandlingen används molekyldynamiksimuleringar (MD) för att karakterisera de specifika förändringar i struktur och dynamik som kan observeras i ett modellmembran under påverkan av interagerande kolhydratmolekyler. Kolhydrater är viktiga för en mängd olika processer i och omkring de komplexa cellmembranen där de både existerar som fria molekyler, men även bundna i form av glykolipider och glykoproteiner. Den i växtriket vanligaste glykolipiden, mono-galaktosyldiacyl-glycerol (MGDG) är en av byggstenarna i tylakoidmembranet. Tylakoidmembranet är i sin tur en del av kloroplasterna, de organeller i växtcellerna där fotsyntesen äger rum, och det har visat sig att MGDG är viktig för deras funktion. Den första artikeln (I) i avhandlingen behandlar en undersökning av membran som utgörs av blandingar av MGDG och fosfolipiden dimyristoylfosfatidylkolin (DMPC). Molekyldynamiksimuleringar användes för att karaktärisera hur membranegenskaperna påverkas vid inblandning av MGDG. Tre system simulerades med 0%, 20% och 45% MGDG vardera. Simuleringarna validerades genom direkta jämförelser mellan beräknade dipolkopplingar och dipolkopplingar från NMR-mätningar (Nuclear Magnetic Resonance, kärnmagnetisk resonans) på liknande system. Beräkningar av parameterar relaterade till membranets struktur och dynamik visade på förändringar vid ökad MGDGkoncentration. Förändringarna anses bero på lateral stress i membranet på grund av den molekylära formen hos MGDG och förändringar i den laterala packningen av molekylerna. I vissa växter och organismer är anrikning av fria sockermolekyler, såsom sackaros och trehalos, känt som en av strategierna för att överleva uttorkning och köld. Att dessa sockerarter har en skyddande effekt på biomembran är sedan länge känt, men mekanismerna på molekylär nivå kan till stor del anses oförklarade. I den andra artikeln (II), undersöks växelverkan mellan trehalos och ett DMPC-membran. Beräkningar av strukturella och dynamiska parametrar användes tillsammans med fria energiberäkningar för att karaktärisera effekten av trehalos på membranegenskaperna. Vi kunde visa att adsorptionen av trehalos på membranytan följer en enkel tvåstegsmodell i enlighet med experimentella studier, och bekräfta vissa av de föreslagna hypoteserna för membran–sockerinteraktioner. I den tredje artikeln (III) valideras simuleringarna med hjälp av dipolkopplingar från vår NMR-undersökning av trelxxix

halos i lösning. Jämförelser mellan beräknade och experimentellt bestämda dipolkopplingar visar att sockerkraftfältet som använts i simuleringarna korrekt beskriver trehalosstrukturen. Beräknade torsionsvinklar i trehalos både från isotropa simuleringar och med ett närvarande DMPC-membran visar att trehalosstrukturen är opåverkad av membranet, vilket till viss del bekräftar de antaganden som görs vid analysen av dipolkopplingar i NMR. Den fjärde artikeln (IV) är ett försök att använda en grovkornig modell för att komma ifrån de problem med långsam dynamik som uppstår när trehalos aggregerar på membranytan. Det visade sig att fortsatt utveckling av kraftfältsparametrarna är nödvändig för att på ett korrekt sätt kunna beskriva interaktionerna mellan trehalos och lipidmembranet. Avslutningsvis undersöktes kurvaturkänsligheten hos trehalos. En trehalosmolekyl simulerades tillsammans med krökta lipidmembran, och det kunde konstateras att trehalos företrädesvis interagerar med membranet i områden med hög negativ kurvatur.

Acknowledgements

First and foremost I would like to thank my supervisor Arnold Maliniak. Without your tremendous support this thesis would never have been written. You believed in me when I first applied for a PhD in Stockholm, and you kept believing. Your door is always open, and I really appreciate the discussions about science and everything else we have had during the years, both in private and at the lunch table. You are a great teacher and supervisor, and the things you have taught me about science and academia will be a great asset in the future, regardless of where I end up. I also want to thank my co-supervisor, Alexander Lyubartsev, who have taught me a lot of useful things about simulations in general, and lipid bilayer simulations in particular. It has been a pleasure to have your knowledge at hand when my simulations did unexpected things. My deepest gratitude goes to Baltzar, with whom I have had the pleasure to work and share office during all my time at the department. I think we have had a great time exchanging knowledge. You taught me everything about data analysis and Fortran, and I kept your computer running and your open terminal windows at a minimum. Göran, Olle and Jakob, thank you for the successful collaboration. I have learned a lot, and it was fun at the same time! Göran, many thanks for your comments on the thesis manuscript, they were really helpful! If you think about carbohydrates and lipids all day, you tend to get hungry. The lunches at the department are always a pleasure. Thanks to Arnold, Baltzar, Jozef, Erik, Amber, Sasha L., Federico and Diana, the discussions at the lunch table never get boring. During my years at the department I have met a lot of nice people and friends that have contributed in one way or another to my work and personal development. Amber, it has been great getting to know you. It was a pleasure sharing teaching and PhD representations with you. Your motivation has really encouraged me to do more, and never give up when things get rough. Erik, thanks for many fruitful discussions about science, work and everything. Martin L, thanks for your input on curvature sensing and lending me your code. The journal club at DBB was a great initiative! Joakim and Sasha M, it was fun going to London with you, and your office has always been the place to go for advice about lipid bilayers, simulations and beer. Zoltan, Sweden is very lxxxi

quiet now, but at least we have public service. Alexander J, you better check the beer cash before I leave. Many thanks to Henrik F, Samrand, Henrik C, Mohsen, Tom, Guido, Alfredo, Cecilia, Erik W, Neda, Christina, Bojan and all of you who have enriched my stay at MMK. I have spent a considerable amount of time as a PhD student representative at the department, but also at the faculty. I would like to thank all people that have engaged in the PhD Council, I really think we made a difference. My thanks also go to the people I have had the pleasure to meet in all the meeting groups where I have been a representative. Special thanks to Gunnar who keep an open mindset and always consider everyone’s opinion. Without my great family and friends I would never have gotten this far. Thank you all for pushing me at the right moments, and taking a step back when needed. *** Mamma, Pappa, Maria, Sofia och Morfar (och Johannes, Ylva, Tora och Leif), tack för att ni alltid finns där! Christina, Bosse, Elise, Anita, Rune, Janne, Reidun, Gunnel och Kjell, det är fantastiskt att ha fått en sån trevlig släkt där man alltid känner sig välkommen! Marcus, Pär och Magnus, titta även jag kan få nåt gjort, det tar bara lite tid! Marcus, tack för korrläsningen! Allt hade tett sig mycket mjöligare utan den. Tobbe, nu är det min tur. Tack för din sista minuten-input, vi får ta en sån där redig pubrunda sen och återkoppla! Jannie och Margit, jag vet inte vad jag hade gjort utan er. Vad som än händer är det alltid skönt att veta att det faktiskt finns saker som är viktigare. Jag älskar er!

Stockholm, 2016-03-21

References

The numbers after each entry denote the pages on which the reference is cited in the thesis.

[1] D. G. S MITH , R. C APPAI , AND K. J. BARNHAM. The redox chemistry of the Alzheimer’s disease amyloid β peptide. Biochim. Biophys. Acta - Biomembr., 1768(8):1976–1990, 2007. 13 [2] A. R AUK. The chemistry of Alzheimer’s disease. Chem. Soc. Rev., 38(9):2698, 2009. 13 [3] J. C ROWE. Anhydrobiosis. Annu. Rev. Physiol., 54(1):579–599, 1992. 13, 57 [4] C. L IU , J. J. G ALLAGHER , K. K. S AKIMOTO , E. M. N ICHOLS , C. J. C HANG , M. C. Y. C HANG , AND P. YANG . Nanowire-bacteria hybrids for unassisted solar carbon dioxide fixation to valueadded chemicals. Nano Lett., 15(5):3634–3639, 2015. 13 [5] M. M. M AY, H.-J. L EWERENZ , D. L ACKNER , F. D IMROTH , AND T. H ANNAPPEL. Efficient direct solar-to-hydrogen conversion by in situ interface transformation of a tandem structure. Nat. Commun., 6:8286, 2015. 13 [6] T. RUNDLÖF, C. L ANDERSJÖ , K. LYCKNERT, A. M ALINIAK , AND G. W IDMALM. NMR investigation of oligosaccharide conformation using dipolar couplings in an aqueous dilute liquid crystalline medium. Magn. Reson. Chem., 36(10):773–776, 1998. 14 [7] M. S TAAF, C. H ÖÖG , B. S TEVENSSON , A. M ALINIAK , AND G. W IDMALM. Conformational investigation of a cyclic enterobacterial common antigen employing NMR spectroscopy and molecular dynamics simulations. Biochemistry, 40(12):3623–3628, 2001. [8] K. LYCKNERT, A. M ALINIAK , AND G. W IDMALM. Analysis of oligosaccharide conformation by NMR spectroscopy utilizing 1 H , 1 H and 1 H , 13 C residual dipolar couplings in a dilute liquid crystalline phase. J. Phys. Chem., 105(A):5119–5122, 2001. [9] B. S TEVENSSON , C. L ANDERSJÖ , G. W IDMALM , AND A. M ALINIAK. Conformational distribution function of a disaccharide in a Liquid Crystalline phase determined using NMR Spectroscopy. J. Am. Chem. Soc., 124(21):5946–5947, 2002. [10] C. L ANDERSJÖ , B. S TEVENSSON , R. E KLUND , J. Ö STERVALL , P. S ÖDERMAN , G. W IDMALM , AND A. M ALINIAK . Molecular conformations of a disaccharide investigated using NMR spectroscopy. J. Biomol. NMR, 35(2):89–101, 2006. 14 [11] J. F. ROBYT. Essentials of Carbohydrate Chemistry. Springer-Verlag, New York, 1998. 15, 16 [12] P. H. R AVEN , R. F. E VERT, AND S. E. E ICHHORN. Biology of Plants. W.H. Freeman and Company, New York, 6 edition, 1999. 17, 18 [13] P. ATKINS AND J. DE PAULA. Physical Chemistry. Oxford University Press, Oxford, 10 edition, 2014. 17, 18

lxxxiii

[14] D. E. VANCE AND J. E. VANCE, editors. Biochemistry of Lipids, Lipoproteins and Membranes. Elsevier B.V., Edmonton, Alberta, Canada, 2008. 17 [15] K. M. S CHMID AND J. B. O HLROGGE. Lipid metabolism in plants. In D. E. VANCE AND J. E. VANCE, editors, Biochem. Lipids, Lipoproteins Membr., chapter 4, pages 97–130. Elsevier B.V., Edmonton, Alberta, Canada, 5 edition, 2008. 18 [16] W. D OWHAN , M. B OGDANOV, AND E. M ILEYKOVSKAYA. Functional roles of lipids in membranes. In D. E. VANCE AND J. E. VANCE, editors, Biochem. Lipids, Lipoproteins Membr., chapter 1, pages 1–37. Elsevier B.V., Edmonton, Alberta, Canada, 5 edition, 2008. 18, 52 [17] C. A. V ERACINI AND D. C ATALANO. Aggregates of amphiphiles in lyotropic liquid crystals. In G. R. L UCKHURST AND C. A. V ERACINI, editors, Mol. Dyn. Liq. Cryst., chapter 20, pages 505–518. Kluwer Academic Publishers, Dordrecht, 1994. 18, 19 [18] J. N. I SRAELACHVILI , D. J. M ITCHELL , AND B. W. N INHAM. Theory of self-assembly of hydrocarbon amphiphiles into micelles and bilayers. J. Chem. Soc. Faraday Trans. 2, 72:1525, 1976. 18 [19] V. V. K UMAR. Complementary molecular shapes and additivity of the packing parameter of lipids. Proc. Natl. Acad. Sci. U. S. A., 88(2):444–448, 1991. 19 ˇ [20] N. K U CERKA , Y. L IU , N. C HU , H. I. P ETRACHE , S. T RISTRAM -NAGLE , AND J. F. NAGLE. Structure of fully hydrated fluid phase DMPC and DLPC lipid bilayers using X-Ray scattering from oriented multilamellar arrays and from unilamellar vesicles. Biophys. J., 88(4):2626– 2637, 2005. 20, 42, 53, 72

[21] R. L. S CHNAAR. Glycolipid-mediated cell–cell recognition in inflammation and nerve regeneration. Archives of Biochemistry and Biophysics, 426(2):163–172, 2004. Highlight issue on glycobiology dedicated to the memory of Victor Ginsburg. 21, 47 [22] P. D ÖRMANN AND C. B ENNING. Galactolipids rule in seed plants. Trends Plant Sci., 7(3):112– 118, 2002. 21, 47 [23] P. J ORDAN , P. F ROMME , H. T. W ITT, O. K LUKAS , W. S AENGER , AND N. K RAUSS. Three-dimensional structure of cyanobacterial photosystem I at 2.5 A resolution. Nature, 411(6840):909–917, 2001. 21, 47 [24] P. JARVIS , P. D ÖRMANN , C. A . P ETO , J. L UTES , C. B ENNING , AND J. C HORY. Galactolipid deficiency and abnormal chloroplast development in the Arabidopsis MGD synthase 1 mutant. Proc. Natl. Acad. Sci. U. S. A., 97(14):8175–8179, 2000. 21 [25] P. D ÖRMANN. Galactolipids in plant membranes. In eLS, pages 1–7. John Wiley & Sons, Ltd, Chichester, UK, 2013. 21, 47, 77 [26] M. L ENTI , C. G ENTILI , A. P IANEZZI , G. M ARCOLONGO , A. L ALLI , R. C ANCEDDA , AND F. D. C ANCEDDA. Monogalactosyldiacylglycerol anti-inflammatory activity on adult articular cartilage. Nat. Prod. Res., 23(8):754–762, 2009. 21 [27] W. B. B ENJAMIN F RANKLIN AND M. FARISH. "Of the stilling of waves by means of oil. Extracted from Sundry Letters between Benjamin Franklin, LL. D. F. R. S. William Brownrigg, M. D. F. R. S. and the Reverend Mr. Farish". Phil. Trans, 64:445–460, 1774. 22 [28] M. E. D ERRICK. Agnes Pockels, 1862-1935. Journal of Chemical Education, 59(12):1030, 1982. 23 [29] L. K. JAMES. Nobel Laureates in Chemistry, 1901-1992. Wiley VCH, 1993. 23 [30] G. B REZESINSKI AND H. M ÖHWALD. Langmuir monolayers to study interactions at model membrane surfaces. Advances in Colloid and Interface Science, 100–102:563–584, 2003. 23

[31] C. P EETLA , A. S TINE , AND V. L ABHASETWAR. Biophysical interactions with model lipid membranes: applications in drug discovery and drug delivery. Molecular Pharmaceutics, 6(5):1264– 1276, 2009. [32] M. E EMAN AND M. D ELEU. From biological membranes to biomimetic model membranes. Biotechnol. Agron. Soc. Environ., 14(4):719–736, 2010. 23 [33] E. A MADO , A. K ERTH , A. B LUME , AND J. K RESSLER. Infrared reflection absorption spectroscopy coupled with Brewster angle microscopy for studying interactions of amphiphilic triblock copolymers with phospholipid monolayers. Langmuir, 24(18):10041–10053, 2008. 23 [34] M. D ELEU , M. PAQUOT, P. JACQUES , P. T HONART, Y. A DRIAENSEN , AND Y. F. D UFRÊNE. Nanometer Scale Organization of Mixed Surfactin/Phosphatidylcholine Monolayers. Biophysical Journal, 77(4):2304–2310, 1999. 23 [35] H. M OTSCHMANN AND H. M ÖHWALD. Langmuir-Blodgett films. In K. H OLMBERG , D. O. S HAH , AND M. J. S CHWUGER, editors, Handbook of Applied Surface and Colloid Chemistry, 2, chapter 5, pages 87–95. Wiley-VCH, 2002. 23 [36] K. E DWARDS AND A. BAEUMNER. Analysis of liposomes. Talanta, 68(5):1432–1441, 2006. 24 [37] N. S TEPANYANTS , G. D. M. J EFFRIES , O. O RWAR , AND A. J ESORKA. Radial sizing of lipid nanotubes using membrane displacement analysis. Nano Lett., 12(3):1372–1378, 2012. 24 [38] M. F OSSATI , B. G OUD , N. B ORGESE , AND J.-B. M ANNEVILLE. An investigation of the effect of membrane curvature on transmembrane-domain dependent protein sorting in lipid bilayers. Cell. Logist., 4:e29087, 2014. 24 [39] R. WATANABE , N. S OGA , T. YAMANAKA , AND H. N OJI. High-throughput formation of lipid bilayer membrane arrays with an asymmetric lipid composition. Sci. Rep., 4:7076, 2014. 24 [40] L. K AM AND S. G. B OXER. Formation of supported lipid bilayer composition arrays by controlled mixing and surface capture. J. Am. Chem. Soc., 122(51):12901–12902, 2000. 24 [41] M.-P. M INGEOT-L ECLERCQ , M. D ELEU , R. B RASSEUR , AND Y. F. D UFRÊNE. Atomic force microscopy of supported lipid bilayers. Nat. Protoc., 3(10):1654–1659, 2008. 24 [42] K. E L K IRAT , S. M ORANDAT, AND Y. F. D UFRÊNE. Nanoscale analysis of supported lipid bilayers using atomic force microscopy. Biochim. Biophys. Acta, 1798(4):750–765, 2010. 24 [43] L. M ÄLER AND A. G RÄSLUND. Artificial membrane models for the study of macromolecular delivery. In M. B ELTING, editor, Methods Mol. Biol. Macromol. Drug Deliv., 480, chapter 9, pages 129–139. Humana Press, New York, 2009. 24 [44] W. Y E , J. L IEBAU , AND L. M ÄLER. New membrane mimetics with galactolipids: lipid properties in fast-tumbling bicelles. J. Phys. Chem. B, 117(4):1044–1050, 2013. 24, 52 [45] P. D. B LOOD AND G. A. VOTH. Direct observation of Bin/amphiphysin/Rvs (BAR) domaininduced membrane curvature by means of molecular dynamics simulations. Proc. Natl. Acad. Sci. U. S. A., 103(41):15068–15072, 2006. 25 [46] W. H UMPHREY, A. DALKE , AND K. S CHULTEN. VMD – Visual Molecular Dynamics. Journal of Molecular Graphics, 14:33–38, 1996. 25, 50, 61, 68 [47] N. M ETROPOLIS , A. W. ROSENBLUTH , M. N. ROSENBLUTH , A. H. T ELLER , AND E. T ELLER. Equation of state calculations by fast computing machines. J. Chem. Phys., 21(6):1087, 1953. 27, 29 [48] B. J. A LDER AND T. E. WAINWRIGHT. Studies in molecular dynamics. General method. J. Chem. Phys., 31(2):459, 1959. 27

[49] J. B. G IBSON , A . N. G OLAND , M. M ILGRAM , AND G. H. V INEYARD. Dynamics of radiation damage. Phys. Rev., 120(4):1229–1253, 1960. 27 [50] A. R AHMAN. Correlations in the motion of atoms in liquid argon. Phys. Rev., 136:A405–A411, 1964. 27 [51] G. E. M OORE. Cramming more components onto integrated circuits. Proc. IEEE, 86(1):82–85, 1965. 27 [52] G. E. M OORE. Progress in digital integrated electronics. IEEE, pages 11–13, 1975. 27 [53] I. T. F OSTER , J. L. T ILSON , A. F. WAGNER , R. L. S HEPARD , R. J. H ARRISON , R. A. K ENDALL , AND R. J. L ITTLEFIELD . Toward high-performance computational chemistry: I. Scalable Fock matrix construction algorithms. J. Comput. Chem., 17(1):109–123, 1996. 29 [54] X. H E , T. Z HU , X. WANG , J. L IU , AND J. Z. H. Z HANG. Fragment quantum mechanical calculation of proteins and its applications. Acc. Chem. Res., 47(9):2748–2757, 2014. 29 [55] A. P. LYUBARTSEV AND A. L. R ABINOVICH. Force field development for lipid membrane simulations. Biochim. Biophys. Acta - Biomembr., 2016. http://dx.doi.org/10.1016/j.bbamem. 2015.12.033. 29 [56] K. A . F ICHTHORN AND W. H. W EINBERG. Theoretical Foundations of Dynamic Monte-Carlo Simulations. J. Chem. Phys., 95(2):1090–1096, 1991. 29 [57] D. F RENKEL AND B. S MIT. Understanding Molecular Simulation: From algorithms to applications. Elsevier Academic Press, San Diego, 2nd edition, 2002. 30 [58] L. V ERLET. Computer "experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Phys. Rev., 159(1):98–103, 1967. 30 [59] R. D. RUTH. A canonical integration technique. IEEE Trans. Nucl. Sci., 30(4):2669–2671, 1983. 30 [60] W. L. J ORGENSEN AND J. T IRADO -R IVES. Potential energy functions for atomic-level simulations of water and organic and biomolecular systems. Proc. Natl. Acad. Sci. U. S. A., 102(19):6665–6670, 2005. 31 [61] W. L. J ORGENSEN , J. C HANDRASEKHAR , J. D. M ADURA , R. W. I MPEY, AND M. L. K LEIN. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys., 79(2):926–935, 1983. 31, 49 [62] H. J. C. B ERENDSEN , J. P. M. P OSTMA , W. F. VAN G UNSTEREN , AND J. H ERMANS. Interaction Models for Water in Relation to Protein Hydration. In B. P ULLMAN, editor, Intermol. Forces, 14, pages 331–342. Springer Netherlands, Dordrecht, 1981. 31, 49 [63] H. J. C. B ERENDSEN , J. R. G RIGERA , AND T. P. S TRAATSMA. The missing term in effective pair potentials. J. Phys. Chem., 91(24):6269–6271, 1987. 31 [64] B. ROUX. Commentary: surface tension of biomembranes. Biophys. J., 71(3):1346–1347, 1996. 32 [65] F. J ÄHNIG. What is the surface tension of a lipid bilayer membrane? Biophys. J., 71(3):1348– 1349, 1996. [66] S. E. F ELLER AND R. W. PASTOR. On simulating lipid bilayers with an applied surface tension: periodic boundary conditions and undulations. Biophys. J., 71(3):1350–1355, 1996. 32 [67] O. B ERGER , O. E DHOLM , AND F. JÄHNIG. Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. Biophys. J., 72(5):2002–2013, 1997. 32, 48, 53

[68] L. D. S CHULER , X. DAURA , AND W. F. VAN G UNSTEREN. An improved GROMOS96 force field for aliphatic hydrocarbons in the condensed phase. J. Comput. Chem., 22(11):1205–1218, 2001. 32, 33 [69] R. S ALOMON -F ERRER , D. A. C ASE , AND R. C. WALKER. An overview of the Amber biomolecular simulation package. Wiley Interdiscip. Rev. Comput. Mol. Sci., 3(2):198–210, 2013. 32 [70] J. P. M. J ÄMBECK AND A. P. LYUBARTSEV. Derivation and systematic validation of a refined all-atom force field for phosphatidylcholine lipids. J. Phys. Chem. B, 116(10):3164–3179, 2012. 32, 60, 64 [71] A. M ACIEJEWSKI , M. PASENKIEWICZ -G IERULA , O. C RAMARIUC , I. VATTULAINEN , AND T. ROG. Refined OPLS all-atom force field for saturated phosphatidylcholine bilayers at full hydration. J. Phys. Chem. B, 118(17):4571–4581, 2014. 32 [72] J. B. K LAUDA , R. M. V ENABLE , J. A. F REITES , J. W. O’C ONNOR , D. J. T OBIAS , C. M ONDRAGON -R AMIREZ , I. VOROBYOV, A. D. M AC K ERELL , AND R. W. PASTOR. Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. J. Phys. Chem. B, 114(23):7830–7843, 2010. 32 [73] S. L EE , A. T RAN , M. A LLSOPP, J. B. L IM , J. H ÉNIN , AND J. B. K LAUDA. CHARMM36 united atom chain model for lipids and surfactants. J. Phys. Chem. B, 118(2):547–556, 2014. [74] C.- J . H ÖGBERG , A. M. N IKITIN , AND A. P. LYUBARTSEV. Modification of the CHARMM force field for DMPC lipid bilayer. J. Comput. Chem., 29(14):2359–2369, 2008. 32 [75] S.-W. C HIU , S. A . PANDIT, H. L. S COTT, AND E. JAKOBSSON. An improved united atom force field for simulation of mixed lipid bilayers. J. Phys. Chem. B, 113(9):2748–2763, 2009. 32 [76] A. K UKOL. Lipid models for united-atom Molecular Dynamics simulations of proteins. J. Chem. Theory Comput., 5(3):615–626, 2009. 59 [77] D. P OGER , W. F. VAN G UNSTEREN , AND A. E. M ARK. A new force field for simulating phosphatidylcholine bilayers. J. Comput. Chem., 31(6):1117–1125, 2010. 32, 59 [78] K. N. K IRSCHNER , A. B. YONGYE , S. M. T SCHAMPEL , J. G ONZÁLEZ -O UTEIRIÑO , C. R. DANIELS , B. L. F OLEY, AND R. J. W OODS. GLYCAM06: A generalizable biomolecular force field. Carbohydrates. J. Comput. Chem., 29(4):622–655, 2008. 32, 60 [79] A. D. M AC K ERELL , E. P. R AMAN , AND O. G UVENCH. CHARMM additive all-atom force field for glycosidic linkages in carbohydrates involving furanoses. J. Phys. Chem. B, 114(40):12981– 12994, 2010. 32 [80] O. G UVENCH , E. R. H ATCHER , R. M. V ENABLE , R. W. PASTOR , AND A. D. M ACKERELL. CHARMM additive all-atom force field for glycosidic linkages between hexopyranoses. J. Chem. Theory Comput., 5(9):2353–2370, 2009. [81] O. G UVENCH , S. S. M ALLAJOSYULA , E. P. R AMAN , E. H ATCHER , K. VANOMMESLAEGHE , T. J. F OSTER , F. W. JAMISON , AND A. D. M AC K ERELL. CHARMM additive all-atom force field for carbohydrate derivatives and its utility in polysaccharide and carbohydrate-protein modeling. J. Chem. Theory Comput., 7(10):3162–3180, 2011. [82] S. S. M ALLAJOSYULA , O. G UVENCH , E. H ATCHER , AND A. D. M AC K ERELL. CHARMM additive all-atom force field for phosphate and sulfate linked to carbohydrates. J. Chem. Theory Comput., 8(2):759–776, 2012. 32 [83] W. DAMM , A. F RONTERA , J. T IRADO -R IVES , AND W. L. J ORGENSEN. OPLS all-atom force field for carbohydrates. J. Comput. Chem., 18(16):1955–1970, 1997. 32, 48

[84] D. KONY, W. DAMM , S. S TOLL , AND W. F. VAN G UNSTEREN. An improved OPLS-AA force field for carbohydrates. J. Comput. Chem., 23(15):1416–1429, 2002. 32 [85] H. S. H ANSEN AND P. H. H ÜNENBERGER. A reoptimized GROMOS force field for hexopyranose-based carbohydrates accounting for the relative free energies of ring conformers, anomers, epimers, hydroxymethyl rotamers, and glycosidic linkage conformers. J. Comput. Chem., 32(6):998–1032, 2011. 32, 59 [86] W. P LAZINSKI , A. L ONARDI , AND P. H. H ÜNENBERGER. Revision of the GROMOS 56A6CARBO force field: Improving the description of ring-conformational equilibria in hexopyranose-based carbohydrates chains. J. Comput. Chem., pages 1–12, 2015. 32 [87] B. L. F OLEY, M. B. T ESSIER , AND R. J. W OODS. Carbohydrate force fields. Wiley Interdiscip. Rev. Comput. Mol. Sci., 2(4):652–697, 2012. 32 [88] M. M. R EIF, P. H. H ÜNENBERGER , AND C. O OSTENBRINK. New interaction parameters for charged amino acid side chains in the GROMOS force field. J. Chem. Theory Comput., 8(10):3705–3723, 2012. 33 [89] S. J. M ARRINK , H. J. R ISSELADA , S. Y EFIMOV, D. P. T IELEMAN , AND A. H. DE V RIES. The MARTINI force field: coarse grained model for biomolecular simulations. J. Phys. Chem. B, 111(27):7812–7824, 2007. 33, 68 [90] S. J. M ARRINK AND D. P. T IELEMAN. Perspective on the Martini model. Chem. Soc. Rev., 42(16):6801–6822, 2013. 33, 71, 77 [91] A. P. LYUBARTSEV AND A. L AAKSONEN. M.DynaMix – a scalable portable parallel MD simulation package for arbitrary molecular mixtures. Comput. Phys. Commun., 128(3):565–589, 2000. 33 [92] J. C. P HILLIPS , R. B RAUN , W. WANG , J. G UMBART, E. TAJKHORSHID , E. V ILLA , C. C HIPOT, R. D. S KEEL , L. K ALÉ , AND K. S CHULTEN. Scalable molecular dynamics with NAMD. J. Comput. Chem., 26(16):1781–1802, 2005. 33 [93] D. VAN D ER S POEL , E. L INDAHL , B. H ESS , G. G ROENHOF, A. E. M ARK , AND H. J. C. B ERENDSEN. GROMACS: Fast, flexible, and free. J. Comput. Chem., 26(16):1701–1718, 2005. 33 [94] B. H ESS , C. K UTZNER , D. VAN DER S POEL , AND E. L INDAHL. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput., 4(3):435–447, 2008. 33 [95] GNU General Public License. http://www.gnu.org/licenses/gpl.html. Accessed: 201602-22. 34 [96] T. DARDEN , D. YORK , AND L. P EDERSEN. Particle mesh Ewald: An Nlog(N) method for Ewald sums in large systems. J. Chem. Phys., 98(12):10089, 1993. 34 [97] R. S PERB , I. G. T IRONI , I. G. T IRONI , R. S PERB , P. E. S MITH , P. E. S MITH , W. F. V. G UN STEREN , AND W. F. V. G UNSTEREN . A generalized reaction field method for molecular dynamics simulations. J. Chem. Phys., 102(13):5451–5459, 1995. 34 [98] H. J. C. B ERENDSEN , J. P. M. P OSTMA , W. F. VAN G UNSTEREN , A. D I N OLA , AND J. R. H AAK. Molecular dynamics with coupling to an external bath. J. Chem. Phys., 81(8):3684, 1984. 35, 36 [99] P. H. H ÜNENBERGER. Thermostat algorithms for Molecular Dynamics simulations. In Adv. Polym. Sci., chapter 173, pages 105–149. Springer-Verlag, Berlin Heidelberg, 2005. 35 [100] S. N OSÉ. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys., 81(1):511–519, 1984. 35

[101] W. H OOVER. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A, 31(3):1695–1697, 1985. 35 [102] G. B USSI , D. D ONADIO , AND M. PARRINELLO. Canonical sampling through velocity rescaling. J. Chem. Phys., 126(1):014101, 2007. 36 [103] M. PARRINELLO AND A. R AHMAN. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys., 52(12):7182, 1981. 36 [104] J. F. NAGLE , R. Z HANG , S. T RISTRAM -NAGLE , W. S UN , H. I. P ETRACHE , AND R. M. S UTER. X-ray structure determination of fully hydrated L alpha phase dipalmitoylphosphatidylcholine bilayers. Biophys. J., 70(3):1419–1431, 1996. 37 [105] N. K UCERKA , S. T RISTRAM -NAGLE , AND J. F. NAGLE. Structure of fully hydrated fluid phase lipid bilayers with monounsaturated chains. J. Membr. Biol., 208(3):193–202, 2005. 37, 53, 73 [106] J. F. NAGLE. Area/lipid of bilayers from NMR. Biophys. J., 64(5):1476–1481, 1993. 37 [107] H. I. P ETRACHE , S. W. D ODD , AND M. F. B ROWN. Area per lipid and acyl length distributions in fluid phosphatidylcholines determined by 2 H NMR spectroscopy. Biophys. J., 79(6):3172– 3192, 2000. 51 [108] A. L EFTIN , T. R. M OLUGU , C. J OB , K. B EYER , AND M. F. B ROWN. Area per lipid and cholesterol interactions in membranes from separated local-field 13 C NMR spectroscopy. Biophys. J., 107(10):2274–2286, 2014. 37 [109] A. TARDIEU , V. L UZZATI , AND F. C. R EMAN. Structure and polymorphism of the hydrocarbon chains of lipids: A study of lecithin-water phases. J. Mol. Biol., 75(4):711–733, 1973. 37 [110] Q. WAHEED AND O. E DHOLM. Undulation contributions to the area compressibility in lipid bilayer simulations. Biophys. J., 97(10):2754–2760, 2009. 37, 41, 43, 72, 73 [111] W. J. A LLEN , J. A. L EMKUL , AND D. R. B EVAN. GridMAT-MD: a grid-based membrane analysis tool for use with molecular dynamics. J. Comput. Chem., 30(12):1952–1958, 2009. 37, 51 [112] W. S HINODA AND S. O KAZAKI. A Voronoi analysis of lipid area fluctuation in a bilayer. J. Chem. Phys., 109(4):1517, 1998. 37 [113] F. B ROCHARD AND J. L ENNON. Frequency spectrum of the flicker phenomenon in erythrocytes. J. Phys., 36(11):1035–1047, 1975. 39 [114] J. FAUCON , M. D. M ITOV, P. M ÉLÉARD , I. B IVAS , AND P. B OTHOREL. Bending elasticity and thermal fluctuations of lipid membranes. Theoretical and experimental requirements. J. Phys., 50(17):2389–2414, 1989. 39 [115] W. R AWICZ , K. C. O LBRICH , T. M C I NTOSH , D. N EEDHAM , AND E. E VANS. Effect of chain length and unsaturation on elasticity of lipid bilayers. Biophys. J., 79(1):328–339, 2000. 39, 72 [116] J. F. NAGLE. Introductory Lecture: Basic quantities in model biomembranes. Faraday Discuss., 161:11–29, 2013. 39 [117] D. C UVELIER , I. D ERÉNYI , P. BASSEREAU , AND P. NASSOY. Coalescence of membrane tethers: experiments, theory, and applications. Biophys. J., 88(4):2714–2726, 2005. 39 [118] W. H ELFRICH. Elastic properties of lipid bilayers: Theory and possible experiments. Zeitschrift fur Naturforsch. - Sect. C J. Biosci., 28(11-12):693–703, 1973. 39 [119] P. B. C ANHAM. The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell. J. Theor. Biol., 26(1):61–81, 1970. 39

[120] E. G. B RANDT, A. R. B RAUN , J. N. S ACHS , J. F. NAGLE , AND O. E DHOLM. Interpretation of fluctuation spectra in lipid bilayer simulations. Biophys. J., 100(9):2104–2111, 2011. 39 [121] G. K HELASHVILI , G. PABST, AND D. H ARRIES. Cholesterol orientation and tilt modulus in DMPC bilayers. J. Phys. Chem. B, 114(22):7524–7534, 2010. 39 [122] G. K HELASHVILI AND D. H ARRIES. How cholesterol tilt modulates the mechanical properties of saturated and unsaturated lipid membranes. J. Phys. Chem. B, 117(8):2411–2421, 2013. [123] G. K HELASHVILI , B. KOLLMITZER , P. H EFTBERGER , G. PABST, AND D. H ARRIES. Calculating the bending modulus for multicomponent lipid membranes in different thermodynamic phases. J. Chem. Theory Comput., 9(9):3866–3871, 2013. 39 [124] S. E. F ELLER AND R. W. PASTOR. Constant surface tension simulations of lipid bilayers: the sensitivity of surface areas and compressibilities. J. Chem. Phys., 111(3), 1999. 40 [125] A. S PAAR AND T. S ALDITT. Short range order of hydrocarbon chains in fluid phospholipid bilayers studied by x-ray diffraction from highly oriented membranes. Biophys. J., 85(3):1576– 1584, 2003. 42 [126] J. S. H UB , T. S ALDITT, M. C. R HEINSTÄDTER , AND B. L. DE G ROOT. Short-range order and collective dynamics of DMPC bilayers: a comparison between molecular dynamics simulations, X-ray, and neutron scattering experiments. Biophys. J., 93(9):3156–3168, 2007. 42 [127] S. T RISTRAM -NAGLE , Y. L IU , J. L EGLEITER , AND J. F. NAGLE. Structure of gel phase DMPC determined by X-Ray diffraction. Biophys. J., 83(6):3324–3335, 2002. 42 [128] D. A XELROD , D. KOPPEL , J. S CHLESSINGER , E. E LSON , AND W. W EBB. Mobility measurement by analysis of fluorescence photobleaching recovery kinetics. Biophys. J., 16(9):1055– 1069, 1976. 42 [129] H. G EERTS , M. D E B RABANDER , R. N UYDENS , S. G EUENS , M. M OEREMANS , J. D E M EY , AND P. H OLLENBECK . Nanovid tracking: a new automatic method for the study of mobility in living cells based on colloidal gold and video microscopy. Biophys. J., 52(5):775–782, 1987. 42 [130] R. N. G HOSH AND W. W. W EBB. Automated detection and tracking of individual and clustered cell surface low density lipoprotein receptor molecules. Biophys. J., 66(5):1301–1318, 1994. 42 [131] A.-L. K UO AND C. G. WADE. Lipid lateral diffusion by pulsed nuclear magnetic resonance. Biochemistry, 18(11):2300–2308, 1979. 42 [132] G. O RÄDD , G. L INDBLOM , AND P. W. W ESTERMAN. Lateral diffusion of cholesterol and dimyristoylphosphatidylcholine in a lipid bilayer measured by pulsed field gradient NMR spectroscopy. Biophys. J., 83(5):2702–2704, 2002. 42, 55 [133] G. P RANAMI AND M. H. L AMM. Estimating error in diffusion coefficients derived from Molecular Dynamics simulations. J. Chem. Theory Comput., 11(10):4586–4592, 2015. 42 [134] H. S UN , S. T IAN , S. Z HOU , Y. L I , D. L I , L. X U , M. S HEN , P. PAN , AND T. H OU. Revealing the favorable dissociation pathway of type II kinase inhibitors via enhanced sampling simulations and two-end-state calculations. Sci. Rep., 5:8457, 2015. 43 [135] D. K ITCHEN , H. D ECORNEZ , J. F URR , AND J. BAJORATH. Docking and scoring in virtual screening for drug discovery: methods and applications. Nat. Rev. Drug Discov., 3(11):935–949, 2004. 43 [136] C. L. W ENNBERG , D. VAN DER S POEL , AND J. S. H UB. Large influence of cholesterol on solute partitioning into lipid membranes. J. Am. Chem. Soc., 134(11):5351–5361, 2012. 43

[137] J. T IAN , A. S ETHI , B. I. S WANSON , B. G OLDSTEIN , AND S. G NANAKARAN. Taste of sugar at the membrane: thermodynamics and kinetics of the interaction of a disaccharide with lipid bilayers. Biophys. J., 104(3):622–632, 2013. 43 [138] G. M. T ORRIE AND J. P. VALLEAU. Nonphysical sampling distributions in Monte Carlo freeenergy estimation: Umbrella sampling. J. Comput. Phys., 23(2):187–199, 1977. 44 [139] S. K UMAR , J. M. ROSENBERG , D. B OUZIDA , R. H. S WENDSEN , AND P. A. KOLLMAN. The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. J. Comput. Chem., 13(8):1011–1021, 1992. 44 [140] A. L AIO AND M. PARRINELLO. Escaping free-energy minima. Proc. Nat. Acad. Sci. USA, 99:12562–12566, 2002. 44 [141] A. BARDUCCI , G. B USSI , AND M. PARRINELLO. Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett., 100(2):020603, 2008. 45 [142] G. A . T RIBELLO , M. B ONOMI , D. B RANDUARDI , C. C AMILLONI , AND G. B USSI. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun., 185(2):604–613, 2014. 45 [143] H. A RONSSON. The galactolipid monogalactosyldiacylglycerol (MGDG) contributes to photosynthesis-related processes in Arabidopsis thaliana. Plant Signal. Behav., 3(12):1093– 1095, 2008. 47, 77 [144] Z.-H. J IANG , A. G EYER , AND R. R. S CHMIDT. The macrolidic glycolipid Calonyctin A, a plant growth regulator: Synthesis, structural assignment, and conformational analysis in micellar solution. Angew. Chemie Int. Ed. English, 34(22):2520–2524, 1995. 47 [145] Z. G ÁLOVÁ AND T. KOŽÁR. Effect of structural variability of the ceramide part on the saccharide-ceramide linkage in model glycolipids studied by Molecular Mechanics and Molecular Dynamics methods. Collect. Czechoslov. Chem. Commun., 61(10):1405–1431, 1996. [146] N. I IDA -TANAKA , K. F UKASE , H. U TSUMI , AND I. I SHIZUKA. Conformational studies on a unique bis-sulfated glycolipid using NMR spectroscopy and molecular dynamics simulations. Eur. J. Biochem., 267(23):6790–6797, 2000. 47 [147] P. R AM , E. K IM , D. S. T HOMSON , K. P. H OWARD , AND J. H. P RESTEGARD. Computer modelling of glycolipids at membrane surfaces. Biophys. J., 63(6):1530–1535, 1992. 47 [148] T. R ÓG , I. VATTULAINEN , AND M. K ARTTUNEN. Modeling glycolipids: take one. Cell. Mol. Biol. Lett., 10(4):625–630, 2005. 47 [149] T. R ÓG , I. VATTULAINEN , A. B UNKER , AND M. K ARTTUNEN. Glycolipid membranes through atomistic simulations: effect of glucose and galactose head groups on lipid bilayer properties. J. Phys. Chem. B, 111(34):10146–10154, 2007. 47 [150] M. T ESSIER , M. D E M ARCO , A. YONGYE , AND R. W OODS. Extension of the GLYCAM06 biomolecular force field to lipids, lipid bilayers and glycolipids. Mol. Simul., 34(4):349–364, 2008. 47, 54 [151] V. C ASTRO , S. V. DVINSKIKH , G. W IDMALM , D. S ANDSTRÖM , AND A. M ALINIAK. NMR studies of membranes composed of glycolipids and phospholipids. Biochim. Biophys. Acta Biomembr., 1768(10):2432–2437, 2007. 47, 50, 52, 54, 64, 65, 75 [152] I. C. P. S MITH. NMR of liquid crystalline lipids in biological membranes. In J. W. E MSLEY, editor, Nuclear Magnetic Resonance of Liquid Crystals, pages 533–566. D. Reidel Publishin Company, Dordrecht, 1985. 48 [153] S. W. C HIU , M. C LARK , V. BALAJI , S. S UBRAMANIAM , H. L. S COTT, AND E. JAKOBSSON. Incorporation of surface tension into molecular dynamics simulation of an interface: a fluid phase lipid bilayer membrane. Biophys. J., 69(4):1230–1245, 1995. 48

[154] J.-P. RYCKAERT AND A. B ELLEMANS. Molecular dynamics of liquid n-butane near its boiling point. Chem. Phys. Lett., 30(1):123–125, 1975. 48 [155] D. A. C ASE , T. A. DARDEN , T. E. C HEATHAM , C. L. S IMMERLING , J. WANG , R. E. D UKE , R. L UO , M. C ROWLEY, R. C. WALKER , W. Z HANG , K. M. M ERZ , B. WANG , S. H AYIK , A. ROITBERG , G. S EABRA , I. KOLOSSVÁRY, K. F. W ONG , F. PAESANI , J. VANICEK , X. W U , S. R. B ROZELL , T. S TEINBRECHER , H. G OHLKE , L. YANG , C. TAN , J. M ONGAN , V. H ORNAK , G. C UI , D. H. M ATHEWS , M. G. S EETIN , C. S AGUI , V. BABIN , AND P. A. KOLLMAN. Amber 11. University of California, San Francisco, 2010. 48 [156] N. C HAKRABARTI , C. N EALE , J. PAYANDEH , E. F. PAI , AND R. P OMÈS. An Iris-like mechanism of pore dilation in the CorA magnesium transport system. Biophys. J., 98(5):784–792, 2010. 48, 60 [157] J. J OYARD AND R. D OUCE. Galactolipid synthesis. In P. K. S TUMPF, editor, The Biochemistry of Plants: Lipids: Structure and Function, 9, chapter 9, pages 215–268. Academic Press Inc. Ltd., London, 1987. 50 [158] C.-J. H ÖGBERG AND A. P. LYUBARTSEV. A Molecular Dynamics investigation of the influence of hydration and temperature on structural and dynamical properties of a dimyristoylphosphatidylcholine bilayer. J. Phys. Chem. B, 110(29):14326–14336, 2006. 50, 52, 53, 55 [159] K. H RISTOVA AND S. H. W HITE. Determination of the hydrocarbon core structure of fluid dioleoylphosphocholine (DOPC) bilayers by x-ray diffraction using specific bromination of the double-bonds: effect of hydration. Biophys. J., 74(5):2419–2433, 1998. 50 [160] S. C. C OSTIGAN , P. J. B OOTH , AND R. H. T EMPLER. Estimations of lipid bilayer geometry in fluid lamellar phases. Biochim. Biophys. Acta, 1468(1-2):41–54, 2000. 51, 53 [161] J. F. NAGLE AND S. T RISTRAM -NAGLE. Structure of lipid bilayers. Biochim. Biophys. Acta Rev. Biomembr., 1469(3):159–195, 2000. 51 [162] J. W OHLERT AND O. E DHOLM. Dynamics in atomistic simulations of phospholipid membranes: Nuclear magnetic resonance relaxation rates and lateral diffusion. J. Chem. Phys., 125(20):204703, 2006. 52 [163] C. B OTTIER , J. G ÉAN , F. A RTZNER , B. D ESBAT, M. P ÉZOLET, A . R ENAULT, D. M ARION , AND V. V IÉ . Galactosyl headgroup interactions control the molecular packing of wheat lipids in Langmuir films and in hydrated liquid-crystalline mesophases. Biochim. Biophys. Acta, 1768(6):1526–1540, 2007. 52 [164] B. D EME , C. C ATAYE , M. A . B LOCK , E. M ARECHAL , AND J. J OUHET. Contribution of galactoglycerolipids to the 3-dimensional architecture of thylakoids. FASEB J., 28(8):3373–3383, 2014. 52, 77 [165] J. S EELIG AND N. WAESPE -S ARCEVIC. Molecular order in cis and trans unsaturated phospholipid bilayers. Biochemistry, 17(16):3310–3315, 1978. 54 [166] P. F. A LMEIDA , W. L. VAZ , AND T. E. T HOMPSON. Lateral diffusion in the liquid phases of dimyristoylphosphatidylcholine/cholesterol lipid bilayers: a free volume analysis. Biochemistry, 31(29):6739–6747, 1992. 55 [167] J. B. K LAUDA , B. R. B ROOKS , AND R. W. PASTOR. Dynamical motions of lipids and a finite size effect in simulations of bilayers. J. Chem. Phys., 125(14):144710, 2006. 55 [168] G. L INDBLOM , G. O RÄDD , L. R ILFORS , AND S. M OREIN. Regulation of lipid composition in Acholeplasma laidlawii and Escherichia coli membranes: NMR studies of lipid lateral diffusion at different growth temperatures. Biochemistry, 41(38):11512–11515, 2002. 55

[169] E. L INDAHL AND O. E DHOLM. Molecular dynamics simulation of NMR relaxation rates and slow dynamics in lipid bilayers. J. Chem. Phys., 115:4938–4950, 2001. 55 [170] C. A NÉZO , A. H. DE V RIES , H.-D. H ÖLTJE , D. P. T IELEMAN , AND S.-J. M ARRINK. Methodological issues in lipid bilayer simulations. J. Phys. Chem. B, 107(35):9424–9433, 2003. 55 [171] J. C ROWE , L. C ROWE , AND D. C HAPMAN. Preservation of membranes in anhydrobiotic organisms: the role of trehalose. Science, 223(4637):701–703, 1984. 57 [172] J. H. C ROWE , L. M. C ROWE , J. F. C ARPENTER , AND C. AURELL W ISTROM. Stabilization of dry phospholipid bilayers and proteins by sugars. Biochem. J., 242(1):1–10, 1987. [173] M. C AFFREY, V. F ONSECA , AND A. C. L EOPOLD. Lipid-sugar interactions: relevance to anhydrous biology. Plant Physiol., 86(3):754–758, 1988. [174] K. L. KOSTER , M. S. W EBB , G. B RYANT, AND D. V. LYNCH. Interactions between soluble sugars and POPC (1-palmitoyl-2-oleoylphosphatidylcholine) during dehydration: vitrification of sugars alters the phase behavior of the phospholipid. Biochim. Biophys. Acta - Biomembr., 1193(1):143–150, 1994. 57 [175] J. H. C ROWE. Anhydrobiosis: An Unsolved Problem with Applications in Human Welfare, chapter Anhydrobio, pages 263–280. Springer International Publishing, Cham, 2015. 57 [176] H. D. A NDERSEN , C. WANG , L. A RLETH , G. H. P ETERS , AND P. W ESTH. Reconciliation of opposing views on membrane–sugar interactions. Proc. Natl. Acad. Sci., 108(5):1874–1878, 2011. 57, 59, 61, 63, 71, 76, 78 [177] M. K ARPLUS AND D. H. A NDERSON. Valence-bond interpretation of electron-coupled nuclear spin interactions; application to methane. J. Chem. Phys., 30(1):6, 1959. 58 [178] M. K ARPLUS. Vicinal proton coupling in nuclear magnetic resonance. J. Am. Chem. Soc., 85(18):2870–2871, 1963. 58 [179] M. H. L EVITT. Spin Dynamics, chapter 20. John Wiley & Sons Ltd, Chichester, 2008. 58 [180] C. O OSTENBRINK , A. V ILLA , A. E. M ARK , AND W. F. VAN G UNSTEREN. A biomolecular force field based on the free enthalpy of hydration and solvation: The GROMOS force-field parameter sets 53A5 and 53A6. J. Comput. Chem., 25(13):1656–1676, 2004. 59 [181] L. R. W INTHER , J. Q VIST, AND B. H ALLE. Hydration and Mobility of Trehalose in Aqueous Solution. J. Phys. Chem. B, 116(30):9196–9207, 2012. 60, 78 [182] M. K. BAKER AND C. F. A BRAMS. Dynamics of lipids, cholesterol, and transmembrane αhelices from microsecond Molecular Dynamics simulations. J. Phys. Chem. B, 118(47):13590– 13600, 2014. 60 [183] K. B. KONOV, N. P. I SAEV, AND S. A. D ZUBA. Low-temperature molecular motions in lipid bilayers in the presence of sugars: insights into cryoprotective mechanisms. J. Phys. Chem. B, 118(43):12478–12485, 2014. 61, 71 [184] S. O. J ÓNSDÓTTIR , S. A. C OOKE , AND E. A. M ACEDO. Modeling and measurements of solidliquid and vapor-liquid equilibria of polyols and carbohydrates in aqueous solution. Carbohydr. Res., 337(17):1563–1571, 2002. 62 [185] K. B. KONOV, D. V. L EONOV, N. P. I SAEV, K. Y. F EDOTOV, V. K. VORONKOVA , AND S. A . D ZUBA. Membrane–sugar interactions probed by pulsed electron paramagnetic resonance of spin labels. J. Phys. Chem. B, 119(32):10261–10266, 2015. 63, 76, 78

[186] G. A. C LARK , J. M. H ENDERSON , C. H EFFERN , B. A KGÜN , J. M AJEWSKI , AND K. Y. C. L EE. Synergistic interactions of sugars/polyols and monovalent salts with phospholipids depend upon sugar/polyol complexity and anion identity. Langmuir, 31(46):12688–12698, 2015. 63 [187] K. L. KOSTER , Y. P. L EI , M. A NDERSON , S. M ARTIN , AND G. B RYANT. Effects of vitrified and nonvitrified sugars on phosphatidylcholine fluid-to-gel phase transitions. Biophys. J., 78(4):1932–1946, 2000. 63 [188] A. B OTAN , F. FAVELA -ROSALES , P. F. J. F UCHS , M. JAVANAINEN , M. K ANDU Cˇ , W. K ULIG , A. L AMBERG , C. L OISON , A. LYUBARTSEV, M. S. M IETTINEN , L. M ONTICELLI , J. M ÄÄTTÄ , O. H. S. O LLILA , M. R ETEGAN , T. R ÓG , H. S ANTUZ , AND J. T YNKKYNEN. Toward atomistic resolution structure of phosphatidylcholine headgroup and glycerol backbone at different ambient conditions. J. Phys. Chem. B, 119(49):15075–15088, 2015. 64 [189] R. Y. D ONG. Nuclear Magnetic Resonance of Liquid Crystals. Springer-Verlag, New York, 1994. 66 [190] C. A. L ÓPEZ , A. J. R ZEPIELA , A. H. DE V RIES , L. D IJKHUIZEN , P. H. H ÜNENBERGER , AND S. J. M ARRINK. Martini coarse-grained force field: extension to carbohydrates. J. Chem. Theory Comput., 5(12):3195–3210, 2009. 68 ˇ [191] N. K U CERKA , M.-P. N IEH , AND J. K ATSARAS. Fluid phase lipid areas and bilayer thicknesses of commonly used phosphatidylcholines as a function of temperature. Biochim. Biophys. Acta Biomembr., 1808(11):2761–2771, 2011. 71

[192] J. C. M ATHAI , S. T RISTRAM -NAGLE , J. F. NAGLE , AND M. L. Z EIDEL. Structural determinants of water permeability through the lipid membrane. J. Gen. Physiol., 131(1):69–76, 2008. 72 ˇ [193] J. PAN , S. T RISTRAM -NAGLE , N. K U CERKA , AND J. F. NAGLE. Temperature dependence of structure, bending rigidity, and bilayer interactions of dioleoylphosphatidylcholine bilayers. Biophys. J., 94(1):117–124, 2008. 72

[194] D. J. M URPHY. The importance of non-planar bilayer regions in photosynthetic membranes and their stabilisation by galactolipids. FEBS Lett., 150(1):19–26, 1982. 77 [195] B.-A. B RÜNING , S. P RÉVOST, R. S TEHLE , R. S TEITZ , P. FALUS , B. FARAGO , AND T. H ELL WEG . Bilayer undulation dynamics in unilamellar phospholipid vesicles: Effect of temperature, cholesterol and trehalose. Biochim. Biophys. Acta - Biomembr., 1838(10):2412–2419, 2014. 78 [196] B. K ENT, T. H UNT, T. A. DARWISH , T. H AUSS , C. J. G ARVEY, AND G. B RYANT. Localization of trehalose in partially hydrated DOPC bilayers: insights into cryoprotective mechanisms. J. R. Soc. Interface, 11:20140069, 2014. 78