UFDC Home  myUFDC Home  Help 



Full Text  
MOLECULAR MODELING OF BIOMEMBRANE DEFORMATIONSTHE 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 11 Membrane diagrams ........................... 2 12 Coat protein assembly and vesicle formation ..... . . 4 13 Potential functions ................ . .... 6 14 Comparision of atomistic and CG structures ..... . . 9 21 Detailed atomic structures and the corresponding CG models . 16 22 Dependence of the size and shape of LPA micelles on \!g2+] ..... 18 23 Dependence of the micelle asphericity on the Mg2+ concentration. .19 24 Bilayer configuration of pure DOPA systems . . ..... 20 25 Temperature induced phase transitions ................ ..21 26 Simulations of temperatureinduced phase transitions . ... 23 27 Molecular geometry characterization for pure LPA systems . 25 28 Molecular geometry characterization in mixed lipid systems . 27 29 Proposed transition states in vesicle fusion .............. ..28 210 Stalk formation in DOPE/LPA system ............. .. 29 31 Comparision of theoretical and experimental vesicle shapes . 33 32 Atomisitic and CG representation of phosphoinositides . ... 35 33 Depictions of PI4P ............... ........ .. 36 34 Angle probability distributions ................ ...... 38 35 Spectral intensity for PI4P/DPPC systems . . 39 36 Bending modulus for PI4P/DPPC bilayer . . ...... 39 37 Definition of rij .................. ........... 43 38 Angle Orientation .................. ......... .. 45 39 Slab orientation for calculation of bilayer surface tension ...... ..47 310 Phases within a bilayer .................. ..... .. 50 311 Patch system with repulsive walls ................ 52 312 Density of pressure of CG systems ................ 54 313 Virial and pressure profile comparison ................ ..55 314 Bond potential comparison .................. ..... 55 315 Density and pressure profiles for patch systems . . ..... 57 316 Tilt and bending deformations .................. ..... 58 317 Construction of H1I phase ............. ..... 59 318 Examination of H1I phase from MD .... . ... 61 319 Tilt in a HII phase ........... . . ... .62 41 Proposed CG models of additional lipids: a) PI(4,5)2. b) DAG. . 65 42 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 DEFORMATIONSTHE 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 phosphatidylinositol4phosphate (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 1la [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 1lb. 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 11: 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 1la. 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 12: 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 12. 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 [79]. 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 nonequilibrium 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 nonbonded interactions are due to van der Waals and electrostatic forces. For computational purposes, only pair potentials are typically considered, taking 3body and higher term interactions becomes very time consuming because they involve summing over triplets (or greater) of molecules [10]. However, the 2body 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 13: Potential functions: a) LennardJones potential b) Electrostatic poten tial The non bonded forces are often specified by a LennardJones 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 cutoff radius is specified to increase the computational efficiency. If two particles are separated by a distance greater than the cutoff radius, it is assumed they have no influence on each other. Special care must be taken when using cutoff 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 cutoff radius. In addition to the LennardJones potential, an electrostatic potential must be accounted for if there are charged particles in the iv1. 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 cutoff 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 cmipl.,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 LennardJones 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 v1, II, to a heat bath using firstorder kinetics. In this weak coupling scheme the v1' 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 Pii 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 firstorder 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 14: 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 46 atoms; shown in Figure 14 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 DOPEDOPA AND LPADOPA 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 nonbilayer intermediates. Consequently, these processes do not occur spontaneouslythey 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 twotailed 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 condensedmatter pir, i. that amphiphilic molecules selfassemble into . . regates whose morphology is intimately linked to the shape of the individual molecules [14]. Specifically, V and Ashaped 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 LPAacyltransferases is ultimately driven by the distinct shapes of LPA and PA molecules. Recently, Kooijman et al. performed experiments in a cellfree 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 31PNMR 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 coarsegrained 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 lengthscales. 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 coarsegrained molecular models is a subject of active ongoing research [? ? ? ]. In this work, we will use the coarsegrained 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 coarsegrained 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 nonbonded beads are modeled by the LennardJones 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 coarsegrained model for water, metal ions, and DOPE lipids proposed in by Marrink et al. [? ]. In addition, we have developed coarsegrained 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 coursegrained representations of some of these lipid molecules are shown in Fig. 21. The structure and CG model for the lipid LPA shown in Fig. 2la consists of 7 CG particles. Each bead is classified with respect to the ','I1 iI published by Marrink et al. [? ]. The glycerolester linkage is modelled by a single bead in LPA with the classification Nda, indicating a nonpolar 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 wormlike micelle, where the experiments predict spherical micelles. The structure and model for dioleoyl phosphatidic acid (DOPA) is shown in Fig. 2lb. 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 NGlyercol o c Oleoyl C Oleoyl Fatty Acid C Fatty Acid SCC (a) (b) Figure 21: Detailed atomic structures and the corresponding coarsegrained mod els of considered lipids: (a) LPA and (b) DOPA. Mapping between different groups of atoms and the coarsegrained 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 coarsegrained (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 coarsegrained 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 selfassembly simulations of the LPA lipids (shown in Fig. 2la), 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 selfassembled into micelles and equilibrated within 800 ns. Figures 22ad 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. 22a). When the Mg2+/LPA ratio is increased to 0.2, five micelles are formed (Fig. 22b). If the Mg2+/LPA ratio is increased even further to 1.0, the larger micelles .i' regate to form a single .. . regate. Figures 22c,d show the top and (a) (b) I _____ I i_ _'_ I (c) (d) Figure 22: 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 disklike configuration. Therefore, as the concentration of the ions increases, the micellar size increases and micelles acquire disklike shapes, which is consistent with experimental 1PNMR measurements of Kooijman et al. [15]. To quantify the shape change of the LPA .,.' regates shown in Fig. 22, we calculated their asphericity, which is defined as y3 2f t 2)2 A 2j (Ri R (2.1) j1= 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 23: Dependence of the micelle asphericity on the Mg2+ concentration. 80 ns of the simulation. Fig. 23 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 ...:iregates shown in Fig. 22. According to the shapestructure 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 disklike 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 24: 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 DOPELPA and DOPEDOPA Mixed System We also investigated phase transitions between the lamellar and inverted hexagonal phases in the mixed DOPELPA and DOPEDOPA systems. The exper imental data [15] summarized in Fig. 25 indicate phase transitions from lamellar to inverted hexagonal phases within the bilayer systems as the temperature of the v1i 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. 25). Upon addition of 10 mol% LPA to DOPE, the phase transition temperature increases to Ts 300 K (O curve of Fig. 25). 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. 25). 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 912 waters per lipid, 100 A Pure DOPE (Exp) 90 9  DOPE/LPA (Exp) \ DOPE/DOPA (Exp) 80  SA 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 25: Temperatureinduced 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 1I. 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. 25). 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 presimulated 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 DOPELPA 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 selfassembly of initially randomly dispersed lipids in water to prepare a lipid bilayer. The selfassembly 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. 26a. 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 DOPELPA system are shown in Fig. 26b and c. At temperatures at or below 265 K the lamellar phase remains stable on the timescale of simulations (4.0pts), see Fig. 26b. 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. 26c. All of the configurations shown in Fig. 26 were generated by copying the simulation cell 3 times in the lateral direction to better illustrate the geometry of the systems. (a) (b) Figure 26: Simulations of temperatureinduced phase transition between lamellar and inverted hexagonal phases in the mixed DOPELPA lipid systems: (a) Initial perturbed configuration for DOPELPA simulations. (b) Final configuration of DOPELPA 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 25 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 DOPEDOPA simulations, a bilayer consisting of 120 lipid molecules (10% DOPA, 90% DOPE) was simulated by selfassembly 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 DOPELPA 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 vI' 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 25 shows that the simulations are in qualitative agreement with the experimental data. The model DOPEDOPA system undergoes a phase transition at a temperature approximately 15 K lower than the phasetransition temperature for DOPELPA 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 DOPEDOPA and DOPELPA systems, but the corresponding simulated values transitions temperatures are above 250 K and 265 K, respectively. As can be seen, the current coarsegrained 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 temperaturedependent 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 temperaturedependent phase transitions, the simulations of the DOPELPA 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. 27a, 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 27: 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. iiif i. the individual measurements of ao, lI, and v are shown in Fig. 27b. 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 DOPELPA 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. 28a. 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. 28b 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[2327]. 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. 29a. 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. 29b; in order to reach this configuration the lipid tails must tilt away from normal to waterlipid 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 28: Molecular Geometry Characterization of DOPE, LPA and DOPA in mixed system: a) DOPELPA i1 I 1i1 packing parameters. b) DOPELPA system measurements of tail volume, tail length and head spacing. c) DOPEDOPA system packing parameters d) DOPEDOPA measurements of tail volume, tail length and head spacing. Temperature (K) 235 240 ~f p  (a) (b) Figure 29: 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. 210ad. 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 210: 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 LennardJones 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 nonuniform 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, 2931]. In this study the effect of introducing a phosphoinositide molecule, phosphatidylinositol4phosphate (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 lipidwater 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 [3638]. 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 31: 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 EulerLagrange 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 31 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 32: 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 32, 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 33a shows how the inositol ring gets pulled down towards the S0) a = 180 ' (a) (b) Figure 33: 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 33b 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 Selfassembly 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 34a 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 bilayerwater 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 pointwise 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 34: 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(lq), these spectral intensity plots for each system are shown in Figure 35ad. The longer wavelength modes, q > 1.0 nm, were used to fit a line by the leastsquares 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 36. 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 35: Spectral intensity for PI4P/DPPC PI4P. c) 10% PI4P. d) 20% PI4P. x 1020 95 085 9 851 T 5 10 15 % PI4P Figure 36: 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 [79, 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 twocomponents: 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 coarsegrain model simulations are described in Section 3.3.2. 3.3.2 Description of Potentials and Forces LennardJones Potential The LennardJones(LJ) potential is a nonbonded 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 cutoff 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 37: Definition of rij where ri is the vector pointing from bead i to bead j as described in Figure 37, 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 2body potentials for electrostatic and bonded interactions described below. Electrostatic Potential The electrostatic potential (ES) is a nonbonded 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 38. The angle potential, Vo, and forces, Fo are described below, F(j) i k F(i) j F(k) Figure 38: 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 2body 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 Vir2Body 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 39. The slab is moved through the zdirection 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 39: 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 IrvingKirkwood 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 IrvingKirkwood contour is used [51]. A generalized expression for the static pressure for an mbody potential can be derived from Eq. 3.16 is described by 3.17. P (Body {Vj~ V :} C6(R 1,// (3.17) 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) Smbody VolAz 0 0 Zi Smbody 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 (IrvingKirkwood) 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 rewritten 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 zdirection 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 310: 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 310 describes the geometry of phases embedded within a bilayer. In this situation pressure would be calculated as a function of the zdirection, 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 2dimensional 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 ydirection 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 310 describes the desired configuration for the line tension calculation in which the aphase 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 311: 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 102. This potential is shifted and cutoff with the same radii and in the same manner as all other LJ potentials as described in Section 3.3.2. Figure 311 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 312a the density profile is plotted, where the units are number of beads per slab. In Figure 312 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 313: 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 314: 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 310. 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 315. 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 1011N and for the 20% system a = 2.0 x 1010N. 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 1013N) and theoretical estimate (1011)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 315: 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 316: Tilt and bending: a)Definition of tilt. b) Nonuniform 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 316a. t N (3.26) nN 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 316b and (a) (b) (c) Figure 317: 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 nonuniform tilt deformation, where tilt changes across the element. Figure 316c 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 317. 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 317b; both situation are energetically unfavorable. To alleviate this l II .11 i. i i. i,y ", as it is know, the construction of the unit cell was reexamined 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 multilamellar configuration formed a stable H11 phase. This H11 was examined to evaluate if the current hexagonal construction was observed. Figure 318a 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 318b illustrates the hexagonal nature of the hydrophobic by plotting the location of all tail beads in the lipids. Figure 318c plots the location of the lipidwater interface, defined as the midpoint of the bond connecting the phosphate bead with the glycerol bead and fitted to a hexagon by leastsquares. It is evident from the Figure 318 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 317c and verified by Figures 318ac, the lipids must tilt away from the normal direction to the lipidwater 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 lipidwater interface location to the center of mass of the tail beads. The tilt along a face of the hexagonal lipidwater 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 318: 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) Lipidwater interface location (defined as the middle of the phosphateglycerol 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 319: 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 iv1' Iin Here we were able to verify the changing molecular geometries in lipid ii1ii 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 DOPADOPE 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 1010N when the concentration of PI4P in the inhomogeneous section reaches 2' . If the inhomogeneous section only contains 5% PI4P no anisotropy in the 2dimension 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 ii 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 phosphatidylinositol4,5bisphosphate (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 41: 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 41. 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 42: 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 42. 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. 182540, 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. 15261533, Mar 1996. [4] J. E. Rothman and F. T. Wieland, "Protein sorting by transport vesicles.," Science, vol. 272, no. 5259, pp. 227234, Apr 1996. [5] L. V. Chernomordik and M. M. Kozlov, "Proteinlipid interplay in fusion and fission of biological membranes.," Ann Rev Biochem, vol. 72, pp. 175207, 2003. [6] K. N. Burger, "Greasing membrane fusion and fission machineries.," Traffic, vol. 1, no. 8, pp. 605613, 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. 42933, Nov 1999. [8] A. Siddhanta and D. Shields, "Secretory vesicle budding from the transGolgi network is mediated by phosphatidic acid levels.," J Biol Chem, vol. 273, no. 29, pp. 179958, 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. 43447, 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. 368490, 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. 23570, 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. 16274, 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. 2097102, Feb 2005. [17] R. Goetz and R. Lipowsky, "Computer simulations of bilayer membranes: Selfassembly and interfacial tension," J Chem Phys, vol. 108, no. 17, pp. 73977409, 1998. [18] R. Faller and S.J. Marrink, "Simulation of domain formation in DLPCDSPC mixed bilayers.," Langmuir, vol. 20, no. 18, pp. 71,,i7693, 2004. [19] B. Smit, ' \I.. ,culardynamics simulations of amphiphilic molecules at a liquidliquid interface," Phys Rev A, vol. 37, pp. 34313433, 1988. [20] B. J. Palmer and J. Liu, "Simulations of micelle selfassembly in surfactant solutions," Langmuir, vol. 12, pp. 746753, 1996. [21] R. P. Rand and N. L. Fuller, "Structural dimensions and their changes in a reentrant hexagonallamellar transition of phospholipids.," B.: ,'1!;1.i J, vol. 66, no. 6, pp. 212738, 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. 663135, 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. 291313, Jan 1999. [24] D. P. Siegel and R. M. Epand, "The mechanism of lamellartoinverted hexagonal phase transitions in phosphatidylethanolamine: implications for membrane fusion mechanisms.," B.:,.l, /. J, vol. 73, no. 6, pp. 30893111, Dec 1997. [25] M. Rappolt, A. Hickel, F. Bringezu, and K. Lohner, i\!. Ii..,iiii of the lamellar/inverse hexagonal phase transition examined by high resolution xray diffraction.," B.:,.'1.7; J, vol. 84, no. 5, pp. 31113122, 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. 882895, 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. 113126, Aug 1983. [28] E. St..ul.., E. J. Prenner, M. Kriechbaum, G. Degovics, R. N. Lewis, R. N. McElhaney, and K. Lohner, "Xray 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. 12, pp. 213230, 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. 4939, Aug 2001. [30] M. Jost, F. Simpson, J. M. Kavran, M. A. Lemmon, and S. L. Schmid, "Phosphatidylinositol4,5bisphosphate is required for endocytic coated vesicle formation.," Curr Biol, vol. 8, no. 25, pp. 1399402, 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. 190206, 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. 92331, 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. 42633, 2000. [34] F. Jlicher and R. Lipowsky, "Domaininduced budding of vesicles.," Phys Rev Lett, vol. 70, no. 19, pp. 29642967, 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. 26702683, Mar 1996. [36] M. Hamm and M. Kozlov, "Elastic energy of tilt and bending of fluid membranes," Eur Phys J E, vol. 3, pp. 323335, 2000. [37] M. Hamm and M. Kozlov, "Tilt model of inverted amphiphilic mesophases," Eur Phys J B, vol. 6, pp. 519528, 1998. [38] D. C. Turner and S. M. Gruner, "Xray diffraction reconstruction of the inverted hexagonal (HII) phase in lipidwater systems.," Biochemistry, vol. 31, no. 5, pp. 13401355, Feb 1992. [39] W. Helfrich, "Elastic properties of lipid bilayers: theory and possible experi ments.," Z Naturforsch [C], vol. 28, no. 11, pp. 693703, 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. 659664, 1990. [41] U. Seifert, K. Berndl, and R. Lipowsky, "Shape transformations of vesicles: Phase diagram for spontaneous curvature and bilayercoupling models.," Phys Rev A, vol. 44, no. 2, pp. 11821202, Jul 1991. [42] E. A. Evans and W. Rawiczl, "Entropydriven tension and bending elasticity in condensedfluid membranes.," Phys Rev Lett, vol. 64, no. 17, pp. 20942097, 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. 612227, 2001. [44] J. Bradshaw, R. Bushby, C. Giles, M. Saunders, and A. Saxena, "The headgroup orientation of dimyristoylphosphatidylinositol4phosphate in mixed lipid bilayers: a neutron diffraction study.," Biochim B:,.!1,.;1 Acta, vol. 1329, no. 1, pp. 12438, Oct 1997. [45] M. V. Berry, "Liquid surfaces," Surf Sci, vol. 1, pp. 291327, 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, PrenticeHall, 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. 231246, 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. 81729, Jun 1950. [51] A. Harasima, "Molecular theory of surface tension," Adv Chem Phys, vol. 1, pp. 203237, 1958. [52] M. Berry, "The molecular mechanism of surface tension," Phys Educ, vol. 6, pp. 7984, 1971. [53] J. Gullingsrud and K. Schulten, "Lipid bilayer pressure profiles and mechanosensitive channel gating.," B.:! ,/,.1 J, vol. 86, no. 6, pp. 34963509, Jun 2004. 71 [54] E. Lindahl and O. Edholm, "Spatial and energeticentropic decomposition of surface tension in lipid bilayers from molecular dynamics simulations," J Chem Phys, vol. 113, no. 9, pp. 388293, 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. 8214, Oct 2003. [56] P. Dn iin. R. Templer, and J. Seddon, "Quantify packing frustration energy in inverse lyotropic mesophases," Langmuir, vol. 13, pp. 351359, 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. 151175, 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. 