<%BANNER%>

Molecular Modeling of Biomembrane Deformations--The Role of Lipids


PAGE 1

MOLECULAR MODELING OF BIOMEMBRANE DEF ORMA TIONS|THE R OLE OF LIPIDS By ERIC R. MA Y A DISSER T A TION PRESENTED TO THE GRADUA TE SCHOOL OF THE UNIVERSITY OF FLORID A IN P AR TIAL FULFILLMENT OF THE REQUIREMENTS F OR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORID A 2006

PAGE 2

Cop yrigh t 2006 b y Eric R. Ma y

PAGE 3

I dedicate this w ork to paren ts Da vid and Linda Ma y who instilled in me the imp ortance of education and encouraged me to undertak e science and engineering. I also dedicate this w ork to m y brother T o dd Ma y who has sho wn me what can b e ac hiev ed through p ersistence and hard w ork. Lastly I dedicate this to m y ancee Jill T o disco who has lo v ed and supp orted me throughout m y academic career and made sacrices so that I could ac hiev e m y goals, for whic h I will alw a ys b e grateful.

PAGE 4

A CKNO WLEDGMENTS I tak e this opp ortunit y thank m y men tor, Dr. A tul Narang, for prop osing this problem to me and in tro ducing me to the eld of biological mo deling. I greatly appreciate the guidance he pro vided me in m y researc h, and also his advice and wisdom on life's other matters. I also w ould lik e to thank Dr. Dmitry Kop elevic h for all the time he has sp en t w orking closely with me o v er the last sev eral y ears. Additionally I w ould also lik e to thank the other mem b ers of m y committee, Dr. Ranganathan Nara y anan and Dr. Gerry Sha w, for their advice and a v ailabilit y I w ould lik e to thank Dr. Karthik Subramanian, Dr. Shakti Gupta, Jason No el, V ed Sharma, Ashish Gupta, Gunjan Mohan and V alere Chen for their assistance, suggestions and the enjo y able scien tic discussions w e'v e had. I w ould also lik e to thank Shirley Kelly and Debbie Sando v al for their assistance throughout m y graduate studies. iv

PAGE 5

T ABLE OF CONTENTS page A CKNO WLEDGMENTS . . . . . . . . . . . . . . iv LIST OF FIGURES . . . . . . . . . . . . . . . . vii ABSTRA CT . . . . . . . . . . . . . . . . . . ix 1 INTR ODUCTION . . . . . . . . . . . . . . . 1 1.1 Sp ecic Aims . . . . . . . . . . . . . . 1 1.2 Bac kground . . . . . . . . . . . . . . . 3 1.2.1 Role of Lipids in Biomem brane Deformations . . . . 3 1.2.2 Molecular Dynamics . . . . . . . . . . . 5 1.3 Broader Impact . . . . . . . . . . . . . . 9 2 PHASE TRANSITIONS IN MIXED LIPID SYSTEMS: A MD STUD Y 11 2.1 In tro duction . . . . . . . . . . . . . . . 11 2.2 Metho ds . . . . . . . . . . . . . . . . 13 2.2.1 Molecular Mo del . . . . . . . . . . . . 13 2.2.2 Sim ulation Details . . . . . . . . . . . 16 2.3 Results . . . . . . . . . . . . . . . . 16 2.3.1 LP A Micelles . . . . . . . . . . . . . 17 2.3.2 Phase Beha vior of Pure DOP A System . . . . . 19 2.3.3 Phase Beha vior of Mixed Lipid Systems . . . . . 20 2.4 Discussion . . . . . . . . . . . . . . . 24 3 DETERMINA TION OF ELASTIC PR OPER TIES F OR BILA YERS . 31 3.1 In tro duction . . . . . . . . . . . . . . . 31 3.2 Bending Mo dulus . . . . . . . . . . . . . 32 3.2.1 Curv ature Elasticit y . . . . . . . . . . . 32 3.2.2 Molecular Mo del . . . . . . . . . . . . 35 3.2.3 Sim ulation Details . . . . . . . . . . . 36 3.2.4 Results . . . . . . . . . . . . . . 37 3.3 Line T ension . . . . . . . . . . . . . . . 40 3.3.1 Pressure Calculation from Molecular Dynamics Sim ulations 40 3.3.2 Description of P oten tials and F orces . . . . . . 42 3.3.3 Calculation of Virial . . . . . . . . . . . 45 3.3.4 Distribution of Virial . . . . . . . . . . 46 3.3.5 Calculation of Line T ension . . . . . . . . . 49 v

PAGE 6

3.3.6 Sim ulations Details . . . . . . . . . . . 51 3.3.7 Results . . . . . . . . . . . . . . 52 3.4 Tilt Deformations . . . . . . . . . . . . . 56 3.4.1 Energy Mo del for Tilt . . . . . . . . . . 56 3.4.2 Nature of the H I I Phase and Evidence of Tilt . . . . 59 4 CONCLUSIONS . . . . . . . . . . . . . . . 63 4.1 Ma jor Conclusions . . . . . . . . . . . . . 63 4.2 F uture Directions . . . . . . . . . . . . . 64 4.2.1 Additional Lipid Sp ecies . . . . . . . . . 64 4.2.2 Con tin uum Scale Mo deling . . . . . . . . . 65 REFERENCES . . . . . . . . . . . . . . . . . 67 BIOGRAPHICAL SKETCH . . . . . . . . . . . . . . 72 vi

PAGE 7

LIST OF FIGURES Figure page 1{1 Mem brane diagrams . . . . . . . . . . . . . 2 1{2 Coat protein assem bly and v esicle formation . . . . . . . 4 1{3 P oten tial functions . . . . . . . . . . . . . . 6 1{4 Comparision of atomistic and CG structures . . . . . . . 9 2{1 Detailed atomic structures and the corresp onding CG mo dels . . 16 2{2 Dep endence of the size and shap e of LP A micelles on [Mg 2+ ] . . 18 2{3 Dep endence of the micelle asphericit y on the Mg 2+ concen tration. . 19 2{4 Bila y er conguration of pure DOP A systems . . . . . . . 20 2{5 T emp erature induced phase transitions . . . . . . . . 21 2{6 Sim ulations of temp erature{induced phase transitions . . . . 23 2{7 Molecular geometry c haracterization for pure LP A systems . . . 25 2{8 Molecular geometry c haracterization in mixed lipid systems . . . 27 2{9 Prop osed transition states in v esicle fusion . . . . . . . 28 2{10 Stalk formation in DOPE/LP A system . . . . . . . . 29 3{1 Comparision of theoretical and exp erimen tal v esicle shap es . . . 33 3{2 A tomisitic and CG represen tation of phosphoinositides . . . . 35 3{3 Depictions of PI4P . . . . . . . . . . . . . . 36 3{4 Angle probabilit y distributions . . . . . . . . . . . 38 3{5 Sp ectral in tensit y for PI4P/DPPC systems . . . . . . . 39 3{6 Bending mo dulus for PI4P/DPPC bila y er . . . . . . . . 39 3{7 Denition of r ij . . . . . . . . . . . . . . . 43 3{8 Angle Orien tation . . . . . . . . . . . . . . 45 3{9 Slab orien tation for calculation of bila y er surface tension . . . 47 vii

PAGE 8

3{10 Phases within a bila y er . . . . . . . . . . . . 50 3{11 P atc h system with repulsiv e w alls . . . . . . . . . . 52 3{12 Densit y of pressure of CG systems . . . . . . . . . . 54 3{13 Virial and pressure prole comparision . . . . . . . . 55 3{14 Bond p oten tial comparison . . . . . . . . . . . . 55 3{15 Densit y and pressure proles for patc h systems . . . . . . 57 3{16 Tilt and b ending deformations . . . . . . . . . . . 58 3{17 Construction of H I I phase . . . . . . . . . . . . 59 3{18 Examination of H I I phase from MD . . . . . . . . . 61 3{19 Tilt in a H I I phase . . . . . . . . . . . . . . 62 4{1 Prop osed CG mo dels of additional lipids: a) PI(4,5) 2 b) D A G. . . 65 4{2 Domains within a mem brane, tak en from ref [ 1 ] . . . . . . 66 viii

PAGE 9

Abstract of Dissertation Presen ted to the Graduate Sc ho ol of the Univ ersit y of Florida in P artial F ulllmen t of the Requiremen ts for the Degree of Do ctor of Philosoph y MOLECULAR MODELING OF BIOMEMBRANE DEF ORMA TIONS|THE R OLE OF LIPIDS By Eric R. Ma y Decem b er 2006 Chair: A tul Narang Ma jor Departmen t: Chemical Engineering Cellular transp ort pro cesses suc h as endo cytosis and exo cytosis in v olv e the budding of v esicles from the donor mem brane (ssion) and their in tegration in to the mem brane of the acceptor compartmen t (fusion). Both ssion and fusion are energetically unfa v orable, since they in v olv e the formation of highly curv ed nonbila y er in termediates. Consequen tly these pro cesses do not o ccur sp on taneously; they are under the strict con trol of sp ecic proteins. In the classical paradigm of ssion, mem brane deformation w as though t to b e driv en en tirely b y in teractions b et w een proteins. Evidence is no w accum ulating that these proteins do not act alone. They op erate in concert with sp ecic mem brane lipids, suc h as phosphoinositides (PI), whic h lo calize in the region of deformation. It is of in terest to understand the exten t to whic h the lo calized phosphoinositides c hange the mec hanical prop erties of biomem branes. In this w ork, a coarse grain molecular mo del is used to conduct molecular dynamics (MD) sim ulations of mixed lipid systems. The initial part of the w ork fo cuses on the v erication of the molecular mo del b y comparing sim ulation data to exp erimen tal phase transition data for mixed lipid systems con taining dioleo yl ix

PAGE 10

phosphatidyle ethanolamine (DOPE), dioleo yl phosphatidic acid (DOP A) and lysophosphatidic acid (LP A). Also, the molecular geometries of the previously men tioned lipids are c haracterized and related the the mec hanisms driving the phase transitions. Additionally MD sim ulations and analyses are p erformed on another mixed lipid system consisting of dipalmito yl phosphatidyl c holine (DPPC) and phosphatidylinositol-4-phosphate (PI4P). In this study t w o elastic parameters of the mem brane are measured from the sim ulations, namely the b ending mo dulus and the co ecien t of line tension b et w een mem brane domains of dieren t comp osition. x

PAGE 11

CHAPTER 1 INTR ODUCTION 1.1 Sp ecic Aims The goal of this researc h is to b etter understand the mec hanisms underlying the deformations of biological mem branes. Biological mem branes are v ery complex structures, consisting of a v ariet y of phospholipid sp ecies, c holesterol, and proteins as sho wn in Figure 1{1 a [ 2 ]. The mem brane is an essen tial comp onen t of the cell. It serv es to enclose the inner con ten ts, while also p erforming a n um b er of other functions. The mem brane selectiv ely allo ws passage of certain ions in to and out of the cell to main tain c hemical gradien ts, whic h are essen tial to the cell's surviv al. Also, the mem brane is em b edded with receptors whic h receiv e c hemical signals from the en vironmen t instructing the cell to p erform sp ecic op erations. Biomem branes, as w ell as b eing comp ositionally complex, are also structurally dynamic, c hanging shap e to p erform a v ariet y of tasks. Euk ary otic cells are capable of endo cytosis and exo cytosis, the pro cess b y whic h material is brough t in to or shipp ed out from the cell, resp ectiv ely or passed b et w een dieren t compartmen ts of the cell. These pro cesses in v olv e the formation of transp ort v esicles, whic h are small pac k ages with an outer lipid mem brane to enclose the con ten ts whic h are transp orted. A sc hematic of the budding and ssion pro cess is sho wn in Figure 1{1 b. Inside of the cell there are sev eral orangelles, some of whic h are enclosed b y a mem brane as w ell, suc h as the Golgi b o dy mito c hondria, endoplasmic reticulum, and n ucleus. These organelles also undergo mem brane shap e transformations. Lipid molecules mak e up the ma jorit y of these mem branes, whic h con tain man y dieren t lipid sp ecies. During these shap e c hanging pro cesses, certain lipid sp ecies will lo calize in and around the region of deformation. 1

PAGE 12

2 (a) (b) Figure 1{1: Mem brane diagrams: a) Example of a biological mem brane [ 2 ]. b) Depiction of budding and ssion of a transp ort v esicle. The h yp othesis of this study is that the alteration of the c hemical comp osition of a biomem brane will lead to a c hange in the mec hanical prop erties of the membrane. This c hange in the mec hanical prop ert y could then lead to a sp on taneous deformation of the mem brane. The sp ecic questions this study will try to answ er are the follo wing: 1. Are certain elastic prop erties aected b y a c hange in the c hemical comp osition of the mem brane? 2. What is the molecular mec hanism driving the alteration of an elastic propert y? 3. Do es a c hange in c hemical comp osition of the mem brane lead to a congurational c hange of the mem brane? With these questions answ ered, it will b e the ultimate goal to relate the c hanges in mec hanical prop erties to the bio c hemistry o ccurring inside the cell when these deformations o ccur. In essence, it is desired to dev elop a complete mo del starting from c hemical reactions inside the bio c hemical path w a ys eecting comp ositional

PAGE 13

3 c hanges inside the mem brane, eecting mec hanical c hanges of the mem brane, leading to a con tin uum lev el mo del whic h can predict global mem brane shap e. 1.2 Bac kground 1.2.1 Role of Lipids in Biomem brane Deformations Lipid molecules are the main comp onen t of biomem branes. These lipids are amphiphilic in nature, ha ving one part whic h has an anit y for h ydrophilic en vironmen ts (head) and one part whic h has an anit y for h ydrophobic en vironmen ts (tail). The tail(s) consists mainly of saturated h ydro carb ons, but unsaturated b onds do commonly exist as w ell. There can b e either a single tail or t w o tails connected to the head. The head and tail are connected through a glycerol group. The head group often con tains a phosphate group, as w ell as other p olar groups, whic h often carry a net c harge. In biological systems the en vironmen t in whic h the cells liv es is aqueous, as is the in terior of the cell. This naturally leads to the bila y er conguration of the lipid molecules to separate the in ternal con ten ts of the cell from the en vironmen t, as sho wn in Figure 1{1 a. In this conguration, the heads are submerged in b oth the in terior of the cell and the exterior en vironmen ts while the tails form a h ydrophobic region at the in terior of the mem brane itself. This in terior region of the mem brane is v ery imp ortan t b ecause it pro vides a region for proteins to reside since they often con tain h ydrophobic regions. The mec hanism b y whic h biomem branes c hange shap e is the fo cus of this w ork. It had originally b een p ostulated that the driving force b ehind the formation of v esicles from biomem branes w as driv en b y protein in teractions. It is b eliev ed that in teractions b et w een mem brane b ound proteins and external proteins kno wn as coat proteins, esp ecially the protein complex clathrin, act to increase curv ature of the mem brane and lead to bud formation [ 3 4 ]. This is an incremen tal pro cess in whic h the addition of more coat proteins con tin uously deform the mem brane un til a bud is formed and it pinc hes o from the original mem brane. The coat then

PAGE 14

4 Figure 1{2: Depiction of coat protein assem bly and the formation of a transp ort v esicle. [ 2 ]. will disso ciate from the v esicle to allo w for ssion with the target destination, see Figure 1{2 The pro cess of fusion and ssion of v esicles from mem branes is an energetically unfa v orable one. The in termediate states can b e highly curv ed, requiring more energy to form than a v ailable from t ypical thermal ructuations of the system [ 5 ]. Ho w ev er, evidence has b een accum ulating that the proteins do not act alone, but act in concert with the lipid molecules themselv es [ 6 ]. Some of the k ey lipids are ones whic h exist in relativ ely lo w concen trations inside the mem brane, but pla y a crucial role in these morphological c hanges. Sp ecically phosphoinositides, phosphatidic acid, lysophosphatidic acid and diacylglycerol ha v e b een implicated as essen tial lipids in the pro cesses of fusion and ssion of transp ort v esicles [ 7 { 9 ]. T o this p oin t, it has b een unclear as to what function these lipids pla y in altering the conguration of biomem branes. It is the goal of this w ork to study the eect dieren t lipids ha v e on bila y ers in the absence of proteins.

PAGE 15

5 1.2.2 Molecular Dynamics Molecular dynamics (MD) is a computational to ol used to study molecular systems b oth in equilibrium and non-equilibrium regimes. In MD, all particles are treated classically where the particles in teract through dieren t b onded and nonb onded p oten tials. The b onded in teractions can consist of b ond length, b ond angle and dihedral p oten tials, whic h are t ypically mo delled b y a harmonic p oten tial. The non-b onded in teractions are due to v an der W aals and electrostatic forces. F or computational purp oses, only pair p oten tials are t ypically considered, taking 3-b o dy and higher term in teractions b ecomes v ery time consuming b ecause they in v olv e summing o v er triplets (or greater) of molecules [ 10 ]. Ho w ev er, the 2-b o dy p oten tials ha v e sho wn to pro duce v ery go o d results in comparison to exp erimen tal measuremen ts. T o determine the particles p osition and v elo cities, Newtons equation's of motion m ust b e solv ed m i d 2 r i dt 2 = F i ; i = 1 :::N (1.1) where i refers to the particle n um b er, running from 1 to N r is the p osition v ector, m is the mass of the particle and F is the force v ector acting on the particle. The force is giv en b y the negativ e deriv ativ e of the p oten tial with resp ect to p osition, F i = dV d r i The algorithm t ypically used to solv e these equations is the V erlet Algorithm [ 10 ]. v ( t + 4 t 2 ) = v ( t 4 t 2 ) + F ( t ) m 4 t (1.2) r ( t + 4 t ) = r ( t ) + v ( t + 4 t 2 ) 4 t (1.3) It can b e seen Eq. 1.2 up dates the v elo cit y at time t + 4 t 2 from the v elo cit y at time t 4 t 2 and the force at time t. The p osition is then up dated, as sho wn in Eq. 1.3 using the v elo cit y at time t + 4 t 2 and the p osition at time t to get the p osition at time t + 4 t

PAGE 16

6 (a) (b) Figure 1{3: P oten tial functions: a) Lennard-Jones p oten tial b) Electrostatic p otential The non b onded forces are often sp ecied b y a Lennard-Jones p oten tial of the form V LJ ( r ij ) = 4 (( =r ij ) 12 ( =r ij ) 6 ) as sho wn in Figure 1{3 a. The parameter giv es the depth of the p oten tial w ell and is the hard sphere radius. Since the p oten tial rattens out at large r a cut-o radius is sp ecied to increase the computational eciency If t w o particles are separated b y a distance greater than the cut-o radius, it is assumed they ha v e no inruence on eac h other. Sp ecial care m ust b e tak en when using cut-o radii to a v oid discon tin uities in the p oten tial functions. T o a v oid this, a shifting function can b e added to the p oten tial in order for it to smo othly approac h a zero slop e at the cut-o radius. In addition to the Lennard-Jones p oten tial, an electrostatic p oten tial m ust b e accoun ted for if there are c harged particles in the system. The form of the electrostatic p oten tial is V E S ( r ij ) = q i q j 4 0 r ij where q i is the c harge on molecule i and 0 is the p ermittivit y of free space, and is sho wn in Figure 1{3 b. Again, a cut-o radius can b e c hosen for the electrostatic p oten tial or a more complex metho d of accoun ting for the long range in teractions kno wn as Ew ald summation ma y b e emplo y ed. In molecular dynamic sim ulations p erio dic b oundary conditions are often used. This is to eliminate the edge eects whic h w ould b e in tro duced b y the presence of w alls. When p erio dic conditions are used, it allo ws molecules to exit

PAGE 17

7 from the sim ulation b o x on one side and then reapp ear, en tering on the opp osite side of the b o x. The molecules will in teract with the nearest p erio dic image of the other molecules; this is kno wn as the minim um image con v en tion. In order for the p erio dic sim ulation to represen t the macroscopic system the length of the sim ulation b o x should b e at least 6 where is the Lennard-Jones parameter represen ting the hard sphere radius of the particles [ 10 ]. It is often desired for the sim ulation system to main tain a constan t temp erature and pressure, so that the results are comparable to exp erimen ts. One metho d to con trol temp erature is through the Berendsen algorithm, whic h couples the system to a heat bath using rst-order kinetics. In this w eak coupling sc heme the system resp onds exp onen tially when it deviates from the bath temp erature, T 0 giv en b y Eq. 1.4 where is the coupling time scale. This is ac hiev ed b y rescaling the v elo cities of all particles b y a factor giv en in Eq. 1.5 T is related to b y Eq. 1.6 where C V is the total heat capacit y of the system, k is the Boltzmann constan t and N d f is the total n um b er of degrees of freedom [ 11 ]. dT dt = T 0 T (1.4) = 1 + 4 t T ( T 0 T 1) (1.5) = 2 C V T = N d f k (1.6) The w eak coupling metho d of temp erature coupling is easily implemen ted, but it do es not represen t a correct statistical mec hanical ensem ble. An alternativ e metho d whic h do es prob e a correct canonical ensem ble, is the extended system approac h, whic h in tro duces a \T emp erature Piston" to con trol the temp erature, thereb y in tro ducing and additional degree of freedom to the system. The equations of motion b ecome mo died as sho wn in Eq. 1.7 where is the heat bath v ariable or friction parameter. This parameter is dynamic and is go v erned b y it's o wn equation of motion Eq. 1.8 Q is dened b y Eq. 1.9 where T is coupling parameter

PAGE 18

8 whic h denes the timescale of oscillation of kinetic energy b et w een the heat bath and the system [ 12 ]. d 2 r i dt 2 = F i m i d r i dt (1.7) d dt = 1 Q ( T T 0 ) (1.8) Q = T T 0 4 2 (1.9) Similar coupling sc hemes exist for pressure coupling as w ell. In pressure coupling, the b o x v ectors are scaled to adjust the v olume and therefore the pressure of the system. This can b een ac hiev ed through w eak coupling where the pressure ob eys rst-order kinetics or through an extend system approac h, where the equations of motion of the particles are mo died and the b o x v ectors m ust ob ey their o wn equation of motion. Pressure coupling can b e done either isotropically or anisotropically whic h allo ws for deformation of the sim ulation b o x from the initial asp ect ratios. When dealing with lipid bila y ers, the c hoice of anisotropic coupling to equal pressure in all directions leads to a zero surface tension of the bila y er, whic h is often desirable. The study of lipid bila y ers using MD, is w ell established, where n umerous studies ha v e sho wn go o d agreemen t with exp erimen tal data regarding the densit y of these systems, diusion co ecien ts and p ermeabilit y [ 13 ]. The molecular mo dels used in MD can b e atomistic, in whic h all molecules, b onds, b ond angles and partial c harges are sp ecied (p ossibly excluding h ydrogens) or coarse grained in whic h m ultiple atoms are dened b y a single in teraction site. There are metho ds of deriving coarse grain p oten tials from, atomistic p oten tials one of whic h is the Iterativ e Boltzmann In v ersion [ ? ]. In this metho d an initial p oten tial is guessed using a kno wn radial distribution function, g ( r ) ; suc h that V 0 ( r ) = k T l n ( g ( r )). Sim ulating the system with this p oten tial will then pro duce a dieren t radial distribution, g 0 ( r ) and the p oten tial is corrected b y adding the term

PAGE 19

9 (a) (b) Figure 1{4: Comparison of atomistic and coarse grain structures a) DPPC b) W ater k T l n ( g 0 ( r ) g ( r ) ) to the initial p oten tial guess. This pro cess can b e iterated sev eral times un til reasonable con v ergence of the radial distribution function is ac hiev ed. T o date sev eral coarse grain mo dels ha v e b een published, including one b y Marrink et al. whic h includes mo dels for sev eral lipid sp ecies [ ? ]. In Marrink's mo del, a coarse grain particle t ypically represen ts 4-6 atoms; sho wn in Figure 1{4 are depictions of the atomistic and coarse grain represen tations of DPPC, a common biological lipid, and w ater. The lab elling of the coarse grain particles corresp onds the the force eld designation of the particles published b y Marrink [ ? ]. The main adv an tage of using a coarse grain mo del is that it reduces the degrees of freedom in the system and also allo ws for a larger time step, whic h mak e this mo del m uc h more computationally ecien t than an atomistic mo del. 1.3 Broader Impact The foundation of the w ork is the adv ancemen t of scien tic kno wledge. The w ork has b een conducted without a direct application in tended. Ho w ev er, it is hop eful that the ndings of this researc h could lead to more directed studies. Biology is a eld whic h often oers qualitativ e and sometimes sp eculativ e explanations of the phenomena. T aking an engineering approac h to biological problems oers

PAGE 20

10 the p ossibilit y of quan titativ e and mec hanistic explanations. The motiv ation b ehind this researc h is understanding the mec hanical implications of observ ed lipid lo calization in biomem branes. F urthermore, the goal of this w ork is to understand ho w c hanges in c hemical comp osition of mem branes aect elastic prop erties and what are the underlying molecular mec hanisms driving these c hanges in the membrane prop erties. These eects can then b e compared with the forces generated inside the cell including pressure, osmotic pressure, p olymerization of lamen ts and protein in teractions to see if c hemical c hanges within the mem brane are pla ying a signican t role in the deformation of mem branes.

PAGE 21

CHAPTER 2 PHASE TRANSITIONS IN MIXED DOPE-DOP A AND LP A-DOP A LIPID SYSTEMS: A MOLECULAR D YNAMICS STUD Y 2.1 In tro duction Cellular transp ort pro cesses suc h as endo cytosis, exo cytosis, and sub cellular trac king in v olv e the budding of v esicles from the donor mem brane (ssion) and their in tegration in to the mem brane of the acceptor compartmen t (fusion). Both ssion and fusion are energetically unfa v orable, since they in v olv e the formation of highly curv ed non-bila y er in termediates. Consequen tly these pro cesses do not o ccur sp on taneously|they are under the strict con trol of sp ecic proteins. In the classical paradigm of ssion, mem brane deformation w as though t to b e driv en en tirely b y in teractions b et w een proteins. According to this mo del, ssion is initiated b y the recruitmen t of certain cytosolic proteins called c o at pr oteins to the mem brane. The step wise assem bly of coat proteins then incremen tally deforms the mem brane in to a spherical shap e, th us pro ducing a v esicle whic h is ultimately released along with its encased coat. Evidence is no w accum ulating that these proteins do not act alone. They op erate in concert with particular mem brane lipids, namely phosphoinositides, diacylglycerol, and phosphatidic acid. One function of the lipids is to recruit proteins in v olv ed in mem brane deformation. It has b een sho wn, for instance, that phospholipase D, an enzyme required for v esicle formation in the Golgi complex, is recruited to the mem brane b y diacylglycerol. More imp ortan tly ho w ev er, lipids app ear to b e directly in v olv ed in deformation of biomem branes. Early exp erimen ts sho w ed that enzymes that acylate lysophosphatidic acid (LP A) to phosphatidic acid (P A) strongly activ ate v esicle formation. Since LP A and P A are oneand 11

PAGE 22

12 t w o-tailed phospholipids, it seems plausible to h yp othesize that LP A is coneshap ed ( 4 ) and P A has the shap e of an in v erted cone ( r ). No w, it is kno wn from condensed-matter ph ysics that amphiphilic molecules self-assem ble in to aggregates whose morphology is in timately link ed to the shap e of the individual molecules [ 14 ]. Sp ecically r and 4 -shap ed molecules form structures with p ositiv e and negativ e curv atures b ecause these geometries reduce b ending stresses within the mem brane. This has led to the h yp othesis that the shap e c hange induced b y LP A-acyltransferases is ultimately driv en b y the distinct shap es of LP A and P A molecules. Recen tly Ko oijman et al. p erformed exp erimen ts in a cell-free system to ascertain the molecular shap e of LP A and P A [ 15 ]. T o this end, systems of LP A, P A, and mixtures of these lipids with v arious other lipids w ere studied in a v ariet y of en vironmen ts in whic h pH, temp erature and div alen t cation concen trations w ere v aried to mimic the conditions in the cytosol and Golgi complex. In the mixed lipid exp erimen ts, LP A or P A w as mixed with either dioleo yl phosphatidylethanolamine (DOPE) or dioleo yl phosphatidylethanolamine (DEPE), and again the pH, temp erature, and ionic conditions w ere altered. Using 31 P-NMR measuremen ts the system phase w as distinguished for a giv en set of conditions. It w as observ ed that a system of a giv en lipid comp osition w ould transition to a dieren t phase b y altering en vironmen tal conditions or temp erature. These phase transitions in v olv e a c hange in curv ature of the of the equilibrium structure, often from a relativ ely rat lamellar phase to a high curv ed in v erted hexagonal phase. LP A and P A are of particular in terest b ecause in the Golgi b o dy the acylation of LP A to form P A has b een sho wn to induce the formation of a v esicle from the Golgi mem brane [ 7 ]. LP A and P A are v ery dieren t from eac h other from a geometrical p ersp ectiv e. The sp on taneous curv ature of LP A w as determined to b e highly p ositiv e, desiring

PAGE 23

13 to form micelles, while P A has a negativ e sp on taneous curv ature, under ph ysiological conditions [ 16 ]. It is b eliev ed that this transformation from LP A to P A induces stress in the mem brane leading to the destabilization of the bila y er and the formation of a bud. The goal of this w ork is to dev elop a mo del for the lipids LP A and P A to b e used in molecular dynamics studies. A coarse grain (CG) mo del dev elop ed b y Marrink et al. is capable of repro ducing exp erimen tally observ ed phase transitions for a phospholipid system consisting of dioleo yl phosphatidylc holine (DOPC) and dioleo yl phosphatidylethanolamine (DOPE) [ ? ]. Using Marrink's force eld parameters, mo dels for P A and LP A w ere constructed and sim ulations w ere conducted and compared to exp erimen tal results of Ko oijman et al. [ 15 ]. In addition, it is desired to c haracterize the molecular geometry of these molecules to gain insigh t in to the mec hanism driving the phase transitions. This w ork is a preliminary study to v erify CG mo dels for the lipids LP A and P A. With the assurance that the mo dels are reliable, further study is planned to fully in v estigate the mec hanism driving mem brane deformation in mixed lipid systems. 2.2 Metho ds 2.2.1 Molecular Mo del In order to accurately ev aluate the elastic (i.e., macroscopic) prop erties using microscopic sim ulations, it is necessary to p erform sim ulations of systems con taining thousands of lipids for h undreds or thousands of nanoseconds. The application of a detailed atomistic mo del to suc h length and time scales is extremely timeconsuming ev en with mo dern computing systems and therefore limits the scop e of the systems that can b e in v estigated within a reasonable time frame. Therefore, sim ulations of lipid systems often emplo y less detailed mo dels, suc h as the Bro wnian dynamics sim ulations [ ? ], dissipativ e particle dynamics [ ? ], and coarse-grained molecular dynamics (CGMD) mo dels [ 17 ? ? ? 18 ]. In this w ork, w e use a CGMD

PAGE 24

14 mo del, whic h allo ws one to p erform sim ulations on near atomic length scales, while reducing the computational time b y orders of magnitude and therefore allo wing exploration of the systems on sucien tly large timeand length-scales. The CGMD mo dels appro ximate small groups of atoms b y a single united atom (b ead). Sev eral suc h mo dels ha v e b een in tro duced and applied to sim ulations of v arious complex molecular systems [ 17 19 ? ? ? ? 20 ]. Dev elopmen t of coarse-grained molecular mo dels is a sub ject of activ e ongoing researc h [ ? ? ? ]. In this w ork, w e will use the coarse-grained mo del prop osed b y Marrink et al. [ ? ]. This mo del has b een sho wn to yield go o d agreemen t with exp erimen ts and atomistically detailed sim ulations for suc h quan tities as densit y and elasticit y of pure lipid bila y ers. Within this mo del, for example, 4 meth yl groups of lipid tails are represen ted b y a single h ydrophobic coarse-grained b ead and 4 w ater molecules are represen ted b y a single spherical p olar b ead. This mo del also pro vides sev eral t yp es of b eads that can b e used to mo del v arious groups of atoms (suc h as glycerol and phosphate) within lipid headgroups. The in teractions b et w een non-b onded b eads are mo deled b y the Lennard-Jones p oten tial, whereas the in teractions of b onded b eads are mo deled b y the harmonic b ond and angle vibration p oten tials. Also, electrostatic in teractions b et w een c harged b eads are tak en in to accoun t. The parameters for these p oten tials as w ell as mo dels for sev eral molecules w ere previously published b y Marrink et al. [ ? ] and the same parameters w ere used in this w ork. F or the non b onded in teractions a cut o distance of 1.2 nm w as used. T o a v oid discon tin uities of the p oten tials at the cut o p oin t the LJ p oten tial is smo othly shifted to zero b et w een 0.9 nm and the cuto p oin t. The electrostatic p oten tial is also shifted to zero b eginning at 0.0 nm, to represen t the eect of ionic screening. In this w ork, w e will use the coarse-grained mo del for w ater, metal ions, and DOPE lipids prop osed in b y Marrink et al. [ ? ]. In addition, w e ha v e dev elop ed

PAGE 25

15 coarse-grained mo dels of the additiv e lipids using the metho dology of Marrink et al. [ ? ] for mapping dieren t functional groups within a molecule to coarsegrained b eads with sp ecic prop erties. The detailed atomistic structures and the corresp onding course-grained represen tations of some of these lipid molecules are sho wn in Fig. 2{1 The structure and CG mo del for the lipid LP A sho wn in Fig. 2{1 a consists of 7 CG particles. Eac h b ead is classied with resp ect to the system published b y Marrink et al. [ ? ]. The glycerol-ester link age is mo delled b y a single b ead in LP A with the classication Nda, indicating a non-p olar particle capable of b oth h ydrogen b onding and donating. Glycerol in double tailed molecules is t ypically mo delled with t w o b eads, ho w ev er for LP A the c hoice to use only a single b ead mo del for glycerol w as based up on a closer represen tation of the true mass of the molecule and the qualitativ e agreemen t of the sim ulations with exp erimen ts using this mo del. An attempt to mo del the glycerol of LP A with 2 b eads, as in DOP A, for a total of 8 b eads in the mo del, pro duced a sim ulation result whic h did not agree with the exp erimen ts. In the absence of ions the selfassem bled structure of the 8 b ead LP A formed a w orm{lik e micelle, where the exp erimen ts predict spherical micelles. The structure and mo del for dioleo yl phosphatidic acid (DOP A) is sho wn in Fig. 2{1 b. The CG mo del consists of 13 CG particles. The glycerol ester link age is mo delled in the t ypical w a y using t w o b eads of t yp e Na, indicating it is non p olar and can act as a h ydrogen acceptor. LP A has the abilit y to b e a h ydrogen donor due to the h ydro xyl group on the glycerol, where in DOP A the o xygen is b ound to the ester link age. Both of these lipids are mo delled with oleo yl fatt y acid tails (18 carb ons, 1 double b ond), the double b ond causing the kink in the c hain. Lastly a mo del for a div alen t cation is in tro duced using a single b ead of t yp e Qda, carrying a c harge of 1.4 e. The c harge is reduced from 2.0 e to accoun t for the h ydration shell surrounding the ion.

PAGE 26

16 (a) (b) Figure 2{1: Detailed atomic structures and the corresp onding coarse-grained mo dels of considered lipids: (a) LP A and (b) DOP A. Mapping b et w een dieren t groups of atoms and the coarse-grained b eads is sho wn for some groups. The t yp es of the lipids (C = h ydrophobic, P = p olar, etc.) and their p oten tial parameters are the same as in Marrink et al. [ ? ] 2.2.2 Sim ulation Details Molecular dynamic (MD) sim ulations w ere conducted, using the GR OMA CS MD sim ulations pac k age [ ? ], v ersion 3.2.1. In all sim ulations anisotropic Berendsen pressure coupling w as used with a reference pressure of 1.0 bar. All systems w ere also coupled to a temp erature bath using Berendsen coupling. The time step used for the sim ulations w as 40 fs. P erio dic b oundary conditions w ere imp osed in all three directions. 2.3 Results In this study molecular dynamic sim ulations w ere used to in v estigate the p olymorphic phase b eha vior of lipid systems. Exp erimen tal results ha v e indicated that lipid systems are sensitiv e to pH, temp erature, ion concen tration and relativ e w ater h ydration lev els [ 15 21 22 ]. W e use the coarse-grained (CG) molecular mo del dev elop ed b y Marrink et al. and discussed ab o v e. Recall that in this mo del four w ater molecules are represen ted b y a single CG particle. Therefore, in what follo ws when referring to the amoun t of w ater in the sim ulation system the amoun t

PAGE 27

17 of \real" w ater molecules will b e giv en. Also, it is kno wn that the dynamics of coarse-grained mo dels is t ypically faster than that of a real system. This is due to the smo othing of the p oten tials. F or the sp ecic mo del considered here, the CG time is 4 times faster than the real time. Therefore, the times rep orted for the sim ulations are giv en as the eectiv e time, whic h is calculated b y scaling the observ ed sim ulation time b y a factor of 4. 2.3.1 LP A Micelles Ko oijman et al. studied the shap e c hange of LP A aggregates in resp onse to increasing concen trations of Mg 2+ (Mg 2+ /LP A = 0, 0.2, and 1.0) [ 15 ]. The exp erimen ts w ere p erformed at 37and a pH of 7.2. They observ ed that in the absence of Mg 2+ LP A forms micellar aggregates. Ho w ev er, as the Mg 2+ /LP A ratio is increased, the shap e of the LP A aggregates c hanges from a micellar conguration to a system of bila y er disks. T o test the abilit y of the CG mo del to repro duce these observ ations, w e conducted self-assem bly sim ulations of the LP A lipids (sho wn in Fig. 2{1 a), w ater and dieren t concen trations of div alen t ions (represen ting Mg 2+ ). The sim ulations w ere p erformed at 37. The sim ulation b o x w as randomly p opulated with 500 molecules of LP A, 192 molecules of w ater p er LP A molecule, and the appropriate amoun t of Mg 2+ (determined b y the Mg 2+ /LP A ratio that w as b eing studied). Initial random disp ersions of LP A and Mg 2+ in w ater quic kly self-assem bled in to micelles and equilibrated within 800 ns. Figures 2{2 a{d sho w the equilibrium congurations ac hiev ed at the end of the sim ulation for the three Mg 2+ /LP A concen tration ratios used in the exp erimen ts. The sim ulations sho w that in the absence of Mg 2+ six small micelles are formed (Fig. 2{2 a). When the Mg 2+ /LP A ratio is increased to 0.2, v e micelles are formed (Fig. 2{2 b). If the Mg 2+ /LP A ratio is increased ev en further to 1.0, the larger micelles aggregate to form a single aggregate. Figures 2{2 c,d sho w the top and

PAGE 28

18 (a) (b) (c) (d) Figure 2{2: Dep endence of the size and shap e of LP A micelles on [Mg 2+ ]: (a) Molar ratio [Mg 2+ ]/[LP A] = 0. (b) Molar Ratio [Mg 2+ ]/[LP A] = 0.2. (c) Molar ratio [Mg 2+ ]/[LP A] = 1.0, side view (d) Molar ratio [Mg 2+ ]/[LP A] = 1.0, top view. In these plots the w ater molecules and ions w ere remo v ed for clarit y side views of this single aggregate. It is clear from these gures that the aggregate has a disk-lik e conguration. Therefore, as the concen tration of the ions increases, the micellar size increases and micelles acquire disk-lik e shap es, whic h is consisten t with exp erimen tal 31 P-NMR measuremen ts of Ko oijman et al. [ 15 ]. T o quan tify the shap e c hange of the LP A aggregates sho wn in Fig. 2{2 w e calculated their asphericit y whic h is dened as A = P 3i
PAGE 29

19 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Mg+2/LPAAsphericity Figure 2{3: Dep endence of the micelle asphericit y on the Mg 2+ concen tration. 80 ns of the sim ulation. Fig. 2{3 sho ws that as the concen tration of Mg 2+ increases, so do es the asphericit y of the lipid aggregates. It is of in terest to understand the molecular basis for the shap e c hanges of the LP A aggregates sho wn in Fig. 2{2 According to the shap e-structure concept of lipid p olymorphism, the shap e of an aggregate is determined largely b y the shap e of the comp onen t molecules. Therefore, w e exp ect that the c hange from spherical to disk-lik e micelles indicates a c hange in the molecular geometry of LP A lipids. 2.3.2 Phase Beha vior of Pure DOP A System Ko oijman et al. also studied the phase b eha vior of pure DOP A systems. It w as observ ed that a DOP A system formed a lamellar phase at neutral pH and a temp erature of 310 K b oth in the absence and equimolar presence of Ca 2+ Selfassem bly sim ulations w ere conducted to matc h this result. 500 DOP A molecules and 17500 CG w aters w ere randomly placed in a sim ulations b o x, for the system with Ca 2+ 500 of CG div alen t cations w ere also added to the system. The systems w ere sim ulated at 310 K for 800ns and in b oth cases a bila y er phase w as observ ed. This indicates the mo del for DOP A and div alen t cation can matc h the exp erimen tal observ ations.

PAGE 30

20 (a) (b) Figure 2{4: Bila y er conguration of pure DOP A system observ ed after 800ns of sim ulation. a) No div alen t ions. b) Div alen t ion concen tration equal to lipid concen tration. 2.3.3 Phase Beha vior of DOPE-LP A and DOPE-DOP A Mixed System W e also in v estigated phase transitions b et w een the lamellar and in v erted hexagonal phases in the mixed DOPE-LP A and DOPE-DOP A systems. The exp erimen tal data [ 15 ] summarized in Fig. 2{5 indicate phase transitions from lamellar to in v erted hexagonal phases within the bila y er systems as the temp erature of the system increases. Ko oijman et al. [ 15 ] observ ed that when the temp erature is increased, pure DOPE undergo es a transition at T 275 K from the lamellar ( L ) phase to the in v erted hexagonal ( H I I ) phase (curv e lab eled 4 in Fig. 2{5 ). Up on addition of 10 mol% LP A to DOPE, the phase transition temp erature increases to T 300 K ( r curv e of Fig. 2{5 ). The addition of 10 mol% DOP A to DOPE results in sligh t increase in the phase transition temp erature to T 280 K ( curv e of Fig. 2{5 ). Therefore, the exp erimen tal data sho w that the c hange of the minor mixture comp onen t from LP A to DOP A has a destabilizing eect on the lamellar phase b y lo w ering the phase transition temp erature b y 20 K. W e p erform sim ulations to c hec k if this transition w as mirrored b y the mo del. Recen tly Marrink et al. [ ? ] studied the phase b eha vior of the CG mo del for pure DOPE. They found, for a h ydration lev el in the range of 9-12 w aters p er lipid,

PAGE 31

21 240 260 280 300 320 340 0 10 20 30 40 50 60 70 80 90 100 % BilayerTemperatuer(K) Pure DOPE (Exp) DOPE/LPA (Exp) DOPE/DOPA (Exp) Pure DOPE (Sim) DOPE/LPA (Sim) DOPE/DOPA (Sim) Figure 2{5: T emp erature-induced phase transition b et w een lamellar and in v erted hexagonal phases in the mixed lipid systems: Comparison of exp erimen tal data [ 15 ] and sim ulations for DOPE/DOP A ( ) DOPE/LP A systems ( r ) and pure DOPE ( 4 ). The exp erimen tal data are sho wn b y op en sym b ols and the sim ulations results are sho wn b y closed sym b ols. The lamellar phase corresp onds to 100% bila y er and the in v erted hexagonal phase corresp onds to 0%. The comp osition of the mixed systems is 90% DOPE and 10% of a minor lipid (LP A or DOP A). The sim ulation data for mixed systems are from the curren t w ork, where as the sim ulation data of the pure DOPE system w as tak en from Marrink et al. [ ? ]. the phase transition temp erature to b e ab o v e 287 K ( N curv e of Fig. 2{5 ). This is quite close to the exp erimen tally observ ed phase transition temp erature. T o determine the phase transition, temp erature sim ulations w ere p erformed at sev eral temp eratures in whic h the initial conguration of the system consisted of four pre-sim ulated bila y ers stac k ed on top of eac h other. In some of the sim ulations, the bila y ers w ould form stalks b et w een the la y ers, whic h w ould then gro w and the ( H I I ) phase w ould sp on taneously form. W e p erformed similar sim ulations to c hec k if the CG mo del could capture the phase b eha vior of the DOPE-LP A system. The initial conditions for our sim ulations w ere tak en to b e lamellar phases prepared follo wing the sim ulation proto col of Marrink et al. [ ? ]. First, w e used self-assem bly of initially randomly

PAGE 32

22 disp ersed lipids in w ater to prepare a lipid bila y er. The self-assem bly sim ulations w ere p erformed in a system with relativ ely lo w lipid concen tration whic h leads to bila y er formation within 200 ns, the system w as allo w ed to run for and additional 600ns to allo w for equilibration. The self assem bly sim ulation w as conducted at a w ater to lipid ratio of 80 and a temp erature of 280 K. W e then deh ydrate the bila y er do wn to a w ater lipid ratio of 15 and stac k 2 iden tical deh ydrated bila y ers on top of eac h other. In an attempt to promote formation of the H I I phase, the bila y ers w ere giv en a p erturbation so they w ould resem ble a sine w a v e in one direction. The direction of the p erturbation w as alternated b et w een the bila y ers so that the p eak of one bila y er w ould coincide with the v alley of the neigh b oring bila y er. This conguration w as then energy minimized using the steep est descen t metho d in GR OMA CS for 50000 steps. This energy minimized conguration w as then copied in the normal direction to obtain a system consisting of 4 bila y ers. This initial conguration is sho wn in Fig. 2{6 a. Eac h bila y er consists of 450 DOPE molecules and 50 LP A molecules. The sim ulations w ere p erformed at v arious temp eratures b et w een 250 K and 300 K. Examples of nal conguration of the sim ulations for the DOPE-LP A system are sho wn in Fig. 2{6 b and c. A t temp eratures at or b elo w 265 K the lamellar phase remains stable on the timescale of sim ulations (4.0 s), see Fig. 2{6 b. On the other hand, at T = 275 K and ab o v e, the same system exhibits an instabilit y whic h leads to formation of a stalk and then ultimately the in v erted hexagonal phase is reac hed. The nal conguration for the sim ulation at 300 K is sho wn in Fig. 2{6 c. All of the congurations sho wn in Fig. 2{6 w ere generated b y cop ying the sim ulation cell 3 times in the lateral direction to b etter illustrate the geometry of the systems.

PAGE 33

23 (a) (b) (c) Figure 2{6: Sim ulations of temp erature-induced phase transition b et w een lamellar and in v erted hexagonal phases in the mixed DOPE-LP A lipid systems: (a) Initial p erturb ed conguration for DOPE-LP A sim ulations. (b) Final conguration of DOPE-LP A system at 250 K (lamellar phase). (c) Final conguration of DOPELP A system at 300 K Lipid headgroup b eads are sho wn b y blac k spheres, tail b eads are sho wn b y gra y spheres, and the w ater b eads are sho wn b y white spheres. Figure 2{5 sho ws a comparison b et w een the exp erimen tal system and the sim ulation results. The mo del system displa ys a phase transition temp erature ab o v e 265 K whereas the exp erimen tal system transitioned around 300 K. F or the DOPE-DOP A sim ulations, a bila y er consisting of 120 lipid molecules (10% DOP A, 90% DOPE) w as sim ulated b y self-assem bly at w ater to lipid ratio of 120 and at 240 K for 80ns. These small bila y ers w ere then copied laterally to form a bila y er consisting of 480 lipids. A p erturb ed bila y er system w as then created in the same manner as for the DOPE-LP A system, again at a reduced w ater to lipid ratio of 15. This p erturb ed system w as then sim ulated for up 4.0 s at v arying temp eratures b et w een 240 K and 310 K. A t a temp erature of 265 K and greater the systems sp on taneously formed the H I I phase. A t temp eratures of 240 K and 250 K the system remained stable in the bila y er phase, but at temp erature at or ab o v e 265 K the systems con v erted in to the H I I phase. Figure 2{5 sho ws that the sim ulations are in qualitativ e agreemen t with the exp erimen tal data. The mo del DOPE-DOP A system undergo es a phase transition at a temp erature appro ximately 15 K lo w er than the phase-transition

PAGE 34

24 temp erature for DOPE-LP A system, compared to an appro ximate 20 K shift for the exp erimen tal system. Ho w ev er, the precise v alues of the phase transition temp erature predicted b y the sim ulations are dieren t from those observ ed in the exp erimen ts. The exp erimen ts indicate phase transition temp eratures of 280 K and 300 K for the DOPE-DOP A and DOPE-LP A systems, but the corresp onding sim ulated v alues transitions temp eratures are ab o v e 250 K and 265 K, resp ectiv ely As can b e seen, the curren t coarse-grained mo del do es not capture the exact phase transition temp erature, it do es capture imp ortan t qualitativ e features of the systems, suc h as the destabilizing eect on the lamellar phase due to the c hange from LP A to DOP A. 2.4 Discussion W e also studied the molecular mec hanism in v olv ed in temp erature-dep enden t phase transitions. It is eviden t that as the temp erature increases, so do the bila y er ructuations, th us facilitating the p ossibilit y of a phase c hange. Ho w ev er, for the hexagonal phase to b e stable, it is necessary that the molecular geometry also c hange. T o in v estigate the molecular shap e c hanges in temp erature-dep enden t phase transitions, the sim ulations of the DOPE-LP A system w ere analyzed. Molecular pac king theory states the microstructure of the lipid aggregate should b e represen tativ e of the mean molecular geometry of the system [ 14 ]. In order to v erify this h yp othesis, the pac king parameter w as measured, P = v =l c a 0 ; where v is the v olume o ccupied b y the h ydro carb on tails, l c is the length of the tails in the direction normal to the w ater lipid in terface and a 0 is the area o ccupied b y the head group. F or the pure LP A system, a transition b et w een sev eral spherical micelles to a single bila y er disk is observ ed up on the addition of an equimolar concen tration of div alen t ions. In Fig. 2{7 a, the pac king parameter is sho wn for the LP A system at the v arying div alen t ion concen tration. A signican t increase in the pac king

PAGE 35

25 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 Divalent Ion/Lipid Molar RatioPacking Parameter (a) 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 Divalent Ion/Lipid Molar RatioMolecular Property Measurements a0 (nm 2 ) Tail Length (nm) Tail Volume (nm 3 ) (b) Figure 2{7: Molecular Geometry Characterization of LP A in a pure LP A system: a) LP A P ac king P arameters. b) LP A measuremen ts of tail v olume, tail length and head spacing. parameter o ccurs up on the increase in ion /lipid ratio from 0.2 to 1.0, whic h coincides with the transition from micelle to disk. T o understand wh y the pac king parameter is c hanging, the individual measuremen ts of a 0 l c and v are sho wn in Fig. 2{7 b. Both a 0 and l c decrease up on addition of the ions, while v remains relativ ely constan t. The v olume of the tails should remain constan t since the temp erature is xed; the c hange in a 0 can b e understo o d b ecause the div alen t ions act to screen the negativ ely c harged LP A headgroups from eac h other, causing an eectiv e decrease in head group size. A t 250 K and 265 K the DOPE-LP A system remained in the lamellar phase, and at 285 K and 300 K the system transitions in to the H I I phase; these state p oin ts w ere used for the pac king parameter analysis. Data w as analysed from the nal 80ns of the sim ulation after whic h the system had formed its nal phase and equilibrated. The pac king parameter for DOPE, LP A and the mean v alue for the system are sho wn in Fig. 2{8 a. Discon tin uous jump is observ ed b et w een the L p oin ts and the H I I p oin ts, indicating that the shap e of the molecules b ecomes more lik e an in v erted cone ( 5 ). In Fig. 2{8 b the measuremen ts of a 0 l c and v are

PAGE 36

26 sho wn. It is observ ed that l c is the parameter whic h is c hanging the most. Since l c is decreasing with increasing temp erature, the tails sample congurations farther a w a y from the equilibrium congurations, essen tially reducing their extension normal to the w ater lipid in terface. The theory states that v alues of the pac king parameter ab o v e 1 lead to in v erted phases, ho w ev er it can b een seen that ev en the lamellar phase p oin ts are ab o v e 1. This is due to the assumption in the calculation of the head group area, a 0 that DOPE and LP A ha v e equal size head group areas. The head group area w as estimated b y calculating the area of the w ater lipid in terface and dividing through b y the n um b er of lipid molecules at that in terface. There w as no accurate w a y to individually measure eac h of the sp ecies a 0 's separately therefore, for the purp ose of the pac king parameter calculation it w as assumed b oth sp ecies had the same a o F rom this analysis it is eviden t that the molecular mec hanism driving the phase transition is the conguration of the tails. The increased temp erature imparts a greater energy on the tails to deviate a w a y from relativ ely straigh t congurations to a more b en t and spla y ed structure. The transition states whic h a lipid systems samples as it transitions from a L to a H I I has b een studied b oth exp erimen tally and theoretically[ 23 { 27 ]. Biologically this transition is of particular in terest for understanding the mec hanism of v esicle fusion and understanding the ho w an timicrobial p eptides destro y bacterial mem branes [ 28 ]. A mec hanism prop osed b y Siegel, prop oses the bila y ers rst form a stalk phase and then the outer monola y ers con tract and form a con tact kno w as the tr ans monola y er con tact (TMC), Fig. 2{9 a. The TMC phase in v olv es formation of h ydrophobic v oid spaces whic h Kozlo vsky et al. states are energetically unfa v orable, and a dieren t in termediate stalk structure is prop osed [ 26 ]. This alternativ e stalk mo del is sho wn in Fig. 2{9 b; in order to reac h this conguration the lipid tails m ust tilt a w a y from normal to w ater{lipid in terface to prev en t the formation of v oid space.

PAGE 37

27 250 255 260 265 270 275 280 285 290 295 300 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 Temperature (K)Packing Parameter (a) 250 260 270 280 290 300 0.4 0.6 0.8 1 1.2 1.4 1.6 Temperature (K)Molecular Property a0 (nm 2 ) Lc (nm) Single Tail Volume (nm 3 ) (b) 235 240 245 250 255 260 265 270 275 280 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Temperature (K)Mean Packing Parameter (c) 235 240 245 250 255 260 265 270 275 280 0.4 0.6 0.8 1 1.2 1.4 1.6 Temperature (K)Molecular Property a0 (nm 2 ) Lc (nm) Single Tail Volume (nm 3 ) (d) Figure 2{8: Molecular Geometry Characterization of DOPE, LP A and DOP A in mixed system: a) DOPE-LP A system pac king parameters. b) DOPE-LP A system measuremen ts of tail v olume, tail length and head spacing. c) DOPE-DOP A system pac king parameters d) DOPE-DOP A measuremen ts of tail v olume, tail length and head spacing.

PAGE 38

28 (a) (b) Figure 2{9: Prop osed transitions in v esicle fusion a) Siegal stalk and TMC mo del [ 23 ]. b) Kozlo vsky and Kozlo v stalk mo del [ 26 ]. T o v erify whic h of these mec hanisms our sim ulations underw en t the a system of 90% DOPE and 10% LP A w as analyzed at 300 K. The formation of a stalk is depicted in Fig. 2{10 a{d. It is clear form examining the sim ulation images that the stalk structure is v ery similar to that prop osed b y Kozlo vsky et al.. There is no evidence of TMC of v oid region, furthermore the w ater p ores are long narro w o v als are sho wn b y Kozlo vsky et al. [ 26 ]. Using a coarse grain mo del for molecular dynamics has pro v ed to b e a qualitativ ely reliable to ol for studying lipid systems. In real system there are man y parameters and a high lev el of complexit y Ho w ev er, in our mo del system the system is greatly simplied, but the sim ulations are capable of matc hing the trends sho wn exp erimen tally The k ey result from these sim ulations is the shift in the transition temp erature, not the temp erature itself. The the transition temp erature is lo w ered b y c hanging the minor comp onen t from LP A to DOP A b oth exp erimen tally and in our system. Therefore there is a range of temp erature, in b et w een the phase transition temp erature, whic h the con v ersion of LP A to DOP A or vice

PAGE 39

29 (a) (b) (c) (d) Figure 2{10: Stalk formation in a DOPE/LP A system: a) 96 ns. b) 112 ns. c) 128 ns. d) 144 ns. The images w ere zo omed in on just 2 of the bila y ers to b etter illustrate the stalk formation. W ater w as remo v ed for clarit y the head groups are colored blac k and the tails are grey v ersa can induce a transformation. This mo del can also b e used to quan tify the molecular geometry of the molecules. The sim ulation results giv e go o d qualitativ e results, but there is a discrepancy b et w een the rep orted exp erimen tal transition temp erature and the observ ed v alues in the sim ulations. There is p oten tial to pro duce b etter quan titativ e agreemen t b y c hanging some parameters in the mo del. Increasing the h ydration of the system will stabilize the bila y er b y increasing the separation of the stac k ed bila y ers, making a transition more dicult, and should shift the transition temp erature to a higher v alue. Also, the Lennard-Jones parameters can b e altered as Marrink et al. did in their study of DOPE and DOPC to mak e the bila y er phase more stable [ ? ]. Ho w ev er in this study w e c hose not to mo dify the force eld to remain consisten t with previous sim ulations.

PAGE 40

30 One limitation of this coarse grain approac h is that the pH of the system cannot b e explicitly mo delled. The c harges on the molecules can b e mo died to rerect the desired pH, but H + ions are to o small to b e mo deled in this approac h. Ko oijman et al. do observ e an eect of c hanging the pH on the phase b eha vior of these lipid systems, without c hanging the net c harge on the lipids.

PAGE 41

CHAPTER 3 MOLECULAR MODELING OF KEY ELASTIC PR OPER TIES F OR INHOMOGENEOUS LIPID BILA YERS 3.1 In tro duction In biology mem branes can b e v ery non-uniform in their structure and comp osition. The ma jor comp onen t of the biological mem brane are lipid molecules, while they also consist of c holesterol and proteins [ 2 ]. One t yp e of molecule whic h is implicated in man y biological shap e c hanging pro cess are phosphoinositides [ 9 29 { 31 ]. In this study the eect of in tro ducing a phosphoinositide molecule, phosphatidylinositol-4-phosphate (PI4P) to a homogeneous dipalmito yl phosphatidyl c holine mem brane is studied. Tw o imp ortan t elastic constan ts will b e measured. The b ending mo dulus is a measure of the amoun t of energy required to c hange the curv ature of a mem brane. The b ending mo dulus has previously b een measured b oth exp erimen tally and from sim ulations for homogeneous systems [ 32 33 ]. In this study molecular dynamics sim ulations of mixed PI4P and DPPC bila y ers are sim ulated and analyzed to determine the eect of PI4P concen tration on the b ending mo dulus. Understanding the eect of concen tration on the b ending mo dulus is imp ortan t since in biological mem branes phosphoinositides will lo calize in one region of the mem brane. This region of lo calization is often the site a mem brane deformation in the form of a protrusion, budding or endo cytosis. A second k ey elastic parameter is the line tension existing b et w een lipid phases. If a tension exists b et w een phases this can act to cause budding or endo cytosis as the tension force acts to minimize the size of the in terface, resulting in an out of plane deformation [ 1 34 35 ]. The line tension existing b et w een a 31

PAGE 42

32 homogeneous DPPC section and an inhomogeneous PI4P/DPPC mem brane section will b e determined, again with the use of molecular dynamics sim ulations. A third elastic parameter will also b e in v estigated. When lipid molecules in a monola y er tilt a w a y from the lipid-w ater in terface normal direction and energy is asso ciated with this t yp e of deformation. The tilt mo dulus c haracterized the energy asso ciated with this t yp e of deformation. In this w ork, it will b e sho wn that tilt deformations in a hexagonal phase observ ed through MD sim ulations are consisten t with theoretical and exp erimen tal predictions [ 36 { 38 ]. 3.2 Bending Mo dulus 3.2.1 Curv ature Elasticit y A con tin uum mo del for describing the shap e of bila y er v esicles w as rst prop osed b y Helfric h [ 39 ]. The Helfric h mo del predicts a quadratic relationship b et w een the mean curv ature, C of a homogeneous bila y er and the free energy of the bila y er p er unit area, f in Eq. 3.1 f = 1 2 ( C C 0 ) 2 + K (3.1) where is the b ending mo dulus, C 0 is the sp on taneous curv ature of the bila y er, is the Gaussian mo dulus and K is the Gaussian curv ature. The mean and Gaussian curv atures are dened b elo w: C = 1 + 2 K = 1 2 where 1 and 2 are the principal curv atures of the bila y er surface. The minim um energy conguration of the bila y er can b e solv ed b y limiting the solutions to spherically symmetric shap es for a giv en v olume, V and surface area, A of the enclosed bila y er. This is a constrained minimization problem, where the ob ject is nd the curv e whic h minimizes the total energy of the bila y er giv en

PAGE 43

33 Figure 3{1: Comparison of theoretical and exp erimen tal v esicle shap es. [ 40 ]. the constan t v olume and surface area restrictions. The total functional is giv en b y Eq. 3.2 where, and P are Lagrange m ultipliers. F b = I f dA + A + P V (3.2) All functions can b e parameterized b y the arc length of the curv e dening the axisymmetric surface. A set of Euler-Lagrange shap e equations can then b e deriv ed and n umerically solv ed to yield the minim um energy shap e giv en a set of constrain ts. The phase diagram has b een mapp ed out b y Seifert et al and sev eral biologically relev an t shap es w ere determined to b e minim um energy congurations [ 40 41 ]. Figure 3{1 sho ws a comparison of exp erimen tal and theoretical v esicles shap es as the temp erature of the system is altered thereb y altering the surface area and v olume of the v esicles. Budding shap es, as w ell as, the biconca v e shap e of the red blo o d cell are observ ed. In solving the shap e equations, the in tegral of the Gaussian curv ature, o v er the surface area for closed v esicles remains constan t and therefore need not b e considered in v ariational problem. This lea v es only the mean curv ature term, and it is desired to kno w the v alue of the b ending mo dulus. This has b een determined

PAGE 44

34 exp erimen tally [ 42 ] and also from molecular dynamic sim ulations [ 33 43 ] and go o d agreemen t has b een sho wn b et w een them. In order to determine from sim ulations the thermal ructuations of an equilibrium bila y er can b e analyzed. If w e assume C 0 = 0 and neglect the Gaussian b ending term the free energy p er unit area is giv en b y f = 1 2 ( h xx + h y y ) 2 = 1 2 ( r 2 ( h ( r )) 2 where h is the heigh t of the bila y er and the subscripts denote deriv ativ es with resp ect to a direction, where the bila y er normal is the z direction. T aking the F ourier transform of the surface, Eq. 3.3 and the in v erse transform, Eq. 3.4 yields and expression for h ( r ) in terms of the w a v e v ector q h ( q ) = 1 p A Z h ( r ) e i q r d r (3.3) h ( r ) = 1 p A X q h ( q ) e i q r (3.4) By taking the F ourier transform of the surface, the b ending mo des b ecome decoupled and the total energy of the surface is giv en b y Eq. 3.5 F = Z f dA = 1 2 X q q 4 j h ( q ) j 2 (3.5) F rom the equipartition theorem, the energy of eac h mo de can b e can b e related to the F ourier co ecien ts sho wn in Eq. 3.6 < j h ( q ) j 2 > = k B T q 4 (3.6) where k B is Boltzmann's constan t and T is the temp erature of the system. The b ending mo dulus can then b e determined from the sim ulation data b y tting < j h ( q ) j 2 > v ersus q 4

PAGE 45

35 (a) (b) Figure 3{2: A tomistic and coarse grain represen tation of phosphoinositides: a) PI4P b) DPPC. 3.2.2 Molecular Mo del DPPC is one of the ma jor comp onen ts of most biological mem branes. PI4P is an imp ortan t molecule in mem brane systems b ecause it has b een implicated to pla y an essen tial role in Golgi mem brane functions and it is a precursor to PI(4 ; 5)P 2 [ 9 ]. F rom a geometric p ersp ectiv e, sim ulating systems whic h con tained PI molecules w as in triguing b ecause of the large head group due to the inositol ring and attac hed phosphate groups. The phosphate group on PI4P carries a net -2 e c harge whic h can cause the eectiv e head size to b e ev en larger if it is in the presence of other negativ ely c harged lipids. The CG mo del and atomistic structure of PI4P and DPPC are sho wn in Figure 3{2 it can b e seen that the inositol ring is mo deled as a single p olar b ead due to it's abilit y to form h ydrogen b onds and its high solubilit y in w ater. Exp erimen tal evidence has sho wn that the orien tation of the inositol ring of PI4P when in the presence of DPPC, is not normal to the bila y er, but b ecomes b en t [ 44 ]. Figure 3{3 a sho ws ho w the inositol ring gets pulled do wn to w ards the

PAGE 46

36 (a) (b) Figure 3{3: Depictions of PI4P a) A tomistic represen tation in the presence of DPPC. b) Head group of the CG mo del sho wn with the imp osed equilibrium angles, no equilibrium v alue is imp osed on 1 plane of the bila y er, rather than extend out in to the aqueous phase. It has b een sp eculated that this structural feature of PI4P is caused b y the electrostatic in teraction b et w een the negativ ely c harge phosphate group of PI4P and the p ositiv ely c harged c holine group of DPPC. Figure 3{3 b sho ws the CG mo del of the head group of PI4P and the equilibrium angles. T o try to mimic the structural features of PI4P no p oten tial has b een imp osed on the angle 1 where the angles 2 and 3 ha v e equilibrium angles of 180 o and 120 o resp ectiv ely Since no angle p oten tial is acting on the 1 angle this will allo w for the top phosphate group to mo v e and adopt a conformation due solely to the in termolecular in teractions and of course the b ond length p oten tial. 3.2.3 Sim ulation Details Self-assem bly sim ulations of DPPC and PI4P w ere conducted at relativ e PI4P concen trations of 0%, 5%, 10% and 20%. This w as done b y p opulating a sim ulations b o x with total of 200 lipid molecules at the prop er relativ e concen trations and 28 w aters p er lipid (7 CG W's), whic h w as the same concen tration used in the

PAGE 47

37 exp erimen tal determination of the head group orien tation b y Bradsha w et al [ 44 ]. The systems w ere then sim ulated for 800ns at a temp erature of 298 K, in all systems a bila y er w as formed. The bila y er system w as then copied laterally four times to create a bila y er system consisting of 800 lipid molecules. The 800 lipid system w as then sim ulated for 1.6 s. The rst 400ns w ere considered an equilibration, and the nal 1.2 s w ere used for in the analysis of the system. In b oth the selfassem bly and bila y er sim ulations Berendsen temp erature coupling and anisotropic pressure coupling (1 bar) w ere used as w ell as p erio dic b oundary conditions in all 3 directions. 3.2.4 Results T o determine the structure of the PI4P head group the 1 angle w as calculated for all systems. The probabilit y distribution of 1 and 2 angles in the PI4P head group for eac h of the systems con taining PI4P is sho wn in Figure 3{4 a and b, resp ectiv ely All of the systems displa y the same trend with a p eak at 2.1 rads (120 o ) for 1 and at = for 2 Since 2 is normal to the bila y er w ater in terface, it can b e determined that the upp er phosphate group is b eing pulled do wn to w ards the bila y er-w ater in terface as indicated b y the 1 angle data. Since the phosphate group is partially lying in the in terface, it acts to increase the head group size of PI4P This data indicates that the mo del is qualitativ ely rerecting the true nature of the PI4P head group in a PI4P/DPPC system. As describ ed in Section 3.2.1 the relationship b et w een the ructuations of the bila y er surface and the b ending mo dulus can b e determined, Eq. 3.6 T o obtain h ( q ) for a surface, a w ell dened surface is required. An imaginary lateral grid w as imp osed and eac h b ead w as assigned to a bin based up on it's lateral p osition ( x; y ). F or eac h bin the a v erage normal ( z ) p osition of all molecules in the bin w as determined. The surface w as dened b y assigning eac h bin cen ter the a v erage z p osition of all molecules in the bin giving a p oin t-wise function z = f ( x; y ),

PAGE 48

38 0 0.5 1 1.5 2 2.5 3 3.5 0 0.005 0.01 0.015 0.02 0.025 0.03 Angle (radians)Probability 5% PI4P 10% PI4P 20% PI4P (a) 0 0.5 1 1.5 2 2.5 3 3.5 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Angle (radians)Probability 5% PI4P 10% PI4P 20% PI4P (b) Figure 3{4: Probabilit y distribution for systems of PI4P/DPPC at concen trations of 5%, 10% and 20% a) 1 angle. b) 2 angle. represen ting the surface. A t eac h p oin t in time, during the analysis, a 2D F ourier transform of the surface w as tak en as describ ed b y Eq. 3.3 Eac h ( h ( q )) 2 w as a v eraged and then log ( < j h ( q ) j 2 > ) w as plotted v ersus log( j q j ), these sp ectral in tensit y plots for eac h system are sho wn in Figure 3{5 a-d. The longer w a v elength mo des, j q j > 1 : 0 nm, w ere used to t a line b y the least-squares metho d. A function can b e t to the data of the form log ( < j h ( q ) j 2 > ) = 4log( j q j ) + C 1 where the in tercept v alue, C 1 is equal to log ( k B T = )(Eq. 3.6 ) allo wing for to b e determined. The uncertain t y in C 1 C 1 w as propagated to giv e the error in as sho wn in Eq. 3.7 = d dC 1 C 1 = ( d ( k b T e C 1 ) dC 1 ) C 1 = C 1 (3.7) The b ending mo dulus and error at v arying concen trations of PI4P is sho wn in Figure 3{6 A t the time of this rep ort there w as only 200ns of data for the pure DPPC system, accoun ting for the larger error in that system. The v alues rep orted here are in the same range as previously rep orted for sim ulated and exp erimen tal systems[ 33 42 ].

PAGE 49

39 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 -4 -3 -2 -1 0 1 2 3 4 Log q(1/nm)Log(fc2) (nm4) (a) 1.5 1 0.5 0 0.5 1 1.5 2 2.5 4 3 2 1 0 1 2 3 4 Log q(1/nm)Log(fc2) (nm4) (b) 1.5 1 0.5 0 0.5 1 1.5 2 2.5 4 3 2 1 0 1 2 3 4 Log q(1/nm)Log(fc2) (nm4) (c) 1.5 1 0.5 0 0.5 1 1.5 2 2.5 4 3 2 1 0 1 2 3 4 Log q(1/nm)Log(fc2) (nm4) (d) Figure 3{5: Sp ectral in tensit y for PI4P/DPPC systems a) Pure DPPC. b) 5% PI4P c) 10% PI4P d) 20% PI4P 0 5 10 15 20 25 8 8.5 9 9.5 x 10 20 % PI4PBending Modulus (J) Figure 3{6: Bending mo dulus of PI4P/DPPC bila y er at v arying concen trations of PI4P

PAGE 50

40 3.3 Line T ension The existence of a tension force existing b et w een lipid phases can b e a driving force for shap e c hange of biomem branes. F rom a molecular prosp ectiv e a tension force b et w een lipid phases can b e generated due to the diering in teraction energies b et w een the molecules of the resp ectiv e phases. The lipid phases ma yb e a pure raft or can b e a region of the mem brane where a sp ecic lipid has lo calized. A k ey mem brane lipid in man y biological pro cesses are phosphoinositides whic h are kno wn to pla y a role b oth in ter and in tracellular trac king [ 7 { 9 31 ]. A t the outer cellular mem brane there is a cytosk eleton whic h pla ys a role of structural stabilit y for the cell as w ell as eecting shap e c hange. A t the outer mem brane the mec hanism of shap e c hange is quite complex in v olving the in teraction b et w een lipids, proteins and the cytosk eleton. Ho w ev er, in tracellular compartmen ts, suc h as the Golgi b o dy and the endoplasmic reticulum do not ha v e a cytosk eleton on their in terior, but are still able to mediate shap e c hanges. Therefore the mec hanism of shap e c hange at the in tracellular mem brane is less complex and is lo cation where this w ork is fo cused. 3.3.1 Pressure Calculation from Molecular Dynamics Sim ulations Pressure in a ruid is made up of t w o-comp onen ts: a kinetic part due to the molecular motions and a static part due the the in termolecular p oten tials [ 45 ]. Pressure is a force acting on the unit v olume, so the kinetic part can b e view ed as the momen tum transfer across an imaginary surface. The static part is due to particles on opp osites sides of the imaginary surface exerting a force on eac h other due to the in termolecular p oten tials. The kinetic part of pressure can b e calculated from the temp erature of the system or from the molecular v elo cities as sho wn in Eq. 3.9 and Eq. 3.10

PAGE 51

41 P = P K + P S (3.8) P K = N k B T V (3.9) P K = 2 V E K = 1 V N X i m i v i n v i (3.10) Where P K is the scalar or h ydrostatic pressure, it is related to the pressure tensor, P K b y P K = tr ace ( P K ) = 3. N is the n um b er of particles in the system, V is the v olume of the system, E K is the kinetic energy tensor, m i is the mass of eac h particles and v i is the v elo cit y v ector of eac h particle. The calculation of the kinetic pressure from molecular dynamics is trivial. If the v elo cities of eac h particle are stored to an output tra jectory le the kinetic pressure can b e calculated easily from Eq. 3.10 The static pressure calculation can b e more complex, it is related to the virial of the system [ 46 ]. P S = 2 V Vir (3.11) Vir = 1 2 N X i r i n F i (3.12) The virial is a 2nd order tensor sho wn in Eq. 3.12 where r i is the v ector describing the particle p osition and F i is the force v ector acting on particle i The complexit y of the virial calculation comes in determining the forces acting on the particles. In theory these forces can b e due to external and in ternal p oten tials, but in this w ork only in ternal forces w ere presen t. The p oten tials are an input to the mo del and the forces on eac h b ead are calculated from the negativ e deriv ativ e of the p oten tial with resp ect to particle p osition, F ( i )= dV dr i The p oten tials presen t in the coarse-grain mo del sim ulations are describ ed in Section 3.3.2

PAGE 52

42 3.3.2 Description of P oten tials and F orces Lennard-Jones P oten tial The Lennard-Jones(LJ) p oten tial is a non-b onded p oten tial whic h represen ts the the long range attractiv e v an der W aals force and a short range repulsiv e part due to o v erlap of the electron clouds. F or the sim ulations conducted in this w ork, the LJ p oten tial is shifted b ey ond a certain radius so the forces go smo othly to zero at the cut-o radius. In the sim ulation describ ed in this w ork, the cut o distance w as r LJ c = 1 : 2 nm and the shift distance w as r LJ s = 0 : 9 nm. The LJ in teractions b et w een nearest b onded neigh b ors are excluded. The p oten tials, V LJ and forces, F LJ ( i ) acting on b ead i for the unshifted and shifted parts of the p oten tial are giv en b y V LJ = 8>>>>>>>>>><>>>>>>>>>>: C 12 r 12 ij C 6 r 6 ij r ij < r LJ s 12 C 12 h 1 12 r 12 ij A 12 3 ( r ij r LJ s ) 3 B 12 4 ( r ij r LJ s ) 4 D 12 i 6 C 6 h 1 6 r 6 ij A 6 3 ( r ij r LJ s ) B 6 4 ( r ij r LJ s ) D 6 i r LJ s r ij < r LJ c 0 r ij r LJ c and F ( i ) LJ = 8>>>>>>>>>><>>>>>>>>>>: 12 C 12 r 13 ij + 6 C 6 r 7 ij r ij r ij r ij < r LJ s 12 C 12 h 1 r 13 ij + A 12 ( r ij r LJ s ) 2 + B 12 ( r ij r LJ s ) 3 i 6 C 6 h 1 r 7 ij + A 6 ( r ij r LJ s ) 2 + B 6 ( r ij r LJ s ) 3 i r ij r ij r LJ s r ij < r LJ c 0 r ij r Lj c

PAGE 53

43 Figure 3{7: Denition of r ij where r ij is the v ector p oin ting from b ead i to b ead j as describ ed in Figure 3{7 r ij is the magnitude of r ij and A = ( + 4) r LJ c ( + 1) r LJ s r +2 LJ c ( r LJ c r LJ s ) 2 B = ( + 3) r LJ c ( + 1) r LJ s r +2 LJ C ( r LJ c r LJ s ) 3 D = 1 r LJ c A 3 ( r LJ c r LJ s ) 3 B 4 ( r LJ c r LJ s ) 4 : The force on b ead j is equal and opp osite to the force on b ead i F ( j ) LJ = F ( i ) LJ This is also true of the 2-b o dy p oten tials for electrostatic and b onded in teractions describ ed b elo w. Electrostatic P oten tial The electrostatic p oten tial (ES) is a non-b onded p oten tial whic h acts b et w een b eads with a net c harge. The p oten tial is shifted smo othly to zero at the cut o radius, r E S c = 1 : 2 nm starting from r E S s = 0 : 0 nm. The cut o and the shifting and of the electrostatic p oten tial is done to mimic the distance dep endan t screening phenomenon [ ? ]. Electrostatic in teractions b et w een nearest b onded neigh b ors are excluded. The electrostatic p oten tial and resultan t force are describ ed b elo w.

PAGE 54

44 V E S 8>><>>: = q i q j 4 o r h 1 r ij A 1 3 ( r ij r E S s ) 3 B 1 4 ( r ij r E S s ) 4 D 1 i r ij < r E S C = 0 r ij r E S c F ( i ) E S 8>><>>: = q i q j 4 0 r h 1 r 2 ij + A 1 ( r ij r E S S ) 2 + B 1 ( r ij r E S S ) 3 i r ij r ij r ij < r E S C = 0 r ij r E S c V E S is the electrostatic p oten tial, F ( i ) E S is the force due to the electrostatic p oten tial acting on b ead i q i is the c harge on b ead i o the p ermittivit y of free space and r is the relativ e dielectric constan t, whic h w as tak en to b e r = 20. The constan ts A 1 B 1 and D 1 are describ e b y the same equations for the constan ts in the LJ p oten tial. Bonded P oten tial The b onded p oten tial, V B is a harmonic p oten tial whic h corrects for deviations from equilibrium b ond lengths. The b ond p oten tial and resultan t force, F B are describ ed b elo w. V B = 1 2 k bond ( r ij r 0 ) 2 F ( i ) B = k bond ( r ij r 0 ) r ij r ij The b onding force constan t for a single b ond is k bond = 1250 k J mol nm 2 and the equilibrium b ond length is tak en to b e r 0 = 0 : 47 nm. Angle P oten tial The angle p oten tial can b e used to constrain three consecutiv ely b onded b eads. The p oten tial is a harmonic function of the cosine of the angle made ab out the cen tral b ead. The angle orien tation and resultan t forces are describ ed in Figure 3{8 The angle p oten tial, V and forces, F are describ ed b elo w,

PAGE 55

45 Figure 3{8: Angle Orien tation V = 1 2 k (cos cos 0 ) 2 F ( i ) = k (cos cos 0 ) ( r j i r j k ) r j i r 3 j i r j k r j k r j i r j k F ( j ) = k (cos cos 0 ) r j i r j k r j i + r j k ( r j i r j k ) r j i r 2 j h + r j k r 2 j i # F ( k ) = k (cos cos 0 ) ( r j i r j k ) r j k r 3 j k r j i r j i r j i r j k # where the angle force constan t for single b onds is k = 25 k J mol and the equilibrium angle is 0 = for. The cosine of the angle can b e determined b y taking the dot pro duct of r j i and r j k cos( ) = r j i r j k r j i r j k [ 47 ]. 3.3.3 Calculation of Virial T o calculate the virial of the system, the particle p ositions and total force on eac h particle m ust b e kno wn as describ ed in Eq. 3.12 F or the 2-b o dy p oten tials the force is a function of separation distance b et w een t w o in teracting particles. It is computationally ecien t to calculate the virial due to eac h in teraction Vir ( i; j ), Eq. 3.13 and then sum o v er all in teractions Eq. 3.14

PAGE 56

46 Vir ( i; j ) = 1 2 [ r i n F ij + r j n F j i ] = 1 2 [ r ij n F ij ] (3.13) Vir 2 Bo dy = 1 2 X i
PAGE 57

47 Figure 3{9: Slab orien tation for calculation of bila y er surface tension P S ( R ; t ) = X i [ r i V ( f r i g )] I C 0 i dl ( R l ) (3.16) The most natural c hoice for the con tour is the Irving-Kirkw o o d con tour whic h is a straigh t line b et w een in teracting particles [ 50 ]. Other con tours ha v e b een prop osed, most notably the Harasima con tour, but in this w ork only the Irving-Kirkw o o d con tour is used [ 51 ]. A generalized expression for the static pressure for an m-b o dy p oten tial can b e deriv ed from Eq. 3.16 is describ ed b y 3.17 P S m Bo dy ( R ; t ) = 1 m X X r j k V m r j l V m I C j l j k ( R l ) dl (3.17) where < j > represen ts all \clusters", whic h means a group of in teracting particles, and k and l are the indices of particles within the cluster. The a v erage pressure in eac h slab is what is desired to b e kno wn, so Eq. 3.17 m ust b e in tegrated o v er the dimensions of the slab, P S m b o dy ( Z i ; t ) = 1 V ol Z Z L x 0 dx Z L y 0 Z Z i + Z Z i dz P S m b o dy ( R ; t ) (3.18)

PAGE 58

48 where Z i is the i th slab, Z is the slab thic kness, and V ol Z is the v olume of the slab. The c hoice of con tour m ust b e dened to ev aluated the in tegral, here w e use a straigh t line con tour (Irving-Kirkw o o d) connecting in teracting particles l = r j k + ( r j l r j k ) ; [0 ; 1] (3.19) dl = ( r j l r j k ) d = r j k j l d (3.20) I C j l j k ( R l ) dl = r j k j l Z 1 0 ( R ( r j k + r j k j l ) d (3.21) where r j k j l = r j l r j k Com bining Eqs. 3.17 3.18 and 3.21 an expression whic h can b e ev aluated for the static pressure in a slab is deriv ed P S ( Z i; t ) = 1 mV ol Z X X r j k V m r j l V m r j k j l f ( z j k ; z j l ; Z i ) where z j k is the z p osition co ordinate of particle k in j cluster, and f ( z j k ; z j l ; Z i ) = Z L x 0 Z L y 0 Z Z i + Z Z i dxdy dz Z 1 0 ( R ( r j k + r j l )) d (3.22) Z L x 0 Z L y 0 Z Z i + Z Z i ( R ( r j k + r j l )) dxdy dz = 8>><>>: 1 R = r j k + r j k j l R V ol Z 0 R 6 = r j k + r j k j l (3.23) so R m ust lie on the r j k j l and b e within slab i for the spatial in tegration to b e nonzero. Since the slab spans the system in the x and y directions only the particle z p ositions are critical in ev aluating the pressure in a slab. Eq. 3.23 can b e re-written

PAGE 59

49 with the use of the Hea viside step function. Z L x 0 Z L y 0 Z Z i + Z Z i ( R ( r j k + r j l )) dxdy dz = ( z j k + ( z j l z j k ) Z i ) ( Z i + Z + z j k + ( z j l z j k )) ( x ) = 8>><>>: 1 x > 0 0 x < 0 f ( z j k ; z j l ; Z i ) = Z 1 0 ( z j k + ( z j l z j k ) Z i ) ( Z i + Z + z j k + ( z j l z j k )) d So dep ending on the lo cation of the in teracting particles with resp ect to slab i lo cation, f ( z j k ; z j l ; Z i ) will tak e on dieren t v alues. If b oth particles lie within slab i f = 1; if b oth particles lie outside of slab i and no part of r j k j l lies within slab i f = 0; if b oth particles do not lie with slab i but part of r j k j l do es pass through slab i then f = r Z i j k j l r j k j l where r Z i j k j l is the p ortion of r j k j l within slab i 3.3.5 Calculation of Line T ension Surface tension is calculated from the dierence b et w een lateral and normal comp onen ts of the pressure tensor [ 52 ]. r = Z + 1 1 ( P N ( z ) P L ( z )) dz The z -direction is normal to the in terface, P N and P L are the normal and lateral pressure, resp ectiv ely This deviation in the isotrop y of the pressure tensor can b e due to b oth densit y deviations and c hanges in conguration energy The in tegration is theoretically tak en from 1 to + 1 since the pressure b ecomes isotropic a w a y from the in terface and do es not con tribute to the in tegration.

PAGE 60

50 Figure 3{10: Phases within a bila y er T o calculate line tension, the pressure m ust b e lo calized to the surface dened b y the mem brane and then the dierences b et w een the comp onen ts of this 2dimensional pressure tensor will result in line tension. Figure 3{10 describ es the geometry of phases em b edded within a bila y er. In this situation pressure w ould b e calculated as a function of the z -direction, with slabs spanning the b o x in the x and y directions. ( t ) = 1 N int Z L Z 0 ( P Z Z 2 D ( z ; t ) P X X 2 D ( z ; t )) dz (3.24) P 2 D ( z ; t ) = Z L y 0 P ( z ; y ; t ) dy is line tension, P 2 D is the 2-dimensional pressure tensor and N int is the n um b er of in terfaces in the system, since the line tension p er in terface is desired. Since P ( z ; t ) is what will b e calculated and it is already a v eraged o v er the y -direction w e can see that P 2 D ( z ; t ) = P ( z ; t ) L y and b y com bining this with Eq. 3.24 an expression for line tension as a function of P ( z ; t ) is deriv ed. ( t ) = L y N int Z L Z 0 ( P Z Z ( z ; t ) P X X ( z ; t )) dz (3.25)

PAGE 61

51 3.3.6 Sim ulations Details The same CG mo dels for PI4P and DPPC are used in sim ulations for the calculation of line tension. The 800 lipid sim ulations generated for the b ending mo dulus calculations as describ ed in Section 3.2.3 w ere used to assem ble the line tension systems, referred to from here forw ard as a patc h system. Figure 3{10 describ es the desired conguration for the line tension calculation in whic h the -phase represen ts a homogeneous pure DPPC section and the -phase represen ts an inhomogeneous section of PI4P and DPPC. The initial bila y er sections w ere tak en from the output of 1.6 s for the 800 lipid systems. A new top ology w as then generated b y placing the inhomogeneous section in b et w een to homogeneous sections. T o do this, the inhomogeneous section had to b e mo died so the median bila y er heigh t w as the same as the homogeneous section. Also, it w as required that b oth systems ha v e iden tical x and y dimensions, so whic h ev er of the t w o systems w as larger w as trimmed to matc h the dimensions of the smaller section. The patc h system w as generated with b oth a 5% PI4P and 20% PI4P inhomogeneous section. These patc h systems w ere then energy minimized for 50000 steps using the steep energy minimization metho d. Additional w ater w as then added to the system to bring the w ater lev el to 22 CG W p er lipid and again the system w as energy minimized for 50000 steps. In order the k eep the sections separated and to prev en t the PI4P molecules from diusing out in to the homogeneous section a repulsiv e w all w as constructed, whic h acted only on the PI4P molecules. The w all spanned the sim ulation b o x in the x and y directions. The \dumm y" atoms whic h made up the w all w ere separated with 0.5 nm spacing in b oth the x and y directions. The z co ordinate lo cation of the w all at high z v alue w as determined b y adding 0.4 nm to the maxim um z co ordinate p osition of all PI4P b eads. Lik ewise the w all at lo w z v alue w as determined b y subtracting 0.4 nm from the minimal z co ordinate p osition of all PI4P b eads. The dumm y atoms w ere k ept at a xed p osition throughout the

PAGE 62

52 Figure 3{11: P atc h system with repulsiv e w alls sim ulations. The in teraction p oten tial b et w een the dumm y atoms and the PI4P b eads w as equiv alen t to the repulsiv e part of the LJ p oten tial for p olar nonp olar in teractions V Dumm y = 8>><>>: C 12 r 12 ij if in teracting with PI4P 0 if in teracting with b ead other than PI4P where C 12 = 8 : 3658 10 2 This p oten tial is shifted and cut-o with the same radii and in the same manner as all other LJ p oten tials as describ ed in Section 3.3.2 Figure 3{11 is a sim ulation snapshot of the patc h system conguration, w ater w as remo v ed for clarit y The blac k b eads extending through the bila y er are the dumm y w all atoms. The section b et w een the w alls is a mixture of PI4P and DPPC, while the outer sections of the bila y er are pure DPPC. These systems w ere then sim ulated for 400ns using Berendsen pressure and temp erature coupling, set to 1 bar anisotropically and 323 K resp ectiv ely The pressure coupling w as turned o after the initial 400ns sim ulation and w as run for an additional 400ns with no pressure coupling. 3.3.7 Results Comparison to A tomistic Results Sev eral studies ha v e b een conducted on the calculation of pressure proles and surface tension of homogeneous bila y ers from atomistic molecular dynamics

PAGE 63

53 sim ulations [ 48 53 54 ]. The p oten tials used in atomistic sim ulations dier from CG sim ulations in man y w a ys. A tomistic sim ulations t ypically ha v e more p otentials, often including a dihedral p oten tial whic h is not included in the CG mo del. In atomistic sim ulations the electrostatics can b e v ery dieren t b ecause of the presence of partial c harges on the ma jorit y of atoms, where in the CG mo del most b eads are mo deled as neutral b eads. Additionally the electrostatics ma y b e not b e cut o in atomistic sim ulations to allo w for long range in teractions. In the CG mo del the LJ p oten tial partially accoun ts for the electrostatic. Lastly the co ecien ts of all of the p oten tials will v ary from CG to atomistic sim ulations making the comparison b et w een virial con tributions dicult. Ho w ev er, if the CG mo del is a reliable mo del the pressure proles b et w een CG and atomistic mo del should b e comparable. The pressure prole for a pure DPPC bila y er w as calculated for the same system used for the b ending mo dulus calculation (800 lipids, 7 CG W/lipid, 298 K, anisotropic pressure coupling). The analysis w as done on 1.2 s of data after 400ns of equilibration. This data w as compared to results obtained for a similar atomistic system study conducted b y Sonne et al [ 48 ]. The system used b y Sonne consisted of 72 DPPC molecules sim ulated at 300 K with a ratio of 28 w aters p er lipid (same concen tration as our system). The pressure prole w as calculated as a function of the bila y er normal direction, the z direction for con v enience. In the analysis done here the system w as divided in to 60 slabs, in the w ork b y Sonne the system w as divided in to 70 slabs [ 48 ]. In Figure 3{12 a the densit y prole is plotted, where the units are n um b er of b eads p er slab. In Figure 3{12 the diagonal comp onen ts of the static pressure tensor are plotted as a function of z p osition. It can b e seen that the normal comp onen t remains constan t through the bila y er.

PAGE 64

54 4 3 2 1 0 1 2 3 4 0 50 100 150 200 250 300 Z position (nm)Number Beads W DPPC (a) 4 3 2 1 0 1 2 3 4 700 600 500 400 300 200 100 Z position (nm)Pressure (bar) Px(z) Py(z) Pz(z) (b) Figure 3{12: Densit y and pressure of CG systems and a function of bila y er p osition: a) Densit y prole. b) Pressure prole. In Figure 3{13 the virial con tributions and pressure proles are compared for a CG and atomistic mo dels. The CG virial con tributions are quite dieren t from the atomistic virial con tributions, esp ecially in the magnitude of the n um b ers. This dierence in magnitude can b e explained b y the dierence in the p oten tials and resultan t forces. T o illustrate the dierence, the b onded p oten tial is examined for b oth the CG and atomistic mo dels. Figure 3{14 plots the b onded p oten tial for b oth mo dels. F or a deviation from the equilibrium b ond length on the order of k B T the forces for eac h mo del w ere calculated and an order of magnitude dierence is observ ed. Ho w ev er, the Pressure proles in Figures 3{13 c and d sho w go o d agreemen t b et w een the atomistic and CG data. F B 8>><>>: 1 : 3 10 3 A tomisitic 7 : 9 10 1 CG

PAGE 65

55 4 3 2 1 0 1 2 3 4 600 500 400 300 200 100 0 100 200 300 400 Z position (nm)PL PN (bar) LJ Bond Angle Electrostaic (a) (b) 4 3 2 1 0 1 2 3 4 300 200 100 0 100 200 300 Z position (nm)PL PN (bar) (c) (d) Figure 3{13: Virial and pressure prole comparison: a) CG virial prole. b) A tomistic virial prole [ 48 ]. c) CG static pressure prole. d) A tomistic pressure prole [ 48 ]. 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 0 5000 10000 15000 Bond Displacement (nm)Energy (kJ/mol)Bond Energy Atomistic CG Figure 3{14: Bond p oten tial comparison

PAGE 66

56 Line T ension for P atc h Bila y er Systems Calculation of line tension w as conducted for patc h systems with inhomogeneous sections con taining 5% and 20% PI4P The analysis w as done on the nal 80ns of data for the sim ulations done without pressure coupling. In the analysis the pressure w as calculated as a function of the z direction, normal to the phase in terface, as describ ed in Figure 3{10 The system w as divided in to 100 slabs to calculate the pressure proles. The pressure dierence w as calculated b et w een the z and x directions (in plane comp onen ts) so that the line tension could b e calculated using Eq. 3.25 The n um b er densit y of PI4P b eads and the pressure dierences for the 5% and 20% patc h systems are sho wn in Figure 3{15 It can b e observ ed that the 5% system sho ws no distinct deviation in the pressures near the phase b oundary Ho w ev er, in the 20% patc h system there is a sharp c hange in the pressure dierence prole at the phase b oundary indicating the presence of a line tension. The line tension w as calculated for b oth systems, for the 5% system = 2 : 3 10 11 N and for the 20% system = 2 : 0 10 10 N. The v alue for the 5% can b e attributed to noise and p o or statistics, ho w ev er the 20% systems sho ws a sharp deviation from isotropic b eha vior near the b oundary The v alue obtained here is higher than rep orted for an exp erimen tal system (9 10 13 N) and theoretical estimate (10 11 )N [ 1 55 ]. 3.4 Tilt Deformations 3.4.1 Energy Mo del for Tilt In the pioneering w ork b y Helfric h, on the description of equilibrium v esicle shap es, he c ho ose to neglect the energy asso ciated with tilt of lipid molecules a w a y from the bila y er normal direction [ 39 ]. In the case of rat bila y ers the lipids on a v erage do not tilt a w a y from the normal direction, ho w ev er other lipid phases including the H I I and in termediate stalk phase (as describ ed in Chapter 2) do require the molecules to tilt and this energy m ust b e accoun ted for when predicting

PAGE 67

57 0 5 10 15 20 25 30 35 40 45 0 5 10 15 20 25 Z position (nm)Number PI4P Beads (a) 0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 35 40 45 50 Z position (nm)Number PI4P Beads (b) 0 5 10 15 20 25 30 35 40 45 400 300 200 100 0 100 200 300 400 Z position (nm)PX X PZ Z (bar) (c) 0 5 10 15 20 25 30 35 400 300 200 100 0 100 200 300 400 Z position (nm)PX X PZ Z (bar) (d) Figure 3{15: Densit y and pressure proles for patc h system:a) 5% PI4P densit y prole. b) 20% PI4P densit y prole. c) 5% PI4P pressure dierence prole. d) 20% PI4P pressure dierence prole

PAGE 68

58 (a) (b) (c) Figure 3{16: Tilt and b ending: a)Denition of tilt. b) Non-uniform tilt deformation on a v olume elemen t of monola y er [ 36 ]. c) Bending deformation on a v olume elemen t of monola y er [ 36 ]. stabilit y of phases. Hamm and Kozlo v in tro duced a energy mo del for tilt and used this to predict the stabilit y of the lamellar and in v erted lipid phases as w ell as geometric prop erties of the in v erted phases [ 37 ]. The denition of tilt and the elastic energy p er unit area p er monola y er are describ ed b y Eq. 3.26 and Eq. 4.2 resp ectiv ely and graphical represen tation of tilt is sho wn in Figure 3{16 a. t = n n N N (3.26) f = 1 2 ( t + J s ) 2 + det ( t ) + 1 2 t 2 (3.27) The tilt v ector, t is p erp endicular to N so if w e tak e N to p oin t in the z direction, t will ha v e comp onen ts in the x and y directions. The tilt tensor is dened b y the deriv ativ es of tilt t = r t where r = @ @ where = x; y t = tr ace ( t )= r t J s is the sp on taneous curv ature of the monola y er, and is the tilt mo dulus. It is in teresting to note that the same elastic parameters, and describ e the energy of b oth tilt and b ending deformations. This is explained in a later pap er b y Hamm and Kozlo v and can b e visualized b y the Figures 3{16 b and

PAGE 69

59 (a) (b) (c) Figure 3{17: Construction of H I I phase: a) Original construction of hexagonal phase unit cell [ 37 ]. b) Multiple original construction cells sho wn to illustrate the h ydrophobic v oid region [ 38 ]. c) Recen t hexagonal phase unit cell construction [ 37 ]. c [ 36 ]. Figure [ 36 ]b describ es a non-uniform tilt deformation, where tilt c hanges across the elemen t. Figure 3{16 c describ es a b ending deformation, it can b e seen that the shear and stretc h of this elemen t a essen tailly essen tially the same in b oth deformations, whic h is allo ws for b oth to b e c haracterized b y the same parameters. 3.4.2 Nature of the H I I Phase and Evidence of Tilt It had long b een assumed that the structure of the H I I phase consisted of a circular h ydrophilic in terface and a hexagonal h ydrophobic in terface b et w een unit cells as sho wn in Figure 3{17 This construction allo ws for a mismatc h b et w een the in terior h ydrophilic b oundary and the outer h ydrophobic b oundary In this construction, the lipid tails m ust either stretc h to ll the corners of outer b oundary or a h ydrophobic v oid region will exist as illustrated in the cen ter of the three cells in Figure 3{17 b; b oth situation are energetically unfa v orable. T o alleviate this \frustration energy", as it is kno w, the construction of the unit cell w as re-examined exp erimen tally and theoretically to construct the curren t description [ 37 38 56 ]. This curren t description allo ws for the corners of the cell to b e lled without stretc hing of the h ydro carb on c hains, but comes at the exp ense of tilt deformations along the faces of the hexagon.

PAGE 70

60 F rom the MD sim ulations describ ed in Chapter 2 the hexagonal phase w as examined. F or a system of 10% LP A and 90% DOPE sim ulated at T = 300 K as describ ed in Section 2.3.3, b eginning from a m ulti-lamellar conguration formed a stable H I I phase. This H I I w as examined to ev aluate if the curren t hexagonal construction w as observ ed. Figure 3{18 a sho ws a sim ulation snapshot from the ab o v e describ ed system after 1.2 s of sim ulation, it w ater is remo v ed to sho w that the w ater p ores ha v e a hexagonal geometry Figure 3{18 b illustrates the hexagonal nature of the h ydrophobic b y plotting the lo cation of all tail b eads in the lipids. Figure 3{18 c plots the lo cation of the lipid-w ater in terface, dened as the mid-p oin t of the b ond connecting the phosphate b ead with the glycerol b ead and tted to a hexagon b y least-squares. It is eviden t from the Figure 3{18 the H I I phase observ ed in our MD sim ulations do rerect the description of Hamm and Kozlo v [ 37 ]. In order for the lipids to accommo date the structure prop osed in Figure 3{17 c and v eried b y Figures 3{18 a-c, the lipids m ust tilt a w a y from the normal direction to the lipid-w ater in terface. A t the cen ter of eac h face of the hexagonal phase the lipid will b e aligned with the normal, but a w a y from the cen ter, the lipids m ust tilt. Hamm and Kozlo v predict the tilt should b e linear, that it will reac h it's maxim um v alues that the outer edge of eac h face, and along the p ore direction tilt should b e zero [ 37 ]. This has b een v eried through analysis of the MD sim ulation of the ab o v e describ ed 10% LP A/DOPE system. T o calculate tilt for eac h lipid the director, n m ust b e suitably dened. In this w ork, n is dened b y the v ector p oin ting from the lipid-w ater in terface lo cation to the cen ter of mass of the tail b eads. The tilt along a face of the hexagonal lipid-w ater in terface, t x and along the p ore direction t y a v eraged o v er 50 data p oin ts in time b eginning after 800ns of sim ulation, with a time step of 1.6ns. Figure a sho ws t x as a function of p osition along the facet of the hexagonal p ore. It can b e seen t x is a linear function with

PAGE 71

61 (a) (b) 1 0 1 2 3 4 5 6 27 28 29 30 31 32 33 (c) Figure 3{18: Examination of H I I from MD sim ulations for a 10% LP A/DOPE system at T = 300 K after 1.2 s : a) Sim ulation snapshot of H I I phase. b) T ail b ead lo cations plotted to illustrate the h ydrophobic b oundary c) Lipid-w ater in terface lo cation (dened as the middle of the phosphate-glycerol b ond) plotted for eac h lipid to illustrate the h ydrophilic b oundary

PAGE 72

62 1.5 1 0.5 0 0.5 1 1.5 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 X position (nm)Tx (a) 0 5 10 15 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 Y position (nm)Ty (b) Figure 3{19: Tilt in a H I I phase :a) t x tilt along the one face of the hexagon. b) t y tilt along the direction of the w ater p ore. extremal v alues the edge of the face. Tilt along the p ore direction is sho wn in Figure b, it is nearly zero and has no sustained gradien t. Both of these results in Figure a and b corresp ond to the theoretical predictions of Hamm and Kozlo v[ 37 ].

PAGE 73

CHAPTER 4 CONCLUSIONS 4.1 Ma jor Conclusions One of the adv an tages of coarse grained mo dels in molecular dynamics is that the systems can b e visualized and measured on near atomic length scales, while relativ ely large systems and long sim ulations can b e sim ulated without b eing extremely computationally exp ensiv e. This abilit y allo ws for insigh t to b e gained on the molecular scale as w ell as the abilit y to compute bulk prop erties for these systems. Here w e w ere able to v erify the c hanging molecular geometries in lipid systems as a function of the en vironmen t parameters. 1. The CG mo del dev elop ed for the lipids LP A and DOP A and the ion Mg +2 exhibit qualitativ e agreemen t with exp erimen tal temp erature and c harge based phase transition data. 2. F or pure LP A, a negativ ely c harged lipid, the head size reduced with increasing cationic concen tration, making the molecule more cylindrical and reducing the curv ature of the system. F or mixed DOP A-DOPE the increase in temp erature of the system caused the tails to spla y further apart, while the head group separation also increased with the increased vibrational energy but to lesser degree. Thereb y it w as observ ed the increase in temp erature for this system causes the molecules to b ecome more conical leading to a system of higher curv ature.This study pro vided a basis for further in v estigation in biological systems with a coarse grained mo del for the lipids LP A and DOP A. The abilit y the qualitativ ely repro duce exp erimen tal phases in b oth pure and mixed lipid systems is justication for extending the w ork in this area. 63

PAGE 74

64 3. The in termediate stalk phase observ ed in phase transition sim ulations of mixed lipid systems v eries the structure prop osed b y Kozlo vsky and Kozlo v [ 26 ]. 4. Increasing the relativ e concen tration of PI4P in a PI4P/DPPC bila y er causes an increase in the b ending mo dulus of the mem brane compared to a pure DPPC bila y er. This indicates that the bila y er b ecomes more dicult to deform with the addition of more PI4P up to 20%. 5. A system in whic h PI4P lo calizes in one region of a DPPC bila y er while neigh b oring regions remain homogeneous, exhibits a p ositiv e line tension on the order of 10 10 N when the concen tration of PI4P in the inhomogeneous section reac hes 20%. If the inhomogeneous section only con tains 5% PI4P no anisotrop y in the 2-dimension pressure tensor is observ ed. 6. The hexagonal phase in these MD sim ulations displa y hexagonal phases at the b oth the h ydrophilic and h ydrophobic b oundaries. The tilt along a side of the hexagon h ydrophilic in terface is a linear increasing function. Both of these results supp ort the theoretical predictions of Hamm and Kozlo v [ 37 ]. 4.2 F uture Directions In this w ork molecular mo deling and analysis tec hniques w ere dev elop ed whic h pro duced results w arran ting further in v estigation. In this w ork only a few lipid systems w ere studied, while it is b eliev ed that other lipid sp ecies could b e k ey molecules in eliciting c hanges in elastic parameters of mem branes. 4.2.1 Additional Lipid Sp ecies The lipid phosphatidylinositol-4,5-bisphosphate (PI(4,5)P 2 ) has b een implicated to pla y a role in stim ulating v esicle formation from lip osomes as w ell as from the outer cellular mem brane [ 29 30 35 ]. PI(4,5)P 2 is structurally similar to PI4P but has an additional phosphate group extending from the inositol ring and carries a net 4 c harge compared to 3 of PI4P at ph ysiological conditions [ 57 ]. It is

PAGE 75

65 Figure 4{1: Prop osed CG mo dels of additional lipids: a) PI(4,5) 2 b) D A G. b eliev e that the increased head group size and increased eectiv e head group size due to c harge w ould p oten tially increase the anisotrop y in the 2 dimension pressure tensor in the line tension study Additionally the lipid diacyl glycerol (D A G) is of in terest, due to it role in Golgi budding and bud ssion from the Golgi net w ork. D A G is similar to the lipid DOP A except that it do es not ha v e the phosphate group extending from the glycerol group. It is desired to conduct similar studies done in this w ork, but using PI(4,5)P 2 and D A G as the minor comp onen t in the lipid systems. Prop osed CG mo dels of the lipids PI(4,5)P 2 and D A G as sho wn in Figure 4{1 4.2.2 Con tin uum Scale Mo deling Ultimately it is desired to use the parameters determined from molecular mo deling and incorp orated these ndings in to a con tin uum scale mo del for predicting equilibrium mem brane congurations, based up on the mem brane c hemical comp osition. W ork has b een done in this area starting with a theoretical description of the b ending energy b y Helfric h and extending b y Lip o wsky for homogeneous mem branes [ 39 41 ]. This w ork w as further extended to incorp orate the existence of dieren t homogeneous lipid phases with the same lipid b y mo difying the Helfric h expression to incorp orate line tension b et w een phases [ 1 34 35 ]. The expression of the b ending and in terface energy is for a system with distinct lipid phases is giv en

PAGE 76

66 Figure 4{2: Domains within a mem brane, tak en from ref [ 1 ] b y Eq. 4.1 F = Z [ ( ) 2 ( C 1 + C 2 + C ( ) 0 ) 2 + ( ) C 1 C 2 ] dA + Z [ ( ) 2 ( C 1 + C 2 + C ( ) 0 ) 2 + ( ) C 1 C 2 ] dA + Z d dl (4.1) where ( ) ( ) C ( ) o are the b ending mo dulus, Gaussian curv ature mo dulus and sp on taneous curv ature of the phase, d is the length of the in terface separating the phases. The geometry of the system is describ ed b y Figure 4{2 A mo del describing the energy of asso ciated with tilt and b ending of lipid structures has b een in tro duced b y Hamm and Kozlo v sho wn in Eq. 4.2 and previously describ ed in Chapter 3 of this w ork [ 36 ]. f = 1 2 ( t + J s ) 2 + det ( t ) + 1 2 t 2 (4.2) It is desired to incorp orate b oth the Julic her and Hamm mo dels to dev elop a mo del whic h accoun ts for b ending, tilt and line tension. F urther w ork m ust b e done to dev elop a metho d to extract the tilt mo dulus from the MD sim ulations [ 35 36 ]. Using these energy mo dels it is desired to compute the equilibrium shap es for v arying concen trations of mixed lipid systems using parameters extracted from molecular sim ulations.

PAGE 77

REFERENCES [1] R. Lip o wsky \Budding of mem branes induced b y in tramem brane domains," J Phys II v ol. 2, no. 10, pp. 1825{40, Oct 1992. [2] B. Alb erts, D. Bra y J. Lewis, M. Ra, K. Rob erts, and J. D. W atson, Mole cular Biolo gy of The Cel l Garland Publishing, 3rd edition, 1994. [3] R. Sc hekman and L. Orci, \Coat proteins and v esicle budding.," Scienc e v ol. 271, no. 5255, pp. 1526{1533, Mar 1996. [4] J. E. Rothman and F. T. Wieland, \Protein sorting b y transp ort v esicles.," Scienc e v ol. 272, no. 5259, pp. 227{234, Apr 1996. [5] L. V. Chernomordik and M. M. Kozlo v, \Protein-lipid in terpla y in fusion and ssion of biological mem branes.," A nn R ev Bio chem v ol. 72, pp. 175{207, 2003. [6] K. N. Burger, \Greasing mem brane fusion and ssion mac hineries.," T r ac v ol. 1, no. 8, pp. 605{613, 2000. [7] R. W eigert, M. G. Silletta, S. Span, G. T uracc hio, C. Cericola, A. Colanzi, S. Senatore, R. Mancini, E. V. P olishc h uk, M. Salmona, F. F acc hiano, K. N. Burger, A. Mirono v, A. Luini, and D. Corda, \CtBP/BARS induces ssion of Golgi mem branes b y acylating lysophosphatidic acid.," Natur e v ol. 402, no. 6760, pp. 429{33, No v 1999. [8] A. Siddhan ta and D. Shields, \Secretory v esicle budding from the trans-Golgi net w ork is mediated b y phosphatidic acid lev els.," J Biol Chem v ol. 273, no. 29, pp. 17995{8, Jul 1998. [9] M. D. Matteis, A. Go di, and D. Corda, \Phosphoinositides and the golgi complex.," Curr Opin Cel l Biol v ol. 14, no. 4, pp. 434{47, Aug 2002. [10] M. Allen and D. Tildesley Computer Simulations of Liquids Oxford Science, 1987. [11] H. J. C. Berendsen, J. P M. P ostma, W. F. v an Gunsteren, A. DiNola, and J. R. Haak, \Molecular dynamics with coupling to an external bath," J Chem Phys v ol. 81, no. 8, pp. 3684{90, 1984. [12] D. v an der Sp o el et al., Gr omacs User Manual version 3.2 www.gromacs.org, 2004. 67

PAGE 78

68 [13] D. P Tieleman, S. J. Marrink, and H. J. Berendsen, \A computer p ersp ectiv e of mem branes: molecular dynamics studies of lipid bila y er systems.," Bio chim Biophys A cta v ol. 1331, no. 3, pp. 235{70, No v 1997. [14] J. N. Israelac h vili, Intermole cular and Surfac e F or c es : With Applic ations to Col loidal and Biolo gic al Systems Academic Press, second edition, 1991. [15] E. E. Ko oijman, V. Ch upin, B. de Kruij, and K. N. J. Burger, \Mo dulation of mem brane curv ature b y phosphatidic acid and lysophosphatidic acid.," T r ac v ol. 4, no. 3, pp. 162{74, Mar 2003. [16] E. E. Ko oijman, V. Ch upin, N. L. F uller, M. M. Kozlo v, B. de Kruij, K. N. J. Burger, and P R. Rand, \Sp on taneous curv ature of phosphatidic acid and lysophosphatidic acid.," Bio chemistry v ol. 44, no. 6, pp. 2097{102, F eb 2005. [17] R. Go etz and R. Lip o wsky \Computer sim ulations of bila y er mem branes: Self-assem bly and in terfacial tension," J Chem Phys v ol. 108, no. 17, pp. 7397{7409, 1998. [18] R. F aller and S.-J. Marrink, \Sim ulation of domain formation in DLPC-DSPC mixed bila y ers.," L angmuir v ol. 20, no. 18, pp. 7686{7693, 2004. [19] B. Smit, \Molecular-dynamics sim ulations of amphiphilic molecules at a liquid-liquid in terface," Phys R ev A v ol. 37, pp. 3431{3433, 1988. [20] B. J. P almer and J. Liu, \Sim ulations of micelle self-assem bly in surfactan t solutions," L angmuir v ol. 12, pp. 746{753, 1996. [21] R. P Rand and N. L. F uller, \Structural dimensions and their c hanges in a reen tran t hexagonal-lamellar transition of phospholipids.," Biophys J v ol. 66, no. 6, pp. 2127{38, Jun 1994. [22] L. Y ang, L. Ding, and H. W. Huang, \New phases of phospholipids and implications to the mem brane fusion problem," Bio chemistry v ol. 42, no. 22, pp. 6631{35, Jun 2003. [23] D. P Siegel, \The mo died stalk mec hanism of lamellar/in v erted phase transitions and its implications for mem brane fusion.," Biophys J v ol. 76, no. 1 Pt 1, pp. 291{313, Jan 1999. [24] D. P Siegel and R. M. Epand, \The mec hanism of lamellar-to-in v erted hexagonal phase transitions in phosphatidylethanolamine: implications for mem brane fusion mec hanisms.," Biophys J v ol. 73, no. 6, pp. 3089{3111, Dec 1997. [25] M. Rapp olt, A. Hic k el, F. Bringezu, and K. Lohner, \Mec hanism of the lamellar/in v erse hexagonal phase transition examined b y high resolution x-ra y diraction.," Biophys J v ol. 84, no. 5, pp. 3111{3122, Ma y 2003.

PAGE 79

69 [26] Y. Kozlo vsky and M. M. Kozlo v, \Stalk mo del of mem brane fusion: solution of energy crisis.," Biophys J v ol. 82, no. 2, pp. 882{895, F eb 2002. [27] S. W. Hui, T. P Stew art, and L. T. Boni, \The nature of lipidic particles and their roles in p olymorphic transitions.," Chem Phys Lipids v ol. 33, no. 2, pp. 113{126, Aug 1983. [28] E. Staudegger, E. J. Prenner, M. Kriec h baum, G. Dego vics, R. N. Lewis, R. N. McElhaney and K. Lohner, \X-ra y studies on the in teraction of the an timicrobial p eptide gramicidin S with microbial lipid extracts: evidence for cubic phase formation.," Bio chim Biophys A cta v ol. 1468, no. 1-2, pp. 213{230, Sep 2000. [29] T. F. Martin, \PI(4,5)P(2) regulation of surface mem brane trac.," Curr Opin Cel l Biol v ol. 13, no. 4, pp. 493{9, Aug 2001. [30] M. Jost, F. Simpson, J. M. Ka vran, M. A. Lemmon, and S. L. Sc hmid, \Phosphatidylinositol-4,5-bisphosphate is required for endo cytic coated v esicle formation.," Curr Biol v ol. 8, no. 25, pp. 1399{402, 1998. [31] T. T ak ena w a and T. Itoh, \Phosphoinositides, k ey molecules for regulation of actin cytosk eletal organization and mem brane trac from the plasma mem brane.," Bio chim Biophys A cta v ol. 1533, no. 3, pp. 190{206, Oct 2001. [32] E. A. Ev ans, \Bending resistance and c hemically induced momen ts in mem brane bila y ers.," Biophys J v ol. 14, no. 12, pp. 923{31, Dec 1974. [33] E. Lindahl and O. Edholm, \Mesoscopic undulations and thic kness ructuations in lipid bila y ers from molecular dynamics sim ulations.," Biophys J v ol. 79, no. 1, pp. 426{33, 2000. [34] F. Jlic her and R. Lip o wsky \Domain-induced budding of v esicles.," Phys R ev L ett v ol. 70, no. 19, pp. 2964{2967, Ma y 1993. [35] F. Jlic her and R. Lip o wsky \Shap e transformations of v esicles with in tramembrane domains.," Phys R ev E v ol. 53, no. 3, pp. 2670{2683, Mar 1996. [36] M. Hamm and M. Kozlo v, \Elastic energy of tilt and b ending of ruid mem branes," Eur Phys J E v ol. 3, pp. 323{335, 2000. [37] M. Hamm and M. Kozlo v, \Tilt mo del of in v erted amphiphilic mesophases," Eur Phys J B v ol. 6, pp. 519{528, 1998. [38] D. C. T urner and S. M. Gruner, \X-ra y diraction reconstruction of the in v erted hexagonal (HI I) phase in lipid-w ater systems.," Bio chemistry v ol. 31, no. 5, pp. 1340{1355, F eb 1992. [39] W. Helfric h, \Elastic prop erties of lipid bila y ers: theory and p ossible exp erimen ts.," Z Naturforsch [C] v ol. 28, no. 11, pp. 693{703, 1973.

PAGE 80

70 [40] K. Berndl, J. Kas, R. Lip o wsky E. Sac kmann, and U. Seifert, \Shap e transformations of gian t v esicles: Extreme sensitivit y to bila y er asymmetry ," Eur ophys L ett v ol. 13, no. 7, pp. 659{664, 1990. [41] U. Seifert, K. Berndl, and R. Lip o wsky \Shap e transformations of v esicles: Phase diagram for sp on taneouscurv ature and bila y er-coupling mo dels.," Phys R ev A v ol. 44, no. 2, pp. 1182{1202, Jul 1991. [42] E. A. Ev ans and W. Ra wiczl, \En trop y-driv en tension and b ending elasticit y in condensed-ruid mem branes.," Phys R ev L ett v ol. 64, no. 17, pp. 2094{2097, Apr 1990. [43] S. J. Marrink and A. E. Mark, \Eect of undulations on surface tension in sim ulated bila y ers," J Phys Chem B v ol. 105, pp. 6122{27, 2001. [44] J. Bradsha w, R. Bush b y C. Giles, M. Saunders, and A. Saxena, \The headgroup orien tation of dim yristo ylphosphatidylinositol-4-phosphate in mixed lipid bila y ers: a neutron diraction study .," Bio chim Biophys A cta v ol. 1329, no. 1, pp. 124{38, Oct 1997. [45] M. V. Berry \Liquid surfaces," Surf Sci v ol. 1, pp. 291{327, 1975. [46] J. Hansen and I. McDonald, The ory of Simple Liquids Academic Press, 1986. [47] L. E. Malv ern, Intr o duction to the Me chanics of a Continuous Me dium Pren tice-Hall, Inc., 1969. [48] J. Sonne, F. Y. Hansen, and G. H. P eters, \Metho dological problems in pressure prole calculations for lipid bila y ers.," J Chem Phys v ol. 122, no. 12, pp. 124903, Mar 2005. [49] P Sc hoeld and J. Henderson, \Statistical mec hanics of inhomogeneous ruids," Pr o c R So c L ond A v ol. 379, pp. 231{246, 1982. [50] J. H. Irving and J. G. Kirkw o o d, \The statistical mec hanical theory of transp ort pro cesses. iv. the equations of h ydro dynamics," J Chem Phys v ol. 18, no. 6, pp. 817{29, Jun 1950. [51] A. Harasima, \Molecular theory of surface tension," A dv Chem Phys v ol. 1, pp. 203{237, 1958. [52] M. Berry \The molecular mec hanism of surface tension," Phys Educ v ol. 6, pp. 79{84, 1971. [53] J. Gullingsrud and K. Sc h ulten, \Lipid bila y er pressure proles and mec hanosensitiv e c hannel gating.," Biophys J v ol. 86, no. 6, pp. 3496{3509, Jun 2004.

PAGE 81

71 [54] E. Lindahl and O. Edholm, \Spatial and energetic-en tropic decomp osition of surface tension in lipid bila y ers from molecular dynamics sim ulations," J Chem Phys v ol. 113, no. 9, pp. 3882{93, Sep 2000. [55] T. Baumgart, S. T. Hess, and W. W. W ebb, \Imaging co existing ruid domains in biomem brane mo dels coupling curv ature and line tension.," Natur e v ol. 425, no. 6960, pp. 821{4, Oct 2003. [56] P Duesing, R. T empler, and J. Seddon, \Quan tify pac king frustration energy in in v erse ly otropic mesophases," L angmuir v ol. 13, pp. 351{359, 1997. [57] S. McLaughlin, J. W ang, A. Gam bhir, and D. Murra y \PIP(2) and proteins: in teractions, organization, and information ro w.," A nnu R ev Biophys Biomol Struct v ol. 31, pp. 151{175, 2002.

PAGE 82

BIOGRAPHICAL SKETCH Eric Ma y w as b orn in New Ha v en, Connecticut, on August 29, 1978. He w as raised in Cheshire, Connecticut and attended Cheshire High Sc ho ol, graduating in 1996. He receiv ed a Bac helor of Science from Buc knell Univ ersit y in c hemical engineering in 2001. During his ten ure at Buc knell he completed summer in ternships with Neurogen Corp oration (Branford, CT) and the Connecticut Departmen t of En vironmen tal Protection (Hartford, CT). In 2001 he joined the Departmen t of Chemical Engineering at the Univ ersit y of Florida and A tul Narang's researc h group. His researc h in terests are fo cused on mathematical mo deling of biological systems. 72


Permanent Link: http://ufdc.ufl.edu/UFE0015693/00001

Material Information

Title: Molecular Modeling of Biomembrane Deformations--The Role of Lipids
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0015693:00001

Permanent Link: http://ufdc.ufl.edu/UFE0015693/00001

Material Information

Title: Molecular Modeling of Biomembrane Deformations--The Role of Lipids
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0015693:00001


This item has the following downloads:


Full Text











MOLECULAR MODELING OF BIOMEMBRANE DEFORMATIONS-THE
ROLE OF LIPIDS















By

ERIC R. MAY


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


2006


































Copyright 2006

by

Eric R. 5\.,-















I dedicate this work to parents David and Linda M.i, who instilled in me the

importance of education and encouraged me to undertake science and engineering.

I also dedicate this work to my brother Todd M_.,v who has shown me what can be

achieved through persistence and hard work. Lastly, I dedicate this to my fiancee

Jill Todisco who has loved and supported me throughout my academic career and

made sacrifices so that I could achieve my goals, for which I will always be grateful.















ACKNOWLEDGMENTS

I take this opportunity thank my mentor, Dr. Atul Narang, for proposing this

problem to me and introducing me to the field of biological modeling. I greatly

appreciate the guidance he provided me in my research, and also his advice and

wisdom on life's other matters. I also would like to thank Dr. Dmitry Kopelevich

for all the time he has spent working closely with me over the last several years.

Additionally, I would also like to thank the other members of my committee,

Dr. Ranganathan Nm,[.i',.iii.,i and Dr. Gerry Shaw, for their advice and availability.

I would like to thank Dr. Karthik Subramanian, Dr. Shakti Gupta, Jason Noel,

Ved Sharma, Ashish Gupta, Gunjan Mohan and Valere Chen for their assistance,

-ii-.-.- -lions and the enj .1'.11 scientific discussions we've had. I would also like

to thank Shirley Kelly and Debbie Sandoval for their assistance throughout my

graduate studies.















TABLE OF CONTENTS


ACKNOW LEDGMENTS .............................

LIST OF FIGURES ................................

A B ST R A C T . . . . . . . . .

1 INTRODUCTION ..............................


1.1 Specific Aims ............
1.2 Background .............
1.2.1 Role of Lipids in Biomembrane
1.2.2 Molecular Dynamics ......
1.3 Broader Impact ............


Deformations


2 PHASE TRANSITIONS IN MIXED LIPID SYSTEMS: A MD STUDY

2.1 Introduction .. ....... ... .. .. .. .. ... .... .
2.2 M ethods . . . . . . . .
2.2.1 M olecular M odel ......................
2.2.2 Sim ulation Details .. ..................
2.3 R results . . . . . . . .
2.3.1 LPA M icelles . . . . . . .
2.3.2 Phase Behavior of Pure DOPA System .. ........
2.3.3 Phase Behavior of Mixed Lipid Systems .. ........
2.4 D discussion . . . . . . . .

3 DETERMINATION OF ELASTIC PROPERTIES FOR BILAYERS


3.1 Introduction . . . . .
3.2 Bending Modulus .............
3.2.1 Curvature Elasticity .. .......
3.2.2 Molecular Model .. .........
3.2.3 Simulation Details .. ........
3.2.4 R results . . . .
3.3 Line Tension .. ..............
3.3.1 Pressure Calculation from Molecular
3.3.2 Description of Potentials and Forces
3.3.3 Calculation of Virial .........
3.3.4 Distribution of Virial .. ......
3.3.5 Calculation of Line Tension .....


Dynamics


Simulations


page

iv









3.3.6 Simulations Details ................ ... 51
3.3.7 Results ............... .......... 52
3.4 Tilt Deformations ............... ........ 56
3.4.1 Energy Model for Tilt ...... . . .. 56
3.4.2 Nature of the H11 Phase and Evidence of Tilt . ... 59

4 CONCLUSIONS ............... ........... 63

4.1 1M., I Conclusions .. ............. ..... .63
4.2 Future Directions ............... ........ 64
4.2.1 Additional Lipid Species ..... .......... 64
4.2.2 Continuum Scale Modeling ..... . .. 65

REFERENCES ................... ............... 67

BIOGRAPHICAL SKETCH .................. ......... 72















LIST OF FIGURES
Figure page

1-1 Membrane diagrams ........................... 2

1-2 Coat protein assembly and vesicle formation ..... . . 4

1-3 Potential functions ................ . .... 6

1-4 Comparision of atomistic and CG structures ..... . . 9

2-1 Detailed atomic structures and the corresponding CG models . 16

2-2 Dependence of the size and shape of LPA micelles on \!g2+] ..... 18

2-3 Dependence of the micelle asphericity on the Mg2+ concentration. .19

2-4 Bilayer configuration of pure DOPA systems . . ..... 20

2-5 Temperature induced phase transitions ................ ..21

2-6 Simulations of temperature-induced phase transitions . ... 23

2-7 Molecular geometry characterization for pure LPA systems . 25

2-8 Molecular geometry characterization in mixed lipid systems . 27

2-9 Proposed transition states in vesicle fusion .............. ..28

2-10 Stalk formation in DOPE/LPA system ............. .. 29

3-1 Comparision of theoretical and experimental vesicle shapes . 33

3-2 Atomisitic and CG representation of phosphoinositides . ... 35

3-3 Depictions of PI4P ............... ........ .. 36

3-4 Angle probability distributions ................ ...... 38

3-5 Spectral intensity for PI4P/DPPC systems . . 39

3-6 Bending modulus for PI4P/DPPC bilayer . . ...... 39

3-7 Definition of rij .................. ........... 43

3-8 Angle Orientation .................. ......... .. 45

3-9 Slab orientation for calculation of bilayer surface tension ...... ..47









3-10 Phases within a bilayer .................. ..... .. 50

3-11 Patch system with repulsive walls ................ 52

3-12 Density of pressure of CG systems ................ 54

3-13 Virial and pressure profile comparison ................ ..55

3-14 Bond potential comparison .................. ..... 55

3-15 Density and pressure profiles for patch systems . . ..... 57

3-16 Tilt and bending deformations .................. ..... 58

3-17 Construction of H1I phase ............. ..... 59

3-18 Examination of H1I phase from MD .... . ... 61

3-19 Tilt in a HII phase ........... . . ... .62

4-1 Proposed CG models of additional lipids: a) PI(4,5)2. b) DAG. . 65

4-2 Domains within a membrane, taken from ref [1] . . ... 66















Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

MOLECULAR MODELING OF BIOMEMBRANE DEFORMATIONS-THE
ROLE OF LIPIDS

By

Eric R. M.,.

December 2006

Chair: Atul Narang
M., i]r Department: Chemical Engineering

Cellular transport processes such as endocytosis and exocytosis involve the

budding of vesicles from the donor membrane (fission) and their integration into

the membrane of the acceptor compartment (fusion). Both fission and fusion are

energetically unfavorable, since they involve the formation of highly curved non-

bilayer intermediates. Consequently, these processes do not occur spontaneously;

they are under the strict control of specific proteins.

In the classical paradigm of fission, membrane deformation was thought to

be driven entirely by interactions between proteins. Evidence is now accumulat-

ing that these proteins do not act alone. They operate in concert with specific

membrane lipids, such as phosphoinositides (PI), which localize in the region

of deformation. It is of interest to understand the extent to which the localized

phosphoinositides change the mechanical properties of biomembranes.

In this work, a coarse grain molecular model is used to conduct molecular

dynamics (\ii)) simulations of mixed lipid systems. The initial part of the work

focuses on the verification of the molecular model by comparing simulation data

to experimental phase transition data for mixed lipid systems containing dioleoyl









phosphatidyle ethanolamine (DOPE), dioleoyl phosphatidic acid (DOPA) and

lysophosphatidic acid (LPA). Also, the molecular geometries of the previously

mentioned lipids are characterized and related the the mechanisms driving the

phase transitions. Additionally, MD simulations and analyses are performed

on another mixed lipid system consisting of dipalmitoyl phosphatidyl choline

(DPPC) and phosphatidylinositol-4-phosphate (PI4P). In this study two elastic

parameters of the membrane are measured from the simulations, namely, the

bending modulus and the coefficient of line tension between membrane domains of

different composition.















CHAPTER 1
INTRODUCTION

1.1 Specific Aims

The goal of this research is to better understand the mechanisms underlying

the deformations of biological membranes. Biological membranes are very complex

structures, consisting of a variety of phospholipid species, cholesterol, and proteins

as shown in Figure 1-la [2]. The membrane is an essential component of the

cell. It serves to enclose the inner contents, while also performing a number of

other functions. The membrane selectively allows passage of certain ions into

and out of the cell to maintain chemical gradients, which are essential to the

cell's survival. Also, the membrane is embedded with receptors which receive

chemical signals from the environment instructing the cell to perform specific

operations. Biomembranes, as well as being compositionally complex, are also

structurally dynamic, changing shape to perform a variety of tasks. Eukaryotic

cells are capable of endocytosis and exocytosis, the process by which material is

brought into or shipped out from the cell, respectively, or passed between different

compartments of the cell. These processes involve the formation of transport

vesicles, which are small packages with an outer lipid membrane to enclose the

contents which are transported. A schematic of the budding and fission process

is shown in Figure 1-lb. Inside of the cell there are several orangelles, some of

which are enclosed by a membrane as well, such as the Golgi body, mitochondria,

endoplasmic reticulum, and nucleus. These organelles also undergo membrane

shape transformations. Lipid molecules make up the majority of these membranes,

which contain many different lipid species. During these shape changing processes,

certain lipid species will localize in and around the region of deformation.





























Figure 1-1: Membrane diagrams: a) Example of a biological membrane [2]. b)
Depiction of budding and fission of a transport vesicle.

The hypothesis of this study is that the alteration of the chemical composition

of a biomembrane will lead to a change in the mechanical properties of the mem-

brane. This change in the mechanical property could then lead to a spontaneous

deformation of the membrane. The specific questions this study will try to answer

are the following:

1. Are certain elastic properties affected by a change in the chemical composi-

tion of the membrane?

2. What is the molecular mechanism driving the alteration of an elastic prop-

erty?

3. Does a change in chemical composition of the membrane lead to a configura-

tional change of the membrane?

With these questions answered, it will be the ultimate goal to relate the changes

in mechanical properties to the biochemistry occurring inside the cell when these

deformations occur. In essence, it is desired to develop a complete model starting

from chemical reactions inside the biochemical pathways effecting compositional









changes inside the membrane, effecting mechanical changes of the membrane,

leading to a continuum level model which can predict global membrane shape.

1.2 Background

1.2.1 Role of Lipids in Biomembrane Deformations

Lipid molecules are the main component of biomembranes. These lipids are

amphiphilic in nature, having one part which has an affinity for hydrophilic envi-

ronments (head) and one part which has an affinity for hydrophobic environments

(tail). The tail(s) consists mainly of saturated hydrocarbons, but unsaturated

bonds do commonly exist as well. There can be either a single tail or two tails

connected to the head. The head and tail are connected through a glycerol group.

The head group often contains a phosphate group, as well as other polar groups,

which often carry a net charge. In biological systems the environment in which

the cells lives is aqueous, as is the interior of the cell. This naturally leads to the

bilayer configuration of the lipid molecules to separate the internal contents of

the cell from the environment, as shown in Figure 1-la. In this configuration, the

heads are submerged in both the interior of the cell and the exterior environments

while the tails form a hydrophobic region at the interior of the membrane itself.

This interior region of the membrane is very important because it provides a region

for proteins to reside since they often contain hydrophobic regions.

The mechanism by which biomembranes change shape is the focus of this

work. It had originally been postulated that the driving force behind the formation

of vesicles from biomembranes was driven by protein interactions. It is believed

that interactions between membrane bound proteins and external proteins known

as coat proteins, especially the protein complex clathrin, act to increase curvature

of the membrane and lead to bud formation [3, 4]. This is an incremental process

in which the addition of more coat proteins continuously deform the membrane

until a bud is formed and it pinches off from the original membrane. The coat then










----N3 .

coated 1
clahuvrt in




receptor adaptin nakel rqns~po'



cargo molecule

COAT ASSEMBLY BUD VESICLE UNCOATING
AND CARGO SELECTION FORMATION FORMATION

Figure 1-2: Depiction of coat protein assembly and the formation of a transport
vesicle. [2].


will dissociate from the vesicle to allow for fission with the target destination, see

Figure 1-2.

The process of fusion and fission of vesicles from membranes is an energetically

unfavorable one. The intermediate states can be highly curved, requiring more

energy to form than available from typical thermal fluctuations of the system [5].

However, evidence has been accumulating that the proteins do not act alone, but

act in concert with the lipid molecules themselves [6]. Some of the key lipids are

ones which exist in relatively low concentrations inside the membrane, but play

a crucial role in these morphological changes. Specifically, phosphoinositides,

phosphatidic acid, lysophosphatidic acid and diacylglycerol have been implicated

as essential lipids in the processes of fusion and fission of transport vesicles [7-9].

To this point, it has been unclear as to what function these lipids play in altering

the configuration of biomembranes. It is the goal of this work to study the effect

different lipids have on bilayers in the absence of proteins.









1.2.2 Molecular Dynamics

Molecular dynamics (\I1)) is a computational tool used to study molecular

;-1I-' I i both in equilibrium and non-equilibrium regimes. In MD, all particles are

treated classically, where the particles interact through different bonded and non-

bonded potentials. The bonded interactions can consist of bond length, bond angle

and dihedral potentials, which are typically modelled by a harmonic potential.

The non-bonded interactions are due to van der Waals and electrostatic forces.

For computational purposes, only pair potentials are typically considered, taking

3-body and higher term interactions becomes very time consuming because they

involve summing over triplets (or greater) of molecules [10]. However, the 2-body

potentials have shown to produce very good results in comparison to experimental

measurements.

To determine the particles position and velocities, Newtons equation's of

motion must be solved
d2ri
mr = Fi, i = 1...N (1.1)

where i refers to the particle number, running from 1 to N, r is the position vector,

m is the mass of the particle and F is the force vector acting on the particle. The

force is given by the negative derivative of the potential with respect to position,

Fi, The algorithm typically used to solve these equations is the Verlet

Algorithm [10].


v(t At t F-(t)At (1.2)
2 2 m

r(t + At) r(t) + v(t + )At (1.3)
2

It can be seen Eq. 1.2 updates the velocity at time t + from the velocity at time
At
t- and the force at time t. The position is then updated, as shown in Eq. 1.3,
Lt
using the velocity at time t + 2 and the position at time t to get the position at

time t + At.













'-- Repulsive vES

0 -------- r



Attractive


(a) (b)

Figure 1-3: Potential functions: a) Lennard-Jones potential b) Electrostatic poten-
tial

The non bonded forces are often specified by a Lennard-Jones potential of the

form VLJ(ry) = 4c(((/riy)12 (a/ri)6) as shown in Figure 1 3a. The parameter

e gives the depth of the potential well and a is the hard sphere radius. Since

the potential flattens out at large r, a cut-off radius is specified to increase the

computational efficiency. If two particles are separated by a distance greater than

the cut-off radius, it is assumed they have no influence on each other. Special care

must be taken when using cut-off radii to avoid discontinuities in the potential

functions. To avoid this, a shifting function can be added to the potential in order

for it to smoothly approach a zero slope at the cut-off radius. In addition to the

Lennard-Jones potential, an electrostatic potential must be accounted for if there

are charged particles in the -iv-1. i, The form of the electrostatic potential is

VEs(Tij) = qi where gq is the charge on molecule i and eo is the permittivity of

free space, and is shown in Figure 1 3b. Again, a cut-off radius can be chosen for

the electrostatic potential or a more complex method of accounting for the long

range interactions known as Ewald summation may be cm-ipl.-,-d.

In molecular dynamic simulations periodic boundary conditions are often

used. This is to eliminate the edge effects which would be introduced by the

presence of walls. When periodic conditions are used, it allows molecules to exit









from the simulation box on one side and then reappear, entering on the opposite

side of the box. The molecules will interact with the nearest periodic image of

the other molecules; this is known as the minimum image convention. In order

for the periodic simulation to represent the macroscopic system the length of the

simulation box should be at least 6a where a is the Lennard-Jones parameter

representing the hard sphere radius of the particles [10].

It is often desired for the simulation system to maintain a constant tempera-

ture and pressure, so that the results are comparable to experiments. One method

to control temperature is through the Berendsen algorithm, which couples the

-v-1, II, to a heat bath using first-order kinetics. In this weak coupling scheme the

-v-1' I, responds exponentially when it deviates from the bath temperature, To,

given by Eq. 1.4, where T is the coupling time scale. This is achieved by rescaling

the velocities of all particles by a factor A given in Eq. 1.5. TT is related to T by

Eq. 1.6, where Cv is the total heat capacity of the system, k is the Boltzmann

constant and Ndf is the total number of degrees of freedom [11].

dT To T
(1.4)
dt T
At To
A 1 1) (1.5)
TT T
S= 2CTT/Ndf k (1.6)


The weak coupling method of temperature coupling is easily implemented,

but it does not represent a correct statistical mechanical ensemble. An alternative

method which does probe a correct canonical ensemble, is the extended system

approach, which introduces a "Temperature Pi-i ii" to control the temperature,

thereby introducing and additional degree of freedom to the system. The equations

of motion become modified as shown in Eq. 1.7, where is the heat bath variable

or friction parameter. This parameter is dynamic and is governed by it's own

equation of motion Eq. 1.8. Q is defined by Eq. 1.9 where TT is coupling parameter









which defines the timescale of oscillation of kinetic energy between the heat bath

and the system [12].

d2ri Fi dri (1.7)
dt2 mn dt
(T To) (1.8)
dt Q

Q (1.9)
472

Similar coupling schemes exist for pressure coupling as well. In pressure

coupling, the box vectors are scaled to adjust the volume and therefore the pressure

of the system. This can been achieved through weak coupling where the pressure

obeys first-order kinetics or through an extend system approach, where the

equations of motion of the particles are modified and the box vectors must obey

their own equation of motion. Pressure coupling can be done either isotropically or

anisotropically, which allows for deformation of the simulation box from the initial

aspect ratios. When dealing with lipid bilayers, the choice of anisotropic coupling

to equal pressure in all directions leads to a zero surface tension of the bilayer,

which is often desirable. The study of lipid bilayers using MD, is well established,

where numerous studies have shown good agreement with experimental data

regarding the density of these systems, diffusion coefficients and permeability [13].

The molecular models used in MD can be atomistic, in which all molecules,

bonds, bond angles and partial charges are specified (possibly excluding hydrogens)

or coarse grained in which multiple atoms are defined by a single interaction site.

There are methods of deriving coarse grain potentials from, atomistic potentials

one of which is the Iterative Boltzmann Inversion [? ]. In this method an initial

potential is guessed using a known radial distribution function, g(r), such that

Vo(r) = -kTln(g(r)). Simulating the system with this potential will then produce a

different radial distribution, go(r) and the potential is corrected by adding the term














C Palmitoyl a
Fatty Acid




C) EC
(a) (b)

Figure 1-4: Comparison of atomistic and coarse grain structures a) DPPC b) Wa-
ter

-kTln(-t') to the initial potential guess. This process can be iterated several

times until reasonable convergence of the radial distribution function is achieved.

To date several coarse grain models have been published, including one by Marrink

et al. which includes models for several lipid species [? ]. In Marrink's model,

a coarse grain particle typically represents 4-6 atoms; shown in Figure 1-4 are

depictions of the atomistic and coarse grain representations of DPPC, a common

biological lipid, and water. The labelling of the coarse grain particles corresponds

the the force field designation of the particles published by Marrink [? ]. The main

advantage of using a coarse grain model is that it reduces the degrees of freedom

in the system and also allows for a larger time step, which make this model much

more computationally efficient than an atomistic model.

1.3 Broader Impact

The foundation of the work is the advancement of scientific knowledge. The

work has been conducted without a direct application intended. However, it is

hopeful that the findings of this research could lead to more directed studies. Biol-

ogy is a field which often offers qualitative and sometimes speculative explanations

of the phenomena. Taking an engineering approach to biological problems offers







10

the possibility of quantitative and mechanistic explanations. The motivation be-

hind this research is understanding the mechanical implications of observed lipid

localization in biomembranes. Furthermore, the goal of this work is to understand

how changes in chemical composition of membranes affect elastic properties and

what are the underlying molecular mechanisms driving these changes in the mem-

brane properties. These effects can then be compared with the forces generated

inside the cell including pressure, osmotic pressure, polymerization of filaments and

protein interactions to see if chemical changes within the membrane are playing a

significant role in the deformation of membranes.















CHAPTER 2
PHASE TRANSITIONS IN MIXED DOPE-DOPA AND LPA-DOPA LIPID
SYSTEMS: A MOLECULAR DYNAMICS STUDY

2.1 Introduction

Cellular transport processes such as endocytosis, exocytosis, and subcellular

trafficking involve the budding of vesicles from the donor membrane (fission) and

their integration into the membrane of the acceptor compartment (fusion). Both

fission and fusion are energetically unfavorable, since they involve the formation of

highly curved non-bilayer intermediates. Consequently, these processes do not occur

spontaneously-they are under the strict control of specific proteins.

In the classical paradigm of fission, membrane deformation was thought to be

driven entirely by interactions between proteins. According to this model, fission

is initiated by the recruitment of certain cytosolic proteins called coat proteins to

the membrane. The stepwise assembly of coat proteins then incrementally deforms

the membrane into a spherical shape, thus producing a vesicle which is ultimately

released along with its encased coat.

Evidence is now accumulating that these proteins do not act alone. They

operate in concert with particular membrane lipids, namely, phosphoinositides,

diacylglycerol, and phosphatidic acid. One function of the lipids is to recruit

proteins involved in membrane deformation. It has been shown, for instance, that

phospholipase D, an enzyme required for vesicle formation in the Golgi complex,

is recruited to the membrane by diacylglycerol. More importantly, however, lipids

appear to be directly involved in deformation of biomembranes. Early experiments

showed that enzymes that acylate lysophosphatidic acid (LPA) to phosphatidic

acid (PA) strongly activate vesicle formation. Since LPA and PA are one- and









two-tailed phospholipids, it seems plausible to hypothesize that LPA is cone-

shaped (A) and PA has the shape of an inverted cone (V). Now, it is known

from condensed-matter pir, -i. that amphiphilic molecules self-assemble into

. .- regates whose morphology is intimately linked to the shape of the individual

molecules [14]. Specifically, V- and A-shaped molecules form structures with

positive and negative curvatures because these geometries reduce bending stresses

within the membrane. This has led to the hypothesis that the shape change

induced by LPA-acyltransferases is ultimately driven by the distinct shapes of LPA

and PA molecules.

Recently, Kooijman et al. performed experiments in a cell-free system to

ascertain the molecular shape of LPA and PA [15]. To this end, systems of LPA,

PA, and mixtures of these lipids with various other lipids were studied in a variety

of environments in which pH, temperature and divalent cation concentrations were

varied to mimic the conditions in the cytosol and Golgi complex. In the mixed lipid

experiments, LPA or PA was mixed with either dioleoyl phosphatidylethanolamine

(DOPE) or dioleoyl phosphatidylethanolamine (DEPE), and again the pH, tem-

perature, and ionic conditions were altered. Using 31P-NMR measurements the

:-v,- I II phase was distinguished for a given set of conditions. It was observed that

a system of a given lipid composition would transition to a different phase by

altering environmental conditions or temperature. These phase transitions involve

a change in curvature of the of the equilibrium structure, often from a relatively

flat lamellar phase to a high curved inverted hexagonal phase. LPA and PA are

of particular interest because in the Golgi body the acylation of LPA to form PA

has been shown to induce the formation of a vesicle from the Golgi membrane [7].

LPA and PA are very different from each other from a geometrical perspective.

The spontaneous curvature of LPA was determined to be highly positive, desiring









to form micelles, while PA has a negative spontaneous curvature, under ]lr, -i-

ological conditions [16]. It is believed that this transformation from LPA to PA

induces stress in the membrane leading to the destabilization of the bilayer and the

formation of a bud.

The goal of this work is to develop a model for the lipids LPA and PA to

be used in molecular dynamics studies. A coarse grain (CG) model developed by

Marrink et al. is capable of reproducing experimentally observed phase transitions

for a phospholipid system consisting of dioleoyl phosphatidylcholine (DOPC)

and dioleoyl phosphatidylethanolamine (DOPE) [? ]. Using Marrink's force

field parameters, models for PA and LPA were constructed and simulations were

conducted and compared to experimental results of Kooijman et al. [15]. In

addition, it is desired to characterize the molecular geometry of these molecules

to gain insight into the mechanism driving the phase transitions. This work is

a preliminary study to verify CG models for the lipids LPA and PA. With the

assurance that the models are reliable, further study is planned to fully investigate

the mechanism driving membrane deformation in mixed lipid systems.

2.2 Methods

2.2.1 Molecular Model

In order to accurately evaluate the elastic (i.e., macroscopic) properties using

microscopic simulations, it is necessary to perform simulations of systems contain-

ing thousands of lipids for hundreds or thousands of nanoseconds. The application

of a detailed atomistic model to such length and time scales is extremely time-

consuming even with modern computing systems and therefore limits the scope of

the systems that can be investigated within a reasonable time frame. Therefore,

simulations of lipid systems often napl-:v less detailed models, such as the Brown-

ian dynamics simulations [? ], dissipative particle dynamics [? ], and coarse-grained

molecular dynamics (CGMD) models [17? ? ? 18]. In this work, we use a CGMD









model, which allows one to perform simulations on near atomic length scales, while

reducing the computational time by orders of magnitude and therefore allowing

exploration of the systems on sufficiently large time- and length-scales. The CGMD

models approximate small groups of atoms by a single united atom (bead). Several

such models have been introduced and applied to simulations of various complex

molecular systems [17, 19? ? ? ? 20]. Development of coarse-grained molecular

models is a subject of active ongoing research [? ? ? ]. In this work, we will use

the coarse-grained model proposed by Marrink et al. [? ]. This model has been

shown to yield good agreement with experiments and atomistically detailed simu-

lations for such quantities as density and elasticity of pure lipid bilayers. Within

this model, for example, 4 methyl groups of lipid tails are represented by a single

hydrophobic coarse-grained bead and 4 water molecules are represented by a single

spherical polar bead. This model also provides several types of beads that can be

used to model various groups of atoms (such as glycerol and phosphate) within

lipid headgroups. The interactions between non-bonded beads are modeled by the

Lennard-Jones potential, whereas the interactions of bonded beads are modeled by

the harmonic bond and angle vibration potentials. Also, electrostatic interactions

between charged beads are taken into account.

The parameters for these potentials as well as models for several molecules

were previously published by Marrink et al. [? ] and the same parameters were

used in this work. For the nonbonded interactions a cut off distance of 1.2 nm was

used. To avoid discontinuities of the potentials at the cut off point the LJ potential

is smoothly shifted to zero between 0.9 nm and the cutoff point. The electrostatic

potential is also shifted to zero beginning at 0.0 nm, to represent the effect of ionic

screening.

In this work, we will use the coarse-grained model for water, metal ions, and

DOPE lipids proposed in by Marrink et al. [? ]. In addition, we have developed









coarse-grained models of the additive lipids using the methodology of Marrink

et al. [? ] for mapping different functional groups within a molecule to coarse-

grained beads with specific properties. The detailed atomistic structures and the

corresponding course-grained representations of some of these lipid molecules

are shown in Fig. 2-1. The structure and CG model for the lipid LPA shown in

Fig. 2-la consists of 7 CG particles. Each bead is classified with respect to the

-','-I1 iI published by Marrink et al. [? ]. The glycerol-ester linkage is modelled

by a single bead in LPA with the classification Nda, indicating a non-polar

particle capable of both hydrogen bonding and donating. Glycerol in double

tailed molecules is typically modelled with two beads, however for LPA the choice

to use only a single bead model for glycerol was based upon a closer representation

of the true mass of the molecule and the qualitative agreement of the simulations

with experiments using this model. An attempt to model the glycerol of LPA with

2 beads, as in DOPA, for a total of 8 beads in the model, produced a simulation

result which did not agree with the experiments. In the absence of ions the self-

assembled structure of the 8 bead LPA formed a worm-like micelle, where the

experiments predict spherical micelles.

The structure and model for dioleoyl phosphatidic acid (DOPA) is shown in

Fig. 2-lb. The CG model consists of 13 CG particles. The glycerol ester linkage is

modelled in the typical way using two beads of type Na, indicating it is non polar

and can act as a hydrogen acceptor. LPA has the ability to be a hydrogen donor

due to the hydroxyl group on the glycerol, where in DOPA the oxygen is bound to

the ester linkage. Both of these lipids are modelled with oleoyl fatty acid tails (18

carbons, 1 double bond), the double bond causing the kink in the chain. Lastly, a

model for a divalent cation is introduced using a single bead of type Qda, carrying

a charge of 1.4 e. The charge is reduced from 2.0 e to account for the hydration

shell surrounding the ion.










Phosphate Phosphate
Ndi Glycerol
aC 'N T -N-Glyercol
o c- Oleoyl
C Oleoyl Fatty Acid
C Fatty Acid
SCC



(a) (b)

Figure 2-1: Detailed atomic structures and the corresponding coarse-grained mod-
els of considered lipids: (a) LPA and (b) DOPA. Mapping between different groups
of atoms and the coarse-grained beads is shown for some groups. The types of the
lipids (C = hydrophobic, P = polar, etc.) and their potential parameters are the
same as in Marrink et al. [? ]


2.2.2 Simulation Details

Molecular dynamic (\Ii)) simulations were conducted, using the GROMACS

MD simulations package [? ], version 3.2.1. In all simulations anisotropic Berendsen

pressure coupling was used with a reference pressure of 1.0 bar. All systems were

also coupled to a temperature bath using Berendsen coupling. The time step used

for the simulations was 40 fs. Periodic boundary conditions were imposed in all

three directions.

2.3 Results

In this study molecular dynamic simulations were used to investigate the

polymorphic phase behavior of lipid systems. Experimental results have indicated

that lipid systems are sensitive to pH, temperature, ion concentration and relative

water hydration levels [15, 21, 22]. We use the coarse-grained (CG) molecular

model developed by Marrink et al. and discussed above. Recall that in this model

four water molecules are represented by a single CG particle. Therefore, in what

follows when referring to the amount of water in the simulation ;-',-li i1i the amount









of "real" water molecules will be given. Also, it is known that the dynamics of

coarse-grained models is typically faster than that of a real system. This is due

to the smoothing of the potentials. For the specific model considered here, the

CG time is 4 times faster than the real time. Therefore, the times reported for

the simulations are given as the effective time, which is calculated by scaling the

observed simulation time by a factor of 4.

2.3.1 LPA Micelles

Kooijman et al. studied the shape change of LPA .,.- -regates in response

to increasing concentrations of Mg2+ (I 2+/LPA = 0, 0.2, and 1.0) [15]. The

experiments were performed at 37C and a pH of 7.2. They observed that in the

absence of Mg2+, LPA forms micellar ..- -'regates. However, as the Mg2+/LPA ratio

is increased, the shape of the LPA .,.-'- regates changes from a micellar configuration

to a system of bilayer disks.

To test the ability of the CG model to reproduce these observations, we

conducted self-assembly simulations of the LPA lipids (shown in Fig. 2-la), water

and different concentrations of divalent ions (representing Mg2+). The simulations

were performed at 37C. The simulation box was randomly populated with 500

molecules of LPA, 192 molecules of water per LPA molecule, and the appropriate

amount of Mg2+ (determined by the Mg2+/LPA ratio that was being studied).

Initial random dispersions of LPA and Mg2+ in water quickly self-assembled into

micelles and equilibrated within 800 ns.

Figures 2-2a-d show the equilibrium configurations achieved at the end of the

simulation for the three Mg2+/LPA concentration ratios used in the experiments.

The simulations show that in the absence of Mg2+, six small micelles are formed

(Fig. 2-2a). When the Mg2+/LPA ratio is increased to 0.2, five micelles are formed

(Fig. 2-2b). If the Mg2+/LPA ratio is increased even further to 1.0, the larger

micelles .-i-'- regate to form a single .. -.- regate. Figures 2-2c,d show the top and


















(a) (b)


I _____ I i_ _-'_ I
(c) (d)

Figure 2-2: Dependence of the size and shape of LPA micelles on \!_22+]: (a) Mo-
lar ratio I\!2+]/[LPA] 0. (b) Molar Ratio [\!2+]/[LPA] 0.2. (c) Molar ratio
\!2+]/[LPA] 1.0, side view (d) Molar ratio [\!2+]/[LPA] 1.0, top view. In
these plots the water molecules and ions were removed for clarity.

side views of this single ..- .- regate. It is clear from these figures that the .i.- egate
has a disk-like configuration. Therefore, as the concentration of the ions increases,
the micellar size increases and micelles acquire disk-like shapes, which is consistent
with experimental 1P-NMR measurements of Kooijman et al. [15].
To quantify the shape change of the LPA .,.-'- regates shown in Fig. 2-2, we
calculated their asphericity, which is defined as

y3 2f t 2)2
A 2j (Ri R (2.1)
j-1= I(R2)2
where R, are the radii of gyration. The asphericity characterizes the extent to
which the shape of an ..- .- regate deviates from that of a sphere. The asphericity of

the .I.- regates was calculated by averaging over the data corresponding to the final













I o
0 03 -
0025 -




0 01 -


0 0 02 04 06 08 1
Mg+21LPA

Figure 2-3: Dependence of the micelle asphericity on the Mg2+ concentration.


80 ns of the simulation. Fig. 2-3 shows that as the concentration of Mg2+ increases,

so does the asphericity of the lipid .I. -. regates.

It is of interest to understand the molecular basis for the shape changes of the

LPA ...:i--regates shown in Fig. 2-2. According to the shape-structure concept of

lipid polymorphism, the shape of an .-. -.regate is determined largely by the shape

of the component molecules. Therefore, we expect that the change from spherical

to disk-like micelles indicates a change in the molecular geometry of LPA lipids.

2.3.2 Phase Behavior of Pure DOPA System

Kooijman et al. also studied the phase behavior of pure DOPA systems. It

was observed that a DOPA system formed a lamellar phase at neutral pH and a

temperature of 310 K both in the absence and equimolar presence of Ca2+. Self-

assembly simulations were conducted to match this result. 500 DOPA molecules

and 17500 CG waters were randomly placed in a simulations box, for the system

with Ca2+, 500 of CG divalent cations were also added to the system. The systems

were simulated at 310 K for 800ns and in both cases a bilayer phase was observed.

This indicates the model for DOPA and divalent cation can match the experimental

observations.





















Figure 2-4: Bilayer configuration of pure DOPA system observed after 800ns of
simulation. a) No divalent ions. b) Divalent ion concentration equal to lipid con-
centration.

2.3.3 Phase Behavior of DOPE-LPA and DOPE-DOPA Mixed System

We also investigated phase transitions between the lamellar and inverted

hexagonal phases in the mixed DOPE-LPA and DOPE-DOPA systems. The exper-

imental data [15] summarized in Fig. 2-5 indicate phase transitions from lamellar

to inverted hexagonal phases within the bilayer systems as the temperature of the

-v-1i ii increases.

Kooijman et al. [15] observed that when the temperature is increased, pure

DOPE undergoes a transition at Te 275 K from the lamellar (L,) phase to the

inverted hexagonal (HII) phase (curve labeled A in Fig. 2-5). Upon addition of

10 mol% LPA to DOPE, the phase transition temperature increases to Ts 300 K

(O curve of Fig. 2-5). The addition of 10 mol% DOPA to DOPE results in slight

increase in the phase transition temperature to Tz 280 K (O curve of Fig. 2-5).

Therefore, the experimental data show that the change of the minor mixture

component from LPA to DOPA has a destabilizing effect on the lamellar phase by

lowering the phase transition temperature by a 20 K. We perform simulations to

check if this transition was mirrored by the model.

Recently, Marrink et al. [? ] studied the phase behavior of the CG model for

pure DOPE. They found, for a hydration level in the range of 9-12 waters per lipid,











100 -A- Pure DOPE (Exp)
90 9 -- DOPE/LPA (Exp)
\ DOPE/DOPA (Exp)
80 -
S-A- Pure DOPE (Sim)
70
60 -I DOPE/LPA (Sim)
60 1 DOPE/DOPA (Sim)
50 i
1 401
II 1
30 'E i
"I

10 A
0-
240 260 280 300 320 340
Temperatuer(K)

Figure 2-5: Temperature-induced phase transition between lamellar and inverted
hexagonal phases in the mixed lipid systems: Comparison of experimental data [15]
and simulations for DOPE/DOPA (0) DOPE/LPA systems (0) and pure DOPE
(A). The experimental data are shown by open symbols and the simulations results
are shown by closed symbols. The lamellar phase corresponds to 100% bilayer and
the inverted hexagonal phase corresponds to (' The composition of the mixed
-1-I-. 1, is 90% DOPE and 10% of a minor lipid (LPA or DOPA). The simulation
data for mixed systems are from the current work, where as the simulation data of
the pure DOPE system was taken from Marrink et al. [? ].


the phase transition temperature to be above 287 K (A curve of Fig. 2-5). This

is quite close to the experimentally observed phase transition temperature. To

determine the phase transition, temperature simulations were performed at several

temperatures in which the initial configuration of the system consisted of four

pre-simulated bilayers stacked on top of each other. In some of the simulations,

the bilayers would form stalks between the layers, which would then grow and the

(H11) phase would spontaneously form.

We performed similar simulations to check if the CG model could capture

the phase behavior of the DOPE-LPA system. The initial conditions for our

simulations were taken to be lamellar phases prepared following the simulation

protocol of Marrink et al. [? ]. First, we used self-assembly of initially randomly









dispersed lipids in water to prepare a lipid bilayer. The self-assembly simulations

were performed in a system with relatively low lipid concentration which leads to

bilayer formation within 200 ns, the system was allowed to run for and additional

600ns to allow for equilibration. The self assembly simulation was conducted at

a water to lipid ratio of 80 and a temperature of 280 K. We then dehydrate the

bilayer down to a water lipid ratio of 15 and stack 2 identical dehydrated bilayers

on top of each other. In an attempt to promote formation of the HII phase, the

bilayers were given a perturbation so they would resemble a sine wave in one

direction. The direction of the perturbation was alternated between the bilayers

so that the peak of one bilayer would coincide with the valley of the neighboring

bilayer. This configuration was then energy minimized using the steepest descent

method in GROMACS for 50000 steps. This energy minimized configuration was

then copied in the normal direction to obtain a system consisting of 4 bilayers.

This initial configuration is shown in Fig. 2-6a. Each bilayer consists of 450 DOPE

molecules and 50 LPA molecules. The simulations were performed at various

temperatures between 250 K and 300 K.

Examples of final configuration of the simulations for the DOPE-LPA system

are shown in Fig. 2-6b and c. At temperatures at or below 265 K the lamellar

phase remains stable on the timescale of simulations (4.0pts), see Fig. 2-6b. On

the other hand, at T = 275 K and above, the same system exhibits an instability

which leads to formation of a stalk and then ultimately the inverted hexagonal

phase is reached. The final configuration for the simulation at 300 K is shown in

Fig. 2-6c. All of the configurations shown in Fig. 2-6 were generated by copying

the simulation cell 3 times in the lateral direction to better illustrate the geometry

of the systems.




















(a) (b)


Figure 2-6: Simulations of temperature-induced phase transition between lamellar
and inverted hexagonal phases in the mixed DOPE-LPA lipid systems: (a) Initial
perturbed configuration for DOPE-LPA simulations. (b) Final configuration of
DOPE-LPA system at 250 K lamellarr phase). (c) Final configuration of DOPE-
LPA system at 300 K Lipid headgroup beads are shown by black spheres, tail
beads are shown by gray spheres, and the water beads are shown by white spheres.

Figure 2-5 shows a comparison between the experimental system and the

simulation results. The model system displays a phase transition temperature

above 265 K whereas the experimental system transitioned around 300 K.

For the DOPE-DOPA simulations, a bilayer consisting of 120 lipid molecules

(10% DOPA, 90% DOPE) was simulated by self-assembly at water to lipid ratio

of 120 and at 240 K for 80ns. These small bilayers were then copied laterally to

form a bilayer consisting of 480 lipids. A perturbed bilayer system was then created

in the same manner as for the DOPE-LPA system, again at a reduced water to

lipid ratio of 15. This perturbed system was then simulated for up 4.0kts at varying

temperatures between 240 K and 310 K. At a temperature of 265 K and greater the

-v-I-' i n spontaneously formed the H11 phase. At temperatures of 240 K and 250

K the system remained stable in the bilayer phase, but at temperature at or above

265 K the systems converted into the H1 phase.

Figure 2-5 shows that the simulations are in qualitative agreement with

the experimental data. The model DOPE-DOPA system undergoes a phase

transition at a temperature approximately 15 K lower than the phase-transition









temperature for DOPE-LPA system, compared to an approximate 20 K shift

for the experimental system. However, the precise values of the phase transition

temperature predicted by the simulations are different from those observed in the

experiments. The experiments indicate phase transition temperatures of 280 K

and 300 K for the DOPE-DOPA and DOPE-LPA systems, but the corresponding

simulated values transitions temperatures are above 250 K and 265 K, respectively.

As can be seen, the current coarse-grained model does not capture the exact

phase transition temperature, it does capture important qualitative features of the

systems, such as the destabilizing effect on the lamellar phase due to the change

from LPA to DOPA.

2.4 Discussion

We also studied the molecular mechanism involved in temperature-dependent

phase transitions. It is evident that as the temperature increases, so do the bilayer

fluctuations, thus facilitating the possibility of a phase change. However, for the

hexagonal phase to be stable, it is necessary that the molecular geometry also

change. To investigate the molecular shape changes in temperature-dependent

phase transitions, the simulations of the DOPE-LPA system were analyzed.

Molecular packing theory states the microstructure of the lipid ..-'--regate should

be representative of the mean molecular geometry of the system [14]. In order to

verify this hypothesis, the packing parameter was measured, P = v/lcao, where

v is the volume occupied by the hydrocarbon tails, l~ is the length of the tails in

the direction normal to the water lipid interface and ao is the area occupied by the

head group.

For the pure LPA system, a transition between several spherical micelles to

a single bilayer disk is observed upon the addition of an equimolar concentration

of divalent ions. In Fig. 2-7a, the packing parameter is shown for the LPA system

at the varying divalent ion concentration. A significant increase in the packing










15
S a0(nm2)
Tall Length (nm)
Tall Volume (nmi)


E i


os 05-*



0 02 04 06 08 1 0 02 04 06 08 1
Divalent lon/Lipid Molar Ratio Divalent lon/Lipid Molar Ratio
(a) (b)

Figure 2-7: Molecular Geometry Characterization of LPA in a pure LPA system:
a) LPA Packing Parameters. b) LPA measurements of tail volume, tail length and
head spacing.


parameter occurs upon the increase in ion /lipid ratio from 0.2 to 1.0, which

coincides with the transition from micelle to disk. To understand why the packing

parameter is (1i. ii-if i. the individual measurements of ao, lI, and v are shown

in Fig. 2-7b. Both ao and l1 decrease upon addition of the ions, while v remains

relatively constant. The volume of the tails should remain constant since the

temperature is fixed; the change in ao can be understood because the divalent ions

act to screen the negatively charged LPA headgroups from each other, causing an

effective decrease in head group size.

At 250 K and 265 K the DOPE-LPA system remained in the lamellar phase,

and at 285 K and 300 K the system transitions into the H11 phase; these state

points were used for the packing parameter analysis. Data was analysed from the

final 80ns of the simulation after which the system had formed its final phase and

equilibrated. The packing parameter for DOPE, LPA and the mean value for the

-v,-1 iI are shown in Fig. 2-8a. Discontinuous jump is observed between the La

points and the H11 points, indicating that the shape of the molecules becomes

more like an inverted cone (V). In Fig. 2-8b the measurements of ao, 1 and v are









shown. It is observed that l1 is the parameter which is changing the most. Since l1

is decreasing with increasing temperature, the tails sample configurations farther

away from the equilibrium configurations, essentially reducing their extension

normal to the water lipid interface. The theory states that values of the packing

parameter above 1 lead to inverted phases, however it can been seen that even the

lamellar phase points are above 1. This is due to the assumption in the calculation

of the head group area, ao, that DOPE and LPA have equal size head group

areas. The head group area was estimated by calculating the area of the water

lipid interface and dividing through by the number of lipid molecules at that

interface. There was no accurate way to individually measure each of the species

ao's separately, therefore, for the purpose of the packing parameter calculation it

was assumed both species had the same ao. From this analysis it is evident that the

molecular mechanism driving the phase transition is the configuration of the tails.

The increased temperature imparts a greater energy on the tails to deviate away

from relatively straight configurations to a more bent and splayed structure.

The transition states which a lipid systems samples as it transitions from a L,

to a H11 has been studied both experimentally and theoretically[23-27]. Biologi-

cally this transition is of particular interest for understanding the mechanism of

vesicle fusion and understanding the how antimicrobial peptides destroy bacterial

membranes [28]. A mechanism proposed by Siegel, proposes the bilayers first form

a stalk phase and then the outer ]i "ii"].r'. is contract and form a contact know as

the trans monolayer contact (TMC), Fig. 2-9a. The TMC phase involves formation

of hydrophobic void spaces which Kozlovsky et al. states are energetically unfavor-

able, and a different intermediate stalk structure is proposed [26]. This alternative

stalk model is shown in Fig. 2-9b; in order to reach this configuration the lipid tails

must tilt away from normal to water-lipid interface to prevent the formation of

void space.







































50 255 260 265 270 275 280 285 290 2e5 3(
Temperature (K)

(a)


* aO (nm2)
Lc (nm)
* Single Tall Volume (nd)


*


B 6


250 260 270 280
Temperature (K)

(b)


290 300


0 0


* *


a (nm2)
Lc (nm)
Single Tall Volume (nnd)











10
245 250 255 260 265 270 275 280
Temperature (K)

(d)


Figure 2-8: Molecular Geometry Characterization of DOPE, LPA and DOPA in
mixed system: a) DOPE-LPA i-1 -I 1i1 packing parameters. b) DOPE-LPA system
measurements of tail volume, tail length and head spacing. c) DOPE-DOPA system
packing parameters d) DOPE-DOPA measurements of tail volume, tail length and
head spacing.


Temperature (K)


235 240











~f-----------



p -








(a) (b)

Figure 2-9: Proposed transitions in vesicle fusion a) Siegal stalk and TMC
model [23]. b) Kozlovsky and Kozlov stalk model [26].


To verify which of these mechanisms our simulations underwent the a system

of 90% DOPE and 10% LPA was analyzed at 300 K. The formation of a stalk is

depicted in Fig. 2-10a-d. It is clear form examining the simulation images that the

stalk structure is very similar to that proposed by Kozlovsky et al.. There is no

evidence of TMC of void region, furthermore the water pores are long narrow ovals

are shown by Kozlovsky et al. [26].

Using a coarse grain model for molecular dynamics has proved to be a qual-

itatively reliable tool for studying lipid systems. In real system there are many

parameters and a high level of complexity. However, in our model system the

-' I i,, is greatly simplified, but the simulations are capable of matching the trends

shown experimentally. The key result from these simulations is the shift in the

transition temperature, not the temperature itself. The the transition temperature

is lowered by changing the minor component from LPA to DOPA both experi-

mentally and in our system. Therefore there is a range of temperature, in between

the phase transition temperature, which the conversion of LPA to DOPA or vice

















(a) (b)


(c) (d)


Figure 2-10: Stalk formation in a DOPE/LPA system: a) 96 ns. b) 112 ns. c) 128
ns. d) 144 ns. The images were zoomed in on just 2 of the bilayers to better il-
lustrate the stalk formation. Water was removed for clarity, the head groups are
colored black and the tails are grey.

versa can induce a transformation. This model can also be used to quantify the

molecular geometry of the molecules.

The simulation results give good qualitative results, but there is a discrepancy

between the reported experimental transition temperature and the observed values

in the simulations. There is potential to produce better quantitative agreement by

changing some parameters in the model. Increasing the hydration of the system

will stabilize the bilayer by increasing the separation of the stacked bilayers,

making a transition more difficult, and should shift the transition temperature to a

higher value. Also, the Lennard-Jones parameters can be altered as Marrink et al.

did in their study of DOPE and DOPC to make the bilayer phase more stable [? ].

However in this study we chose not to modify the force field to remain consistent

with previous simulations.







30

One limitation of this coarse grain approach is that the pH of the system

cannot be explicitly modelled. The charges on the molecules can be modified to

reflect the desired pH, but H+ ions are too small to be modeled in this approach.

Kooijman et al. do observe an effect of changing the pH on the phase behavior of

these lipid systems, without changing the net charge on the lipids.















CHAPTER 3
MOLECULAR MODELING OF KEY ELASTIC PROPERTIES FOR
INHOMOGENEOUS LIPID BILAYERS

3.1 Introduction

In biology, membranes can be very non-uniform in their structure and com-

position. The major component of the biological membrane are lipid molecules,

while they also consist of cholesterol and proteins [2]. One type of molecule

which is implicated in many biological shape changing process are phosphoinosi-

tides [9, 29-31]. In this study the effect of introducing a phosphoinositide molecule,

phosphatidylinositol-4-phosphate (PI4P) to a homogeneous dipalmitoyl phos-

phatidyl choline membrane is studied. Two important elastic constants will be

measured.

The bending modulus is a measure of the amount of energy required to change

the curvature of a membrane. The bending modulus has previously been measured

both experimentally and from simulations for homogeneous systems [32, 33]. In

this study, molecular dynamics simulations of mixed PI4P and DPPC bilayers

are simulated and analyzed to determine the effect of PI4P concentration on

the bending modulus. Understanding the effect of concentration on the bending

modulus is important since in biological membranes phosphoinositides will localize

in one region of the membrane. This region of localization is often the site a

membrane deformation in the form of a protrusion, budding or endocytosis.

A second key elastic parameter is the line tension existing between lipid

phases. If a tension exists between phases this can act to cause budding or en-

docytosis as the tension force acts to minimize the size of the interface, resulting

in an out of plane deformation [1, 34, 35]. The line tension existing between a









homogeneous DPPC section and an inhomogeneous PI4P/DPPC membrane section

will be determined, again with the use of molecular dynamics simulations.

A third elastic parameter will also be investigated. When lipid molecules in a

monolayer tilt away from the lipid-water interface normal direction and energy is

associated with this type of deformation. The tilt modulus characterized the energy

associated with this type of deformation. In this work, it will be shown that tilt

deformations in a hexagonal phase observed through MD simulations are consistent

with theoretical and experimental predictions [36-38].

3.2 Bending Modulus

3.2.1 Curvature Elasticity

A continuum model for describing the shape of bilayer vesicles was first

proposed by Helfrich [39]. The Helfrich model predicts a quadratic relationship

between the mean curvature, C, of a homogeneous bilayer and the free energy of

the bilayer per unit area, f, in Eq. 3.1,



f = (C C)2 + K (3.1)

where K is the bending modulus, Co is the spontaneous curvature of the bilayer, K

is the Gaussian modulus and K is the Gaussian curvature. The mean and Gaussian

curvatures are defined below:


C = \ + K2

K = K,12


where Ki and K2 are the principal curvatures of the bilayer surface.

The minimum energy configuration of the bilayer can be solved by limiting

the solutions to spherically symmetric shapes for a given volume, V, and surface

area, A, of the enclosed bilayer. This is a constrained minimization problem, where

the object is find the curve which minimizes the total energy of the bilayer given















ii m







Figure 3-1: Comparison of theoretical and experimental vesicle shapes. [40].

the constant volume and surface area restrictions. The total functional is given by

Eq. 3.2 where, E, and P are Lagrange multipliers.


Fb = fdA + ZA + PV (3.2)

All functions can be parameterized by the arc length of the curve defining the

axisymmetric surface. A set of Euler-Lagrange shape equations can then be

derived and numerically solved to yield the minimum energy shape given a set

of constraints. The phase diagram has been mapped out by Seifert et al. and

several biologically relevant shapes were determined to be minimum energy

configurations [40, 41]. Figure 3-1 shows a comparison of experimental and

theoretical vesicles shapes as the temperature of the system is altered thereby

altering the surface area and volume of the vesicles. Budding shapes, as well as, the

biconcave shape of the red blood cell are observed.

In solving the shape equations, the integral of the Gaussian curvature, over

the surface area for closed vesicles remains constant and therefore need not be

considered in variational problem. This leaves only the mean curvature term, and

it is desired to know the value of the bending modulus. This has been determined









experimentally [42] and also from molecular dynamic simulations [33, 43] and good
agreement has been shown between them. In order to determine K from simulations

the thermal fluctuations of an equilibrium bilayer can be analyzed. If we assume

Co = 0 and neglect the Gaussian bending term the free energy per unit area is

given by
1 1
f = K(h + hy,)2 (V2 (h(r))2
2 2
where h is the height of the bilayer and the subscripts denote derivatives with

respect to a direction, where the bilayer normal is the z direction. Taking the

Fourier transform of the surface, Eq. 3.3, and the inverse transform, Eq. 3.4 yields

and expression for h(r) in terms of the wave vector q.



h(q) h(r)-iq.rdr (3.3)

h(r) = h(q)eiq.r (3.4)


By taking the Fourier transform of the surface, the bending modes become decou-

pled and the total energy of the surface is given by Eq. 3.5.

F fdA q4h(q) 2 (3.5)
q

From the equipartition theorem, the energy of each mode can be can be related to

the Fourier coefficients shown in Eq. 3.6


kBT
< Ih(q) 2 > k (3.6)
KCq4
where kB is Boltzmann's constant and T is the temperature of the system. The

bending modulus can then be determined from the simulation data by fitting

< lh(q) 2 > versus q4











Phosphate


Choline


X- Palmitl. _. .- -Glycerol
C C amityl ,Na Na
FatFattycAcid id
S= C C Palmitoyl
c Fatty Acid
C) (c) C C
S(c) C C (



(a) (b)

Figure 3-2: Atomistic and coarse grain representation of phosphoinositides: a)
PI4P. b) DPPC.

3.2.2 Molecular Model

DPPC is one of the major components of most biological membranes. PI4P is

an important molecule in membrane systems because it has been implicated to play

an essential role in Golgi membrane functions and it is a precursor to PI(4, 5)P2 [9].

From a geometric perspective, simulating systems which contained PI molecules

was intriguing because of the large head group due to the inositol ring and attached

phosphate groups. The phosphate group on PI4P carries a net -2e charge which

can cause the effective head size to be even larger if it is in the presence of other

negatively charged lipids. The CG model and atomistic structure of PI4P and

DPPC are shown in Figure 3-2, it can be seen that the inositol ring is modeled as a

single polar bead due to it's ability to form hydrogen bonds and its high solubility

in water.

Experimental evidence has shown that the orientation of the inositol ring of

PI4P, when in the presence of DPPC, is not normal to the bilayer, but becomes

bent [44]. Figure 3-3a shows how the inositol ring gets pulled down towards the





















S0) a = 180 '



(a) (b)

Figure 3-3: Depictions of PI4P a) Atomistic representation in the presence of
DPPC. b) Head group of the CG model shown with the imposed equilibrium an-
gles, no equilibrium value is imposed on al.


plane of the bilayer, rather than extend out into the aqueous phase. It has been

speculated that this structural feature of PI4P is caused by the electrostatic

interaction between the negatively charge phosphate group of PI4P and the

positively charged choline group of DPPC. Figure 3-3b shows the CG model of

the head group of PI4P and the equilibrium angles. To try to mimic the structural

features of PI4P no potential has been imposed on the angle ci, where the angles

a2 and ca have equilibrium angles of 180Oand 1200, respectively. Since no angle

potential is acting on the ac angle this will allow for the top phosphate group to

move and adopt a conformation due solely to the intermolecular interactions and of

course the bond length potential.

3.2.3 Simulation Details

Self-assembly simulations of DPPC and PI4P were conducted at relative PI4P

concentrations of 0' ., 5' .; 10% and 2' i` This was done by populating a simula-

tions box with total of 200 lipid molecules at the proper relative concentrations

and 28 waters per lipid (7 CG W's), which was the same concentration used in the









experimental determination of the head group orientation by Bradshaw et al. [44].

The systems were then simulated for 800ns at a temperature of 298 K, in all sys-

tems a bilayer was formed. The bilayer system was then copied laterally four times

to create a bilayer system consisting of 800 lipid molecules. The 800 lipid system

was then simulated for 1.6ps. The first 400ns were considered an equilibration,

and the final 1.2pts were used for in the analysis of the system. In both the self-

assembly and bilayer simulations Berendsen temperature coupling and anisotropic

pressure coupling (1 bar) were used as well as periodic boundary conditions in all 3

directions.

3.2.4 Results

To determine the structure of the PI4P head group the al angle was calculated

for all systems. The probability distribution of al and a2 angles in the PI4P head

group for each of the systems containing PI4P is shown in Figure 3-4a and b,

respectively. All of the systems display the same trend with a peak at 0 W2.1 rads

(1200) for al and at 0 = 7 for a2. Since a2 is normal to the bilayer water interface,

it can be determined that the upper phosphate group is being pulled down towards

the bilayer-water interface as indicated by the al angle data. Since the phosphate

group is partially lying in the interface, it acts to increase the head group size of

PI4P. This data indicates that the model is qualitatively reflecting the true nature

of the PI4P head group in a PI4P/DPPC system.

As described in Section 3.2.1 the relationship between the fluctuations of the

bilayer surface and the bending modulus can be determined, Eq. 3.6. To obtain

h(q) for a surface, a well defined surface is required. An imaginary lateral grid

was imposed and each bead was assigned to a bin based upon it's lateral position

(x, y). For each bin the average normal (z) position of all molecules in the bin
was determined. The surface was defined by assigning each bin center the average

z position of all molecules in the bin giving a point-wise function z = f(x, y),











003 008
-5% P14P 007- -5% P14P
0025 ---10% P14P ---10% P14P
002 ....... 20% P14P 006 ....... 20% PI4P
0 05
S0015 004
2 2
6 003-
001
002-
0005
S005001

0 05 1 15 2 25 3 35 0 05 1 15 2 25 3 35
Angle (radians) Angle (radians)
(a) (b)

Figure 3-4: Probability distribution for systems of PI4P/DPPC at concentrations
of 5' .; 10% and 20% a) alangle. b) a2angle.


representing the surface. At each point in time, during the analysis, a 2D Fourier

transform of the surface was taken as described by Eq. 3.3. Each (h(q))2 was

averaged and then log(< lh(q) 2 >) was plotted versus log(l|q), these spectral

intensity plots for each system are shown in Figure 3-5a-d. The longer wavelength

modes, |q| > 1.0 nm, were used to fit a line by the least-squares method. A function

can be fit to the data of the form


log(< |h(q) 2 >) = -4log(lq) + C1


where the intercept value, C1, is equal to log(kBT/K)(Eq. 3.6) allowing for cto be

determined. The uncertainty in C1, acwas propagated to give the error in V as

shown in Eq. 3.7.
d(kbTJ)
aK = d"rC1 ( )ac1 = -KTac (3.7)
dC1 d1

The bending modulus and error at varying concentrations of PI4P is shown

in Figure 3-6. At the time of this report there was only 200ns of data for the pure

DPPC system, accounting for the larger error in that system. The values reported

here are in the same range as previously reported for simulated and experimental

-,1- I. 1 -[33, 42].





























-05 0 05 1 15 2 25
Log q(1/nm)


-05 0 05 1 15 2 25
Log q(1/nm)

(b)


-1 -05 0 05 1 15 2 25
Log q(1/nm)


-1 -05 0 05
Log q(1/nm:


1 15 2 25


Figure 3-5: Spectral intensity for PI4P/DPPC
PI4P. c) 10% PI4P. d) 20% PI4P.



x 10-20



95


085
9-

851 T


5 10 15
% PI4P


Figure 3-6: Bending modulus of PI4P/DPPC bilayer
PI4P


systems a) Pure DPPC. b) 5%


20 25


at varying concentrations of









3.3 Line Tension

The existence of a tension force existing between lipid phases can be a driving

force for shape change of biomembranes. From a molecular prospective a tension

force between lipid phases can be generated due to the differing interaction energies

between the molecules of the respective phases. The lipid phases maybe a pure

raft or can be a region of the membrane where a specific lipid has localized. A key

membrane lipid in many biological processes are phosphoinositides which are known

to play a role both inter and intracellular trafficking [7-9, 31]. At the outer cellular

membrane there is a cytoskeleton which plays a role of structural stability for the

cell as well as effecting shape change. At the outer membrane the mechanism of

shape change is quite complex involving the interaction between lipids, proteins

and the cytoskeleton. However, intracellular compartments, such as the Golgi body

and the endoplasmic reticulum do not have a cytoskeleton on their interior, but

are still able to mediate shape changes. Therefore the mechanism of shape change

at the intracellular membrane is less complex and is location where this work is

focused.

3.3.1 Pressure Calculation from Molecular Dynamics Simulations

Pressure in a fluid is made up of two-components: a kinetic part due to the

molecular motions and a static part due the the intermolecular potentials [45].

Pressure is a force acting on the unit volume, so the kinetic part can be viewed

as the momentum transfer across an imaginary surface. The static part is due to

particles on opposites sides of the imaginary surface exerting a force on each other

due to the intermolecular potentials. The kinetic part of pressure can be calculated

from the temperature of the system or from the molecular velocities as shown in

Eq. 3.9 and Eq. 3.10.












P PK+ Ps (3.8)

PK (3.9)

2 1
PK = EK MiVi 0 Vi (3.10)


Where PK is the scalar or hydrostatic pressure, it is related to the pressure

tensor, PK, by PK = trace(PK)/3. N is the number of particles in the system,

V is the volume of the system, EK is the kinetic energy tensor, mi is the mass of

each particles and vi is the velocity vector of each particle. The calculation of the

kinetic pressure from molecular dynamics is trivial. If the velocities of each particle

are stored to an output trajectory file the kinetic pressure can be calculated easily

from Eq. 3.10.

The static pressure calculation can be more complex, it is related to the virial

of the system [46].

2
Ps Vir (3.11)
N
Vir- 2 r F (3.12)
i

The virial is a 2nd order tensor shown in Eq. 3.12, where ri is the vector describing

the particle position and Fi is the force vector acting on particle i. The complexity

of the virial calculation comes in determining the forces acting on the particles.

In theory these forces can be due to external and internal potentials, but in this

work only internal forces were present. The potentials are an input to the model

and the forces on each bead are calculated from the negative derivative of the

potential with respect to particle position, F(i) The potentials present in the

coarse-grain model simulations are described in Section 3.3.2.










3.3.2 Description of Potentials and Forces

Lennard-Jones Potential

The Lennard-Jones(LJ) potential is a non-bonded potential which represents

the the long range attractive van der Waals force and a short range repulsive part

due to overlap of the electron clouds. For the simulations conducted in this work,

the LJ potential is shifted beyond a certain radius so the forces go smoothly to zero

at the cut-off radius. In the simulation described in this work, the cut off distance

was TLJ, = 1.2 nm and the shift distance was rLJ, = 0.9 nm. The LJ interactions

between nearest bonded neighbors are excluded. The potentials, VLJ, and forces,

FLJ(i) acting on bead i for the unshifted and shifted parts of the potential are

given by


C12 C6
r12 r6


12C12 [I

-6C6 1
VLJ
01


rij < TLJ,


A12 (T
3 (rij
3


FLJ )3

TLJ)


B12 (ri_ FLJ)4 D12]

B6 (r -LJ) -D6
4


FLJ, < ij < rLJ,


rij > TrLJ


and


-12C12 6C6 rij
ri13 +rj rij

12C12 A +- A12(rj rLJj2 + B12(Ti
F(i)LJ =
-6C6 + A6(rij FLJ 2 + B6(r

0


rij < rLJ,


'LJ, )3]

LJ )3] ^j


rLJ, < rij < rLJ,


rij > rLFj























Figure 3-7: Definition of rij


where ri is the vector pointing from bead i to bead j as described in Figure 3-7,

rij is the magnitude of ri, and

Ae (a + 4)rLJ, (a + 1)rLJ,
A a+2 29
FJi(F LJC bLJe)
B (a + 3)rLJ, (a + 1)rLJ,
FJ f(cLJc rLJ,)3
1 A B
D a- 0(LJC 'rLJ,) rLJ, rLJ.
arc, 3 4

The force on bead j is equal and opposite to the force on bead i, F(j)LJ =

-F(i)LJ. This is also true of the 2-body potentials for electrostatic and bonded

interactions described below.

Electrostatic Potential

The electrostatic potential (ES) is a non-bonded potential which acts between

beads with a net charge. The potential is shifted smoothly to zero at the cut off

radius, TESc = 1.2 nm starting from rES, = 0.0 nm. The cut off and the shifting

and of the electrostatic potential is done to mimic the distance dependant screening

phenomenon [? ]. Electrostatic interactions between nearest bonded neighbors are

excluded. The electrostatic potential and resultant force are described below.











i_ qij 1 A, i3 43 4 (rij < rES c
VES 0 1ir7 3' rEVS -3 rtij rESJ D rD j < KrES0

= 0 rij > ESC

F (i) Es = [ + A(r^ -rESS)2 B ru -VESS)] 3 Vi K VES0
= 0 rij > TESC

VES is the electrostatic potential, F(i)Es is the force due to the electrostatic

potential acting on bead i, qi is the charge on bead i, co the permittivity of free

space and c, is the relative dielectric constant, which was taken to be c, = 20. The

constants A1, Bland D1 are describe by the same equations for the constants in the

LJ potential.

Bonded Potential

The bonded potential, VB, is a harmonic potential which corrects for devia-

tions from equilibrium bond lengths. The bond potential and resultant force, FB,

are described below.


1
VB = I bond Tij ro)2



F(i)B = bond(rij ) Tri
rij
The bonding force constant for a single bond is kbond = 1250 -kJ and the

equilibrium bond length is taken to be ro = 0.47 nm.

Angle Potential

The angle potential can be used to constrain three consecutively bonded

beads. The potential is a harmonic function of the cosine of the angle made about

the central bead. The angle orientation and resultant forces are described in

Figure 3-8. The angle potential, Vo, and forces, Fo are described below,










F(j)
i k



F(i) j F(k)




Figure 3-8: Angle Orientation




Vo = ko(cosO cosOo)2
2
F(i)o ko(cos 0 cos 0) (rji 3rik)rji _rk
I r iFjk Fjifjk]
ko(cos 0 cos 00) rr r (jk)
F(j) ri + rk (rji" rk) 2 + 2-
jifjk [h ji

F(k)o ko(cos0 cos0) (rji rjk)rjk rji
rjkfrji Fjijkfk

where the angle force constant for single bonds is ko 25 k and the equilibrium

angle is0o = r for. The cosine of the angle can be determined by taking the dot

product of rji and rjk, cos(6) rjrjk [47].
rjirjk
3.3.3 Calculation of Virial

To calculate the virial of the system, the particle positions and total force on

each particle must be known as described in Eq. 3.12. For the 2-body potentials

the force is a function of separation distance between two interacting particles. It

is computationally efficient to calculate the virial due to each interaction Vir(i, j),

Eq. 3.13 and then sum over all interactions Eq. 3.14.











1 1
Vir(i, j) [ri 0 FF + rj 0 Fi] [rij Fu-] (3.13)
2 2
Vir2-Body Y > [rij 0 Fij] (3.14)
i
Fij is the force on i due to j, and Fi, is equal and opposite to Fij from the

conservation of linear momentum. The calculation of the virial due to the angle

potential can be achieved by summing over all angles in the system, Eq. 3.15.



Viro > [r, 0 F, + rj 0 Fj + rk 0 Fk] (3.15)
Angles

3.3.4 Distribution of Virial

In this work the pressure profile, or the spatial variation in the pressure

tensor, is of interest. This variation in components of the pressure tensor is the

source of surface tension and line tension. To determine pressure as a function of

space the pressure must be determined locally within the system. To do this the

system is section off into slabs and the pressure is calculated in each slab. The

slab is taken to have the same dimensions in as the simulation box in two of the

dimensions, but in the direction normal to the interface of interest the dimension

will be much smaller. The slab orientation for the calculation of a bilayer surface

tension is described by Figure 3-9. The slab is moved through the z-direction and

the pressure is calculated at each slab location so P(z) is obtained. To calculate

P(z)K Eq. 3.10 is used in which the sum is taken over only those molecules which

lie within the slab.

The calculation of P(z)s is not well defined since the location of the force

acting between particles is ambiguous [48]. The static pressure at a point R can

be written in a general form by Eq 3.16, where Coi is any contour connecting the

particle position, ri, with and arbitrary point, Ro [49].









Z



1 Bilayer



/_ Slab
Y





Figure 3-9: Slab orientation for calculation of bilayer surface tension



PS(R, t) =- [VV({rY})] j dl(R- 1) (3.16)

The most natural choice for the contour is the Irving-Kirkwood contour which is a

straight line between interacting particles [50]. Other contours have been proposed,

most notably the Harasima contour, but in this work only the Irving-Kirkwood

contour is used [51]. A generalized expression for the static pressure for an m-body

potential can be derived from Eq. 3.16 is described by 3.17.




P (Body {Vj~ V :} C6(R 1,// (3.17)
klPk

where < j > represents all "clusters", which means a group of interacting particles,

and k and I are the indices of particles within the cluster. The average pressure in

each slab is what is desired to be known, so Eq. 3.17 must be integrated over the

dimensions of the slab,



1 "Lxz oLY Zi+AZ
pab -(Z, t) = ] dx f dz P, bd (R, t) (3.18)
Sm-body VolAz 0 0 Zi Sm-body









where Z, is the ith slab, AZ is the slab thickness, and VolAz is the volume of the
slab. The choice of contour must be defined to evaluated the integral, here we use a
straight line contour (Irving-Kirkwood) connecting interacting particles


(3.19)

(3.20)


13 r + A(r r A [0, 1]


S6(R 1li' -r 6(R (rk + Arjkjl)dA (3.

where rjkl = rj rjk Combining Eqs. 3.17, 3.18 and 3.21, an expression which
can be evaluated for the static pressure in a slab is derived


P (Zi,t)


1)


mVolAz {V, m } r f (zk zi, Z, )



where Zjkis the z position coordinate of particle k in j cluster, and
Lx r'Ly Zi+AZ Yl
f(zj, z, Z) = dxdydz 6(R (rj + ris))dA (3.22)
Jo tJO Jo


f Lx LY zf+Az
0 0 1 zi


6(R (rjk + rj))dxdydz


1 R = rj + ArJkj Re VolAz

0 R / rj + Arijkj


(3.23)


so R must lie on the rJkJ and be within slab i for the spatial integration to be non-
zero. Since the slab spans the system in the x and y directions only the particle z
positions are critical in evaluating the pressure in a slab. Eq. 3.23 can be re-written









with the use of the Heaviside step function.

rL, rLy, rZi+AZ
/ LY 6(R (r, + r,))dxdydz =(zk + A(zi zj) Zj)x
0 JO 0 Zi

E(Z, + AZ + z, + A(zi zik))

1 x>O
O(x) = x
0 x < 0

f(z, Zh, Z ) j- (z + A(z A z ) Zi)x

E(Z' + AZ + zkj + A(zi Zjk))dA

So depending on the location of the interacting particles with respect to slab i

location, f(zjk, z,, Z,) will take on different values. If both particles lie within slab

i, f = 1; if both particles lie outside of slab i and no part of rjkji lies within slab

i, f = 0; if both particles do not lie with slab i, but part of rjkjI does pass through

slab i, then
Zi.
SJkJl

where rz7 is the portion of riJikwithin slab i.

3.3.5 Calculation of Line Tension

Surface tension is calculated from the difference between lateral and normal

components of the pressure tensor [52].



7- (PN(Z) PL(z))dz
-00
The z-direction is normal to the interface, PN and PL are the normal and lateral

pressure, respectively. This deviation in the isotropy of the pressure tensor can be

due to both density deviations and changes in configuration energy. The integration

is theoretically taken from -oo to +oo since the pressure becomes isotropic away

from the interface and does not contribute to the integration.









Y
X Slab


Bilayer

'^aT~w Tao


I .Z 1I, Z

Figure 3-10: Phases within a bilayer

To calculate line tension, the pressure must be localized to the surface defined

by the membrane and then the differences between the components of this 2-

dimensional pressure tensor will result in line tension. Figure 3-10 describes the

geometry of phases embedded within a bilayer.

In this situation pressure would be calculated as a function of the z-direction,

with slabs spanning the box in the x and y directions.



1 yLZ
a(t) (PZZ, PXX(z, t))dz (3.24)
Nint Jo

P2D(z, t) P(z, y, t)dy

a is line tension, P2Dis the 2-dimensional pressure tensor and Nint is the number of

interfaces in the system, since the line tension per interface is desired. Since P(z, t)

is what will be calculated and it is already averaged over the y-direction we can see

that P2D(z, t) = P(z, t)Ly and by combining this with Eq. 3.24 an expression for

line tension as a function of P(z, t) is derived.


L Lz
(t) = / (PZZ( t) PXX(z, t))dz (3.25)
int o0


a P a









3.3.6 Simulations Details

The same CG models for PI4P and DPPC are used in simulations for the

calculation of line tension. The 800 lipid simulations generated for the bending

modulus calculations as described in Section 3.2.3 were used to assemble the line

tension systems, referred to from here forward as a patch system. Figure 3-10

describes the desired configuration for the line tension calculation in which the

a-phase represents a homogeneous pure DPPC section and the /-phase represents

an inhomogeneous section of PI4P and DPPC. The initial bilayer sections were

taken from the output of 1.6 ps for the 800 lipid systems. A new topology was

then generated by placing the inhomogeneous section in between to homogeneous

sections. To do this, the inhomogeneous section had to be modified so the median

bilayer height was the same as the homogeneous section. Also, it was required that

both systems have identical x and y dimensions, so which ever of the two systems

was larger was trimmed to match the dimensions of the smaller section. The patch

-v,- I in was generated with both a 5% PI4P and 20% PI4P inhomogeneous section.

These patch systems were then energy minimized for 50000 steps using the steep

energy minimization method. Additional water was then added to the system

to bring the water level to 22 CG W per lipid and again the system was energy

minimized for 50000 steps. In order the keep the sections separated and to prevent

the PI4P molecules from diffusing out into the homogeneous section a repulsive

wall was constructed, which acted only on the PI4P molecules. The wall spanned

the simulation box in the x and y directions. The "dummy" atoms which made up

the wall were separated with 0.5 nm spacing in both the x and y directions. The z

coordinate location of the wall at high z value was determined by adding 0.4 nm to

the maximum z coordinate position of all PI4P beads. Likewise the wall at low z

value was determined by subtracting 0.4 nm from the minimal coordinate position

of all PI4P beads. The dummy atoms were kept at a fixed position throughout the





















Figure 3-11: Patch system with repulsive walls


simulations. The interaction potential between the dummy atoms and the PI4P

beads was equivalent to the repulsive part of the LJ potential for polar nonpolar

interactions


C12 if interacting with PI4P
1Dummy -
0 if interacting with bead other than PI4P

where C12 = 8.3658 x 10-2. This potential is shifted and cut-off with the same radii

and in the same manner as all other LJ potentials as described in Section 3.3.2.

Figure 3-11 is a simulation snapshot of the patch system configuration, water

was removed for clarity. The black beads extending through the bilayer are the

dummy wall atoms. The section between the walls is a mixture of PI4P and DPPC,

while the outer sections of the bilayer are pure DPPC. These systems were then

simulated for 400ns using Berendsen pressure and temperature coupling, set to

1 bar anisotropically and 323 K respectively. The pressure coupling was turned

off after the initial 400ns simulation and was run for an additional 400ns with no

pressure coupling.

3.3.7 Results

Comparison to Atomistic Results

Several studies have been conducted on the calculation of pressure profiles

and surface tension of homogeneous bilayers from atomistic molecular dynamics









simulations [48, 53, 54]. The potentials used in atomistic simulations differ from

CG simulations in many ways. Atomistic simulations typically have more poten-

tials, often including a dihedral potential which is not included in the CG model.

In atomistic simulations the electrostatics can be very different because of the

presence of partial charges on the majority of atoms, where in the CG model most

beads are modeled as neutral beads. Additionally, the electrostatics may be not

be cut off in atomistic simulations to allow for long range interactions. In the CG

model the LJ potential partially accounts for the electrostatic. Lastly the coeffi-

cients of all of the potentials will vary from CG to atomistic simulations making

the comparison between virial contributions difficult. However, if the CG model is

a reliable model the pressure profiles between CG and atomistic model should be

comparable.

The pressure profile for a pure DPPC bilayer was calculated for the same

system used for the bending modulus calculation (800 lipids, 7 CG W/lipid, 298 K,

anisotropic pressure coupling). The analysis was done on 1.2ps of data after 400ns

of equilibration. This data was compared to results obtained for a similar atomistic

system study conducted by Sonne et al. [48]. The system used by Sonne consisted

of 72 DPPC molecules simulated at 300 K with a ratio of 28 waters per lipid (same

concentration as our system). The pressure profile was calculated as a function of

the bilayer normal direction, the z direction for convenience. In the analysis done

here the system was divided into 60 slabs, in the work by Sonne the system was

divided into 70 slabs [48]. In Figure 3-12a the density profile is plotted, where the

units are number of beads per slab. In Figure 3-12 the diagonal components of the

static pressure tensor are plotted as a function of z position. It can be seen that

the normal component remains constant through the bilayer.












200 -300-
DPPC Epp z


50 -400 -

100 -500 A






S 3 2 1 0 1 2 3 4 -3 -2 1 0 1
Position (nm) Z position (nm)
(a) (b)

Figure 3 12: Density and pressure of CG systems and a function of bilayer posi-
tion: a) Density profile. b) Pressure profile.


In Figure 3 13 the virial contributions and pressure profiles are compared for a

CG and atomistic models. The CG virial contributions are quite different from the

atomistic virial contributions, especially in the magnitude of the numbers.

This difference in magnitude can be explained by the difference in the po-

tentials and resultant forces. To illustrate the difference, the bonded potential is

examined for both the CG and atomistic models. Figure 3 14 plots the bonded

potential for both models. For a deviation from the equilibrium bond length on the

order of k T the forces for each model were calculated and an order of magnitude

difference is observed. However, the Pressure profiles in Figures 3 13c and d show

good agreement between the atomistic and CG data.




1.3 x 103 Atomisitic
FB
7.9 x 101 CG














WAO, ...


-30 -20 -10 0 10 20 30
rAi

(b)



P P
IK
-H
IK SE(IK)











20 0 20 40
[LA]

(d)


Figure 3-13: Virial and pressure profile comparison:
istic virial profile [48]. c) CG static pressure profile.
file [48].


Bond Energy


Atomistic
CG


5000


-05 -04 -03 -02 -01 0 01 02 03 04
Bond Displacement (nm)


Figure 3-14: Bond potential comparison


4000

2000




-2000


a-


a) CG virial profile. b) Atom-
d) Atomistic pressure pro-









Line Tension for Patch Bilayer Systems

Calculation of line tension was conducted for patch systems with inhomoge-

neous sections containing 5% and 20% PI4P. The analysis was done on the final

80ns of data for the simulations done without pressure coupling. In the analysis

the pressure was calculated as a function of the z direction, normal to the phase

interface, as described in Figure 3-10. The system was divided into 100 slabs to

calculate the pressure profiles. The pressure difference was calculated between the z

and x directions (in plane components) so that the line tension could be calculated

using Eq. 3.25. The number density of PI4P beads and the pressure differences for

the 5% and 20% patch systems are shown in Figure 3-15.

It can be observed that the 5% system shows no distinct deviation in the

pressures near the phase boundary. However, in the 20% patch system there is a

sharp change in the pressure difference profile at the phase boundary indicating the

presence of a line tension. The line tension was calculated for both systems, for the

5% system a = -2.3 x 10-11N and for the 20% system a = 2.0 x 10-10N. The value

for the 5% can be attributed to noise and poor statistics, however the 20% systems

shows a sharp deviation from isotropic behavior near the boundary. The value

obtained here is higher than reported for an experimental system (9 x 10-13N) and

theoretical estimate (10-11)N [1, 55].

3.4 Tilt Deformations

3.4.1 Energy Model for Tilt

In the pioneering work by Helfrich, on the description of equilibrium vesicle

shapes, he choose to neglect the energy associated with tilt of lipid molecules away

from the bilayer normal direction [39]. In the case of flat bilayers the lipids on

average do not tilt away from the normal direction, however other lipid phases

including the H11 and intermediate stalk phase (as described in Chapter 2) do

require the molecules to tilt and this energy must be accounted for when predicting








































/V\^\/"VMAA\


Z position (nm)


Z position (nm)
(d)


Figure 3-15: Density and pressure profiles for patch system:a) 5% PI4P density
profile. b) 20% PI4P density profile. c) 5% PI4P pressure difference profile. d) 20%
PI4P pressure difference profile


--111














N n



N

t
(a) (b) (c)

Figure 3-16: Tilt and bending: a)Definition of tilt. b) Non-uniform tilt deforma-
tion on a volume element of ]ii",i'1.,', r [36]. c) Bending deformation on a volume
element of ni, 'ii'1.'.w r [36].


stability of phases. Hamm and Kozlov introduced a energy model for tilt and

used this to predict the stability of the lamellar and inverted lipid phases as well

as geometric properties of the inverted phases [37]. The definition of tilt and the

elastic energy per unit area per monolayer are described by Eq. 3.26 and Eq. 4.2,

respectively and graphical representation of tilt is shown in Figure 3-16a.



t N (3.26)
n-N
1 1
f 2(tl + J,)2 + kdet(t) + ot2 (3.27)

The tilt vector, t, is perpendicular to N, so if we take N to point in the

z direction, t will have components in the x and y directions. The tilt tensor

is defined by the derivatives of tilt t Vt, where V = where = y.

t trace(t=) V t, Jjis the spontaneous curvature of the ni..i ii1...''. r, ,,,nl1.,,i the

tilt modulus. It is interesting to note that the same elastic parameters, K and K,

describe the energy of both tilt and bending deformations. This is explained in a

later paper by Hamm and Kozlov and can be visualized by the Figures 3-16b and





















(a) (b) (c)

Figure 3-17: Construction of H11 phase: a) Original construction of hexagonal
phase unit cell [37]. b) Multiple original construction cells shown to illustrate the
hydrophobic void region [38]. c) Recent hexagonal phase unit cell construction [37].

c [36]. Figure I ;i.]1, describes a non-uniform tilt deformation, where tilt changes

across the element. Figure 3-16c describes a bending deformation, it can be seen

that the shear and stretch of this element a essentaillyessentially the same in both

deformations, which is allows for both to be characterized by the same parameters.

3.4.2 Nature of the H11 Phase and Evidence of Tilt

It had long been assumed that the structure of the H11 phase consisted of

a circular hydrophilic interface and a hexagonal hydrophobic interface between

unit cells as shown in Figure 3-17. This construction allows for a mismatch

between the interior hydrophilic boundary and the outer hydrophobic boundary.

In this construction, the lipid tails must either stretch to fill the corners of outer

boundary or a hydrophobic void region will exist as illustrated in the center of

the three cells in Figure 3-17b; both situation are energetically unfavorable. To

alleviate this l I-I .11- i. i i. i,y -", as it is know, the construction of the unit

cell was re-examined experimentally and theoretically to construct the current

description [37, 38, 56]. This current description allows for the corners of the cell to

be filled without stretching of the hydrocarbon chains, but comes at the expense of

tilt deformations along the faces of the hexagon.









From the MD simulations described in Chapter 2 the hexagonal phase was

examined. For a system of 10% LPA and 90% DOPE simulated at T = 300 K, as

described in Section 2.3.3, beginning from a multi-lamellar configuration formed

a stable H11 phase. This H11 was examined to evaluate if the current hexagonal

construction was observed. Figure 3-18a shows a simulation snapshot from the

above described system after 1.2ps of simulation, it water is removed to show

that the water pores have a hexagonal geometry. Figure 3-18b illustrates the

hexagonal nature of the hydrophobic by plotting the location of all tail beads in

the lipids. Figure 3-18c plots the location of the lipid-water interface, defined as

the mid-point of the bond connecting the phosphate bead with the glycerol bead

and fitted to a hexagon by least-squares. It is evident from the Figure 3-18 the

H11 phase observed in our MD simulations do reflect the description of Hamm and

Kozlov [37].

In order for the lipids to accommodate the structure proposed in Figure 3-17c

and verified by Figures 3-18a-c, the lipids must tilt away from the normal direction

to the lipid-water interface. At the center of each face of the hexagonal phase the

lipid will be aligned with the normal, but away from the center, the lipids must

tilt. Hamm and Kozlov predict the tilt should be linear, that it will reach it's

maximum values that the outer edge of each face, and along the pore direction tilt

should be zero [37]. This has been verified through analysis of the MD simulation

of the above described 10% LPA/DOPE system. To calculate tilt for each lipid

the director, n, must be suitably defined. In this work, n, is defined by the vector

pointing from the lipid-water interface location to the center of mass of the tail

beads. The tilt along a face of the hexagonal lipid-water interface, t,, and along

the pore direction ty averaged over 50 data points in time beginning after 800ns of

simulation, with a time step of 1.6ns. Figure a shows tx as a function of position

along the facet of the hexagonal pore. It can be seen t, is a linear function with








61

























32 6.O



30 O0L








(a) (b) (c)

Figure 3-18: Examination of HI, from MD simulations for a 10% LPA/DOPE sys-
tem at T = 300 Kafter 1.2/ps. a) Simulation snapshot of Hjjphase. b) Tail bead
locations plotted to illustrate the hydrophobic boundary. c) Lipid-water interface
location (defined as the middle of the phosphate-glycerol bond) plotted for each
lipid to illustrate the hydrophilic boundary.













08 08

06 06

04 04

02 02

0 0

-02 -02

-04 -04

-06 -06

-0 -1 -05 0 05 1 5 0 5 10 15
X position (nm) Y position (nm)

(a) (b)


Figure 3-19: Tilt in a HI, phase :a) t1, tilt along the one face of the hexagon. b)

ty, tilt along the direction of the water pore.



extremal values the edge of the face. Tilt along the pore direction is shown in


Figure b, it is nearly zero and has no sustained gradient. Both of these results in


Figure a and b correspond to the theoretical predictions of Hamm and Kozlov[37].















CHAPTER 4
CONCLUSIONS

4.1 1M. i 'r Conclusions

One of the advantages of coarse grained models in molecular dynamics is

that the systems can be visualized and measured on near atomic length scales,

while relatively large systems and long simulations can be simulated without being

extremely computationally expensive. This ability allows for insight to be gained

on the molecular scale as well as the ability to compute bulk properties for these

-iv-1' Iin- Here we were able to verify the changing molecular geometries in lipid

--ii-1ii In as a function of the environment parameters.

1. The CG model developed for the lipids LPA and DOPA and the ion

Mg+2exhibit qualitative agreement with experimental temperature and

charge based phase transition data.

2. For pure LPA, a negatively charged lipid, the head size reduced with in-

creasing cationic concentration, making the molecule more cylindrical and

reducing the curvature of the;-, -1I I. For mixed DOPA-DOPE the increase

in temperature of the system caused the tails to splay further apart, while the

head group separation also increased with the increased vibrational energy

but to lesser degree. Thereby it was observed the increase in temperature for

this system causes the molecules to become more conical leading to a system

of higher curvature.This study provided a basis for further investigation in

biological systems with a coarse grained model for the lipids LPA and DOPA.

The ability the qualitatively reproduce experimental phases in both pure and

mixed lipid systems is justification for extending the work in this area.









3. The intermediate stalk phase observed in phase transition simulations

of mixed lipid systems verifies the structure proposed by Kozlovsky and

Kozlov [26].

4. Increasing the relative concentration of PI4P in a PI4P/DPPC bilayer causes

an increase in the bending modulus of the membrane compared to a pure

DPPC bilayer. This indicates that the bilayer becomes more difficult to

deform with the addition of more PI4P up to 21' .

5. A system in which PI4P localizes in one region of a DPPC bilayer while

neighboring regions remain homogeneous, exhibits a positive line tension on

the order of 10-10N when the concentration of PI4P in the inhomogeneous

section reaches 2' .- If the inhomogeneous section only contains 5% PI4P no

anisotropy in the 2-dimension pressure tensor is observed.

6. The hexagonal phase in these MD simulations display hexagonal phases at

the both the hydrophilic and hydrophobic boundaries. The tilt along a side

of the hexagon hydrophilic interface is a linear increasing function. Both of

these results support the theoretical predictions of Hamm and Kozlov [37].

4.2 Future Directions

In this work molecular modeling and analysis techniques were developed which

produced results warranting further investigation. In this work only a few lipid

-i-i -I I were studied, while it is believed that other lipid species could be key

molecules in eliciting changes in elastic parameters of membranes.

4.2.1 Additional Lipid Species

The lipid phosphatidylinositol-4,5-bisphosphate (PI(4,5)P2) has been impli-

cated to play a role in stimulating vesicle formation from liposomes as well as from

the outer cellular membrane [29, 30, 35]. PI(4,5)P2is structurally similar to PI4P,

but has an additional phosphate group extending from the inositol ring and carries

a net -4 charge compared to -3 of PI4P at pi- -i. .lgical conditions [57]. It is











Phosphate


Na -- Ia I
o c Palmitoyl GI, c roI
S Fatty Acid C C Oleoyl
C Fatty Acid
o c ( ( c

Cc

Figure 4-1: Proposed CG models of additional lipids: a) PI(4,5)2. b) DAG.


believe that the increased head group size and increased effective head group size

due to charge would potentially increase the anisotropy in the 2 dimension pressure

tensor in the line tension study. Additionally the lipid diacyl glycerol (DAG) is of

interest, due to it role in Golgi budding and bud fission from the Golgi network.

DAG is similar to the lipid DOPA except that it does not have the phosphate

group extending from the glycerol group. It is desired to conduct similar studies

done in this work, but using PI(4,5)P2 and DAG as the minor component in the

lipid systems. Proposed CG models of the lipids PI(4,5)P2and DAG as shown in

Figure 4-1.

4.2.2 Continuum Scale Modeling

Ultimately, it is desired to use the parameters determined from molecular mod-

eling and incorporated these findings into a continuum scale model for predicting

equilibrium membrane configurations, based upon the membrane chemical com-

position. Work has been done in this area starting with a theoretical description

of the bending energy by Helfrich and extending by Lipowsky for homogeneous

membranes [39, 41]. This work was further extended to incorporate the existence

of different homogeneous lipid phases with the same lipid by modifying the Helfrich

expression to incorporate line tension between phases [1, 34, 35]. The expression of

the bending and interface energy is for a system with distinct lipid phases is given



















Figure 4-2: Domains within a membrane, taken from ref [1]


by Eq. 4.1

F [~-(C1 + C2 + C()2 + )CC2]CA+

2 (Cl + C2 + C))2 + rC1C2]dA+

o dl (4.1)

where K(, (a), CG() are the bending modulus, Gaussian curvature modulus

and spontaneous curvature of the a phase, da is the length of the interface

separating the phases. The geometry of the system is described by Figure 4-2.

A model describing the energy of associated with tilt and bending of lipid

structures has been introduced by Hamm and Kozlov shown in Eq. 4.2 and

previously described in Chapter 3 of this work [36].

1 1 t( t2
f (t| + J,)2 + idet(t)+ 0t2 (4.2)

It is desired to incorporate both the Julicher and Hamm models to develop a

model which accounts for bending, tilt and line tension. Further work must be done

to develop a method to extract the tilt modulus from the MD simulations [35, 36].

Using these energy models it is desired to compute the equilibrium shapes for

varying concentrations of mixed lipid systems using parameters extracted from

molecular simulations.















REFERENCES


[1] R. Lipowsky, "Budding of membranes induced by intramembrane domains," J
Phys II, vol. 2, no. 10, pp. 1825-40, Oct 1992.

[2] B. Alberts, D. Bray, J. Lewis, M. Raff, K. Roberts, and J. D. Watson, Molecu-
lar Biology of The Cell, Garland Publishing, 3rd edition, 1994.

[3] R. Schekman and L. Orci, "Coat proteins and vesicle budding.," Science, vol.
271, no. 5255, pp. 1526-1533, Mar 1996.

[4] J. E. Rothman and F. T. Wieland, "Protein sorting by transport vesicles.,"
Science, vol. 272, no. 5259, pp. 227-234, Apr 1996.

[5] L. V. Chernomordik and M. M. Kozlov, "Protein-lipid interplay in fusion and
fission of biological membranes.," Ann Rev Biochem, vol. 72, pp. 175-207,
2003.

[6] K. N. Burger, "Greasing membrane fusion and fission machineries.," Traffic,
vol. 1, no. 8, pp. 605-613, 2000.

[7] R. Weigert, M. G. Silletta, S. Span, G. Turacchio, C. Cericola, A. Colanzi,
S. Senatore, R. Mancini, E. V. Polishchuk, M. Salmona, F. Facchiano, K. N.
Burger, A. Mironov, A. Luini, and D. Corda, "CtBP/BARS induces fission of
Golgi membranes by acylating lysophosphatidic acid.," Nature, vol. 402, no.
6760, pp. 429-33, Nov 1999.

[8] A. Siddhanta and D. Shields, "Secretory vesicle budding from the trans-Golgi
network is mediated by phosphatidic acid levels.," J Biol Chem, vol. 273, no.
29, pp. 17995-8, Jul 1998.

[9] M. D. Matteis, A. Godi, and D. Corda, "Phosphoinositides and the golgi
complex.," Curr Opin Cell Biol, vol. 14, no. 4, pp. 434-47, Aug 2002.

[10] M. Allen and D. Tildesley, Computer Simulations of Liquids, Oxford Science,
1987.

[11] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and
J. R. Haak, i\! I 'cular dynamics with coupling to an external bath," J Chem
Phys, vol. 81, no. 8, pp. 3684-90, 1984.

[12] D. van der Spoel et al., Gromacs User Manual version 3.2, www.gromacs.org,
2004.









[13] D. P. Tieleman, S. J. Marrink, and H. J. Berendsen, "A computer perspective
of membranes: molecular dynamics studies of lipid bilayer systems.," Biochim
B.:, !,.;l; Acta, vol. 1331, no. 3, pp. 235-70, Nov 1997.

[14] J. N. Israelachvili, Intermolecular and Surface Forces : With Applications to
Colloidal and Biological SplAt- ii -. Academic Press, second edition, 1991.

[15] E. E. Kooijman, V. Chupin, B. de Kruijff, and K. N. J. Burger, "Modulation
of membrane curvature by phosphatidic acid and lysophosphatidic acid.,"
Traffic, vol. 4, no. 3, pp. 162-74, Mar 2003.

[16] E. E. Kooijman, V. Chupin, N. L. Fuller, M. M. Kozlov, B. de Kruijff, K. N. J.
Burger, and P. R. Rand, "Spontaneous curvature of phosphatidic acid and
lysophosphatidic acid.," Biochemistry, vol. 44, no. 6, pp. 2097-102, Feb 2005.

[17] R. Goetz and R. Lipowsky, "Computer simulations of bilayer membranes:
Self-assembly and interfacial tension," J Chem Phys, vol. 108, no. 17, pp.
7397-7409, 1998.

[18] R. Faller and S.-J. Marrink, "Simulation of domain formation in DLPC-DSPC
mixed bilayers.," Langmuir, vol. 20, no. 18, pp. 71,,i-7693, 2004.

[19] B. Smit, '- \I.. ,cular-dynamics simulations of amphiphilic molecules at a
liquid-liquid interface," Phys Rev A, vol. 37, pp. 3431-3433, 1988.

[20] B. J. Palmer and J. Liu, "Simulations of micelle self-assembly in surfactant
solutions," Langmuir, vol. 12, pp. 746-753, 1996.

[21] R. P. Rand and N. L. Fuller, "Structural dimensions and their changes in a
reentrant hexagonal-lamellar transition of phospholipids.," B.: ,'1!;1.i J, vol. 66,
no. 6, pp. 2127-38, Jun 1994.

[22] L. Yang, L. Ding, and H. W. Huang, "New phases of phospholipids and
implications to the membrane fusion problem," Biochemistry, vol. 42, no. 22,
pp. 6631-35, Jun 2003.

[23] D. P. Siegel, "The modified stalk mechanism of lamellar/inverted phase
transitions and its implications for membrane fusion.," B.:,!';l ,- J, vol. 76, no.
1 Pt 1, pp. 291-313, Jan 1999.

[24] D. P. Siegel and R. M. Epand, "The mechanism of lamellar-to-inverted
hexagonal phase transitions in phosphatidylethanolamine: implications for
membrane fusion mechanisms.," B.:,.l, /.- J, vol. 73, no. 6, pp. 3089-3111, Dec
1997.

[25] M. Rappolt, A. Hickel, F. Bringezu, and K. Lohner, i\!. Ii..,i-iii of the
lamellar/inverse hexagonal phase transition examined by high resolution x-ray
diffraction.," B.:,.'1.7;- J, vol. 84, no. 5, pp. 3111-3122, M.I.- 2003.









[26] Y. Kozlovsky and M. M. Kozlov, "Stalk model of membrane fusion: solution
of energy crisis.," B.:',..1!.;1 J, vol. 82, no. 2, pp. 882-895, Feb 2002.

[27] S. W. Hui, T. P. Stewart, and L. T. Boni, "The nature of lipidic particles and
their roles in polymorphic transitions.," Chem Phys Lipids, vol. 33, no. 2, pp.
113-126, Aug 1983.

[28] E. St..u-l-.-.--,- E. J. Prenner, M. Kriechbaum, G. Degovics, R. N. Lewis,
R. N. McElhaney, and K. Lohner, "X-ray studies on the interaction of the
antimicrobial peptide gramicidin S with microbial lipid extracts: evidence
for cubic phase formation.," Biochim B.:, ,,,..l Acta, vol. 1468, no. 1-2, pp.
213-230, Sep 2000.

[29] T. F. Martin, "PI(4,5)P(2) regulation of surface membrane traffic.," Curr
Opin Cell Biol, vol. 13, no. 4, pp. 493-9, Aug 2001.

[30] M. Jost, F. Simpson, J. M. Kavran, M. A. Lemmon, and S. L. Schmid,
"Phosphatidylinositol-4,5-bisphosphate is required for endocytic coated vesicle
formation.," Curr Biol, vol. 8, no. 25, pp. 1399-402, 1998.

[31] T. Takenawa and T. Itoh, "Phosphoinositides, key molecules for regulation
of actin cytoskeletal organization and membrane traffic from the plasma
membrane.," Biochim B.:,.! ;1..- Acta, vol. 1533, no. 3, pp. 190-206, Oct 2001.

[32] E. A. Evans, "Bending resistance and chemically induced moments in
membrane bilayers.," B.:i,'1;I J, vol. 14, no. 12, pp. 923-31, Dec 1974.

[33] E. Lindahl and O. Edholm, -'!. --..- ,pic undulations and thickness fluctua-
tions in lipid bilayers from molecular dynamics simulations.," B.: ,.1!;1i. J, vol.
79, no. 1, pp. 426-33, 2000.

[34] F. Jlicher and R. Lipowsky, "Domain-induced budding of vesicles.," Phys Rev
Lett, vol. 70, no. 19, pp. 2964-2967, M.v 1993.

[35] F. Jlicher and R. Lipowsky, "Shape transformations of vesicles with intramem-
brane domains.," Phys Rev E, vol. 53, no. 3, pp. 2670-2683, Mar 1996.

[36] M. Hamm and M. Kozlov, "Elastic energy of tilt and bending of fluid
membranes," Eur Phys J E, vol. 3, pp. 323-335, 2000.

[37] M. Hamm and M. Kozlov, "Tilt model of inverted amphiphilic mesophases,"
Eur Phys J B, vol. 6, pp. 519-528, 1998.

[38] D. C. Turner and S. M. Gruner, "X-ray diffraction reconstruction of the
inverted hexagonal (HII) phase in lipid-water systems.," Biochemistry, vol. 31,
no. 5, pp. 1340-1355, Feb 1992.

[39] W. Helfrich, "Elastic properties of lipid bilayers: theory and possible experi-
ments.," Z Naturforsch [C], vol. 28, no. 11, pp. 693-703, 1973.









[40] K. Berndl, J. Kas, R. Lipowsky, E. Sackmann, and U. Seifert, "Shape
transformations of giant vesicles: Extreme sensitivity to bilayer ..i- .mmetry,"
E,.../,,~.; i Lett, vol. 13, no. 7, pp. 659-664, 1990.

[41] U. Seifert, K. Berndl, and R. Lipowsky, "Shape transformations of vesicles:
Phase diagram for spontaneous- curvature and bilayer-coupling models.," Phys
Rev A, vol. 44, no. 2, pp. 1182-1202, Jul 1991.

[42] E. A. Evans and W. Rawiczl, "Entropy-driven tension and bending elasticity
in condensed-fluid membranes.," Phys Rev Lett, vol. 64, no. 17, pp. 2094-2097,
Apr 1990.

[43] S. J. Marrink and A. E. Mark, "Effect of undulations on surface tension in
simulated bilayers," J Phys Chem B, vol. 105, pp. 6122-27, 2001.

[44] J. Bradshaw, R. Bushby, C. Giles, M. Saunders, and A. Saxena, "The
headgroup orientation of dimyristoylphosphatidylinositol-4-phosphate in mixed
lipid bilayers: a neutron diffraction study.," Biochim B:,.!-1,.;1- Acta, vol. 1329,
no. 1, pp. 124-38, Oct 1997.

[45] M. V. Berry, "Liquid surfaces," Surf Sci, vol. 1, pp. 291-327, 1975.

[46] J. Hansen and I. McDonald, Ti, i.,y of Simple Liquids, Academic Press, 1986.

[47] L. E. Malvern, Introduction to the Mechanics of a Continuous Medium,
Prentice-Hall, Inc., 1969.

[48] J. Sonne, F. Y. Hansen, and G. H. Peters, '\!I 1I,,dological problems in
pressure profile calculations for lipid bilayers.," J Chem Phys, vol. 122, no. 12,
pp. 124903, Mar 2005.

[49] P. Schofield and J. Henderson, "Statistical mechanics of inhomogeneous
fluids," Proc R Soc Lond A, vol. 379, pp. 231-246, 1982.

[50] J. H. Irving and J. G. Kirkwood, "The statistical mechanical theory of
transport processes. iv. the equations of hydrodynamics," J Chem Phys, vol.
18, no. 6, pp. 817-29, Jun 1950.

[51] A. Harasima, "Molecular theory of surface tension," Adv Chem Phys, vol. 1,
pp. 203-237, 1958.

[52] M. Berry, "The molecular mechanism of surface tension," Phys Educ, vol. 6,
pp. 79-84, 1971.

[53] J. Gullingsrud and K. Schulten, "Lipid bilayer pressure profiles and
mechanosensitive channel gating.," B.:-! ,/,.1 J, vol. 86, no. 6, pp. 3496-3509,
Jun 2004.







71

[54] E. Lindahl and O. Edholm, "Spatial and energetic-entropic decomposition of
surface tension in lipid bilayers from molecular dynamics simulations," J Chem
Phys, vol. 113, no. 9, pp. 3882-93, Sep 2000.

[55] T. Baumgart, S. T. Hess, and W. W. Webb, Iiii.,-iig coexisting fluid domains
in biomembrane models coupling curvature and line tension.," Nature, vol.
425, no. 6960, pp. 821-4, Oct 2003.

[56] P. Dn iin. R. Templer, and J. Seddon, "Quantify packing frustration energy
in inverse lyotropic mesophases," Langmuir, vol. 13, pp. 351-359, 1997.

[57] S. McLaughlin, J. Wang, A. Gambhir, and D. Murray, "PIP(2) and proteins:
interactions, organization, and information flow.," Annu Rev B.:1/,.l;1- Biomol
Struct, vol. 31, pp. 151-175, 2002.















BIOGRAPHICAL SKETCH

Eric _\.iv was born in New Haven, Connecticut, on August 29, 1978. He was

raised in Cheshire, Connecticut and attended Cheshire High School, graduating

in 1996. He received a Bachelor of Science from Bucknell University in chemical

engineering in 2001. During his tenure at Bucknell he completed summer intern-

ships with Neurogen Corporation (Branford, CT) and the Connecticut Department

of Environmental Protection (Hartford, CT). In 2001 he joined the Department

of Chemical Engineering at the University of Florida and Atul Narang's research

group. His research interests are focused on mathematical modeling of biological

systems.