Applied Science, Faculty of
Electrical and Computer Engineering, Department of
DSpace
UBCV
Duan, Yan
Master of Applied Science - MASc
University of British Columbia
This work presents a study of cholesterol's effects on archaebacterial cell membranes using coarse grain molecular dynamic simulations. As a major component in biological membranes, cholesterol is closely related to the dynamics of lipids and biomechanical properties of the membrane. A coarse grained molecular dynamics (CGMD) model is constructed to study the membrane properties. The CMGD model provides insights into the diffusion dynamics of lipids, membrane thickness, line tension, and surface tension as a function of cholesterol content. The CGMD simulation results are validated using experimental measurements from a tethered archaebacterial bilayer lipid membrane. The membrane is tethered to an inert gold bioelectronic interface, which allows the experimental measurements to be performed using standard laboratory equipments. A fractional order macroscopic model is introduced to link microscopic simulation results with macroscopic experimental measurements. To ensure the bioelectronic interface does not affect the membrane dynamics and biomechanics, it is shown that variations in the position dependent water density are negligible near the surface of the membrane. Furthermore, the Percus-Yevick equation is used to confirm that harsh repulsive forces play a negligible role in the long range dynamics of the water density profile.
Study of Cholesterol in TetheredMembrane Using Coarse-GrainedMolecular Dynamics SimulationsbyYan DuanB.Sc in Engineering Physics, The University of Alberta, 2014A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Electrical and Computer Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)September 2015c Yan Duan 2015AbstractThis thesis presents a study of cholesterol's e ects on archaebacterial cellmembranes using coarse grain molecular dynamic simulations. As a majorcomponent in biological membranes, cholesterol is closely related to the dy-namics of lipids and biomechanical properties of the membrane. A coarsegrained molecular dynamics (CGMD) model is constructed to study themembrane properties. The CMGD model provides insights into the dif-fusion dynamics of lipids, membrane thickness, line tension, and surfacetension corresponding to di erent cholesterol levels. The CGMD simulationresults are validated using experimental measurements from a tethered ar-chaebacterial bilayer lipid membrane. The membrane is tethered onto aninert gold bioelctronic interface, which allows the experimental measure-ments to be performed using standard laboratory equipments. A fractionalorder macroscopic model is introduced to link microscopic simulation resultswith macroscopic experimental measurements. To ensure the bioelectronicinterface does not a ect the membrane dynamics and biomechanics, it isshown that variations in the position dependent water density are negligiblenear the surface of the membrane. Furthermore, the Percus-Yevick equationis used to con rm that harsh repulsive forces play a negligible role in thelong range dynamics of the water density pro le.iiPrefaceThis thesis consists of three main chapters. The rst chapter gives the mo-tivation of the thesis and ve levels of abstraction for cell membrane model-ing. The second chapter provides a comprehensive introduction on necessaryconcepts and theories to understand molecular dynamics, and its simulationschemes and tools. The third chapter is largely based on work conduct-ed in statistical signal processing lab at University of British Columbia, byDr.William Hoiles, Prof.Vikram Krishnamurthy and Yan Duan. I was re-sponsible for all the coarse-grained simulation, derivation and applicationof Percus-Yevick equation for the water density pro le, numerical solutionof Volterra di erential-integral equation, and parameter estimations for thefractional order macroscopic model.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation: the Biosensor . . . . . . . . . . . . . . . . . . . . 11.2 Cell Membrane System . . . . . . . . . . . . . . . . . . . . . 21.3 Bioelectronic Interface . . . . . . . . . . . . . . . . . . . . . . 31.4 Abstractions in Bio-System Modeling . . . . . . . . . . . . . 41.4.1 Ab-initio Molecular Dynamics . . . . . . . . . . . . . 51.4.2 Molecular Dynamics . . . . . . . . . . . . . . . . . . . 51.4.3 Coarse-Grained Molecular Dynamics . . . . . . . . . 61.4.4 Continuum Theory . . . . . . . . . . . . . . . . . . . 61.4.5 Macroscopic Model . . . . . . . . . . . . . . . . . . . 61.5 Thesis Contributions . . . . . . . . . . . . . . . . . . . . . . 72 Fundamental Molecular Dynamics . . . . . . . . . . . . . . . 92.1 Basic Algorithm and Force Fields . . . . . . . . . . . . . . . 92.1.1 Global MD Algorithm . . . . . . . . . . . . . . . . . . 92.1.2 Force Fields . . . . . . . . . . . . . . . . . . . . . . . 112.2 GROMACS: the MD Simulation Engine . . . . . . . . . . . . 132.2.1 Initialize GROMACS Files and Programs . . . . . . . 132.2.2 Energy Minimization . . . . . . . . . . . . . . . . . . 15ivTable of Contents2.2.3 Running the Production Simulation . . . . . . . . . . 162.2.4 Analysis of the Simulation Results . . . . . . . . . . . 182.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 Study of Cholesterol in Tethered Membranes . . . . . . . . 193.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 193.2 CGMD Simulation Setup and Membrane Formation . . . . . 223.2.1 CGMD Model . . . . . . . . . . . . . . . . . . . . . . 223.2.2 Coarse-Grained Molecular Dynamics Simulation Pro-tocol . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.2.3 Formation of Membrane in Experiment . . . . . . . . 253.3 Laterial Di usion Dynamics of Lipids and Cholesterol . . . . 263.3.1 Numerical Solution of Volterra Di erential-Integral E-quation . . . . . . . . . . . . . . . . . . . . . . . . . . 273.3.2 CGMD Simulation Results on Di usion Dynamics . . 313.4 Water Density Pro le at Bioelectronic Interface . . . . . . . 333.4.1 Derivation of Percus-Yevick Equation . . . . . . . . . 343.4.2 Match the Percus-Yevick Equation with CGMD Re-sults . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.5 Archaebacterial Membrane Biomechanics . . . . . . . . . . . 383.5.1 Archaebacterial Membrane Thickness . . . . . . . . . 383.5.2 Archaebacterial Membrane Line Tension and SurfaceTension . . . . . . . . . . . . . . . . . . . . . . . . . . 403.6 Fractional Order Macroscopic Model . . . . . . . . . . . . . . 423.6.1 Fractional Order Macroscopic Model and ParameterEstimation . . . . . . . . . . . . . . . . . . . . . . . . 433.6.2 Experimental Results Utilizing the Macroscopic Mod-el . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484 Summary and Future Work . . . . . . . . . . . . . . . . . . . 51Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .A Matlab Code to Calculate the Percus Yevick Pair Distribu-tion Function for Spherical Particles . . . . . . . . . . . . . . 60vList of Tables3.1 Lipid and Cholesterol Di usion (nm2/s) . . . . . . . . . . . 323.2 Biomechanic Parameters of Membrane . . . . . . . . . . . . . 393.3 Macroscopic Model Parameters for Archaebacterial Membrane 48viList of Figures1.1 Overview of biosensor . . . . . . . . . . . . . . . . . . . . . . 21.2 Schematic of natural cell membrane . . . . . . . . . . . . . . 31.3 Schematic of ve levels of abstraction for membrane models . 52.1 Some MARTINI mapping examples . . . . . . . . . . . . . . . 123.1 A brief summary of recent related works on study of choles-terol in membrane using MD/CGMD simulation methods . . 203.2 Ball structures of Dphpc, GDPE lpids, and Cholesterol . . . . 213.3 Coarse graiend molecular dynamics structure . . . . . . . . . 243.4 Computed mean-square displacement . . . . . . . . . . . . . . 313.5 The computed water density computed from the CGMD sim-ulation results, and analytical results from the Percus-Yevickequation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.6 Ribbon structure of archaebacterial membrane . . . . . . . . 413.7 Fractional order macroscopic model of the tethered archae-bacterial membrane . . . . . . . . . . . . . . . . . . . . . . . 453.8 Parameter estimation using least squares method . . . . . . . 473.9 Experimentally measured (gray dots) and numerically pre-dicted current response for tethered membranes containingcholesterol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49viiAcknowledgementsI would like to rst thank my supervisor Prof.Vikram Krishnamurthy atUBC. Not only he o ers me valuable directions and support for my graduatestudy, but also his great passion and insights on science and engineeringencourage me to make more contributions to our community.I am truly honored to work with Dr.William Hoiles, as a beginner inthis eld, I got valuable training from him. His keen insights, comprehensiveknowledge and kind patience showed me the characteristics that an excellentscholar would possess.I am very grateful for my family members. Their encouragements andloving supports are so important, that I could have been a totally di erentperson without them.Finally I would like to express my appreciation to my dear friends andcolleagues, it's always good time to talk to them and learn from them.viiiDedicationTo my family: Shaoan, Defang, Bing, my cousins and my grandparents.ixChapter 1Introduction1.1 Motivation: the BiosensorIn 1956, Dr.Leland C. Clark conducted his famous experiment on oxygen de-tection using platinum, where he placed the enzyme Glucose Oxidase (GOD)close to the surface of platinum, thus to trap oxygen against the electrodesusing a piece of dialysis membrane [13]. It was the rst time that elemen-t recognition was achieved by using bio-materials. Since then, substantialresearch has emerged to investigate the element recognition process involv-ing bio-materials, which has further led to the developments of biologicalsensors, or the simply-called "biosensors".Commercialized biosensors are already in use in areas such as healthcaremonitoring, industrial processing and environmental pollution control [50].Even though this is the case, further improvements on biosensors to achievehigher sensitivity of element detection, wider scope of applications, and low-er production costs for biosensors are strongly demanded by markets [51].This necessitates a thorough understanding of element recognition processinvolved in biosensors.A standard biosensor consists of a recognition component called bio-receptor, and a corresponding transducer component. A bioreceptor inter-acts with the desired chemical analytes and the interaction is detected by thetransducer. The transducer produces a measurable electrical signal, which issent to an electronic system composed of an ampli er and a display circuit.The biosensor structure is summarized in Fig. 1.1. The class of biorecep-tor includes a variety of biomolecules such as membrane proteins, enzymes,antibodies and nucleic acids. An electrochemical transducer is normally ametal electrode, while other types such as optical and thermal transducersare also in use [20]. More details on choices of transducers are introducedin Sec. 1.3.Since many of the element recognition processes occur at the cell mem-brane, investigation of cell membrane properties has played an importantrole in biosensor research. Furthermore, due to the diculty of using livinganimal cell membranes, arti cial lipid membranes have become the main11.2. Cell Membrane Systemsource to obtain insights into the recognition processes occurring in theproximity of the cell membrane. In this thesis, a study on synthetic archae-bacterial membrane will be presented.Figure 1.1: Overview of a biosensor, consisting of a bioreceptor (membraneproteins) and a transducer (metal electrode/bioelectronic interface), whichare connected to an ampli er and a display circuit1.2 Cell Membrane SystemThe cell membrane is composed of three primary components: membranelipids, macromolecules, and the cytoskeletal laments [44]. Fig. 1.2 summa-rizes the membrane structure.Membrane lipids are a group of biological compounds that possess am-phiphilic property, which means that each lipid molecule has one water-soluble end and one non-water-soluble end. The amphiphilic property forceslipid molecules to group together to form a double-layer structure in waterenvironment, so a cell membrane is often referred as "lipid bilayer" in manycontexts. Di erent components of membrane lipids have di erent e ectson cell shape, permeability and organization of macromolecules [44]. Mem-brane lipids contain a wide variety of biological compounds, the most com-mon ones include phospholipids and sterols such as cholesterol. For arti cialmembranes, various synthetic phospholipid derivatives can be used [16]. Inthis work, a synthetic archaebacterial membrane composed of zwittrion-ic C20 diphytanylether-glycero-phosphatidylcholine lipid (DphPC) and C20diphytanyl-diglyceride ether lipid (GDPE) is investigated.Macromolecules refer to a group of large (more than 1000 atoms) moleculesincluding nucleic acids, carbonhydrates and proteins. In the case of cellmembrane, macromolecules usually refer to proteins. The membrane pro-21.3. Bioelectronic Interfaceteins can be classi ed into two functional classes: transport protein andmembrane receptor [44]. As one type of bioreceptors, membrane recep-tor proteins transmit electrical signals between intracecullar and extracel-lular environments. Transport proteins, on the other hand, aid in trans-membrane movements of ions and large molecules by forming channels/-pores. Furthermore, cytoskeletal laments provide the physical structuralsupport for the membrane.This thesis studies the properties of lipids, such as di usion dynamics,membrane thickness and surface/line tension. Therefore, a simple membranemodel (archaebacterial membrane) that only includes membrane lipids isconstructed.Figure 1.2: Schematic of a natural biological membrane. The extracellularuid represents the contents outside the cell, and the cytosol is the interi-or of the cell with the membrane separating the two domains. Note thatbiological membranes are composed of thousands of di erent components(macromolecules, lipids, chemical species)1.3 Bioelectronic InterfaceElectrical instrumentation is required to measure the electrical signals gener-ated from arti cial membrane. Thus to perform measurements for a biosen-sor, we need to introduce a transducer, that is, a bioelectronic interface toconnect the biological system to the electrical instrumentation.We consider a tethered membrane system connected to a bioelectronic31.4. Abstractions in Bio-System Modelinginterface, which is the standard set up for an electrochemical biosensor. Aparticularly useful bioelectronic interface comprises of inert gold electrodes.They have advantages over redox active electrodes such as biopolymers,because redox active electrodes will tend to force the tethers to dissociatefrom the electrode surface, thus to destroy the membrane. Also, redoxactive electrodes will release metal ions into solution, which can interferewith the electrophysiological response of proteins and peptides [14]. Theinert gold electrode, however, capacitively couples the electronic domain tothe physiological domain without the issues associated with redox electrodes.However, with the presence of gold interface, the so-called "di usionlimited e ect" need to be considered in modeling of tethered membranesystem. The di usion limited e ect refers to the fact that a fast reaction willbe triggered if membrane components di use and contact the gold surface,such e ect at the interface can cause the charge transfer of bio-molecules,which are highly undesirable [53]. The problem of di usion limited e ect,and its solution will be discussed with more details in Sec. 3.6 of this thesis.1.4 Abstractions in Bio-System ModelingAs mentioned above, a bioelectronic interface should connect the tetheredmembrane system and electric instrumentation. However, the problem is,how to make interpretations of electrical measurements thus to explain mem-brane characteristics? The investigation on membrane system is to studymicroscopic properties such as di usion dynamics and membrane thickness,but the data that can be measured is macroscopic, such as voltage and cur-rent. Thus a key to the development of novel membrane biosensor is anaccurate mathematical model of the cell membrane, which is able to inter-pret the macroscopic measurements and reects corresponding microscopicmembrane properties correctly. Such a model must link the microscop-ic dynamics of water, lipids, peptides to experimental measurements at amacroscopic time and length scale.There are ve levels of modeling abstraction of mathematical models thathave been used for modeling dynamics of cell membrane. For these abstrac-tions, more the details considered for the system, more is the computationalpower and time required for the simulation, and vice versa. A schematicdiagram for the ve levels of abstraction for bio-systems is summarized inFig. 1.3.41.4. Abstractions in Bio-System ModelingFigure 1.3: Schematic diagram illustrating the length and timescale achiev-able by the atomistic to macroscopic simulation methods1.4.1 Ab-initio Molecular DynamicsThe "Ab-initio" molecular dynamics, or the "from the beginning" moleculardynamics model, is the model that takes into account the most comprehen-sive details for a system. Ab-initio includes both classical Newton's physicsand quantum mechanical Schrodinger's equations for the dynamics of tar-get system, which might include particles such as water, ions, membranelipids, proteins, and peptides. This is the model that provides the mostdetailed description of a real system. However, due to these (sometimesunnecessarily) detailed equations, an excessive computational power is re-quired. Thus Ab-initio molecular dynamics could only attain membranelength scale of nanometer and simulation time in an order of femtosecond,which is too small compared to the data that can be obtained from any ex-perimental measurement [39]. Therefore such method is not quite popularin bio-system modeling so far.1.4.2 Molecular DynamicsMolecular dynamics (MD) is a simpli cation of the above "Ab-initio" molec-ular dynamics model. The quantum mechanical Schrodinger's equationsconsidered in "Ab-initio" molecular dynamics are ignored in this model,that is, the matrix representation of semi-empirical potential from quan-tum mechanics is not used in MD, and only empirical potentials associatedwith chemical bonds, bond angles and non-bonded forces are considered intoMD [42]. More details of molecular dynamics would be discussed in the nextchapter. It's also worth noting that, although only Newtonian equations areevaluated for each time step of MD simulation, still merely a length orderof nanometer and time order of nanosecond can be achieved by using MD.51.4. Abstractions in Bio-System Modeling1.4.3 Coarse-Grained Molecular DynamicsCoarse-Grained Molecular Dynamics (CGMD) is a further simpli cationbased on molecular dynamics. By grouping certain atoms together intocoarse-grained beads, with the bead-to-bead interactions empirically pa-rameterized, membrane dynamics can still be evaluated using Newton's e-quations of motion [47]. Such simpli cation allows CGMD simulations toachieve simulation time scale of microseconds with a system size of tens ofnanometers, which are good enough to match most real membrane dynamicsthat can be measured by experiments. The main results presented in thisthesis would be based on CGMD simulations.1.4.4 Continuum TheoryDespite their di erent levels of abstraction, the above three molecular dy-namics models are considering discrete particles of a system. Thus if thenumber of particles in target system is large, the eciency of computationcan be limited. Therefore, a more simpli ed class of models is to treat thediscrete entities as continuous densities, which represent the space-time aver-age of the microscopic motion of particles. The most well-known continuummodel should be the Poisson-Nernst-Planck system of equations for di usionprocess of ion transport, which combines the Poisson equation from electro-statics, and the Nernst-Planck equation for di usion [60]. Such a signi cantsimpli cation allows continuum models to achieve simulation time scale ofthe the order of microseconds, with a system size of micrometers.1.4.5 Macroscopic ModelAs its name suggests, a macroscopic model describes a bio-system systemby using macroscopic parameters, which can be obtained from experimentalmeasurements. Normally these parameters can be de ned to be any entity,and there's even no speci c physical interpretation required for them: aslong as these estimated parameters can t the experimental measurements,the macroscopic model can be regarded to be successful. In this thesis, afractional macroscopic model is introduced to estimate biological param-eters, such as tethered membrane conductance and capacitance based onexperimental measurements. Due to the free choice of parameters for themacroscopic models, there is no limit to the attainable scales of length andtime.61.5. Thesis Contributions1.5 Thesis ContributionsThis thesis is composed of four chapters, which include an introductorychapter for fundamental molecular dynamics, and a chapter on the study ofcholesterol in tethered membrane using CGMD simulations.Chapter 2 provides an introduction to the fundamental newtonian physic-s involved in molecular dynamics algorithms, along with common force elds(potential functions) including MARTINI force eld used in CGMD simula-tions. Then a four-step MD simulation scheme is introduced along with thesimulation engine GROMACS. Details on con guration of GROMACS lesand programs are provided, with an introduction to common terminologiesand algorithms used in MD simulations.Chapter 3 presents the core research work in this thesis. The e ectsof cholesterol on tethered archaebacterial membrane system are studied byusing Coarse-Grained molecular dynamics simulations. The archaebacterialmembrane consists of 70% DphPc and 30% GDPE lipids, which are teth-ered onto a gold bioelectronic interface. By following the simulation schemeintroduced in Chapter 2, the CGMD simulation based on MARTINI force eld is used to study properties of membrane lipids, as cholesterol level variesfrom 0% to 50%.The following results are presented regarding membrane properties af-fected by cholesterol: In order to match the simulation model with real membrane used inexperiments, a CGMD model for a tethered archaebacterial membraneis introduced. The CGMD simulation protocol is provided followingthe four-step scheme introduced in Chapter 2. To make reliable measurements in experiments, it is crucial to gureout the region that is not a ected by the bioelectronic interface, whichcan change the di usion dynamics of membrane and water. Thus theposition-dependent water density pro le at the bioelectronic interfaceis studied. An analytical solution of Percus-Yevick equation is ob-tained, to show that the water density pro le doesn't vary much as faras 4nm from the interface, and short-range interactions in Lennard-Jones potential contribute little to a ect water density pro le. Lipids undergo three regimes of di usion: ballistic, subdifusion andFickian di usion, which are predicted by mode-coupling theory. Al-though the three di usion dynamics can be obtained by Volterra integral-di erential equation, however, due to the excessive computational pow-71.5. Thesis Contributionser required by its numerical solution, CGMD simulation results areused to analyze di usion dynamics of membrane instead. It is shownthat cholesterol concentration does not a ect the transition time be-tween the subdi usion and Fickian di usion, and cholesterol concen-tration only a ects the di usion coecient of lipids in Fickian di usionregime. As cholesterol concentration increases, the lipid di usion co-ecients decrease. The e ects that cholesterol have on biomechanics properties of mem-brane are also presented by using CGMD simulation results. It isshown that as the cholesterol concentration increases from 0% to 50%,the membrane thickness increases, the membrane line tension increas-es, while the membrane surface tension increases up to 30% cholesterolin membrane then decreases. To link the macroscopic experimental measurements with microscop-ic simulation results, and to account for the di usion limited e ectsdue to the bioelectronic interface presented in the tethered membranesystem, a fractional order macroscopic model is introduced. The teth-ered membrane is modeled by a fractional order RC circuit, whoseparameters are obtained by a least-square estimator.Finally, Chapter 4 briey summarizes the key results, and commentson the suggestions for future research related to the work presented in thisthesis.8Chapter 2Fundamental MolecularDynamics2.1 Basic Algorithm and Force FieldsCurrently there are two classes of simulation techniques in computationalchemistry and biology: molecular dynamics (MD) and Monte Carlo simu-lation [22]. Monte Carlo molecular modeling method starts from an initialmicrostate, then keeps moving to the desired state according to the desiredensemble's Boltzman probability distribution [22]. In compasiron, MD is amore universal technique, especially in non-equilibrium ensemble and anal-ysis of dynamic scenarios. Ever since it was rst proposed by Alder andWainwright in the late 1950's [1] to study the interactions of hard spheres,MD has been used extensively for studying the structural and dynamicalproperties of molecules.2.1.1 Global MD AlgorithmAn MD algorithm starts by assigning initial positions and velocities to allparticles within the system being considered. The algorithm looks for aglobal optimal value of energy potential function for the initial system. Astep called "energy minimization" guarantees that the initial positions andvelocities of all particles can be accepted (i.e states of atoms are within therange pre-de ned by the simulation engine). Finally, for each time step, aset of classical Newtonian equations of motion is solved for ithparticle of asystem, which includes N interacting particles [28]:mi@2ri@t2= Fidridt= vividt=dvidt(2.1)92.1. Basic Algorithm and Force Fieldswhere miis the mass of the ithparticle, t the time, rithe vector of parti-cle relative to the origin, and Fithe interactive forces applied on the par-ticle. The interactive forces are the derivatives of the potential functionE(r1; r2; :::rN):Fi=@E@ri(2.2)where the potential function is given as a combination of bonded and non-bonded interactions between the particles. A potential function can beshown as follows:E(r1; :::; rN) =Xbondska(l  l0)2+Xbondskb(  0)2+Xtorsionkc[1 + cos(n!  )]+Xatompairs4ij[(ffijrij)12 (ffijrij)6]+Xatompairskqiqjrij(2.3)The rst and second term are the covalently bonding energies, which areinduced by deformations of bond length l and bond angle , respectively;the third term describes the energy due to rotation around the chemicalbond, the rst three terms represent the bonded potential. The last twoterms describe the non-bonded interactions between all atom pairs, whichinclude Lennard-Jones potential and the Coulomb electrostatic potential,respectively. ka, kb, kcare all constants related to speci c atom interactions, is the permittivity constant and k is the Coulomb constant. r is thedistance between atoms, while ff is the distance at which inter-atom potentialis zero.By solving the newtonian equations 2.1, the particle trajectories whichare represented by coordinates ri=1;2;:::;Ncan be obtained. Based on theinitial conditions: potential interaction E, initial positions r and velocitiesv , with temperature and pressure con gured to be desired values, the sim-ulation keeps updating the coordinates of the atoms by solving Newtonianequations, thus the coordinates of particles are generated at speci c timesteps.102.1. Basic Algorithm and Force FieldsTo summarize, a molecular dynamics simulation involves the followingthree steps:1. Input initial conditions, which include pre-de ned potential functions(force elds), initial positions and velocities of particles2. Compute forces by solving equation 2.23. Update the system states (positions, velocities and energies of all par-ticles) by solving equation 2.12.1.2 Force FieldsAs described in equation 2.3, force eld (or equivalently, potential func-tion) can be fully described by bonded and non-bonded interactions. Beforestarting any molecular dynamics simulation, we need to pre-de ne the pa-rameters that describe the potential function properly, thus to compute theforces and update the system trajectories. Di erent force elds are madeand used for di erent purposes. The most comprehensive force eld is the so-called "all-atom" force eld{that is, all atoms are considered to be assignedinteraction parameters [69]. However, limited by computational power andsimulation time, mainstream force elds usually ignore some unimportantatoms such as non-polar hydrogens, which contribute insigni cant interac-tions to a system. In this section, one common widely used force eld isintroduced, which will be used in the simulations in Chapter 3 of this thesis.Coarse Grained Force Fields: MARTINIThe signi cance of molecular dynamics simulations is based on the fact that,mechanical and chemical processes at atomistic level are able to reect realstates of system at macroscopic level. Therefore, an MD simulation needsto match a real experiment at a reasonable scale of space and time. Thismeans that, an ideal MD simulation should not only take into account max-imum details for topology (nano-meter scale) of the target system, but alsoproceed with long enough time that is comparable to the real experiment,which is at time scale of at least micro-second. However, a bio-system forstudy can contain thousands of particles, thus to keep essential topology ofthe system, simulation time is limited to nano-seconds due to restriction ofcomputational power.A compromising solution for such problem is to use coarse-graining{thatis, to select atoms that share common chemical/biological characteristics,112.1. Basic Algorithm and Force Fieldsthen group them together and treat each group as a single bead. Clearly,such approximation will reduce degrees of freedoms of the original system,and ignore many inside interactions between the atoms. It is expected thatsuch approximation might induce errors of simulation results. However, itturns out that, with clever choices of grouped atoms, coarse-grained modelsare able to generate simulation results that are in excellent agreements withexperiments, while the time scale of simulation can be signi cantly extendedto micro-seconds.The MARTINI force eld [47] is a popular coarse-grained molecular dy-namics force eld, which is suitable for biomolecular systems. It follows asimple four-to-one mapping, that is, four selected atoms are grouped andtreated as one coarse-grained bead, and only polar, non-polar, apolar andcharged interactions are de ned in the force eld. Some MARTINI mappingexamples are visualized in Fig 2.1 [47].Figure 2.1: Martini mapping examples of selected molecules: (A) Stan-dard water particle representing four water molecules, (B) Polarizable wa-ter molecule with embedded charges, (C) DMPC lipid, (D) Polysaccharidefragment, (E) Peptide, (F) DNA fragment, (G) Polystyrene fragment, (H)Fullerene molecule.122.2. GROMACS: the MD Simulation Engine2.2 GROMACS: the MD Simulation EngineAmong currently available simulation packages of molecular dynamics, GRO-MACS (GROningen MAchine for Chemical Simulations)[6] is known for be-ing lightweight, fast, free and open source. As an ecient engine to performmolecular dynamics simulations and energy minimization, GROMACS hasvarious applications in computational biology and molecular modeling[6].To accomplish a molecular dynamics simulation, GROMACS follows a four-step scheme to obtain the simulation results. In the following subsections,each of the four steps will be introduced in details.2.2.1 Initialize GROMACS Files and ProgramsGROMACS programs perform MD simulation steps by reading and writing les in several speci c formats [69], which are listed below: xxx.mdp: the main parameter input le that records parameters andconditions for energy minimization, position restraints, and main MDsimulation step xxx.pdb: short for the Protein DataBank le format,this input lecontains molecular structure and coordinate information of the parti-cles xxx.gro: similar to .pdb le, this input le contains molecular structurein Gromos 86 format, which indexes particles according to their types xxx.itp: an input topology le that de nes the system particles' topolo-gies, such as charge, mass, and radius of particles xxx.top: the input topology le, which de nes the force elds, particletypes and number of particles, its particles speci cations are based onthe .itp le xxx.tpr: the input le for simulation steps, it contains starting coor-dinates and velocities of the system particles xxx.xtc: the output le generated by production run, it includes tra-jectory information of all system particles for each time step xxx.edr: the output le generated by production run, it contains en-ergy information at the end of simulation time132.2. GROMACS: the MD Simulation EngineFurthermore, being a Unix-based program package, GROMACS com-mand follows standard Unix/Linux command protocol, the most commonfunctional routines are: mdrun: the main computational chemistry engine to perform molecu-lar dynamics simulation and energy minimization editconf: a program that convert one le format to another grompp: before any production run, this program could take .top and.mdp le as inputs and generate the input le .tpr for the productionrun genbox: add solvants (water in most cases) into the coordinates lelike .gro les g msd: it takes .xtc le as input, and compute mean square displace-ment from the trajectory le, thus to compute di usion constant ofspeci c particlesInitial ConditionsBefore any MD simulation, the topology le (.itp, .top) has to be properlyedited, thus to load the force eld, which de nes the interactions betweenparticles of the system. Then three parameters need to be initialized: coor-dinates and velocities of the particles, along with the pre-de ned box size.These parameters will be included in .gro le. Setting initial velocities isoptional, if velocities are not initialized manually, then GROMACS wouldassign initial velocities to the particles following MaxwellBoltzmann distri-bution in statistical mechanics [69], given temperature T :p(vi) =rmi2BTexp(miv2i2BT) (2.4)where B is the Boltzman constant.Periodic Boundary ConditionsOne possible problem in molecular dynamics simulation is the "edge e ec-t": the simulated particles in a system might move out of the pre-de nedsimulation box. To solve such problem, periodic boundary conditions areintroduced into GROMACS, that is, the original single box is replaced bya box-array, which contains multiple translated copy of the same unit. For142.2. GROMACS: the MD Simulation Enginethe convenience of de ning various systems, GROMACS allows di erent s-tandard shapes for cell units (while for eciency, rectangular cells are themost popular).It is also worth noting that, to solve the problem that there are multi-ple particle images inside translated cells, GROMACS follows the so-called"minimum image" convention, that is, only the nearest particle image is con-sidered for short-range non-bonded interactions. Speci cally, a cut-o radiusRcutis applied to truncate short range non-bonded interactions, where Rcutis less than half the minimum box vector [69]:Rcut<12min(jjajj; jjbjj; jjcjj) (2.5)where a; b; c are the vectors de ning the simulation box. The cut-o radiusis set in the .mdp le.Solvating the SystemFor the simulations of bio-systems, such as cell membrane, the whole systemmust be immersed into water to imitate a real bio-system. The program toadd solvent water is "genbox". After solvating the system, topology les(.top) need to be edited to record the added water molecules. The numberof added water molecules can be obtained by checking the output .gro legenerated by genbox.2.2.2 Energy MinimizationIf initial state of the system (i.e initial velocities and positions of particles)is out of acceptable range, the force induced by the interactions will betoo large for a simulation to start, due to a large value of potential energy.Thus a step of potential energy minimization is required to make sure thesimulation can be started properly, as the system state will be close toequilibrium if its potential energy is approximately minimized.This is a typical optimization problem that searches for the global min-imum point of the potential function given by (2.3). There are numerousoptimization algorithms (conjugate gradient, L-BFGS) available for energyminimization. In GROMACS, the most common one is the steepest descentmethod [69]. It chooses a step in the direction of the negative gradient,which guarantees the step to be descending, thus if the step size is proper,the global minimum point would be reached after enough steps of simulation.To make sure all simulation steps start with a close-to-equilibrium state,energy minimization must be done whenever a system is changed (i.e after152.2. GROMACS: the MD Simulation Enginea new system le is created or after solvents are added to a system). Theparameters involved in energy minimization are de ned in a speci c .mdp le, which is basically same as the .mdp le used in production run, exceptthat the "integrator" needs to be set to be "steep".2.2.3 Running the Production SimulationOnce the system is well-equilibrated via energy minimization, it is ready forthe production run of MD simulation. The program "grompp" records all theparameters in .mdp le to set the simulation in the desired environment, anexample production run .mdp le containing common parameters is shownas follows:integrator =mdnsteps =500dt =0:002nstlist =10rlist =0:9coulombtype =pmercoulomb =0:9vdw  type =cut offrvdw =0:9tcoupl =Berendsentc grps =protein non proteintau t =0:1 0:1ref  t =298 298Pcoupl =Berendsentau p =1:0 1:0compressibility =5e 5 5e 5ref  p =1:0(2.6)This .mdp le initializes the simulation environment, such that equationsof motion are integrated using the so-called leap-frog integrator ("integrator =md"), which is the default integrator in GROMACS, in 500 time steps with162.2. GROMACS: the MD Simulation Engine0.002 ps for each step. "nstlist = 10" means the system is updated for ev-ery 10 steps, rlist is the cut-o distance (in nm), beyond which short-rangenon-bonded interactions for a certain particle would be ignored. Similarly,"rcoulomb" and "rvdw" set the cut-o distance for coulomb force and vander waals interactions, with units in nm. "coulomtype" de nes the summa-tion method used to calculate total energy of long range electrostatics, formore details, reader can refer to [69]. The rest of above parameters involvea concept of temperature and pressure coupling in GROMACS, which isintroduced as follows.Groups: Temperature and Pressure CouplingGenerally there are four types of ensembles for simulation of bio-systems: NVE: Constant number of particles (N), system volume (V) and energy(E) NVT: Constant number of particles (N), system volume (V) and tem-perature (T) NPT: Constant number of particles (N), system pressure (P) and tem-perature (T) NAPxyT: Constant number of particles (N), system area (A), systempressure (P) along x-y direction, and temperature (T)Since isothermal and isobaric simulations (NPT) are the most relevantto experimentl data, the temperature and pressure should be controlledproperly in simulations. However, due to a result of integration errors andheating e ects from interactive forces, temperature and pressure tend to driftslightly inside the simulated system. To solve such problems, GROMACSallows user to de ne groups to control the temperature and pressure.For example, in the .mdp le presented in 2.6, two particle groups "pro-tein" and "non-protein" are pre-de ned, thus temperature-coupling andpressure-coupling algorithms can be applied to control the desired tempera-ture/pressure. The most common coupling algorithm is the Berendsen weakcoupling [69], it corrects the deviation of temperature T and pressure P withrates according to:dTdt=T0 TfitdPd=P0 Tfip(2.7)where T0, P0are the reference temperature and pressure, fitand fipare pre-de ned constants for correction, with units in ps.172.3. Summary2.2.4 Analysis of the Simulation ResultsThe production run will output the trajectory le (.xtc) to record all thesystem state in each time step, and energy le (.edr) that contains all the en-ergy terms at the end of simulation. With these information as a function oftime, various system properties such as di usion coecient, thickness, sur-face tension and line tension of membrane can be obtained. Visualizationtools such as pymol and VMD [31] can be used to visualize simulation result-s. Although output les can be processed and analyzed by any user-de nedprogram, GROMACS o ers various programs to analyze the simulation re-sults. For example, g msd can obtain mean square displacement of particles,which is used to calculate di usion coecients. g energy can extract poten-tial energy information, thus to calculate line tension and surface tension ofa membrane system.2.3 SummaryIn this chapter, the global MD algorithm is explained using classical Newto-nian physics equations. A coarse grained MARTINI force eld is introducedin details. The MD simulation engine GROMACS is discussed by introduc-ing its commonly used les and programs, along with details of con guration.Then a four-step scheme of MD simulation is introduced: Set up GROMACS les and programs with proper initializations Minimize potential energy to ensure the system start at a near-equilibriumstate The production run simulation with con gured .mdp le Analysis of the simulation results using various toolsAs shown in next chapter, although CGMD simulation applies quite sim-pli ed physics to particles of a bio-system, by following the above four stepsin GROMACS, the obtained simulation results are in excellent agreementswith experimental measurements.18Chapter 3Study of Cholesterol inTethered Membranes3.1 IntroductionAs a major sterol component in most eukaryotic membranes, Cholesterol(C27H45OH) is an important membrane components, which can regulatemembrane properties, such as lipid di usion and membrane stability [34, 49].Although signi cant research work has focused on the e ects that choles-terol has on eukaryotic membranes, little attention has been paid to howcholesterol a ects the properties of archaebacterial membranes. The aim ofthis chapter is to present a study of cholesterol on a multicomponent syn-thetic archaebacterial membrane, where CGMD simulation makes crucialcontribution to justify the experimental results.From coarse-grained molecular dynamics simulation, it can be shownthat archaebacterial membrane properties including electroporation, di u-sion dynamics, and biomechanics, are all a ected by the concentration ofcholesterol present in cell membrane. The experimental data also con rmsthe simulation results. Given the unique structure of archaebacterial lipid-s (i.e. the hydrocarbon chains are ctionalized with methyl groups), thecholesterol content has a noticeably di erent e ect on the membrane prop-erties compared to that of eukaryotic and prokaryotic membranes. Theresults provide novel insights into how cholesterol a ects archaebacterialmembrane properties, which are of use for the design of synthetic biomimet-ic membranes [57].Related WorksStudies of cholesterol span over the past several decades, and have involvedboth experimental measurements and molecular dynamics techniques. Fora review the reader can be referred to [34, 49] and citations therein. Dueto the microscopic scale of length and time involved in membrane research,atomistic-level simulations are necessary for predicting the dynamic proper-193.1. Introductionties of membranes, which are dicult to be measured experimentally. Sig-ni cant insights have been gained based on the use of molecular dynamicssimulations. In [5, 10, 17, 52] it is shown that increasing the cholesterolcontent in POPC, DOPC, DSPC, and DPPC membranes will result in anincrease in membrane stability, that is, an increasing physical resistance tomembrane defects. While increasing the cholesterol in DOPC membranesreduces the membrane line tension [7]. In DPPC membranes, increasing theconcentration of cholesterol will cause an increase in membrane thickness,and a decrease in the lateral di usion of lipids [5, 25].Recently coarse-grained molecular dynamics (CGMD) has been appliedto study the atomistic e ects cholesterol has on POPC, DOPC, and DPPCmembranes [15]. The CGMD results show that for cholesterol concentra-tions between 0% to 40% the membrane thickness decreases, however for50% the membrane thickness increases:possibly as a result of the interdigi-tation between lipid tails resulting from the free space under cholesterol [15].So far a substantial amount of work has focused on lipids containing PO,DO, DS, and DP lipid tails. However, an important question is, what e ectsdoes cholesterol have on lipids with phytanyl tails, which are typically foundin acrhaebacterial membranes? By using experimental measurements fromblack lipid membranes, Uitert et. al. [68] show that low concentrations ofcholesterol increase the membrane stability, however above 20% the mem-brane stability begins to decrease. Several questions remain of how choles-terol a ects the dynamics and biomechanics of archaebacterial membranescomposed of phytanyl lipids.Recent related works are summarized in Fig. 3.1.Figure 3.1: A brief summary of recent related works on study of cholesterolin membrane using MD/CGMD simulation methods203.1. IntroductionSimulation Setup and Key ResultsIn this thesis a CGMD model based on the MARTINI force eld [45, 46] isconstructed to study the di usion dynamics and biomechanics of archaebac-terial membranes containing cholesterol, whose concentrations range from0% to 50%. The archaebacterial membrane is composed of 70% DphPC and30% GDPE lipids (both of which contain phytanyl tails). Furthermore, tosimulate a tethered archaebacterial membrane, two bioelectronic gold sur-face attaching membrane are included in the CGMD model. In Fig. 3.2, thebrief ball structures of Dphpc, GDPE and cholesterol are provided, we canobserve that compared with Dphpc and GDPE lipids, cholesterol is smallerin size.Figure 3.2: Ball structures of Dphpc, GDPE lpids, and CholesterolFirst, the CGMD simulation results are used to compute the membraneproperties such as: di usion dynamics of lipids, membrane thickness, surfacetension, and line tension of the archaebacterial membrane. In [19], it wasshown that lipids undergo three primary regimes of di usion: the ballistic,subdi usion, and Fickian di usion. By using the results from mode-couplingtheory and CGMD simulations, it will be shown that the transition time be-tween the subdi usion and Fickian di usion regimes is not dependent onthe concentration of cholesterol in the membrane. Therefore, the choles-terol only a ects the di usion coecient of lipids in the Fickian di usionregime. To validate the CGMD simulation results, we use experimentalmeasurements from tethered archaebacterial membranes containing di er-ent concentrations of cholesterol.Second, to ensure that the tethers and bioelectronic interface would nota ect the membrane response, we must be able to compute the positiondependent density of water at the bioelectronic interface. The fact thatdensity has negligible variations at the membrane surface suggests that thebioelctronic interface does not impact the dynamics: this result is validatedusing di usion measurements of the lipids in the distal and proximal layer213.2. CGMD Simulation Setup and Membrane Formationof the archaebacterial membrane. Additionally, by using approximationsto the Yvan-Born-Green integral equation [72], we will show that harshrepulsive forces play a negligible role in the long range dynamics of theposition dependent density pro le of water at the bioelectronic interface.Finally, in order to link the microscopic molecular dynamics simulationresults and macroscopic experimental results, a fractional order macroscopicmodel is introduced. Experimental measurements are performed by mea-suring the current response of the tethered archaebacterial membranes to agiven voltage excitation.3.2 CGMD Simulation Setup and MembraneFormationTo investigate the e ects of Cholesterol on archaebacterial membranes, aCGMD model of an archaebacterial membrane containing cholesterol is con-structed. By using the MARTINI force eld [45, 46] in simulation, thedynamics of lipid di usion, membrane thickness, surface tension, and linetension will be studied. The CGMDmodel is validated by experimental mea-surements from a tethered archaebacterial lipid membrane, while a fractionalorder macroscopic model is to be introduced, thus to allow experimental re-sults to be compared with simulation results. The CGMD model, along withthe experimental measurements, provides key insights into how cholesteroland the bioelectronic interface impact the dynamics and biomechanics ofarchaebacterial membranes.3.2.1 CGMD ModelTo model the tethered archaebacterial membrane, we use the MARTINIforce eld [45, 46], which is a popular CGMD force eld designed for biomolec-ular systems. The main set-up of the MARTINI force- eld is to map ap-proximately four heavy non-hydrogen atoms into one coarse-grained bead.Normally, each bead has a mass of 72 amu. For example, in the CGMDmodel four water atoms are represented by a single coarse-grained bead. Tokeep the model as simple as possible, only four main interactions are de- ned: polar (P), non-polar (N), apolar (C) and charged (Q), each of themhas several subtypes. Q and N types have four subtypes, Qda, Qd, Qa, Q0and Nda, Nd, Na, N0, which mean they have di erent hydrogen-bondingcapabilities of the atom group: da = donor or acceptor, d = donor, a =acceptor, 0 = no hydrogen bonding. On the other hand, P and C types223.2. CGMD Simulation Setup and Membrane Formationhave ve subtypes, P1, P2, P3, P4, P5and C1, C2, C3, C4, C5, where thesubscripts 1-5 denote their increasing polar anity [45, 46].The CGMD model is constructed to imitate the essential dynamics ofthe tethered archaebacterial membrane system, which is composed of lipids,cholesterol, a gold bioelectronic interface, tethers and spacers. Given the factthat tether density of the tethered archaebacterial membrane is only aboutless than 1% (i.e. for every 268 lipid in the proximal layer one is tethered),the contribution of the tethers and spacers is regarded to be negligible, andtherefore the tethering is not included in the CGMD model. The mappingof the lipids and gold surface into the MARTINI force eld is provided asfollows:Lipids: The molecular components include the zwittrionic C20 diphytanyl-ether-glycero-phosphatidylcholine lipid (DphPC), C20 diphytanyl-diglycerideether lipid (GDPE), cholesterol, and the gold surface. The tethered archae-bacterial membrane we consider is composed of a 30% GDPE to 70% DphPCratio of lipids. The lipid ratio is identical to that used for the experimentalmeasurements. The phosphatidylcholine headgroup of the DphPC lipid isrepresented by two beads: the positive choline by the Q0bead, and thenegative phosphate by the Qabead. The ether-glycol is represented by aSNabead, and each of the phytanyl tails by four C1beads.The phytanyland either glycerol moieties of GDPE are represented by the same mappingas for the DphPC, however the hydroxyl headgroup of GDPE is represent-ed by a P4bead. In total the DphPC lipid is composed of 12 beads, andthe GDPE lipid by 11 beads. The cholesterol is represented by 8 beads asde ned in [15].Gold Surface: The gold surface is composed of a square lattice withcustom Pfbeads. The distance between adjacent beads is 0.3 nm. Theinteraction of the Pfbead is designed to reduce the e ects of excess adsorp-tion to the surface. The interaction between Pfand P4is 1/3 the valuebetween P4and P4, and the interaction between Pfand other bead typesis  12% of the MARTINI value between P4and respective bead types.The following interactions are excluded: interaction between Pfbeads, andbetween the C5beads of the tethers and spacers, and P4and Qobeads ofthe lipids. Note that a similar interaction is used in [43] to represent thegold surface in the MARTINI force- eld.The complete CGMD simulation structure is provided in Fig.3.3 for ref-erence. The tethering reservoir is selected to have a height of 4 nm to matchthe experimentally measured tethering reservoir thickness from [23].233.2. CGMD Simulation Setup and Membrane FormationBioelectronic InterfaceBulk ElectrolyteTethering Reservoirhm4 nmFigure 3.3: Coarse grained molecular dynamics structure of 0% tetheredDphPC membrane with hmdenoting the membrane thickness. Lipid tailsare represented by the green beads, Qabead is displayed in blue, the Qobead in orange, P4hydroxyl headgroup of GDPE by the red bead, the SNabead as pink, and the water beads as a translucent blue. The gold surface isindicated by the gold planes. Cholesterol is represented by the black beads.The coloring scheme of the axis is red for x, green for y, and blue for z.3.2.2 Coarse-Grained Molecular Dynamics SimulationProtocolThe molecular dynamics simulations were performed using GROMACS [24]version 4.6.2 with the MARTINI force eld [45, 46]. The interactions of theCGMD beads are de ned by the Lennard-Jones (LJ) potential, and harmon-ic potentials (constraints=none) are utilized for bond and angle interactions.A shift function is added to the Coulombic force (coulombtype=shift) to s-moothly and continuously decay to zero from 0 nm (rcoulomb-switch) to1.2 nm (rcoulomb). The LJ interactions were treated likewise except thatthe shift function was turned on between 0.9 nm (rvdw-switch) and 1.4 nm243.2. CGMD Simulation Setup and Membrane Formation(rvdw). The grid-type neighbour searching algorithm is utilized for the sim-ulation, that is, atoms in the neighbouring grid were updated every 10 timesteps. The equations of motion are integrated using the leapfrog algorithmwith a time step of 0.02 ps. Periodic boundary conditions are implementedin xy-dimensions (Fig. 3.3). Simulations are performed in the NAPzT en-semble using a temperature of 350 K to match that used in [43] for similarmembrane structures. The temperature is held constant using a velocityrescaling algorithm [9] with a time constant of 0.5 ps. Furthermore, Berend-sen pressure coupling was applied with semi-isotropic type. The lipid andwater molecules are coupled separately for temperature and pressure con-trol. The gold surface is modeled using the walls option in GROMACS.Note that CGMD simulation times are reported as e ective time, that is,four times the actual simulation time. The e ective time is introduced toaccount for the speed-up in the CGMD model [45, 46].All the systems studied here were rst energy minimized using the steep-est descent method in GROMACS. A 50 ns equilibration run is performedprior to the production run. Production runs are performed for a simulationtime of up to 1 s. Visualization of the CGMD results are reported usingPyMOL.3.2.3 Formation of Membrane in ExperimentThe tethered archaebacterial membrane is constructed using the solvent-exchange formation process presented in [28,29]. The tethered membrane issupported on a polycarbonate slide containing a 100 nm thick sputtered goldelectrodes each with dimensions 0:73 mm. The formation of the tetheredmembrane proceeds in two steps.First, benzyl disul de tetra-ethyleneglycol and benzyl disul de tetra-ethyleneglycol are xed to the surface of the gold electrode. Speci cal-ly, an ethanolic solution containing 370 M of 1% benzyl disul de tetra-ethyleneglycol and 99% benzyl disul de tetra-ethyleneglycol is exposed tothe gold surface for 30 min, then the surface is ushed with ethanol and airdried for approximately 2 min.The second stage involves the formation of the tethered archaebacterialmembrane. 8 L of 3 mM ethanolic solution containing a mixture of 70%DphPC (zwitterionic C20 diphytanylether-glycerophosphatidylcholine) and30% GDPE (C20 diphytanyl-diglyceride ether) lipids, and the rest choles-terol is brought into contact with the gold surface from the rst step. Thissolution is incubated for 2 min at 20C in allowing the formation of thetethered archaebacterial membrane. Proceeding the 2 min incubation, 300253.3. Laterial Di usion Dynamics of Lipids and CholesterolL of phosphate bu ered saline solution at a pH of 7.2 is ushed throughthe chamber.The tethered archaebacterial membrane is equilibrated for 30 min priorto performing any experimental measurements. The quality of the tetheredarchaebacterial membrane is measured continuously using an SDx tetheredmembranes tethaPodTM swept frequency impedance reader operating atfrequencies of 1000, 500, 200,100,40,20,10,5,2,1,0.5,0.1 Hz and an excitationpotential of 20 mV (SDx Tethered Membranes, Roseville, Sydney).3.3 Laterial Di usion Dynamics of Lipids andCholesterolThe di usion dynamics of a homogeneous medium like water can be de-scribed using standard Fickian di usion, where the mean-square displace-ment (MSD) is proportional to time:MSD = h(x(t) x0)2i = 4Dt (3.1)Here xois the initial position, x(t) is the current position at time t, D is thedi usion coecient and hi is the ensemble average, the ensemble is takenover all particles.However, as a result of the polymeric structure of lipids, the dynamic-s of lipids are more complex than simple liquids like water. Recently theso-called Mode-coupling theory (MCT), which is originally used for investi-gating dynamics of glass-forming liquids, has been applied to describe thedi usion dynamics of lipids in membranes [19]. In this section, we provide amethod to model the di usion dynamics of lipids, to see if the lipids are inthe ballistic, subdi usion, or Fickian di usion regimes, by using the resultsfrom CGMD simulationsAs predicted by MCT theory, the time evolution of the mean-squaredisplacement of lipids is given by a generalized form of Fickian di usion:h(x(t) xo)2i / t (3.2)where is the power-law exponent. For standard Fickian di usion = 1,in which case the proportionality constant in Eq. 3.2 is related to the lipiddi usion coecient D by ht Einstein relation D = h(x(t) xo)2i=4t.The MCT theory [12] predicts that at the femtosecond timescale, thelipids are in a ballistic regime where = 2. As time evolves, the lipids enterthe subdi usive region, with < 1, as a result of local caging e ects from263.3. Laterial Di usion Dynamics of Lipids and Cholesterolneighboring aggregated lipids. Then at the nanosecond timescale, the lipidsenter the Fickian di usion regime with = 1, allowing for the di usioncoecient to be computed as in Eq. 3.1. An estimate of the power-law coef- cient in Eq. 3.2 can be computed from the CGMD simulation trajectoriesusing the following relation: (t) =@ ln(h(x(t) xo)2i)@ ln(t): (3.3)To model the di usion dynamics of the lipids, we consider the dynamicsto satisfy the following linear Volterra integro-di erential equation, whichdescribes the time dependence of MSD in lipid bilayers [19]:@h(x(t))2i@t+tZ0M(t s)h(x(s))2ids = 4(kBTmL)t;M(t) =(t)fi3+Bet=fi1(1 + (t=fi2) ); (3.4)for xo= 0.where M(t) is the memory kernel, kBis Bolztmann's constant, T istemperature, mLis the mass of the lipid, and fi1; fi2; fi3; and B are modelparameters. fi3is the transition time at which the MSD transitions fromthe ballistic region to subdi usion region, fi2the onset of the subdi usionregion, and fi1the transition from subdi usion to Fickian di usion region.The di usion coecient D can be computed from the results of Eq. 3.4 inthe Fickian di usion region, where for t ! 1 the di usion coecient isrelated to the MSD by h(x(t))2i  4Dt, and can be evaluated using:D = (kBTmL)h1ZoM(t)dti1: (3.5)Notice that Eq. 3.5 provides another method to compute the di usion coef- cient D compared to the standard Eienstein relation given below Eq. 3.2.3.3.1 Numerical Solution of Volterra Di erential-IntegralEquationAlthough to our best knowledge, there is not any ecient method to solveEqn. 3.4 analytically, it is still possible to obtain an approximate numerical273.3. Laterial Di usion Dynamics of Lipids and Cholesterolsolution, where numerical di erentiation and numerical integration can becalculated by using nite di erence method:f0(t) f(t+ h) f(t)h(3.6)where the subdivided interval of integration h is represented by:h =tN aN(3.7)here N is the number of subdivided intervals chosen and tNis the last timestep. Similarly, tiis the ithtime step. Then by using the trapezoidal rule:Ztn0f(x)  h(12(f(x1) + 2[f(x2) + f(x3) + :::f(xn1)] +12fn) (3.8)and simpson's rule:Zt2n0f(x)dx h3[f0+ 4(f1+ f3+ :::+ f2n1)+2(f2+ f4+ f2n2) + f2n](3.9)the integration term in volterra integral-di erential equation 3.4 can beapproximated by plugging in trapezoidal rule and simpson's rule for oddand even time steps, respectively, then clearly we have the approximatedintegrals for each time step:Zt10M(t s)h(x(s))2ids h[12M(t1 s1)h(x(s1))2i]Zt20M(t s)h(x(s))2ids h3[M(t2 s0)h(x(s0))2ids+4M(t2 s1)h(x(s1))2ids+M(t2 s2)h(x(s2))2ids]283.3. Laterial Di usion Dynamics of Lipids and CholesterolZt30M(t s)h(x(s))2ids h[12M(t3 s0)h(x(s0))2ids+M(t3 s1)h(x(s1))2ids+M(t2 s2)h(x(s2))2ids+12M(t3 s3)h(x(s3))2ids]     ZtN10M(t s)h(x(s))2ids h[12M(tN1 s0)h(x(s0))2]+M(tN1 s1)h(x(s1))2]ids+ :::M(tN1 sN2)h(x(sN2))2]ids+12M(tN1 sN1)h(x(sN1))2]idsZtN0M(t s)h(x(s))2ids h3[M(tN s0)h(x(s0))2ids+4[M(tN s1)h(x(sN))2ids+ :::M(tN sN1)h(x(sN1))2ids]2[M(tN s1)h(x(sN))2ids+ :::M(tN sN2)h(x(sN2))2ids]+M(t2 s2)h(x(s2))2]ids(3.10)So if we re-write the above approximated representations in a more com-pact form, that is, if we renameM(tisj) to beMij, and rename h(x(sj))2]ito be xj, then the numerical volterra integral-di erential equation 3.4 canbe written as:x2 x02h=4kBTmLt1h2[M10x0+M11x1]x3 x12h=4kBTmLt2h3[M20x0+ 4M21x1+M22x2]     xN xN22h=4kBTmLtN1h2[MN1;0x0+2[MN1;1x1+ :::MN1;N2xN2] +MN1;N1xN1]293.3. Laterial Di usion Dynamics of Lipids and CholesterolxN xN12h=4kBTmLtNh3[MN0x0+ 4[MN1x1+ :::MN;N1xN1]+2[MN2x2+ :::MN;N2xN2] +MN;NxN](3.11)Therefore the Eqn 3.11 can be re-written in a matrix form AX = B, whereX = [x1; x2; :::xN]0, andA =266664h2k111 0 : : : 008h23k21+2h23k22 1 1 : : : 0: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :0 : : : AN1;N21 00 0 : : : AN;N11377775(3.12)withAN1;N2= 2[h2MN1;1+ :::h2MN1;N2] + h2MN1;N1andAN;N1=8h23[MN1+MN3+:::MN;N1]+4h23[MN2+MN4+:::MN;N2]+2h23MNNFurthermore,B =266666666666648hkBTmLt1 h2M10x1+ 18hkBTmLt22h23M20x28hkBTmLt3 h2M10x38hkBTmLtN1 h2MN1;0xN18hkBTmLtN2h23MN0xN37777777777775(3.13)Now we are able to solve the matrix equation AX = B and obtain the mean-square displacements as a function of time. Thus the ballistic, subdi usiveand Fickian region can be observed. However, as mentioned above, the bal-listic region and subdi usive region will occur at timescale of femtosecond,or 106nanosecond. This means that if we require a 1 ms computation, thedimension of matrix A is at least 109 109. Although A is a sparse matrix,still the required computational power and time for solving this equation isconsidered to be excessive. Further work of matrix sparsi cation is expectedto be done regarding this part. Since in this thesis, we are only consideringthe e ects of cholesterol on Fickian di usion dynamics of membrane, CGMDsimulation results are used to analyze the di usion dynamics of lipids.303.3. Laterial Di usion Dynamics of Lipids and Cholesterol3.3.2 CGMD Simulation Results on Di usion DynamicsIn this subsection, the results of CGMD simulations are utilized to gaininsights into the concentration e ects, that cholesterol has on the di usiondynamics of lipids in the archaebacterial membrane.The key question addressed in this section is to con rm that, if the lipiddynamics in CGMD model indeed shows a short time ballistic regime, anextended subdi usion regime, and a Fickian di usion regime, as predictedby the MCT theory. Furthermore, if these regimes are present, then whatis the transition time between each regime. Knowledge of these transitiontimes is crucial, as based on Eq. 3.1, the di usion coecients can only becomputed using CGMD trajectories in the Fickian di usion regime.10−1 100 101 102 10310−1100101102103Time [ns]MSD[nm2]t [ns]〈(x(t))2〉[nm2]Subdiffusion Fickianβ ≈ 0.5 β ≈ 1Figure 3.4: Computed mean-square displacement for the DphPC, GDPE,and cholesterol for the 0% and 50% cholesterol membranes with de ne inEq. 3.3. Notice that the di usion dynamics are in the subdi usion regime(  0:5) for t  3 ns, and for t  20 ns the di usion dynamics are inthe Fickian regime (  1) . This is in agreement with the mode-couplingtheory for exible macromolecules [11, 19].Fig. 3.4 presents the numerically computed mean-square displacementof DphPC, GDPE, and cholesterol from the CGMD simulation results forcholesterol concentrations of 0% and 50%, which are generated from the313.3. Laterial Di usion Dynamics of Lipids and Cholesteroltrajectory .xtc le of simulation. From Fig.3.4, we see that for t  3 nsthe molecules di use in the subdi usion regime, and for t  20 ns thedi usion of the molecules is in the Fickian di usion regime. This is inagreement with the results predicted using mode-coupling theory for exiblemacromolecules [11, 19]. Given the time step of the CGMD simulation is 20fs, the ballistic region > 1 is not observed in Fig. 3.4 for any of the lipidsor cholesterol. This is expected as the ballistic region is typically observedfor t  10 fs [19].Surprisingly, the transition times between the subdi usion and Fickiandi usion regime is not dependent on the cholesterol content. This suggeststhat the concentration of cholesterol presented contributed negligibly to thecaging e ect, that is, only the Fickian di usion dynamics are strongly de-pendent on the concentration of cholesterol present. The results in [18, 25]for SOPC, SLPC, SAPC, SDPC, and DPPC suggest that as the cholesterolcontent increases there should be a decrease in the di usion coecient D.Table 3.1: Lipid and Cholesterol Di usion (nm2/s)DphPC GDPE Cholesterol0% 69.63.2 57.42.7 -10% 69.00.4 56.32.4 106.540.920% 49.31.3 44.42.6 92.533.730% 50.03.8 48.51.2 86.829.140% 54.510 67.57.3 55.18.650% 39.51.6 31.31.8 46.78.2To gain insight into how the concentration of cholesterol a ects the d-i usion dynamics in the Fickian regime, we compute the di usion coe-cient D for DphPC, GDPE, and cholesterol for archaebacterial membranescontaining 0% to 50% cholesterol. The results are provided in Table 3.1.The di usion coecient of DphPC for 0% cholesterol is in excellent agree-ment with the experimentally measured di usion coecient of 18:1  5:6nm2/s [4]. Furthermore, the numerically computed di usion coecientof cholesterol are in excellent agreement with the experimentally measureddi usion coecient of cholesterol which are in the range of 10 nm2/s to100 nm2/s [56, 63]. The cholesterol has a higher di usion coecient thanDphPC and GPDE, and the di usion of cholesterol monotonically decreasesas the concentration of cholesterol increases. If we consider only the massand size of cholesterol, it is expected that the lower mass and size of choles-323.4. Water Density Pro le at Bioelectronic Interfaceterol compared to DphPC and GDPE will allow cholesterol to have a largerdi usion coecient. Another contributing factor is that the headgroup ofDphPC and GDPE both have a larger dipole moment than that of choles-terol which will also reduce the di usion coecient of the lipids comparedto the cholesterol [25].An interesting question is, why a "slow" lipid system could be evenslowed down by faster-moving cholesterol, which has higher di usion co-ecients than GDPE and Dphpc molecules? Such phenomenon can bewell explained by the free volume theory [67], which is used to predict aliquid-ordered-liquid-disordered (lo ld) coexistence region in phosphatidyl-choline/cholesterol mixture [55]. If the environment temperature is abovethe phase transition temperature, the membrane lipids are in the so-calleddisordered phase (ld), where the acyl chains of the phospholipid moleculesare in a disordered state that contains a high fraction of gauche conformer-s [2]. When cholesterol are added into the membrane, cholesterol rst spanover the hydrocarbon cores of the bilayer, if concentration of cholesterol in-creases, free-volum models [54, 55] predict that excessive cholesterol wouldpack tightly into the lipids, thus the cavities and deects from the membraneare lled by cholesterols, and more regions of membrane are converted in-to the ordered state (lo), which has much less free volume, and therefore,smaller di usion coecients.3.4 Water Density Pro le at BioelectronicInterfaceBefore conducting actual experimental measurements for the tethered mem-brane, a crucial problem is to gure out the region where the properties oftethered membrane are not a ected by bioelectronic interface. Thus it be-comes necessary that water density pro le nearby the gold surface should beinvestigated: the membrane should be set in the region where water densitypro le is not a ected much by the bioelectronic interface.As the dynamics of water in proximity to the bioelectronic interface havevery di erent characteristics from that in the bulk region [3, 21, 76], the ther-modynamic and structural properties of water molecules at the interface area ected by two primary factors: a smaller number of neighboring moleculeinteractions, and a change in the potential energy of the uid as a result ofinteractions with the surface. The density pro le of water near the interfacetypically consists of oscillations that are similar to sinusoidal waves , withperiod close to the mean thickness of each water layer in proximity to the333.4. Water Density Pro le at Bioelectronic Interfaceinterface [21, 76]. This observation suggests that the oscillations occur at asimilar length scale to the molecular diameter of the water molecules.In our simulation the interface is modeled using a lattice of coarse-grainedbeads which interact with the water beads via a Lennard-Jones potential.Then a key question we would like to answer is: does the density pro lefrom the CGMD model match that from a hard-sphere uid at a hard-wall interface? In this section we will show that, approximations to theYvan-Born-Green integral equation [72], which lead to an analytical resultrepresented by the so-called "Percus-Yevick Equation", will be able to e-valuate the density pro le of a hard-sphere uid at a hard-wall interface.This density pro le can then be compared with the density pro le from theCGMD simulation, thus to evaluate the e ects of harsh repulsive forces thatCGMD water beads have on the water density.Now let's provide a derivation of the Percus-Yevick equation rst.3.4.1 Derivation of Percus-Yevick EquationFirst we consider a monoatomic uid with density , which is interactingwith pair potential u, and an external potential ffi that represents the con-tribution of interactions from the interface. Then the mean force acting onparticle 1 by a particle 2 can be represented by the gradient of potentialr1u(r1; r2). The total mean force on particle 1 by the other particles andinterface is given by the Yvan-Born-Green integral equation [72]:kBTr1ln((r1)) =r1ffi(r1)Zr1u(r1; r2)(r2jr1)dr2(3.14)where kBis the Boltzmann's constant, T is the temperature, and (r2jr1)the conditional singlet density (i.e. the conditional density at r2given aparticle is xed at r1).In (3.14) the so-called correlation function is represented as function ofthe external eld ffi. Note that the right hand side of (3.14) is the averageforce acting on a particle xed at r1. To estimate (r), approximations to(3.14) are to be made to only account for long-range interactions. The mainidea is to construct an external reference potential ffiR, which only accountsfor long range interactions of the particles and the interface, and pair inter-action that only accounts for short range interactions, which is denoted asthe reduced system. The only requirement is that the density pro le of thefull system and reduced system must be equal. Using this method, Weekset. al. [72] propose the approximation that the singlet density for both the343.4. Water Density Pro le at Bioelectronic Interfacefull system and reduced system must be similar at short range. The externalpotential ffi and reference potential ffiRare then related by [21, 72]:ffiR(r1) = ffi(r1) +Z[R(r2) B]u1(r1; r2)dr2(3.15)with Bthe bulk density, u1(r1; r2) the attractive part of the pair potential,and Rthe density in the external reference potential. Note that for agiven (r), which can be computed from the CGMD simulation, an e ectivereference potential ffiRcan be evaluated by solving (3.15) self-consistently.To solve for (r) analytically we must determine an expression for ffiRin(3.15). It is reasonable (we will see why very soon) to approximate ffiRby ffiR0where we have neglected the harsh repulsive forces at the interfaceleaving only the attractive forces in the full system. Then, for a hard-wall,ffi =1 for r  1 and 0 otherwise, the density (r) is given by [71]:(r1) = B+ BZ[(r1) B]c(r12)dr2;(r1) = 0 for r1 1 (3.16)with c() the direct correlation function of the uid, r12= r1r2. For a hard-sphere, the direct correlation function is given by a cubic polynomial whichis dependent on the radius of the hard-sphere, denoted by R, and the bulkwater density B[73]. For a hard-sphere and hard-wall interface Eq. 3.16 isgiven by the Percus-Yevick equation which can be solved analytically usingthe method in [70], which is shown as follows.The close form of Percus-Yevick equation is derived directly from Ornstein-Zernike equation, which simply states that total inuence h of particle i onparticle j is a superposition of both direct and indirect inuence:h(rij) = c(rij) + BZdrkc(rik)[g(rkj) 1] = g(r) 1 (3.17)where g(r) is the "scaled" density of particles, or the pair distribution func-tion; the direct correlation function c(rij) describes the direct inuence ofparticle i on particle j; and second term of the equation, which is the indi-rect correlation function, suggests that the indirect inuence of particle i onparticle j is an integrated result of particle i acting on a reference particlek, which in turn has inuence on j.By assuming the case of hard sphere potential:u(r) =(1 if r < b0 if r  b:(3.18)353.4. Water Density Pro le at Bioelectronic Interfacewhere b is the radius of particle, after de ning function y(r) to be:h(r) c(r) = y(r) 1;withy(r) =(c(r) if r < bg(r) if r  b:(3.19)it's straightforward to conclude thatc(r) =(y(r) if r < b0 if r  b:(3.20)Then by setting particle j to be the original point, rename ri= r, rj= r0,we can express Ornstein-Zernike equation in terms of y(r) [65]:y(r) = 1 + BZr0<bdr0y(r0)[1 g(r  r0)] (3.21)compared with Eq. 3.16, clearly they are equivalent if we set (r)=B=g(r).Wertheim had obtained a close-form analytical solution for the P-Y e-quation. It was suggested that the solution for direct correlation function isactually a cubic polynomial: c(r) = + (r=b) + (r=b)3(3.22)with =(1 + 2f)2(1 f)4; = 6f(1 + f=2)2(1 f)4; =f(1 + f)22(1 f)4; f =Bb36for x  1, where f is de ned as the fractional volume of the particle.As the indirect correlation function in the Ornstein-Zernike equation is aconvolution integral, then in Fourier Transform domain, it's proportional tothe product of fourier transforms of h(r) and c(r), we could take advantageof this property.H(p) =1(2)3Z11dreiprh(r)C(p) =1(2)3Z11dreiprc(r) (3.23)363.4. Water Density Pro le at Bioelectronic InterfaceSolving the Ornstein-Zernike equation in Fourier Transform domain givesus:H(p) =C(p)1 B(2)3C(p)Then by doing inverse fourier transform:h(r) =Z11dpeiptH(p) = 4Z10dp(sinprpr)p2H(p) (3.24)Finally, the pair distribution function g(r) is given by g(r) = h(r) + 1 [65].3.4.2 Match the Percus-Yevick Equation with CGMDResultsFig. 3.5 shows the water density pro le from the CGMD simulation, alongwith the density pro le predicted by the Percus-Yevick equation (Eq. 3.16).Notice that the two density pro les in Fig. 3.5 are in excellent agreement.Remember that the solution of Percus-Yevick equation 3.24 does not involveshort-range interactions, this suggests that harsh short-range interactions inthe LJ potential do not play a signi cant role in the density pro le. Theradius of the hard-sphere in the Percus-Yevick equation is b = 0:495 nm,which is in agreement with the potential well of the LJ potential of theCGMD water of 0:47 nm. As seen in Fig. 3.5, the distance between the localminima/maxima is approximately b, as expected from the discussion in [21,76]. At 4 nm from the bioelectronic interface, there is negligible variation inthe density pro le, which suggests that the lipid di usion coecients shouldnot be e ected by the interface. To validate this claim, we computed thedi usion coecients in the proximal and distal layers and found them to bein excellent agreement. For example, for the 20% cholesterol membrane thedi usion coecients of the DphPC lipids in the proximal and distal layersare 49:1 1:2 and 49:3 1:3 nm2/s respectively.Given the fact that bioelectronic interface is in close proximity (i.e. ap-proximately 4 nm) to the membrane surface, we would like to investigatethe e ect that the bioelectronic interface have on the water dynamics anddi usion dynamics of lipids. Given the distance between the bioelectron-ic interface and membrane surface is 4 nm, with the investigations madeabove, we are able to show that the density variations are negligible at adistance of 4 nm from the interface. Indeed, we have successfully estimat-ed that the negligible importance of harsh short range interactions of thewater beads on the density pro le. Therefore the bioelectronic interface con-tributes negligibly to the di usion dynamics of the lipid membrane. Using373.5. Archaebacterial Membrane Biomechanics0 5 10 15 190500100015002000z [nm]ρ(z)[kg/m3]Membrane SurfaceCGMDPercus-YevickFigure 3.5: The computed water density computed from the CGMD sim-ulation results, and analytical results from the Percus-Yevick equation (E-q. 3.16) with R = 0:49 nm and B= 1000 kg/m3. The bioelectronic interfaceis located at r = 0 nm and r = 19:8 nm.the Percus-Yevick equation (Eq. 3.16), it is shown that the long-range inter-actions of water dominate the density pro le of the LJ uid in the CGMDmodel.3.5 Archaebacterial Membrane BiomechanicsHow does the concentration of cholesterol a ect the archaebacterial mem-brane biomechanics? By using the results from CGMD simulations, westudy how cholesterol content a ects membrane thickness hm, line tension, and surface tension ff. Recall that the fractional order macroscopic mod-el is dependent on hm; ; and ff, which allows experimental measurementsfrom the tethered archaebacterial membrane to be utilized to validate theresults from the CGMD model. The computed hm; ff; and for the 0% to50% cholesterol membranes is provided in Table 3.2.3.5.1 Archaebacterial Membrane ThicknessTable 3.2 provides the numerically computed membrane thickness for % c-holesterol from 0% to 50%. The membrane thickness for the 0% cholesterol isin agreement with the experimentally measured thickness for DphPC basedtethered membranes, which do not contain cholesterol [23]. The position383.5. Archaebacterial Membrane BiomechanicsTable 3.2: Biomechanic Parameters of Membrane0% 10% 20% 30% 40% 50%Membrane Thickness hm(nm)PO43.8 4.09 4.13 4.17 3.92 3.61OH 3.36 3.65 3.69 3.76 3.81 3.98ROH - 2.92 3.01 3.06 3.09 3.17Surface Tension (mN/m)ff 57.8 73.5 86.1 97.1 68.9 43.9Line Tension (pN) 50.8 51.8 54.8 60.5 69.4 70.1The membrane thickness is computed from the distance betweenthe molecule in the distal to the molecule in the proxyl layer of themembrane.of cholesterol's hydroxyl group (ROH) is always less than the associatedhead groups (OH and PO4) of the DphPC and GDPE. This is in agreementwith the results from molecular dynamics simulation of DPPC moleculesfor 11% and 50% cholesterol [62]. As expected, as the cholesterol contentincreases from 0% to 30%, there is a decrease in membrane thickness as aresult of the cholesterol forming a complex with the hydrocarbon tails. Asimilar e ect has been observed for cholesterol in DPPC, DMPC, SOPC,and POPC membranes (6, 45, 46). Interestingly for archaebacterial mem-branes containing 40% and 50% cholesterol there is a decrease in membranethickness. This results suggests that the phytanyl chains in the DphPC andGDPE form a complex with the cholesterol causing the membrane thicknessto decrease. Note that this is similar to the condensing e ect suggestedfor DPPC, DMPC, SOPC, and POPC membranes [17, 32, 66] in which thecholesterol forms a complex with the hydrocarbon tails. However, in thecase of lipids with a phytanyl tail, there is a decrease in thickness not anincrease.To validate the numerically computed membrane thickness, we use ex-perimental measurements from tethered archaebacterial membranes. Recallthat the capacitance of the tethered membrane Cmis dependent on the di-electric permittivity "m, thickness hm, and surface area of the membraneAm. For a constant "mand surface area Am xed to be 2.1 mm2, as hmdecreases the associated capacitance of the membrane must increase. FromTable 3.3,as the concentration of cholesterol from 0% to 30% increases the393.5. Archaebacterial Membrane Biomechanicsassociated capacitance of the membrane decreases, this suggests that themembrane thickness is increasing. From 40% to 50% cholesterol concentra-tion, the tethered membrane capacitance increases, this suggests a decreasein membrane thickness. These results validate the numerically computed ar-chaebacterial membrane thickness provided in Table 3.2 for concentrationsof cholesterol from 0% to 50%.3.5.2 Archaebacterial Membrane Line Tension and SurfaceTensionBefore introducing the macroscopic model, methods of computing the linetension and surface tension ff based on CGMD simulation results areprovided, because the calculation of membrane conductance will dependon the surface tension and line tension. The line tension is de ned asthe energy cost per unit length at the boundary, where the hydrocarbonlipid tails and water are splited. The surface tension ff is de ned as theenergy required to increase the surface area of the membrane by a unit area.Therefore, the defect density (instability) of the membrane increases as theratio =ff increases.To compute the line tension of the membrane, we use the procedureprovided in [35]. The line tension can be computed from the ribbon likestructure (Fig.3.6) using =12DLxLyPxx+ Pyy2 PzzE(3.25)with Pxx; Pyy; Pzzthe diagonal elements of the pressure tensor, Lxand Lythe simulation size in the x and y directions respectively, and h   i denotingthe ensemble average over time.To construct the lipid structure in (Fig. 3.6), an intact bilayer containing320 lipids, with a 70% DphPC and 30% GDPE composition is used. Thehydrophilic interior of the bilayer is initially adjacent to the x and z dimen-sions of the simulation cell. The simulation cell is then expanded in the xdirection from 14 nm to 16 nm, and in the y direction from 10 nm to 13nm to ensure the membrane forms an edge. Initially a 50 ns equilibrationrun was performed to allow the edge to form, this was followed by a 250 nsproduction from from which (3.25) can be estimated. Simulation are per-formed in a NPxyLzT ensemble at a temperature of 320 K. Temperature iskept constant using the velocity rescaling algorithm [9] with a time constantof 0.5 ps. Pressure is coupled semi-isotropically using the weak couplingscheme [9] with a time constant of 3 ps, compressibility of 0.3 nm2/nN, anda reference pressure of 100 kN/m2.403.5. Archaebacterial Membrane BiomechanicsFigure 3.6: Ribbon structure of archaebacterial membrane. The bead col-oring is identical to that used in Fig. 3.3. The coloring scheme of the axisis red for x, blue for y, and green for z. Note that this axis is only used forcomputing the line tension.The surface tension of the membrane is computed using [35]:ff =12DLzPzzPxx+ Pyy2E(3.26)with the parameters de ned below (3.25). The evaluation of (3.26) is per-formed in the NAPxyT ensemble using a total production run of 250 ns.From Table 3.2 the cholesterol content both a ects the surface tension andline tension of the membrane. The computed values for the surface tensionin Table 3.2 are in excellent agreement with the experimental results pro-vided in [75], and simulation results from [36, 37, 74], whose work focus onsimilar DphPC based membranes. From molecular dynamics simulations ofDPPC membranes, as the percentage of cholesterol increases there is a de-crease in the surface area of the membrane [25]. Since the surface tension isde ned as the energy required to increase the surface area of the membraneby a unit area, due to the interactions between cholesterol and lipids, if theconcentration of cholesterol increases the surface tension of the membraneis expected to increase. From Table 3.2, for 0% to 30% cholesterol the sur-face tension ff increases with increasing cholesterol content. However, it is413.6. Fractional Order Macroscopic Modelunexpected that the surface tension decreases for 40% and 50% cholesterolcontent. Furthermore, from Table 3.2 we see that as the line tension ofthe membrane monotonically increases as the cholesterol content increases.An interesting question is, why the line tension of the lipid/cholesterolmixture increases monotonically when concentration of cholesterol increas-es, while surface tension increased till 30% cholesterol, then dropped downafterwards? The process of raft formation over membrane surface is ableto provide an explanation. As it has been well studied, when cholesterol isintroduced into saturated phospholipids, a raft-like domain could be formedand enhanced as the concentration of cholesterol increases [64]. Since surfacetension is de ned as the energy cost per unit area associated with decreasingthe membrane area, when small clusters of cholesterol formed, a graduallyordered state of phospholipid/cholesterol mixture is stabilized, thus surfacetension increases. However, if the concentration of cholesterol is very high(above 30%), small clusters of cholesterol are expected to collide with eachother to form larger ones, thus the whole system is again de-stabilized bylarge clusters and surface tension decreases. Furthermore, for the line ten-sion case, Equation . 3.25 showed that line tension represents a force alongone direction (in our case, the zz-direction in Fig. 3.6), thus when more raftstructures form on the membrane, more ordered states are introduced alongzz-direction, thus line tension increases monotonically.To validate the computed surface tension ff and line tension in Table 3.2we use experimental results from tethered archaebacterial membranes. Re-call that the population of membrane defects increases as the ration =ffincreases. The population of membrane defects is given by the equilibriummembrane conductance Gm, which can be measured experimentally fromthe fractional order macroscopic model. From the experimentally measuredmembrane conductance Gmin Table 3.3, for 0% to 30% cholesterol the ra-tio =ff is expected to decrease, and from 40% to 50% =ff is expected toincrease. This result is in agreement with the numerically computed ff and in Table 3.2.3.6 Fractional Order Macroscopic ModelIn this section the fractional order macroscopic model of the tethered archae-bacterial membrane is to be provided. Experimental measurements from thetethered membrane and macroscopic model are used to validate the resultsfrom the CGMD model. The tethered membrane system studied in thisthesis consists of three di erent regions: the electrical double layers at the423.6. Fractional Order Macroscopic Modelgold electrodes, the bulk electrolyte reservoir, and the tethered membrane.The electrical double layers are composed of a bound region of ions, anda di usive region, which can be modeled using an overall capacitance Ctdlonone side of bioelectronic interface, and Cbdlon the other side. Furthermore,the bulk electrolyte reservoir can be modeled by a pure ohmic resistanceRe. Moreover, since a tethered membrane can be regarded as a uniformlypolarized structure, it could be modeled by another capacitance Cm. Addi-tionally, the so-called "electroporation" phenomenon, which states the e ectthat numerous pores can form on surface of membrane when an excitationpotential Vs(t) is presented, can be modeled by a membrane conductanceGm.The macroscopic model is based on the equivalent circuit model statedabove, however, a fractional order operator must be included to accountfor the di usion limited processes present at the bioelectronic interface ofthe tethered membrane. A di usion limited process refers to the e ect thatmembrane molecules at the interface will di use and contact the bioelectron-ic interface, thus to trigger a fast reaction of charge transfer of bio-molecules.An overview of the fractional order macroscopic model is presented in Fig.3.7.3.6.1 Fractional Order Macroscopic Model and ParameterEstimationIn this subsection, a fractional order macroscopic model of the tetheredarchaebacterial membrane is provided. This model estimates the membraneconductance Gmand capacitance Cm, from the experimentally measuredcurrent response of the tethered membrane. The membrane conductanceGmis dependent on the di usion dynamics of lipids, line tension , and thesurface tension ff; while the membrane capacitance Cmis dependent on thethickness of the membrane hm.How do we obtain the necessary parameters to calculate membrane con-ductance and inductance? Recall that the di usion dynamics, line tension,surface tension, and membrane thickness can be computed from the CGMDsimulation results. Fractional order operators are utilized in the macroscopicmodel, as the gold surface (the bioelectronic interface) of the tethered mem-brane may contain di usion-limited charge transfer, quasi-reversible chargetransfer, and ionic adsorption dynamics. These double-layer charging e ectscan be modeled using fractional order operators [8].The tethered archaebacterial membrane is composed of three distinc-t regions: the bioelectronic interface at the gold electrodes, the tethered433.6. Fractional Order Macroscopic Modelmembrane, and the bulk electrolyte solution. The membrane is assumed tobe uniformly polarizable and contain aqueous pores as a result of randomthermal uctuations. This allows the tethered membrane to be modeledby a capacitance Cmin parallel with the tethered membrane conductanceGm[23, 26, 27]. The membrane conductance is time dependent as a result ofvariations in the concentration of conducting aqueous pores in the tetheredarchaebacterial membrane. The bulk electrolyte solution is assumed to bepurely ohmic with a resistance Re.Due to the di usion limited process, there exists an electrical doublelayer [30] at the bioelectronic interface, which can be modeled using a ca-pacitor if di usion-limited charge transfer, quasireversible charge transfer,and ionic adsorption dynamics are not present. If these double-layer charg-ing e ects are present then the bioelectronic interface can be modeled usinga constant-phase-element composed of a capacitance and the fractional or-der operator p. If p < 1 then a di usion-limited process is present, and ifp = 1 then a di usion-limited process is not present. An excitation potentialVsis applied across the two electrodes of the tethered membrane and thecurrent response I is measured. The fractional order macroscopic model ofthe tethered membrane is given by (Fig. 3.7):dVmdt= (1CmRe+GmCm)Vm1CmReVdl+1CmReVs;dpVdldtp= 1CdlReVm1CdlReVdl+1CdlReVs; (3.27)I(t) =1ReVs Vm Vdl; (3.28)where Cdlis the total capacitance of Ctdland Cbdlin series with p in Eq. 3.27the fractional order operator, Vmthe transmembrane potential, and Vdlthedouble-layer charging potential.The tethered membrane conductance Gmis modeled using asymptoticapproximations to the Smoluchowski-Einstein equation [27]. The equationsgoverning the dynamics of Gmare provided as follows:443.6. Fractional Order Macroscopic ModelGm=bN(t)cXi=1Gp(ri); (3.29)dridt= DkBT@W@rifor ri2 f1; 2; : : : ; bN(t)cg;dNdt= e(VmVep)21NNoeq(VmVep)2;W (r; Vm) = 2r  ffr2+ (Cr)4+Wes(r; Vm)WmIn Eq. 3.29, is the pore creation rate coecient, Vepis the characteristicvoltage of electroporation (i.e. the voltage at which the e ects of electropo-ration are non-negligible), Nois the equilibrium pore density at Vm= 0, andq = (rm=r)2is the squared ratio of the minimum energy radius rmat Vm= 0with rthe minimum energy radius of hydrophilic pores [29, 41, 48, 61]. Wis the hydrophobic aqueous pore energy and consists of four energy terms:the pore edge energy , the membrane surface tension ff, the electrostaticinteraction between lipid heads, and the transmembrane potential energycontribution Wes. The parameters Gpand Wesare provided in [27].For a membrane containing negligible defects the parameters Cm; Cdl; p; Reare constant with Gmdependent on the transmembrane potential.Cbdl+VsGmCm+VmReCtdlI(t)Figure 3.7: Fractional order macroscopic model of the tethered archae-bacterial membrane. The circuit parameters are de ned in the TetheredArchaebacterial Membrane section.453.6. Fractional Order Macroscopic ModelFor a transmembrane potential below 50 mV, Gmcan be regarded asconstant and represents the equilibrium number of aqueous pores in thetethered membrane. To estimate the equilibrium conductance and parame-ters in Eq. 3.27, impedance measurements from the tethered membrane canbe used. For a sinusoidal drive potential Vs(t) = Vosin(2ft) with frequencyf and magnitude Vobelow 50 mV, the impedance, denoted by Z(f) of thetethered membrane is given by:Z(f) = Re+1Gm+ j2fCm+1(j2f)pCdl: (3.30)In Eq. 3.30, j denotes the imaginary part of complex numberp1. Giventhe di erence of the experimentally measured impedance and the numeri-cally predicted impedance from Eq. 3.30, the parameters in Eq. 3.27 can beevaluated by using a least-square estimator. For instance, the constrainedoptimization function "fmincon" in Matlab is able to estimate the param-eters. For the case of 0% cholesterol, by using the experimental measure-ments of impedance, fmincon gives membrane capacitance Cm= 34:4nF ,membrane conductance Gm= 0:48S, the fractional order p = 0:95. Asshown in Fig 3.8, the parameters obtained from least-square estimator arein excellent agreement with experimental data.For transmembrane potentials above 50 mV, the membrane conductanceis dependent on parameters such as the line tension and surface tension ofthe membrane. To model the dynamics of Gmand the double-layer charginge ects of the tethered membrane, Eq. 3.27 is required to be coupled with theSmoluchowski-Einstein equation 3.29, which governs the current responseI of the tethered membrane given a drive potential Vs.3.6.2 Experimental Results Utilizing the MacroscopicModelIn this subsection, experimental measurements from the tethered archae-bacterial membrane are provided, such that the CGMD simulation resultscan be validated. By using the fractional order macroscopic model, theexperimental measurements can be used to calculate the membrane conduc-tance Gmand membrane capacitance Cm. Gmis dependent on the di usioncoecient D, surface tension ff, line tension . Cmis dependent on themembrane thickness hmand area Am. Prior to all experimental measure-ments, the membrane integrity and static parameters in the fractional ordermacroscopic model are evaluated using impedance measurements. The pa-rameters are summarized in Table 3.3. The fractional order operator p is463.6. Fractional Order Macroscopic Model10−1 100 101 102 103 104103104105106107108FrequencyMagnetude of ImpedanceFigure 3.8: The data tting of estimated parameters in Eq 3.27 (red), whichare obtained from experimental impedance measurements (black dot) withVmbelow 50mV, for the membrane with 0% cholesterol.in the range of 0.95 to 0.98, this con rms that a di usion-limited processis present at the bioelectronic gold interface of the tethered archaebacterialmembrane. The associated capacitance Cdlis in the range of 120 nF to 180nF. The membrane capacitance Cm, conductance Gm, and characteristicvoltage of electroporation Vepare provided in Table 1.Fig. 3.9 provides the experimentally measured and numerically predict-ed current response of the tethered archaebacterial membrane for di erentmol % concentrations of cholesterol. The computed current is in excellent a-greement with the experimentally measured current response of the tetheredarchaebacterial membrane. This con rms the experimental measurementsand macroscopic model to estimate important biological parameters, whichcan be used to validate results from the CGMD model.Electroporation refers to a biophysical method that uses an electricalpulse to create temporary pores in cell membranes. The characteristicvoltage of electroporation, denoted by Vep, is the threshold voltage un-der which the rst pores appear on the membrane surface. An interest-ing question is, how does the concentration of cholesterol a ect Vepin thetethered archaebacterial membrane? From Table 3.3 we see that as the473.7. Summaryconcentration of cholesterol increases there is an increase in Vep. This isin agreement with the results that typically seen for membranes containingPOPC [10], DPPC [25], DphPC [68], phosphatidylcholine, phosphatidylser-ine and phosphatidyl-glycerol lipid bilayers [38], and egg yolk phosphatidyl-choline bilayer [40], in which an increase in cholesterol increases the stabilityof the membrane (i.e. increases Vepand decreases Gm).Furthermore, what e ects cholesterol has on the equilibrium membraneconductance Gm? For 0% to 30% cholesterol concentrations the membraneconductance decreases for increasing % of cholesterol. At an atomistic level,one possible explanation for the decrease in Gmis the packing of the hydro-carbon chains in the DphPC and GDPE lipids. For DphPC membranes thephytanyl tails form a tightly packed network with neighbouring hydrocarbonchains being interdigitated [33, 58, 59]. The addition of a low percentage ofcholesterol (below 30 %) causes the DphPC and GDPE molecules to con-dense the phospholipid network resulting in a decrease in Gm. Interestingly,as the concentration of cholesterol increases above 30%, there is an increasein Gm. This suggests that large concentrations of cholesterol in the archae-bacterial membrane cause the disentanglement of the phytanyl chains inthe DphPC and GDPE lipids which introduce membrane defects into themembrane [68].Table 3.3: Macroscopic Model Parameters for Archaebacterial Membrane% Cholesterol p CmGmVep0 0.95 34.4 nF 0.48S 270 mV10 0.95 32.4 nF 0.26S 290 mV20 0.95 31.4 nF 0.25S 300 mV30 0.95 31 nF 0.24S 330 mV40 0.95-0.98 35 nF 1.00S 345 mV50 0.95-0.98 41 nF 1.11S 350 mVGmis the equilibrium membrane conductance, Cmthe membrane capaci-tance, Vepthe characteristic voltage of electroporation, and p the fractionalorder parameter.3.7 SummaryIn order to gain insights into the di usion dynamics and biomechanics of ar-chaebacterial membranes containing di erent concentrations of cholesterol,483.7. Summary0 2 4 6 8 10−4−202468Time [ms]Current[µA] 0%10%20%30%40%50%Figure 3.9: Experimentally measured (gray dots) and numerically predict-ed current response for tethered membranes containing cholesterol. Theparameters for the numerical predictions are given in Table 3.3. The exci-tation Vsis given by a linearly increasing voltage with rate of 100 V/s for 5ms, then a linearly decreasing voltage with a rate of -100 V/s for 5 ms.a CGMD model based on the MARTINI force eld is constructed for a syn-thetic archaebacterial membrane, which is composed of 70% DphPC and30% GDPE lipids. By using the CGMD model and experimental measure-ments, several key insights are provided.The di usion dynamics of cell membrane lipids is investigated using tra-jectory results of CGMD simulation. It is shown that the transition timebetween the subdi usion and Fickian di usion regimes for lipids is not de-pendent on the concentration of cholesterol in the membrane. The con-centration of cholesterol only a ects the di usion coecient of lipids in theFickian di usion regime{for increasing cholesterol content, the lipid di usioncoecients decrease.To ensure that the bioelectronic interface do not e ect the membraneresponse, we computed the position dependent density of water at the bio-electronic interface. It is illustrated that the bioelectronic interface does not493.7. Summarycontribute to the dynamics and biomechanics of the archaebacterial mem-brane. Additionally, by using approximations to the Yvan-Born-Green inte-gral equation, the analytical solution of Percus-Yevick solution is matchedwith water density pro le from CGMD simulation. It is shown that harshrepulsive forces play a negligible role in the long range dynamics of theposition dependent density pro le of water at the bioelectronic interface.Biomechanic properties of membrane lipids with cholesterol are checkedbased on CGMD simulation results. For concentrations of cholesterol be-low 30% the membrane surface tension increases with increasing cholesterolcontent, however beyond 30% the surface tension decreases. For increas-ing concentrations of cholesterol from 0% to 50%, the line tension of themembrane increases. These results are validated using experimental mea-surements from tethered archaebacterial membranes, when a fractional ordermacroscopic model is used to link the simulation results with experimentalmeasurements.50Chapter 4Summary and Future WorkA study of cholesterol in tethered membrane system using coarse-grainedmolecular dynamics simulations was presented in this thesis. By discussingmodi ed properties of synthetic archaebacterial membrane with di erentconcentrations of cholesterol, this research work presented a balance of accu-racy and eciency of CGMD simulations, along with various mathematicaltechniques that are used to investigate bio-systems at atomistic level.Chapter 1 gave a brief description of biosensor, which motivates the re-search on membrane system tethered onto a bioelectronic interface. Sincethis study is based on simulations to predict membrane properties, ve lev-els of modeling abstraction in bio-system modeling were explained. Amongthese abstractions, coarse-grained molecular dynamics simulation was cho-sen for this work, due to its fair accuracy and high eciency of computation.Chapter 2 presented an introduction to concepts and theories involvedin molecular dynamics. The MD simulation engine GROMACS was intro-duced, along with an explanation on the MARTINI coarse grained force eld. Then a four-step scheme of CGMD simulation using GROMACS wasdiscussed in details, which includes: Setting up speci c formatted les and programs in GROMACS withproper con gurations for the initial state of simulated system. Minimizing the potential energy for the states whenever the system les are changed to ensure the simulation starts at a state that is closeenough to equilibrium. Executing the production run simulation after proper energy mini-mization as long as the xxx.mdp le is con gured to include parame-ters for a desired ensemble. Using the recorded coordinate and energy information as a function oftime to gain insights into characteristics of the bio-system.Chapter 3 presented the core research work in this thesis. A CGMDmod-el based on MARTINI force eld was constructed to study e ects of choles-terol on a synthetic tethered archaebacterial membrane, which is composed51Chapter 4. Summary and Future Workof 70% DphPC and 30% GDPE lipids. The di usion dynamics of lipids,thickness of membrane, surface tension, and line tension of membrane werestudied corresponding to di erent concentrations of cholesterol added to themembrane. By using the CGMD model and experimental measurements,several key insights were provided.In order to ensure that the bioelectronic surface does not a ect mem-brane response, position dependent water density pro le at the two bioelec-tronic interface was investigated. The analytical solution of Percus-Yevickequation obtained as an approximation from Yvan-Born-Green integral e-quation was discussed. It was shown that at about 4 nm from the interface,the water density pro le doesn't vary much and harsh short-range interac-tions in Lennard-Jones potential do not a ect water density pro le.Furthermore, it was observed that the concentrations of cholesterol donot a ect the transition time between the sub-di usion region and Fickiandi usion, which are predicted by MCT theory and can be computed fromVolterra integro-di erential equation. The concentration of cholesterol onlya ects the di usion coecient of lipids in the Fickian di usion regime - forincreasing cholesterol levels, the lipid di usion coecients decrease.Finally, experimental measurements were used to validate the simulationresults using a fractional order macroscopic model. It was con rmed that,as the concentrations of cholesterol increase from 0% to 50% - the di usioncoecients of tethered archaebacterial membrane decrease, the membranethickness and line tension increases, while surface tension increases up to30% cholesterol then decreases.Potential future work includes obtaining numerical solution of Volterraintegral-di erential equation up to femto-second scale, and using resonanceenergy transfer [12] to obtain experimental data for the di usion coecientof membrane components to further validate the simulation results.52Bibliography[1] BJ Alder and TEf Wainwright. Phase transition for a hard spheresystem. The Journal of Chemical Physics, 27(5):1208, 1957.[2] PFF Almeida and WLC Vaz. Lateral di usion in membranes. Handbookof biological physics, 1:305{357, 1995.[3] A. Archer and R. Evans. Relationship between local molecular eld the-ory and density functional theory for non-uniform liquids. The Journalof chemical physics, 138(1):014502, 2013.[4] T. Baba, H. Minamikawa, M. Hato, and T. Handa. Hydrationand molecular motions in synthetic phytanyl-chained glycolipid vesi-cle membranes. Biophysical journal, 81(6):3377{3386, 2001.[5] D. Bennett, J. MacCallum, and P. Tieleman. Thermodynamic analysisof the e ect of cholesterol on dipalmitoylphosphatidylcholine lipid mem-branes. Journal of the American Chemical Society, 131(5):1972{1978,2009.[6] Herman JC Berendsen, David van der Spoel, and Rudi van Drunen.Gromacs: A message-passing parallel molecular dynamics implementa-tion. Computer Physics Communications, 91(1):43{56, 1995.[7] M. Berkowitz. Detailed molecular dynamics simulations of model bi-ological membranes containing cholesterol. Biochimica et BiophysicaActa (BBA)-Biomembranes, 1788(1):86{96, 2009.[8] P. Bresslo and J. Newby. Stochastic models of intracellular transport.Review Modern Physics, 85:135{196, 2013.[9] G. Bussi, D. Donadio, and M. Parrinello. Canonical sampling throughvelocity rescaling. J. Chem. Phys., 126:014101, 2007.[10] M. Casciola, D. Bonhenry, M. Liberti, F. Apollonio, and M. Tarek. Amolecular dynamic study of cholesterol rich lipid membranes: compar-ison of electroporation protocols. Bioelectrochemistry, 100:11{17, 2014.53Bibliography[11] S. Chong, M. Aichele, H. Meyer, M. Fuchs, and journal= Baschnagel, J.Structural and conformational dynamics of supercooled polymer melts:Insights from rst-principles theory and simulations.[12] S. Chong and M. Fuchs. Mode-coupling theory for structural andconformational dynamics of polymer melts. Physical review letters,88(18):185702, 2002.[13] Jr Leland C Clark. Electrochemical device for chemical analysis,November 17 1959. US Patent 2,913,386.[14] Charles G Cran eld, Bruce A Cornell, Stephan L Grage, Paul Duck-worth, Sonia Carne, Anne S Ulrich, and Boris Martinac. Transientpotential gradients and impedance measures of tethered bilayer lipidmembranes: pore-forming peptide insertion and the e ect of electropo-ration. Biophysical journal, 106(1):182{189, 2014.[15] M. Daily, B. Olsen, P. Schlesinger, D. Ory, and N. Baker. Improvedcoarse-grained modeling of cholesterol-containing lipid bilayers. Journalof chemical theory and computation, 10(5):2137{2150, 2014.[16] Nathalie Fa, Sebastien Ronkart, Andre Schanck, Magali Deleu, An-thoula Gaigneaux, Erik Goormaghtigh, and M-P Mingeot-Leclercq. Ef-fect of the antibiotic azithromycin on thermotropic behavior of dopc ordppc bilayers. Chemistry and physics of lipids, 144(1):108{116, 2006.[17] L. Fernndez, G. Marshall, F. Sagus, and R. Reigada. Struc-tural and kinetic molecular dynamics study of electroporation incholesterol-containing bilayers. The Journal of Physical Chemistry B,114(20):6855{6865, 2010.[18] A. Filippov, G. Ordd, and G. Lindblom. Domain formation in modelmembranes studied by Pulsed-Field Gradient-NMR: The role of lipidPolyunsaturation. Biophysical Journal, 93(9):3182 { 3190, 2007.[19] E. Flenner, J. Das, M. Rheinstadter, and I. Kosztin. Subdi usion andlateral di usion coecient of lipid atoms and molecules in phospholipidbilayers. Physical Review E, 79(1):011907, 2009.[20] Inga Gustafson. Phospholipid membranes in biosensor applications:Stability, activity and kinetics of reconstituted proteins and glycolipidsin supported membranes. 2004.54Bibliography[21] J. Hansen and I. McDonald. Theory of simple liquids. Elsevier, 1990.[22] W Keith Hastings. Monte carlo sampling methods using markov chainsand their applications. Biometrika, 57(1):97{109, 1970.[23] F. Heinrich, T. Ng, D. Vanderah, P. Shekhar, M. Mihailescu, H. Nanda,and M. Losche. A new lipid anchor for sparsely tethered bilayer lipidmembranes. Langmuir, 25(7):4219{4229, 2009.[24] B. Hess, C. Kutzner, D. Spoel, and E. Lindahl. Gromacs 4: algorithmsfor highly ecient, load-balanced, and scalable molecular simulation.J. Comput. Chem, 4:435{447, 2008.[25] C. Hofsa, E. Lindahl, and O. Edholm. Molecular dynamics simu-lations of phospholipid bilayers with cholesterol. Biophysical journal,84(4):2192{2206, 2003.[26] W. Hoiles, V. Krishnamurthy, and B. Cornell. Modelling the bioelec-tronic interface in engineered tethered membranes: from biosensing toelectroporation. IEEE Transactions on Biomedical Circuits and Sys-tems, PP(99):1{13, 2014.[27] W. Hoiles, V. Krishnamurthy, C. Cran eld, and B. Cornell. An en-gineered membrane to measure electroporation: E ect of tethers andbioelectronic interface. Biophysical journal, 107(6):1339{1351, 2014.[28] Wm G Hoover. Molecular dynamics. In Molecular Dynamics, volume258, 1986.[29] Q. Hu and R. Joshi. Transmembrane voltage analyses in spheroidal cellsin response to an intense ultrashort electrical pulse. Physical Review E,79:011901, 2009.[30] Z. Hughes and T. Walsh. Structure of the electrical double layer ataqueous gold and silver interfaces for saline solutions. Journal of Colloidand Interface Science, 436(0):99 { 110, 2014.[31] William Humphrey, Andrew Dalke, and Klaus Schulten. Vmd: visualmolecular dynamics. Journal of molecular graphics, 14(1):33{38, 1996.[32] W. Hung, M. Lee, F. Chen, and H. Huang. The condensing e ectof cholesterol in lipid bilayers. Biophysical Journal, 92(11):3960{3967,2007.55Bibliography[33] T. Husslein, D. Newns, P. Pattnaik, Q. Zhong, P. Moore, and M. Klein.Constant pressure and temperature molecular-dynamics simulation ofthe hydrated diphytanolphosphatidylcholine lipid bilayer. The Journalof chemical physics, 109(7):2826{2832, 1998.[34] E. Ikonen. Cellular cholesterol tracking and compartmentalization.Nature Reviews Molecular Cell Biology, 9(2):125{138, 2008.[35] F. Jiang, Y. Bouret, and J. Kindt. Molecular dynamics simulations ofthe lipid bilayer edge. Biophys. J., 87:182{192, 2004.[36] Y. Jiang and J. Kindt. Simulations of edge behaviour in a mixed bilayer:Flctuation analysis. J. Chem. Phys., 126:045105{9, 2007.[37] J. Joannis, F. Jiang, and J. Kindt. Coarse-grained model simulationsof mixed-lipid systems: Composition and line tension of a stabilizedbilayer edge. Langmuir, 22:998{1005, 2006.[38] S. Kakorin, U. Brinkmann, and E. Neumann. Cholesterol reduces mem-brane electroporation and electric deformation of small bilayer vesicles.Biophysical chemistry, 117(2):155{171, 2005.[39] Barbara Kirchner, Philipp J di Dio, and Jurg Hutter. Real-world pre-dictions from ab initio molecular dynamics simulations. In MultiscaleMolecular Methods in Applied Chemistry, pages 109{153. Springer,2012.[40] S. Koronkiewicz and S. Kalinowski. Inuence of cholesterol on elec-troporation of bilayer lipid membranes: chronopotentiometric stud-ies. Biochimica et Biophysica Acta (BBA)-Biomembranes, 1661(2):196{203, 2004.[41] W. Krassowska and P. Filev. Modeling electroporation in a single cell.Biophysical Journal, 92(2):404 { 417, 2007.[42] Raphael D Levine. Molecular reaction dynamics. Cambridge UniversityPress, 2005.[43] C. Liu and R. Faller. Conformational, dynamical. and tensional studyof tethered bilayer lipid membranes in coarse-grained molecular simu-lations. Langmuir, 28:15907{15915, 2012.[44] Harvey F Lodish, Arnold Berk, S Lawrence Zipursky, Paul Matsudaira,David Baltimore, James Darnell, et al. Molecular cell biology, volume 4.Citeseer, 2000.56Bibliography[45] S. Marrink, H. Risselada, D. Ye mov, D. Tieleman, and A. Vries. Themartini force eld: Coarse grained model for biomolecular simulations.J. Phys.Chem.B, 111:7812{7824, 2007.[46] S. Marrink and P. Tieleman. Perspective on the MARTINI model.Chem. Soc. Rev., 42:6801{6822, 2013.[47] Siewert J Marrink, H Jelger Risselada, Serge Ye mov, D Peter Tiele-man, and Alex H De Vries. The martini force eld: coarse grainedmodel for biomolecular simulations. The Journal of Physical Chem-istry B, 111(27):7812{7824, 2007.[48] S. Movahed and D. Li. A theoretical study of single-cell electroporationin a microchannel. Journal of Membrane Biology, 246:151{160, 2013.[49] H. Ohvo-Rekila, B. Ramstedt, P. Leppimaki, and P. Slotte. Choles-terol interactions with phospholipids in membranes. Progress in lipidresearch, 41(1):66{97, 2002.[50] I Palchetti and M Mascini. Biosensor technology: a brief history. InSensors and Microsystems, pages 15{23. Springer, 2010.[51] JJ Pancrazio, JP Whelan, DA Borkholder, W Ma, and DA Stenger. De-velopment and application of cell-based biosensors. Annals of biomedi-cal engineering, 27(6):697{711, 1999.[52] R. Reigada. Electroporation of heterogeneous lipid membranes.Biochimica et Biophysica Acta (BBA)-Biomembranes, 1838(3):814{821, 2014.[53] WR Reinmuth. Theory of di usion limited charge-transfer processes inelectroanalytical techniques. Analytical Chemistry, 34(11):1446{1454,1962.[54] Mantripragada B Sankaram and Thomas E Thompson. Modulation ofphospholipid acyl chain order by cholesterol. a solid-state deuteriumnuclear magnetic resonance study. Biochemistry, 29(47):10676{10684,1990.[55] Mantripragada B Sankaram and Thomas E Thompson. Cholesterol-induced uid-phase immiscibility in membranes. Proceedings of theNational Academy of Sciences, 88(19):8686{8690, 1991.57Bibliography[56] H. Scheidt, D. Huster, and K. Gawrisch. Di usion of cholesterol andits precursors in lipid membranes studied by 1h pulsed eld gradientmagic angle spinning NMR. Biophysical Journal, 89(4):2504 { 2512,2005.[57] Y. Shen, P. Saboe, I. Sines, M. Erbakan, and M. Kumar. Biomimet-ic membranes: a review. Journal of Membrane Science, 454:359{381,2014.[58] K. Shinoda, W. Shinoda, T. Baba, and M. Mikami. Comparative molec-ular dynamics study of ether-and ester-linked phospholipid bilayers.The Journal of chemical physics, 121(19):9648{9654, 2004.[59] W. Shinoda, M. Mikami, T. Baba, and M. Hato. Molecular dynamicsstudy on the e ect of chain branching on the physical properties oflipid bilayers: structural stability. The Journal of Physical ChemistryB, 107(50):14030{14035, 2003.[60] A Singer, D Gillespie, J Norbury, and RS Eisenberg. Singular per-turbation analysis of the steady-state poisson{nernst{planck system:Applications to ion channels. European journal of applied mathematics,19(05):541{560, 2008.[61] K. Smith, J. Neu, and W. Krassowska. Model of creation and evolutionof stable electropores for DNA delivery. Biophysical Journal, 86(5):2813{ 2826, 2004.[62] A. Smondyrev and M. Berkowitz. Structure of dipalmitoylphosphatidyl-choline/cholesterol bilayer at low and high cholesterol concentrations:Molecular dynamics simulation. Biophysical Journal, 77(4):2075 { 2089,1999.[63] P. Stroeve and I. Miller. Lateral di usion of cholesterol in monolayers.Biochimica et Biophysica Acta (BBA) - Biomembranes, 401(2):157 {167, 1975.[64] Seyed R Tabaei, Joshua A Jackman, Seong-Oh Kim, Bo Liedberg,Wolfgang Knoll, Atul N Parikh, and Nam-Joon Cho. Formation ofcholesterol-rich supported membranes using solvent-assisted lipid self-assembly. Langmuir, 30(44):13345{13352, 2014.[65] Leung Tsang, Jin Au Kong, Kung-Hau Ding, and Chi On Ao. Front-matter and Index. Wiley Online Library, 2001.58[66] K. Tu, M. Klein, and D. Tobias. Constant-pressure molecular dynamicsinvestigation of cholesterol e ects in a dipalmitoylphosphatidylcholinebilayer. Biophysical Journal, 75(5):2147 { 2156, 1998.[67] David Turnbull and Morrel H Cohen. On the free-volume model of theliquid-glass transition. The journal of chemical physics, 52(6):3038{3041, 1970.[68] I. Uitert, S. Gac, and A. Berg. The inuence of di erent membrane com-ponents on the electrical stability of bilayer lipid membranes. Biochim-ica et Biophysica Acta (BBA)-Biomembranes, 1798(1):21{31, 2010.[69] David van der Spoel, Erik Lindahl, Berk Hess, AR Van Buuren, E Apol,PJ Meulenho , DP Tieleman, ALTM Sijbers, KA Feenstra, Rudi vanDrunen, et al. Gromacs user manual version 3.3. 2008.[70] E. Waisman, D. Henderson, and J. Lebowitz. Solution of the meanspherical approximation for the density pro le of a hard-sphere uidnear a wall. Molecular Physics, 32(5):1373{1381, 1976.[71] J. Weeks, K. Katsov, and K. Vollmayr. Roles of repulsive and attractiveforces in determining the structure of nonuniform liquids: Generalizedmean eld theory. Physical review letters, 81(20):4400, 1998.[72] J. Weeks, K. Vollmayr, and K. Katsov. Intermolecular forces and thestructure of uniform and nonuniform uids. Physica A: Statistical Me-chanics and its Applications, 244(1):461{475, 1997.[73] M. Wertheim. Exact solution of the percus-yevick integral equation forhard spheres. Physical Review Letters, 10(8):321{323, 1963.[74] J. Wohlert, W. Otter, O. Edholm, and W. Briels. Free energy of a trans-membrane pore calculated from atomistic molecular dynamics simula-tions. The Journal of chemical physics, 124(15):154905, 2006.[75] D. Zhelev and D. Needham. Tension-stabilized pores in giant vesicles:determination of pore size and pore line tension. Biochim. Biophys.Acta, 1147:89{104, 1993.[76] S. Zhou and A. Jamnik. Structure of inhomogeneous lennard-jones uidnear the critical region and close to the vapor-liquid coexistence curve:Monte Carlo and density-functional theory studies. Physical Review E,73(1):011202, 2006.59Appendix AMatlab Code to Calculatethe Percus Yevick PairDistribution Function forSpherical ParticlesMatlab Code1 % This code i s o r i g i n a l l y cont r ibuted by K.H. Ding(2008) and ed i t ed by2 % Yan Duan (2015) at UBC to c a l c u l a t e the PercusYevick pa i r d i s t r i b u t i o n3 % func t i on f o r s p h e r i c a l p a r t i c l e s45 c l c ;6 c l e a r a l l ;7 % I n i t i a l i z a t i o n o f parameters from Eq . ( 3 . 5 ) Eq .( 3 . 1 1 )8 n0=1;9 % Set the diameter o f s p h e r i c a l p a r t i c l e s , un i t s in nm10 dia =0.99;11 s=1;12 tau=9;13 tauc=(2 s q r t (2 ) ) /6 ;14 x i c =(3 s q r t (2 )4) /2 ;15 rm=20;16 nk=1000;17 dk=pi /(rm dia ) ;18 dr=rm dia /(nk+1) ;19 hk=ze ro s (2nk+2 ,1) ;20 vo l=pi  dia ^3/6 ;60Matlab Code21 % Frac t i ona l Volume22 fv=n0 vo l ;23 f=fv ;24 % ca l c u l a t e the parameters o f c ( r ) in Eq . ( 3 . 9 )25 denom=1 f ;26 nu=tau+f /denom ;27 ep s i l o n=f (1+ f /2) /(3denom^2) ;28 lamd=6(nus q r t (nu^2 ep s i l o n ) ) / f ;29 mu=lamd f denom ;30 i f mu > 1+2 f31 break ;32 end3334 alpha=f /denom ;35 beta=1lamd f+3alpha ;36 gamma=3lamddenom ;37 % Open the f i l e to wr i t e c ( r )38 fpk=fopen ( ' pys f . dat ' , 'w+' ) ;3940 pk=0;41 pp1=alpha(4 lamd+3alpha )+1;42 pp2=0;43 pyhk=(1/(pp1^2+pp2^2)1)/n0 ;44 pys f=1+pyhkn0 ;45 f p r i n t f ( fpk , '%6u %14.9 f nn ' , pk , pys f ) ;4647 % ca l c u l a t e h( r ) in terms o f c ( r )48 f o r i k =1:nk49 pk=ik dk ;50 x=pk dia /2 ;51 snx=s i n (x ) ;52 csx=cos (x ) ;53 ps ix=snx/x ;54 phix=3(snxx csx ) /(x^3) ;55 pp1=alpha ( beta phix+gamma ps ix )+csx ;56 pp2=alpha xphix+snx ;57 pyhk=(1/(pp1^2+pp2^2)1)/n0 ;58 pys f=1+pyhkn0 ;59 f p r i n t f ( fpk , '%6u %14.9 f nn ' , pk , pys f ) ;60 hk ( ik+1)=pkpyhk ;61Matlab Code61 end62 f c l o s e ( fpk ) ;6364 % f o u r i e r trans form65 hw=2 f f t ( hk ) ;66 hr=imag (hw( 2 : nk+1) ) ;6768 f p r=fopen ( ' pypdf . dat ' , 'w+' ) ;69 % reve r s e f o u r i e r trans form and wr i t e h( r )70 f o r i r =1:nk71 r=i r dr/ dia ;72 i f r >= 173 g=1+hr ( i r ) /(rm4 pi  r  dia ^2) ;74 f p r i n t f ( fpr , ' %14.9 f %14.9 f nn ' , r , g ) ;75 end76 end77 f c l o s e ( f p r ) ;78 % plo t pa i r d i s t r i b u t i o n func t i on79 load pypdf . dat a s c i i ;80 load pys f . dat a s c i i ;81 % sc a l e the d i s t r i b u t i o n p l o t82 f o r i =[1 : l ength ( pypdf ( : , 2 ) ) ]83 i f pypdf ( i , 1 ) +41.55>39 && pypdf ( i , 1 ) +41.55<4184 pypdf ( i , 2 )=pypdf ( i , 2 ) /1 .5+0 .31 ;85 end86 end87 % Plot the a n a l y t i c a l s o l u t i o n o f PercusYevickequat ion88 f i g u r e (1 ) ;89 x = 19.7 ones (1 ,2300) ;90 y = 1 : 2300 ;91 p lo t (x , y , ' l i n ew id th ' , 2 ) ;92 hold on ;93 x an a l y t i c a l=(pypdf ( 8 : 8 00 , 1 ) +41.55) 0 . 4755 ;94 y an a l y t i c a l=(pypdf ( 8 : 8 00 , 2 ) 0.45) 1800 ;95 p lo t ( x ana l y t i c a l , y ana l y t i c a l , ' r ' , ' l i n ew id th ' , 2 ) ;96 ax i s ( [ 1 2 , 2 0 , 0 . 0 , 2 7 0 0 ] ) ;9798 x l ab e l ( ' r (nm) ' ) ;99 y l ab e l ( ' n rho ( r ) ' ) ;62Matlab Code100 hold on ;101 % Import the CGMD s imula t i on r e s u l t s102 A=importdata ( ' dens i ty . xvg ' ) ;103104 x CGMD=A( : , 1 ) ;105 y CGMD=A( : , 2 ) ;106 p lo t (A( : , 1 ) ,A( : , 2 ) , ' ' , ' l i n ew id th ' , 2 ) ;107108 l egend ( ' Ana ly t i c a l ' , 'CGMD' ) ;109 f i l l P a g e ( gcf , ' pape r s i z e ' , [ 5 3 ] , ' margins ' , [ 0 0 0 0 ] ) ;110 x l ab e l ( ' p o s i t i o n z [nm] ' ) ;111 y l ab e l ( ' water dens i ty ' ) ;112 pr in t deps epsFig63
Thesis/Dissertation
eng
Electrical and Computer Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivs 2.5 Canada
Graduate
Study of cholesterol in tethered membrane using coarse grained molecular dynamics simulations
Text