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 eects 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 dierent 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 aect 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 conrm that harsh repulsive forces play a negligible role in thelong range dynamics of the water density prole.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 prole, numerical solutionof Volterra dierential-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 Diusion Dynamics of Lipids and Cholesterol . . . . 263.3.1 Numerical Solution of Volterra Dierential-Integral E-quation . . . . . . . . . . . . . . . . . . . . . . . . . . 273.3.2 CGMD Simulation Results on Diusion Dynamics . . 313.4 Water Density Prole 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 Diusion (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 oers 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 dierentperson 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 amplier 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, articial 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 amplier 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. Dierent components of membrane lipids have dierent eectson 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 articialmembranes, 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 classied 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 diusion 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 dierent components(macromolecules, lipids, chemical species)1.3 Bioelectronic InterfaceElectrical instrumentation is required to measure the electrical signals gener-ated from articial 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 "diusionlimited eect" need to be considered in modeling of tethered membranesystem. The diusion limited eect refers to the fact that a fast reaction willbe triggered if membrane components diuse and contact the gold surface,such eect at the interface can cause the charge transfer of bio-molecules,which are highly undesirable [53]. The problem of diusion limited eect,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 diusion 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 simplication 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 simplicationbased 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 simplication 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 dierent 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 simplied 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 diusionprocess of ion transport, which combines the Poisson equation from electro-statics, and the Nernst-Planck equation for diusion [60]. Such a signicantsimplication 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 dened to be any entity,and there's even no specic 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 conguration 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 eectsof 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 forceeld 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 aected by the bioelectronic interface, whichcan change the diusion dynamics of membrane and water. Thus theposition-dependent water density prole at the bioelectronic interfaceis studied. An analytical solution of Percus-Yevick equation is ob-tained, to show that the water density prole doesn't vary much as faras 4nm from the interface, and short-range interactions in Lennard-Jones potential contribute little to aect water density prole. Lipids undergo three regimes of diusion: ballistic, subdifusion andFickian diusion, which are predicted by mode-coupling theory. Al-though the three diusion dynamics can be obtained by Volterra integral-dierential equation, however, due to the excessive computational pow-71.5. Thesis Contributionser required by its numerical solution, CGMD simulation results areused to analyze diusion dynamics of membrane instead. It is shownthat cholesterol concentration does not aect the transition time be-tween the subdiusion and Fickian diusion, and cholesterol concen-tration only aects the diusion coecient of lipids in Fickian diusionregime. As cholesterol concentration increases, the lipid diusion co-ecients decrease. The eects 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 diusion limited eectsdue 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-dened 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 specic 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 congured 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 specic timesteps.102.1. Basic Algorithm and Force FieldsTo summarize, a molecular dynamics simulation involves the followingthree steps:1. Input initial conditions, which include pre-dened 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-dene the pa-rameters that describe the potential function properly, thus to compute theforces and update the system trajectories. Dierent force elds are madeand used for dierent 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 insignicant 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 signicance 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 signicantly 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 dened 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 writingles in several specic 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 denes the system particles' topolo-gies, such as charge, mass, and radius of particles xxx.top: the input topology le, which denes the force elds, particletypes and number of particles, its particles specications 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 diusion constant ofspecic particlesInitial ConditionsBefore any MD simulation, the topology le (.itp, .top) has to be properlyedited, thus to load the force eld, which denes the interactions betweenparticles of the system. Then three parameters need to be initialized: coor-dinates and velocities of the particles, along with the pre-dened 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 eec-t": the simulated particles in a system might move out of the pre-denedsimulation 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 dening various systems, GROMACS allows dierent 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. Specically, 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 dening 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 dened in a specic .mdple, 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" denes 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 eects from interactive forces, temperature and pressure tend to driftslightly inside the simulated system. To solve such problems, GROMACSallows user to dene 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-dened, 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-dened 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 diusion 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-denedprogram, GROMACS oers various programs to analyze the simulation re-sults. For example, g msd can obtain mean square displacement of particles,which is used to calculate diusion 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 conguration.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 congured .mdp le Analysis of the simulation results using various toolsAs shown in next chapter, although CGMD simulation applies quite sim-plied 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 diusion and membrane stability [34, 49].Although signicant research work has focused on the eects that choles-terol has on eukaryotic membranes, little attention has been paid to howcholesterol aects 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, diu-sion dynamics, and biomechanics, are all aected by the concentration ofcholesterol present in cell membrane. The experimental data also conrmsthe 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 dierent eect on the membrane prop-erties compared to that of eukaryotic and prokaryotic membranes. Theresults provide novel insights into how cholesterol aects 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-nicant 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 diusion of lipids [5, 25].Recently coarse-grained molecular dynamics (CGMD) has been appliedto study the atomistic eects 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 eectsdoes 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 aects 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 diusion 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: diusion dynamics of lipids, membrane thickness, surfacetension, and line tension of the archaebacterial membrane. In [19], it wasshown that lipids undergo three primary regimes of diusion: the ballistic,subdiusion, and Fickian diusion. By using the results from mode-couplingtheory and CGMD simulations, it will be shown that the transition time be-tween the subdiusion and Fickian diusion regimes is not dependent onthe concentration of cholesterol in the membrane. Therefore, the choles-terol only aects the diusion coecient of lipids in the Fickian diusionregime. To validate the CGMD simulation results, we use experimentalmeasurements from tethered archaebacterial membranes containing dier-ent concentrations of cholesterol.Second, to ensure that the tethers and bioelectronic interface would notaect 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 diusion 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 prole 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 eects 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 diusion, 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 dierent 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 asdened 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 eects 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 dened 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 eective time, that is,four times the actual simulation time. The eective 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 disulde tetra-ethyleneglycol and benzyl disulde tetra-ethyleneglycol are xed to the surface of the gold electrode. Specical-ly, an ethanolic solution containing 370 M of 1% benzyl disulde tetra-ethyleneglycol and 99% benzyl disulde 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 Diusion Dynamics of Lipids and CholesterolL of phosphate buered 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 Diusion Dynamics of Lipids andCholesterolThe diusion dynamics of a homogeneous medium like water can be de-scribed using standard Fickian diusion, 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 thediusion 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 thediusion dynamics of lipids in membranes [19]. In this section, we provide amethod to model the diusion dynamics of lipids, to see if the lipids are inthe ballistic, subdiusion, or Fickian diusion 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 diusion:h(x(t) xo)2i / t(3.2)where is the power-law exponent. For standard Fickian diusion = 1,in which case the proportionality constant in Eq. 3.2 is related to the lipiddiusion 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 subdiusive region, with < 1, as a result of local caging eects from263.3. Laterial Diusion Dynamics of Lipids and Cholesterolneighboring aggregated lipids. Then at the nanosecond timescale, the lipidsenter the Fickian diusion regime with = 1, allowing for the diusioncoecient 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 diusion dynamics of the lipids, we consider the dynamicsto satisfy the following linear Volterra integro-dierential 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 subdiusion region, fi2the onset of the subdiusionregion, and fi1the transition from subdiusion to Fickian diusion region.The diusion coecient D can be computed from the results of Eq. 3.4 inthe Fickian diusion region, where for t ! 1 the diusion 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 diusion coef-cient D compared to the standard Eienstein relation given below Eq. 3.2.3.3.1 Numerical Solution of Volterra Dierential-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 Diusion Dynamics of Lipids and Cholesterolsolution, where numerical dierentiation and numerical integration can becalculated by using nite dierence 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-dierential 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 Diusion 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-dierential 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 Diusion 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, subdiusiveand Fickian region can be observed. However, as mentioned above, the bal-listic region and subdiusive 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 sparsication is expectedto be done regarding this part. Since in this thesis, we are only consideringthe eects of cholesterol on Fickian diusion dynamics of membrane, CGMDsimulation results are used to analyze the diusion dynamics of lipids.303.3. Laterial Diusion Dynamics of Lipids and Cholesterol3.3.2 CGMD Simulation Results on Diusion DynamicsIn this subsection, the results of CGMD simulations are utilized to gaininsights into the concentration eects, that cholesterol has on the diusiondynamics of lipids in the archaebacterial membrane.The key question addressed in this section is to conrm that, if the lipiddynamics in CGMD model indeed shows a short time ballistic regime, anextended subdiusion regime, and a Fickian diusion 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 diusion coecients can only becomputed using CGMD trajectories in the Fickian diusion 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 dene inEq. 3.3. Notice that the diusion dynamics are in the subdiusion regime( 0:5) for t 3 ns, and for t 20 ns the diusion 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 Diusion Dynamics of Lipids and Cholesteroltrajectory .xtc le of simulation. From Fig.3.4, we see that for t 3 nsthe molecules diuse in the subdiusion regime, and for t 20 ns thediusion of the molecules is in the Fickian diusion 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 subdiusion and Fickiandiusion regime is not dependent on the cholesterol content. This suggeststhat the concentration of cholesterol presented contributed negligibly to thecaging eect, that is, only the Fickian diusion 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 diusion coecient D.Table 3.1: Lipid and Cholesterol Diusion (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 aects the d-iusion dynamics in the Fickian regime, we compute the diusion coe-cient D for DphPC, GDPE, and cholesterol for archaebacterial membranescontaining 0% to 50% cholesterol. The results are provided in Table 3.1.The diusion coecient of DphPC for 0% cholesterol is in excellent agree-ment with the experimentally measured diusion coecient of 18:1 5:6nm2/s [4]. Furthermore, the numerically computed diusion coecientof cholesterol are in excellent agreement with the experimentally measureddiusion coecient of cholesterol which are in the range of 10 nm2/s to100 nm2/s [56, 63]. The cholesterol has a higher diusion coecient thanDphPC and GPDE, and the diusion 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 Prole at Bioelectronic Interfaceterol compared to DphPC and GDPE will allow cholesterol to have a largerdiusion 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 diusion 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 diusion 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 diusion coecients.3.4 Water Density Prole 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 aected by bioelectronic interface. Thus it be-comes necessary that water density prole nearby the gold surface should beinvestigated: the membrane should be set in the region where water densityprole is not aected much by the bioelectronic interface.As the dynamics of water in proximity to the bioelectronic interface havevery dierent characteristics from that in the bulk region [3, 21, 76], the ther-modynamic and structural properties of water molecules at the interface areaected 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 prole 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 Prole 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 prolefrom 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 prole of a hard-sphere uid at a hard-wall interface.This density prole can then be compared with the density prole from theCGMD simulation, thus to evaluate the eects 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 prole 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 Prole 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 eectivereference 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 Prole at Bioelectronic Interfacewhere b is the radius of particle, after dening 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 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