University of Southern California Dissertations and Theses

Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations

(USC Thesis Other)

Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content SHOCK-INDUCED PORATION, CHOLESTEROL FLIP-FLOP AND SMALL INTERFERING RNA TRANSFECTION IN A PHOSPHOLIPID MEMBRANE: MULTIMILLION ATOM, MICROSECOND MOLECULAR DYNAMICS SIMULATIONS By Amit Choubey _____________________________________________________________________ A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PHYSICS) May 2014 Copyright 2014 Amit Choubey ii Dedicated to my family iii Acknowledgments I would like to express sincere gratitude towards my advisor Prof. Rajiv Kalia. He has been very supportive and a constant source of encouragement for past five years. I am very thankful to Prof. Aiichiro Nakano, he has been a great source of motivation and has helped me wholeheartedly in every possible way. I also thank Prof. Priya Vashishta for always providing me with sound advice and keen insights. I would also like to thank my committee members Prof. Noah Malmstadt and Prof. Moh El-Naggar. I am very grateful to the USC Provost’s office for providing me fellowship for past five years. I must thank Dr. Ken-ichi Nomura for helping me in the research projects. I would also like to thank my colleagues Dr. Moe Vedadi, Van Ngo, Dr. Adarsh Shekhar, Ying Li, Pankaj Rajak, Dr. Zaoshi Yuan, Dr. Manaschai Kunaseth, Dr. Richard Seymour, and Dr. Weiwei Mou for all their help and support during my graduate studies. I take this opportunity to thank my friends Dr. Siddhartha Santra, Dr. Saptarshi Basu, Rahul Gupta, Avik Halder, Prashant Kumar and Mayukh Banik. They have been very motivating and always supported me through tough times. Finally, I would like to express my very sincere gratitude towards my parents Mr. Sushil Kumar Choubey and Mrs. Mansha Choubey, and my sister Mrs. Sweta Mishra for their unconditional love and support. I owe all my achievements to them and will be forever indebted to them. iv Table of Contents Acknowledgments: page iii Abstract: page 1 Chapter 1: Introduction: page 3 Chapter 2: Molecular Dynamics: page 20 Chapter 3: Structure and Dynamics of Shock-Induced Nanobubble Collapse in Water: page 35 Chapter 4: Poration of Lipid Bilayers by Shock-Induced Nanobubble Collapse: page 50 Chapter 5: Cholesterol Translocation in a Phospholipid Membrane: page 65 Chapter 6: Small Interfering RNA Induces Liquid-to-Ripple Phase Transformation in a Phospholipid Membrane: page 90 Chapter 7: Conclusions: page 109 References: page 112 1 Abstract Biological cell membranes provide mechanical stability to cells and understanding their structure, dynamics and mechanics are important biophysics problems. Experiments coupled with computational methods such as molecular dynamics (MD) have provided insight into the physics of membranes. We use long-time and large- scale MD simulations to study the structure, dynamics and mechanical behavior of membranes. We investigate shock-induced collapse of nanobubbles in water using MD simulations based on a reactive force field. We observe a focused jet at the onset of bubble shrinkage and a secondary shock wave upon bubble collapse. The jet length scales linearly with the nanobubble radius, as observed in experiments on micron-to-millimeter size bubbles. Shock induces dramatic structural changes, including an ice-VII-like structural motif at a particle velocity of 1 km/s. The incipient ice VII formation and the calculated Hugoniot curve are in good agreement with experimental results. We also investigate molecular mechanisms of poration in lipid bilayers due to shock-induced collapse of nanobubbles. Our multimillion-atom MD simulations reveal that the jet impact generates shear flow of water on bilayer leaflets and pressure gradients across them. This transiently enhances the bilayer permeability by creating nanopores through which water molecules translocate rapidly across the bilayer. Effects of nanobubble size and temperature on the porosity of lipid bilayers are examined. The second research project focuses on cholesterol (CHOL) dynamics in phospholipid bilayers. Several experimental and computational studies have been performed on lipid bilayers consisting of dipalmitoylphosphatidylcholine (DPPC) and 2 CHOL molecules. CHOL interleaflet transport (flip-flop) plays an important role in interleaflet coupling and determining CHOL flip-flop rate has been elusive. Various studies report that the rate ranges between milliseconds to seconds. We calculate CHOL flip-flop rates by performing a 15 µs all-atom MD simulation of a DPPC–CHOL bilayer. We find that the CHOL flip-flop rates are on the sub microsecond timescale. These results are verified by performing various independent parallel replica (PR) simulations. Our PR simulations provide significant boost in sampling of the flip-flop events. We observe that the CHOL flip-flop can induce membrane order, regulate membrane- bending energy, and facilitate membrane relaxation. The rapid flip-flop rates reported here have important implications for the role of CHOL in mechanical properties of cell membranes, formation of domains, and maintaining CHOL concentration asymmetry in plasma membrane. Our PR approach can reach submillisecond time scales and bridge the gap between MD simulations and Nuclear Magnetic Resonance (NMR) experiments on CHOL flip-flop dynamics in membranes. The last project deals with transfection barriers encountered by a bare small interfering RNA (siRNA) in a phospholipid bilayer. SiRNA molecules play a pivotal role in therapeutic applications. A key limitation to the widespread implementation of siRNA- based therapeutics is the difficulty of delivering siRNA-based drugs to cells. We have examined structural and mechanical barriers to siRNA passage across a phospholipid bilayer using all-atom MD simulations. We find that the electrostatic interaction between the anionic siRNA and head groups of phospholipid molecules induces a phase transformation from the liquid crystalline to ripple phase. Steered MD simulations reveal that the siRNA transfection through the ripple phase requires a force of ~ 1.5 nN. 3 1 Introduction 1.1 Cells Cells are the building blocks of life; they are basic structural, functional and biological units of living organisms that can replicate independently. The dimensions of cells range from 1 to 100 micrometers. Organisms can be unicellular (bacteria) or multicellular (e.g. humans contain 100 trillion cells). The cells are classified into prokaryotes and eukaryotes. Prokaryotic cells are usually single-celled organism devoid of a nucleus, and are 1 - 5 micrometer in size. Eukaryotic cells are much more complex, bigger in size (~10 to 100 micrometers) and contain a nucleus. Eukaryotic cells can be either single-celled or part of multicellular organisms. All prokaryotic and eukaryotic cells have membranes, which surround the cells and regulate permeability. Cell membranes will be described in the next section. Inside the membrane, a gel-like cytoplasm exists. Most cellular activities such as cell division, glycolysis occur via the organelles, which are suspended in the cytoplasm. Some major organelles are the mitochondria, the endoplasmic reticulum, the Golgi apparatus, vacuoles and lysosomes. The mitochondria are responsible for generating cell’s supply of adenosine triphosphate (ATP). The endoplasmic reticulum has the major function of producing proteins, lipids and other vital molecules, and it also participates in folding and transporting proteins. The Golgi apparatus modifies, sorts and packages macromolecules for cell secretion or use within a cell. Vacuoles and lysosomes play an important role in the cell’s waste disposal system. Most cells also possess DNA and RNA. DNA is the genetic material for the development and functioning of organisms. RNA is used for 4 genetic information transport and generation of proteins. The nucleus is the cell compartment, which stores all the genetic material. Figure 1.1 shows the structure of an animal cell which has a nucleus, cytoplasm and several organelles. Figure 1.1: A schematic diagram of an animal cell accompanied by electron micrographs of its organelles. Adapted from Biochemistry, 2 nd edition, Donald Voet and Judith G. Voet. 1.2 Lipid Bilayers Cell membranes are responsible for separating the contents of a cell from the outside world. The complexity of cell membranes is shown schematically in figure 1.2. The cell membrane has lipids, cholesterol, proteins, carbohydrates and sugars. Apart from providing mechanical support to cells they actively participate in several transport mechanisms such as endocytosis, exocytosis and transmembrane protein channels. 5 Figure 1.2: A schematic diagram of a plasma membrane. Integral proteins (orange) are embedded in a bilayer composed of phospholipids (blue spheres with two wiggly tails; shown, for clarity, in much more greater proportion than they have in biological membranes) and cholesterol (yellow). The carbohydrates components of glycoproteins (yellow beaded chains) and glycolipids (green beaded chains) occur only on the external face of the membrane. Adapted from Biochemistry, 2 nd edition, Donald Voet and Judith G. Voet. The lipid bilayer is a major component of a cell membrane. The lipid bilayer consists of two leaflets of lipid molecules arranged such that the core of the bilayer is hydrophobic and the outer boundaries are hydrophilic. In this thesis our focus will be on phospholipid bilayer consisting of Dipalmitoylphosphatidylcholine (DPPC) molecules. The structure of a DPPC molecule is shown in figure 1.3(a). The molecule has a polar zwitterionic head group and a nonpolar hydrocarbon region. Cholesterol (CHOL) molecules are also commonly found in cell membranes. They play an important role in regulating the membrane rigidity and intracellular transport. In this thesis we also model our membranes with both DPPC and CHOL molecules. Figure 1.3(b) shows the structure of a CHOL molecule. It has a hydroxyl hydrophilic head group, a rigid carbon ring and a floppy hydrocarbon tail. 6 Figure 1.3: (a) Structure of a DPPC molecule. The hydrophilic dipole group (NC 3 + -PO 4 - ) is highlighted in blue and grey. The DPPC also has two hydrocarbon tails, which are hydrophobic. (b) Structure of a CHOL molecule. The cholesterol molecule has a hydrophilic hydroxyl (OH) head group, a rigid carbon ring and a hydrocarbon tail group. Lipids are amphiphilic molecules. In the presence of water, lipids self assemble in such a way as to minimize the contact between hydrophobic tails and water and maximize the contact between the hydrophilic head groups and water. Depending on the relative concentration of lipids and water, several structures – micelles, inverted micelles, bilayer and vesicles – can form; see figure 1.4. In this thesis we will only focus on lipid bilayers (figure 1.4(c)) composed solely of DPPC molecules or a mixture of DPPC and CHOL molecules. The DPPC bilayer exists in crystalline (subgel), gel, liquid ordered and liquid disordered phases depending on the temperature of the system; see fig. 1.5. The crystalline phase is an orthorhombic lattice [1]. In the gel phase the lipids are hexagonally packed and the tails show conformational ordering. In the liquid ordered phase the hexagonal structure is lost but the tails still retain significant conformational order. In the liquid disordered phase the lipid tails are conformationally disordered. 7 Figure 1.4: Cross sections of some of the structures formed by lipids in water: (a) micelle; (b) inverted micelles; (c) bilayer; (d) bilayer vesicle. Adapted from http://dspace.jorum.ac.uk/xmlui/bitstream/handle/10949/1024/Items/T356_3_section2.ht ml Figure 1.5 also shows that the DPPC membrane forms a ripple phase when the temperature is decreased from 40° C to 32° C. X-ray studies indicate that the ripple phase in a DPPC bilayer consists of major (M) and minor (m) regions separated by kinks [2-4]; see figure 1.6. Experimental values for thicknesses of major and minor domains are consistent with the thickness of the gel and liquid crystalline phases, respectively [5]. Fourier transform infrared and 13 C-NMR measurements show a high degree of tail stretching of DPPC molecules in the ripple phase [6, 7]. 8 Figure 1.5: Molecular volume (open circles) and heat capacity (solid line) vs. temperature for DPPC bilayers in excess water. The figure is adapted from ref. [8]. Figure 1.6: Ripple phase structure. The ripple vector is in the horizontal direction; the vertical is the stacking direction. This figure is adapted from ref. [5]. 9 Figure 1.7: Temperature-composition phase diagram of a DPPC-CHOL bilayer. This figure is adapted from ref. [9]. In a bilayer of DPPC and CHOL molecules, the ring structure of CHOL molecules orients perpendicular to the membrane. CHOL ring lies around the DPPC tail region and arrests the movement of tails. Thus the presence of CHOL molecules enhances ordering of the DPPC tails and causes condensation by decreasing the area per lipid. CHOL molecules also increase the mechanical stability and bending rigidity of the membrane (see chapter 5). 10 Figure 1.7 shows the phase diagram of a DPPC-CHOL bilayer [9]. The DPPC- CHOL bilayer exists in three phases: solid-ordered (s o ) or gel, liquid-ordered (l o ) and liquid-disordered (l d ). The figure shows that the phases depend not only on the temperature but also on CHOL concentration in the bilayer. At low CHOL concentrations, the membrane exists mainly in the s o or l d phases. When CHOL concentration is increased, both s o and l d transform to the l o phase where the lipids have no long-range translational order but their chains have higher order and the area per lipid is lower than in the l d phase. In reality cell membranes exhibit highly complex phase behavior because of the presence of a large number of different lipid species. It is proposed that the lipids in cell membranes may form rafts, which can move within the fluid bilayer [10]. When these rafts move inside the cell, they selectively attach to or detach from proteins. Lipid rafts are small (< 200 nm), rich in cholesterol and sphingolipids, and may play a role in intracellular signaling. Domains formed by tertiary mixtures composed of saturated lipids, unsaturated lipids and cholesterol molecules resemble rafts in cell. In vivo observation of rafts is difficult but raft-like mixtures in model membranes have been visualized experimentally [11-15]. Risselada et al. [16] have used coarse-grained molecular dynamics (MD) simulations to study such domains at the molecular level. They show that a tertiary mixture of saturated phosphatidylcholine (PC)/unsaturated PC/cholesterol spontaneously phase separates into a liquid-ordered and a liquid-disordered phase with structural and dynamic properties closely matching experimental data. These simulations also reveal features of domains and the existence of a small surface tension between the domains and 11 a short-ranged line repulsion between the domain boundaries. Such simulations help us understand lipid rafts and cell membrane processes such as self-assembly of functional protein complexes. 1.3 Molecular motions in membranes The timescale of molecular motions in lipid bilayers ranges between picoseconds and hours. Some of the most important dynamical processes in membranes are: (1) Molecular flip-flop from one membrane leaflet to another; (2) lateral diffusion of molecules in the membrane plane; (3) rotational motions; (4) hydrocarbon chain dynamics; and (5) molecular movements perpendicular to the bilayer plane which cause membrane undulations. The timescale for lipid flip-flop dynamics is reported to be 0.5 to 50 hours [17]. Various studies show that the CHOL flip-flop dynamics in membranes is significantly faster than the flip-flop dynamics of lipid molecules. Experimental studies indicate that the CHOL flip-flop timescale is between microseconds to seconds [18, 19]. The reason for faster CHOL flip-flop is that these molecules are smaller and can move more easily than lipid molecules. The discrepancy in estimating CHOL flip-flop rate persists in literature because external perturbations in experimental techniques tend to influence flip- flop rates [18]. Bruckner et al.[20] demonstrated that CHOL stabilizes changes in membrane shape by flipping between the bilayer leaflets. The movement of CHOL molecules decreases the surface area of relatively compressed leaflet while increasing the surface area of the relatively expanded leaflet. This mechanism is illustrated in fig 1.8. They report that the CHOL interleaflet diffusion time is of the order of milliseconds. They also 12 suggest that the evolutionary reason for the presence of CHOL molecules in membranes is that these molecules relax stresses in membranes. Figure 1.8: Deformation of a pure phospholipid membrane compresses the inner leaflet and expands the outer leaflet. CHOL molecules rapidly flip from the compressed to the expanded leaflet, resulting in stress relaxation and dissipation of input energy [20]. The lateral diffusion coefficient (D) of lipid molecules in a bilayer is studied by fluorescence recovery after photobleaching (FRAP), single particle tracking (SPT) and NMR techniques [21-24]. In the disordered fluid phase, D varies between 10 −7 cm 2 /s to 10 −8 cm 2 /s [22, 24]; and in the gel phase D lies between 10 −11 cm 2 /s and 10 −16 cm 2 /s. Rotational motion of lipids is measured by NMR and fluorescence techniques. The time scale for rotational motion ranges between 10 picoseconds and 1 nanosecond. NMR has also been used to determine the dynamics of hydrocarbon chains. The size and dynamics of membrane undulation have been correlated to the elastic properties. The length and time scales of undulations range from 10 nanometers to 10 micrometer and sub nanosecond to several milliseconds, respectively [25]. 13 1.4 RNA Interference by Small Interfering RNA In 2006, Fire and Mello were awarded the Nobel Prize "for their discovery of RNA interference - gene silencing by double-stranded RNA". RNA interference (RNAi) is a pathway whereby a double-stranded small interfering RNA (siRNA) targets and destroys complementary mRNA in eukaryotic cells [26]. Nowadays, synthetic siRNA can be directly introduced into a cell. Once in the cytoplasm, siRNA binds with proteins to form an RNA-induced silencing complex (RISC) in which a protein called Argonaute 2 (AGO2) unwinds the double-stranded siRNA. Subsequently, the sense (or passenger) strand of siRNA is cleaved and the RISC with the antisense (or guide) strand of the siRNA is thus activated to seek out and cleave an mRNA that is complementary to the antisense strand. A single activated RISC can destroy many mRNAs, one after another, thereby silencing genes. Figure 1.9 shows the RNAi mechanism schematically. Potentially effective strategies based on siRNA treatment have been investigated for many types of cancers [27]. There are currently over a dozen RNA-based anti-cancer drugs in clinical trials and many more in early stages of development [27]. A key limitation to the widespread implementation of siRNA techniques is the difficulty of delivering siRNA-based drugs to cells. A double-stranded siRNA has 21-25 nucleotides on each strand and 2 nucleotide overhangs at the 3'- ends to enhance its stability, facilitate its inclusion into RISC, and reduce the immune responses. Because of its large molecular weight (~ 13 kDa) and anionic charge, siRNA cannot passively penetrate the cell membrane. Chapter 6 describes the molecular dynamics study of mechanical barriers to siRNA transfection in a DPPC membrane. 14 Figure 1.9: Long double-stranded RNA (dsRNA) is introduced into the cytoplasm, where it is cleaved into small interfering RNA (siRNA) by the enzyme Dicer. Alternatively, siRNA can be introduced directly into the cell. The siRNA is then incorporated into the RNA-induced silencing complex (RISC), resulting in the cleavage of the sense strand of RNA by argonaute 2 (AGO2). The activated RISC–siRNA complex seeks out, binds to and degrades complementary mRNA, which leads to the silencing of the target gene. The activated RISC–siRNA complex can then be recycled for the destruction of identical mRNA targets. This figure is adapted from ref. [27]. 15 Parisien et al.[28] have developed a machine-learning method for searching protein partners for a specific RNA. The interactions between RNA proteins and a specific RNA are very important for RNA functioning. The method depends on scoring protein’s binding potential for an RNA structure. This computational approach provides a powerful complement to experiments and MD simulations in discovering new RNA proteins. 1.5 Drug Delivery Methods Our focus is to explore methods and barriers to drug delivery into individual cells. The cell membrane acts as a major barrier to delivery of drugs into cells. We have studied how drug molecules (siRNA) interact with DPPC membranes at the atomistic level. Two drug delivery methods, electroporation and sonoporation, are described. We also briefly discuss liposome and nanoparticle based gene therapy methods here. Electroporation refers to reversible permeabilization in membranes induced by an external electric field across the cell membrane. The transient permeabilized state can be used to drive molecules such as ions, drugs, dyes, RNA and DNA into the cells. The membrane is subjected to few hundred millivolt of transmembrane potential to cause electroporation. Tarek and others have performed MD simulations of electroporation and provided detailed picture of the electroporation phenomenon [29, 30]. Simulations reveal that prior to electroporation water molecules organize linearly like wires and penetrate the hydrophobic core of the bilayer. The lipid head groups stabilize these pores by migrating from the membrane-water interface to the middle of the bilayer. Switching off the transmembrane potential allows for complete resealing and reconstruction of the 16 bilayer. Using MD studies Breton has shown that an electric pulse facilitates delivery of charged siRNA to the cytoplasm [29]. Sonoporation technique refers to the rupture of cell membranes by acoustic waves. Acoustic waves in tandem with bubble collapse drastically enhance sonoporation [31]. Ohl et al. [31] have performed sonoporation experiments in the presence of cavitation bubbles. They demonstrate that cavitation bubbles and not the shock wave itself causes drug delivery by the following mechanism: the bubble collapse close to the substrate causes a jetting flow towards the substrates and the impact spreads the jet radially along the substrate. The resulting radial flow exerts a strong shear stress on the cell and causes cell detachment. Cells along the border of the detached area are permanently damaged whereas cells further away show repairable poration. We have performed MD simulations to study the effect of shock waves in the presence of nanobubbles on a DPPC lipid bilayer. The simulations reveal atomistic mechanism underlying the formation of a nanojet and poration of a DPPC bilayer. These simulations are described in details in chapter 4. Liposomes (see fig 1.4) made of neutral and cationic lipids have garnered a lot of interest as vectors for gene therapy. Liposomes are easy to prepare, they are biocompatible and can be loaded with a variety of drugs, DNA and diagnostic agents [32]. Also, liposome based vectors do not have severe restrictions on the payload size and have no proteins to trigger the immune system. The challenge is to achieve high efficiency in delivery using liposomes. The transfection efficiency of the lipid-DNA complexes depends on the lipid/DNA charge ratio [33]. The efficiency increases as the charge ratio increases up to a 17 saturation value. The transfection efficiency also depends on the membrane charge density per unit area. Safinya et al. have shown that the membrane charge density is a universal parameter for lipid-DNA complexes and the optimal charge density for transfection efficiency is 17.0±0.1 ×10 –3 (e/Å 2 ) [34]. Nanoparticles (NPs) containing drugs have also been studied for targeted drug delivery [35]. The size and surface characteristics of nanoparticles can be easily manipulated to achieve targeted delivery. For example, site-specific (selective) targeting can be achieved by attaching ligands to the NP surface or by using a magnetic field to guide magnetic nanoparticles. Also nanoparticles can control the release of a drug which can improve therapeutic efficacy and reduce side effects. Nanoparticles based on polymers, solid-metal, albumin, and gelatin have been tested for drug delivery. These delivery systems can be designed to have the drugs absorbed or conjugated onto the nanoparticle surface, encapsulated inside the polymer/lipid nanoparticle or dissolved within the nanoparticle matrix. The efficiency of these drug delivery systems depends on the size, shape and surface characteristics of the nanoparticles. Nanoparticle carrying drugs can be targeted to release drug in specific pH conditions [36-39]. This allows the drugs to be delivered directly to a tumor site. Tumors are more acidic than normal human cells. Tumors have pH around 6.8 whereas normal tissue has a pH of 7.4. Thus vectors that deliver drugs by degrading themselves in acidic tumor environments but do not release drugs in neutral or basic environments are good- targeted drug delivery systems. Also tumors are generally at higher temperature than the rest of the body. Thus drug delivery systems can also have temperature based specificity [37]. 18 Quantum mechanical simulation approaches have been developed to include the effect of solution pH in computer simulations, but these methods are computationally very expensive. Constant pH MD simulations [40] can address the limitations of QM methods. In these simulations an explicit solvent is used with the λ-dynamics. The Hamiltonian of the system is given by H(λ)=(1−λ) ˆ H 0 +λ ˆ H 1 +H env + m 2  λ 2 +U * (λ), (1.1) where λ is the titration coordinate. 𝐻 ! and 𝐻 ! are λ-dependent reactant and product Hamiltonians, 𝐻 !"# is a λ-independent term, the fourth term is the kinetic energy associated with λ, and the last term is a λ-dependent potential which biases the system to physically relevant protonated state (λ = 0) and deprotonated (λ = 1) states. The dynamics of the λ variable is governed by the potential energy function, V(λ)=(1−λ) ˆ V 0 +λ ˆ V 1 +V env +U * (λ) . (1.2) In this scheme the average value of λ is constrained to be 0 or 1. This means that the system spends most of the time in the protonation or deprotonation state. The transition time between the two states is short and the transition frequency can be controlled in this scheme. In chapter 2 we describe the MD simulation methodology including the interatomic force fields, and the integration algorithm for equations of motion of atoms. The chapter also covers isobaric and isothermal simulations. In chapter 3 the MD simulation of shock wave propagation and nanobubble collapse in water are described [41]. In chapter 4 we describe 30 million atoms MD simulation of poration of lipid bilayers caused by the collapse of a nanobubble [42]. In chapter 5 we discuss 19 microsecond MD simulations of cholesterol interleaflet transport (flip-flop) in lipid membranes [43]. These simulations shed light on the stress relaxation mechanism of membranes. In chapter 6 we describe MD simulation of how a siRNA induces phase transformation from liquid crystalline to ripples phase in a phospholipid bilayer. Concluding remarks are presented in chapter 7. 20 2 Molecular Dynamics 2.1 Overview Molecular dynamics (MD) is a simulation technique that captures the motion of atoms over a wide range of length (~ 1 Å → 100 nm) and time scales (~ 10 -15 → 10 -6 s). The atoms interact with each other and Newton’s equations of motion are numerically solved to obtain phase-space trajectories of the atoms. MD is widely used in chemical physics, material science and biophysics. The major focus of this work is the study of biomembranes using the MD simulation technique. The interatomic force fields and the MD methodology are described in sections 2.2 and 2.3, respectively. In section 2.4 we will describe how isothermal and isobaric MD simulations are performed. Finally, in section 2.5 we will discuss how physical quantities are calculated from phase-space trajectories. MD simulations are usually carried out in the microcanonical ensemble where the number of particles (N), volume (V) and energy (E) are fixed. MD calculations are also performed in isothermal (fixed N, V and temperature T) or isobaric (constant N, pressure P and T) conditions to mirror the experimental conditions. In a MD simulation, the number of atoms and size of the MD cell are specified according to the density of the system. To perform a MD simulation, it is also necessary to specify initial positions and velocities of all the atoms in the system. Subsequently, atomic coordinates and velocities are obtained numerically by integrating the equations of motion. Many integration algorithms exist. We have used the velocity-Verlet algorithm. The schematic in figure 2.1 gives an overview of the MD technique. 21 Figure 2.1: Gives an overview of the MD technique. We specify the initial conditions, compute forces, and update the coordinates and velocities. 22 2.2 Interatomic Force Fields In all our simulations we use Dipalmitoylphosphatidylcholine (DPPC) molecules to model the lipid bilayer. The force field for DPPC molecules is taken from CHARMM 36 or the united atom force field. For water we use either TIP3P (in CHARMM 36) or SPC (when using united atom force field) model. The potential functions for DPPC and water molecules consist of bonded and non-bonded interactions as described below. 2.2.1 Bonded Interactions These interatomic interactions include bond stretching (2-body), bond angles (3- body), and dihedral angles (4-body) interactions between atoms. The bond stretching contribution to potential energy function is given by: V b = 1 2 k ij b (r ij −b ij ) 2 , (2.1) where b ij is the equilibrium bond length and 𝑘 !" ! is a force constant. The bond-angle term between a triplet of atoms i - j - k is also represented by a harmonic potential: V a (θ ijk )= 1 2 k ijk θ (θ ijk −θ ijk 0 ) 2 . (2.2) The dihedral interaction is usually represented by: V d (ϕ ijkl )=k ϕ (1+cos(nϕ ijkl −ϕ s )) . (2.3) Here ϕ ijkl is the angle between planes formed by atoms (i,j,k) and (j,k,l), and n is an integer related to the multiplicity. The parameters 𝑘 !" ! ,𝑏 !" ,𝑘 !"# ! ,𝜃 !"# ! ,𝑘 ! ,𝑛,𝜑 ! depend on atom types and they are given in the force field tables. The bonded terms (2.1) – (2.3) are the most commonly used functional forms in the bioforce fields. Other functional forms for force-field exist 23 and have been applied in GROMACS [44]; see GROMACS manual for details. 2.2.2 Non-bonded Interactions All non-bonded interactions are pairwise additive and central symmetric i.e. V(  r 1 ,...  r N )= V ij i<j ∑ (  r ij ) , (2.4)  r i is the coordinate of the i th atom and  r ij is the vector from atom i to atom j. The force  F i is given by  F i =− dV ij (r ij ) dr ij j ∑  r ij r ij =−  F j . (2.5) The non-bonded interactions contain steric, dispersion, and Coulomb terms. The steric and dispersion terms are combined in the Lennard-Jones interaction (V LJ ): V LJ (r ij )= C ij (12) r ij 12 − C ij (6) r ij 6 . (2.6) The parameters 𝐶 !" (!") and 𝐶 !" (!) depend on pairs of atom types. These parameters are given in CHARMM 36 or united atom force fields. Since the Lenard-Jones interaction decays rapidly with distance, it is a common practice to use a cut-off distance r c in the calculation and shift the potential energy to avoid discontinuities at r c : V SF (r ij )= V(r ij )−V(r c )− dV(r ij ) dr ij " # $ $ % & ' ' r ij =r c (r ij −r c ) r ij ≤r c 0 r ij >r c ) * + + , + + . (2.7) This form of the shifted potential ensures that forces vanish smoothly at the cutoff distance. In the Coulomb interaction, 24 V c = q i q j 4πε 0 ε r r ij i<j ∑ , (2.8) the charges (q i ) are specified in the force field parameters. The Coloumb interaction decays slowly and the summation in Eq. (2.8) is only conditionally convergent. V c is usually calculated by the particle mesh Ewald (PME) method or by an approximate approach such as the reaction force field. In the reaction field method, the Coulomb interaction for homogeneous systems is modified by assuming a constant dielectric environment with a dielectric constant of ε rf beyond the cut-off r c . The interaction is given by: V crf = q i q j 4πε 0 ε r r ij 1+ ε rf −ε r 2ε rf +ε r r ij 3 r c 3 " # $ $ % & ' ' − q i q j 4πε 0 ε r r c 3ε rf 2ε rf +ε r . (2.9) In the PME method the total electrostatic energy of N particles and their periodic images is given by V = 1 8πε 0 q i q j r ij,  n j N ∑ i N ∑ n z * ∑ n y ∑ n x ∑ , (2.10) (n x ,n y ,n z ) =  n is the box index vector, and the star indicates that i = j terms are omitted when (n x ,n y ,n z ) = (0,0,0). The distance r ij,n is the real distance between the charges and not the minimum-image. The Ewald summation is used to convert the slowly-converging sum into two rapidly converging terms plus a constant: V =V dir +V rec +V 0 , (2.11) V dir = 1 8πε 0 q i q j erfc(βr ij,  n ) r ij,  n n z * ∑ n y ∑ n x ∑ i,j N ∑ , (2.12) V rec = 1 8π 2 ε 0 V q i q j i,j N ∑ exp(−(π  m β) 2 +2πi  m⋅(  r i −  r j ))  m 2 m z * ∑ m y ∑ m x ∑ , (2.13) 25 V 0 =− β 4π 1.5 ε 0 q i 2 i N ∑ , (2.14) where the parameter β determines the relative weight of the direct and reciprocal sums and  m = (m x ,m y ,m z). In this way one can use a short cut-off (of the order of 1 nm) in the direct space and also a short cut-off in the reciprocal space (e.g. 10 wave vectors in each direction) to compute (2.12) and (2.13). Unfortunately, the computational cost of the reciprocal part of the sum increases as 𝑁 ! ! and it is therefore not suitable for large systems. Tom Darden [45] proposed the Particle-mesh Ewald to improve the performance of the reciprocal sum. Instead of directly summing wave vectors, the charges are assigned to a grid. The implementation uses cardinal B-spline interpolation, which is referred to as smooth PME (SPME). The grid is then Fourier transformed with a 3D FFT algorithm and the reciprocal energy term is obtained by a single sum over the grid in k-space. The potential at the grid points is calculated by inverse transformation, and by using the interpolation factors for forces on atoms. The PME algorithm scales as 𝑁𝑙𝑜𝑔(𝑁), and is substantially faster than the Ewald summation on system containing N > 10 5 atoms. The simulations are only as accurate as the parameters used in interaction potentials. Thus parameterization of force fields for biomolecules is very crucial and is still an active field of research [46-48]. In our simulations we primarily use two types of force fields: CHARMM or Berger force field. In CHARMM[47], all atoms are described explicitly. In the Berger force field [49], hydrogen atoms are described together with aliphatic carbons. This kind of description is called a united atom model. Even more coarse-grained models exist where few atoms are lumped together and referred to as a bead. For example, in the MARTINI model approximately four heavy atoms are 26 described by a single bead [50]. In all force fields the physical properties (e.g. density, heat of vaporization, etc.) are first calculated from ab initio methods and then refined by comparing the calculated values with the available experimental data [51]. The parameters are adjusted such that the force field reproduces the experimental data. Ab initio or spectroscopy methods provide the parameters for bond stretching and bond angle terms. The partial atomic charges and the parameters in the Lennard-Jones interactions are harder to determine. Partial charges are sometimes obtained from ab initio calculations and sometimes chosen empirically. Once the partial charges, bond stretching and bond angle terms are determined, experimental data are used to fit the parameters in the Lennard-Jones and dihedral interactions. 2.3 MD Methodology 2.3.1 Equations of motion In MD the following equations of motion are solved numerically: d 2  r i dt 2 =  F i m i ;  F i =− ∂V ∂  r i ; v i = d  r i dt . (2.15) Here  r i ,  v i ,  F i are the coordinate, velocity and force vector of the i th atom. In our simulations the bond vibrations are replaced by holonomic constraints, which allows us to use a bigger time step (1 fs) in the integration of equations of motion. We use the Linear Constraint Solver (LINCS) method to maintain the constraints during the simulation. In LINCS the Newton’s equations of motion are solved under the bond constraints using the Lagrange multiplier method. In the matrix form, the equations of motion become: 27 d 2  r dt 2 =  M −1  f =−  M −1 ∂V ∂ r , (2.16) where for a N-particle system  r and  f are 3N position and force vectors, respectively.  M is a 3N×3N diagonal matrix containing the masses of the particles. Let us denote the time-independent holonomic bond-length constraints by: g i ( r)= 0 i=1 ,...,K , (2.17) where the system has K such constraints. Implementing the constraints in the equation of motion using Lagrange multipliers we get: −  M d 2  r dt 2 = ∂ ∂ r (V −  λ⋅  g) , (2.18) where  λ is a vector of Lagrange parameters. The above equation can be simplified into the following form: −  M d 2  r dt 2 +  B T  λ+  f =0 , (2.19) where  B is a K×3N matrix whose elements are B hi = ∂g h ∂r i . (2.20) Since the constraints are time-independent, the first and second derivatives of constraints with respect to time are equal to zero. Thus the constrained equations of motion become: d 2  r dt 2 =(  I −  T  B)  M −1  f −  T d  B dt d r dt , (2.21) where  T =  M −1  B T (  B  M −1  B T ) −1 . The integration algorithm and implementation of these equations of motion are described in ref. [52]. 28 2.3.2 Velocity-Verlet Integration Algorithm The velocity-Verlet algorithm is used to integrate the equation of motion. The algorithm has three steps:  v(t+ Δt 2 )=  v(t)+ Δt 2m  F(t), (2.22)  r(t+Δt)=  r(t)+Δt  v(t+ Δt 2 ) , (2.23)  v(t+Δt)=  v(t+ Δt 2 )+ Δt 2m  F(t+Δt) . (2.24) The velocity-Verlet algorithm is chosen because it is time-reversible and preserves phase- space volume. The velocity-Verlet algorithm is stable in that it prevents any drift in the total energy. 2.3.3 Boundary Conditions To avoid finite-size effects, we use periodic boundary conditions. The calculations are done by imagining that the MD box is replicated in all the three directions, so that the atoms close to the boundary of the MD box interact with atoms in the replicated image box. Also the box is wrapped in the sense that any atom leaving the box renters the box from the other side. Periodic boundary conditions are used together with the minimum-image convention in which an atom i interacts with another atom j or an image of the atom depending on whether the j th atom or its image is closer to the i th atom. The minimum- image convention introduces a constraint in the system, i.e., the cut-off radius for non- bonded interactions cannot exceed half the shortest box length. The minimum image convention is implemented as follows (shown for x-direction only, similarly for other two 29 directions): x ij = x ij + L x 2 , if x ij ≤− L x 2 , (2.25) x ij = x ij − L x 2 , if x ij ≥ L x 2 . (2.26) The non-bonded interactions are calculated between atomic pairs i and j for which r ij < r c . A linked list is constructed to compute pairwise non-bonded interactions. This list is updated every (~ 10) time steps. From the linked list, a neighbor list is constructed to further expedite the calculation of non-bonded interactions. For details see GROMACS manual. 2.3.4 Parallel Implementation of MD Parallel computers consisting of thousands of CPU and GPU processors are commonly used to perform MD simulations. These multiple processors increase the simulation speed and also enable simulation of large systems (10 6 - 10 9 atoms). In a parallel MD scheme the computation load is shared among the processors. The parallel scheme involves interprocessor communication, which increases with an increase in the number of processors. At some point the communication overhead exceeds the computation cost. Thus for each simulation, there exists an optimal number of processors that maximizes the simulation speed and minimizes the communication overhead. The most common parallelization scheme is the spatial domain decomposition. In this scheme, a spatial domain is assigned to each processor and the processors are logically arranged according to the topology of the physical system. Each processor integrates the equations of motion for the particles residing in its domain. The details of 30 the domain decomposition can be found in the GROMACS 4 paper [44]. The parallel implementation involves Message Passing Interface (MPI), which is a library for interprocessor communication in a distributed-memory parallel programming environment. Our simulations are performed on 100 to 1000 processors and the speed of these simulations is about 20 – 100 ns a day. 2.4 Isothermal and Isobaric MD Simulations Berendsen’s algorithm is commonly used to control the temperature of the system. This algorithm corrects the temperature deviation in the system temperature using the following equation: dT dt = T 0 −T τ , (2.27) where T o is the reference temperature and T is the temperature of the system at time t. τ is the strength of the temperature coupling and usually taken to be ~ 0.5 ps. The integration time step is ~ 1 fs. The Berendsen thermostat does not generate a proper canonical ensemble. The errors in ensemble averages scale as 1/N and for very large systems the averages will not be significantly affected (except the distribution of the kinetic energy and fluctuation-dependent properties such as the heat capacity). In Berendsen algorithm the velocities of each atom is scaled every n TC steps by a factor λ: λ= 1+ n TC Δt τ T T 0 T −1 # $ % & ' ( ) * + , - . 1/2 . (2.28) The parameter τ T is close, but not exactly equal, to the time constant τ. In practice, the ratio τ/τ T ranges from 1 (gas) to 2 (harmonic solid) to 3 (water). The Berendsen scheme does not simulate the canonical ensemble. Nose and 31 Hoover introduced a temperature coupling scheme to simulate systems in the canonical ensemble. The Hamiltonian of the system is extended by introducing a thermal reservoir and a friction term. The frictional force is proportional to the product of each particle’s velocity and a friction parameter, ξ, which has a momentum (p ξ ). The equation of motion for p ξ is: dp ξ dt =(T−T 0 ) , (2.29) where T and T o are the current temperature and reference temperature, respectively. The equation of motion for each particle is: d 2  r i dt 2 =  F i m i − p ξ Q d  r i dt , (2.30) where Q is a coupling parameter. Here the conserved quantity is the extended Hamiltonian: H = p i 2 2m i i=1 N ∑ +V(r 1 ,....,r N )+ p ξ 2 2Q +N f KTξ , (2.31) where N f is the total number of degrees of freedom. We use the Parrinello-Rahman pressure coupling to perform simulations in the NPT ensemble [53]. Here, the MD box vectors represented by a matrix b (it is a diagonal matrix for our rectangular simulation box with the box lengths along the diagonal) obey the following equation of motion: d  b 2 dt 2 =V  W −1  b '−1 (  P−  P ref ), (2.32) where V is the volume of the box, and  W is a matrix related to the strength of pressure coupling. The matrices  P and  P ref are the current and reference pressures, respectively. 32 The equation of motion for each particle is: d 2  r i dt 2 =  F i m i −  M d  r i dt , (2.33)  M =  b −1  b d  b ' dt + d  b dt  b ' " # $ % & '  b '−1 . (2.34) See GROMACS manual for details. 2.5 Calculation of Physical Quantities The potential energy is calculated by summing the potential energy contributions from bonded and non-bonded interactions. The kinetic energy is calculated from: E kin = 1 2 m i v i 2 i=1 N ∑ , (2.35) where <…> represents time averaging. The temperature of the system is given by: T = 2E kin N df k , (2.36) where k is the Boltzmann’s constant and N df is the number of degrees of freedom, N df =3N−N c −N com . (2.37) Here N is the total number of atoms in the system, N c is the total number of constraints imposed on the system and N com = 3 because the center of mass velocities are the constants of motion and are usually set to zero. The radial distribution function (RDF) g(r) can be easily calculated using the following formula: g(r)= V N 2 δ(r−r ij ) j≠i ∑ i ∑ , (2.38) 33 where r ij is the distance between atoms i and j. RDF is calculated by generating a histogram with bins of length δr. The bin (k+1) stores the number of pairs with separation between kδr and (k+1)δr. <…> again represents time averaging. The self-diffusion coefficient (D) of atoms in d dimensions is calculated by using the following relation. D= 1 2d lim t→∞  r(t 0 +t)−  r(t 0 ) [ ] 2 t . (2.39) This relation holds only when the observation time (t) goes to infinity. <…> is an average over all the atoms and time origins t 0 . The pressure tensor  P is calculated from the Virial:  P= 1 V m i  v i ⊗  v i +  r ij ⊗  F ij i<j ∑ i N ∑ # $ % % & ' ( ( . (2.40) Here m i is the mass of the i th atom,  v i is its velocity, and and are the position vector and force between atoms i and j, respectively. V is the volume of the simulation box. For lipid bilayers, the stress π(z) is defined as: π(z)=(P xx +P yy ) 2−P zz , (2.41) where P xx , P yy , and P zz are the diagonal components of the pressure tensor  P . For a system consisting of particles with n-body potentials U n , the local pressure tensor  P(  r) is defined as [54]: P αβ (  r)= m i i ∑ v i α v i β δ(  r−  r i )− 1 n n ∑ ∇ jk α U n −∇ jl α U n ( ) k,l ∑ j ∑ dl β  ∫ δ  r−  l ( ) , (2.42) where the integral is along a contour from particle j l to particle j k , and <j> stands for summation over all n clusters in the system. <k,l> describes the summation over all pairs € r ij € F ij 34 of particles within a given cluster n; m i , v i , and  r i refer to the mass, velocity, and location of atom i, respectively; and α and β are Cartesian components. Using this definition of the local pressure tensor, we can calculate the pressure tensor for a volume element V as follows:  P V =  P(  r)d  r V ∫ . (2.43) Here the integral is over the element volume V. Thus the pressure tensor for volume V is P V αβ = 1 V m i i∈V ∑ v i α v i β + 1 nV n ∑ ∇ jk α U n −∇ jl α U n ( ) k,l ∑ j ∑ r j l j k β N f V  r j l + λ N  r j l j k % & ' ( ) * , (2.44) where f V  r ( ) =1, if  r∈V and zero otherwise. Each vector  r j l j k =  r j k −  r j l is divided into N parts and the contribution of a given part λ is added only if the contour goes through V. The contour used is the Irving-Kirkwood contour [54], which is a straight line from particle j l to particle j k . The surface tension of a lipid bilayer between –z 0 and z 0 is given by γ = π(z)dz −z 0 z 0 ∫ . (2.45) The bending modulus of the bilayer is the first moment of the stress profile c 0 κ = z−z 1 ( ) −z 0 z 0 ∫ π z ( ) dz , (2.46) where c 0 denotes the spontaneous curvature with respect to a pivotal plane z 1 . 35 3 Structure and Dynamics of Shock-Induced Nanobubble Collapse in Water 3.1 Introduction When a bubble interacts with a shock wave, it collapses because the surface tension cannot provide enough restoring force. Experimental studies on micron-size bubbles reveal that the impact of a pulsed shock on the proximal side of the bubble shrinks and accelerates it in the direction of the shock propagation. Furthermore, the liquid around the collapsing bubble forms a jet which creates a protrusion and secondary water hammer shock wave when it hits the distal side of the bubble and breaks up [55]. The behavior of a collapsing bubble depends on its location relative to objects embedded in the fluid as well as the shock wave characteristics (amplitude and pulse width) and the initial size of the bubble. Shock-induced bubble collapse phenomena have been studied in the context of a single bubble [56-58] as well as multiple bubbles [59] near a rigid boundary and in confined environments (e.g., microfluidic or lab-on-a-chip systems) [60, 61]. Jetting and secondary water hammer shocks can cause significant damage in materials. This problem is encountered in the disintegration of blades of ship propellers, pipelines and pump blades [62]. In medicine, however, collapsing bubbles have found useful applications such as extracorporeal shock wave lithotripsy [63] and targeted drug delivery [64]. In the so-called sonoporation approach, the collapse of microbubbles generates liquid jets and radial spreading flows [65] that can make the cell membrane 36 transiently permeable to molecular entry. This has potential applications in gene therapy and anticancer drug delivery [66, 67]. In this chapter, we discuss molecular dynamics (MD) simulations of shock- induced collapse of nanobubbles in water. The MD approach is well suited to study this problem because it can provide direct information about the structure and dynamics of nanobubble collapse over microscopic spatio-temporal scales. Contained in this microscopic information are some of the subtle but salient features of bubble collapse that experiments or continuum simulations may not be able to capture. 3.2 Simulation Setup The MD simulations described here are based on quantum-mechanically informed reactive force fields (ReaxFF), which can accurately describe bond breaking/formation and chemical reactions in the system [68]. In ReaxFF, charge transfer between atoms changes dynamically with changes in the atomic configuration and atomic charges are obtained by minimizing the electrostatic energy at every MD step [69]. The ReaxFF also includes three- and four-body interactions. The ReaxFF potential for water was trained against a wide range of quantum mechanical data, including binding energies for n-water clusters (n = 2-30), proton transfer reactions in H 3 O + -H 2 O and OH-H 2 O dimers, concerted proton transfer reactions in [H 2 O] n=2-6 clusters and dissociation energies for H-H, O-H, O-O and O=O bonds. The parameters were validated for bulk water by comparing the density, cohesive energy, self-diffusion coefficients and O-H, O-O and H-H radial distribution functions with experiments. 37 Figure 3.1: Schematic of a simulation cell. The grey plate is the momentum mirror. Shocks are applied in the x direction using a momentum mirror positioned at the left end of the box. At time t = 0 all atoms are given an initial particle velocity, u p , towards the momentum mirror and whenever an atom crosses the y-z plane of the mirror, the x- component of its velocity is reversed. Using a scalable fast reactive force field algorithm (F-ReaxFF) [70, 71], we have performed MD simulations of water subjected to planar shock with and without a nanobubble. Initial dimensions of the MD box are 38.5 nm, 18.6 nm, and 16.6 nm in the x, y, and z directions, respectively, and the system contains a million atoms. We equilibrate the system in the (N, V, T) ensemble i.e. constant number of atoms, constant volume, and constant temperature with periodic boundary condition (PBC) in all directions. The initial mass density and temperature are 0.98 g/cc and 300 K, respectively. After equilibration, we insert nanobubbles of different diameters (D = 6, 8, 38 and 10 nm) at the center of the MD cell by removing 90 % of water molecules from spheres of corresponding diameters. Figure 3.1 shows the simulation setup. These equilibrated systems are then subjected to planar shocks with particle velocities u p = 1.0, 2.5, 3.0, or 3.5 km/s. Shocks are applied in the x direction using a momentum mirror positioned at the left end of the box. At time t = 0 all atoms are given an initial particle velocity, u p , towards the momentum mirror and whenever an atom crosses the y-z plane of the mirror, the x-component of its velocity is reversed. In shock simulations we insert a 2-nm-thick vacuum layer at the end of the MD cell in the x direction (shock direction), turn off the thermostat coupling, and apply PBCs in lateral directions to minimize surface effects normal to the shock direction. The distance between a nanobubble and its nearest mirror image (≥ 6.6 nm) is much larger than the cut-off length (1 nm) of ReaxFF, which rules out any interaction between the bubble and its images. The integration time step in these MD simulations is 0.1 femtosecond. 39 Figure 3.2: Experimental (blue crosses) and MD (red circles) Hugoniot compression curves for particle velocities between 0.4 and 3 km/s. The inset shows the simulation cell and the momentum mirror (grey plate). 3.3 Results We first performed shock simulations without the nanobubble to validate the reactive force field of water. Figure 3.2 shows the setup for this simulation in the inset, and MD (red circles) and experimental (blue crosses) [72] results for the shock velocity, u s , as a function of the particle velocity, u p , i.e. the Hugoniot curve. Shock velocity is obtained from the difference in the shock-front boundaries at two time frames. In each frame, the abrupt change in the density identifies the location of the shock front. The 40 simulation results for the Hugoniot curve are in good agreement with the experimental results [72]. Figure 3.3: Oxygen-oxygen pair-correlation functions, g(r), for ambient and shocked water with and without a nanobubble as a function of distance at u p = 3 km/s. The bubble diameter is 10 nm. Shock produces significant structural changes in water. Figure 3.3 shows that the oxygen-oxygen radial distribution function (RDF) for water under planar shock is significantly different from that for ambient water. The maxima and minima in the RDF for shocked water (u p = 3 km/s) are sharper and shifted with respect to those in ambient water. The most interesting structural changes occur at u p = 1 km/s. In the primary shock region, we find molecular clusters in which a water molecule has 8 nearest neighbors in a body-centered cubic lattice configuration, indicating the nucleation of ice VII; see Fig. 41 3.4. In this ice VII structure, the central oxygen atom is connected to four of its nearest- neighbor oxygen atoms through hydrogen bonds [73, 74]. We do not find an ice-like structure at other particle velocities (u p = 2.5, 3.0 and 3.5 km/s), which is in agreement with experimental observations of ice VII only for particle velocities between 0.75 and 2 km/s [72]. Figure 3.4: Snapshot of a water molecule cluster that forms an ice-VII-like structure in the compressed region at u p = 1 km/s. Red and white spheres represent oxygen and hydrogen atoms, respectively. Dotted lines indicate the directions of nearest neighbor oxygen atoms (within 3Å from the central oxygen atom). Figures 3.5(a) shows how the bubble volume (normalized by the initial bubble volume) and shape of the largest bubble change with time as a result of shock compression. The insets show oxygen atoms of water molecules at the periphery of the largest nanobubble (D = 10 nm) before the shock wave strikes the nanobubble (t = 1.7 ps) and just before the nanobubble collapses (t = 3.3 ps). The volume versus time data 42 correspond to bubble diameters D = 6, 8 and 10 nm and particle velocity u p = 3 km/s. The bubble collapse time τ is estimated to be 1.1, 1.4 and 1.7 ps for D = 6, 8 and 10 nm, respectively. According to the Rayleigh formula (𝜏= 0.45𝐷 ! !" where ρ is the mass density and ΔP is the pressure difference across the bubble surface), τ for the three bubble sizes are 0.8, 1.1 and 1.4 ps. The differences between the MD results and Rayleigh formula arise for the following reasons: (1) in Rayleigh collapse it is assumed that the bubble is within a uniform fluid, whereas the pressure and density around the nanobubble become non-uniform when the shock front reaches the nanobubble; and (2) unlike MD simulations, Rayleigh equation does not include viscosity and surface tension effects arising from interatomic interactions. 43 Figure 3.5: (a) Normalized bubble volume, Ω(t), versus time, t, for bubbles of initial sizes D = 6, 8 and 10 nm. Insets show snapshots of the largest nanobubble before the shock wave strikes (t = 1.7 ps) and before the nanobubble collapses (t = 3.3 ps). (b) and (c) show velocity vector fields in a cross section of the system just before (t = 3.3 ps) and soon after (t = 3.9 ps) the largest nanobubble collapses. The velocities are measured relative to the reference frame of momentum mirror. In (b) the region inside the white dash lines is enlarged to show a focused jet in the inset; in (c) the lateral flow indicates the onset of jet disintegration. Figures 3.5(b) and 3.5(c) show velocity profiles of water molecules around a shrinking nanobubble (D = 10 nm) at u p = 3 km/s. At t = 3.3 ps the nanobubble has shrunk significantly but not collapsed completely under shock compression. The snapshot in Fig. 3.5(c) is taken immediately after the nanobubble collapses (t = 3.9 ps). Here, 44 contours of the magnitude of velocity are color-coded and each arrow represents the direction and the arrow color the magnitude of the average molecular velocity (averaged over all molecules in a voxel of length 0.5 nm) at that position. Regions inside the white dash lines are enlarged in the insets, and the white arrows indicate the direction of shock propagation. We observe that water molecules around the top and bottom of the nanobubble change directions and their average molecular velocities point towards the center of the bubble soon after the shock wave reaches the proximal side of the bubble. This focusing feature of high-speed molecules from the onset of shrinkage to complete collapse of a nanobubble is akin to the jetting phenomenon observed experimentally in microbubble collapse. The difference between our simulation and experiment [55] is that nanojets do not elongate as much as microjets beyond the distal side of collapsed bubbles. The formation of nanojets has also been reported in an MD simulation of pressurized fluid injection through a nanoscale nozzle [75]. Figure 3.5(c) shows molecular motion around the collapsed nanobubble and lateral flow at the periphery of the collapsed region. The lateral flow is smaller than that observed in the vicinity of a solid surface or a lipid bilayer [31, 65]. It should be pointed out that not all nanobubbles collapse under shock. Experimentally, surface gas nanobubbles are stable under tensile and compression pressures of 6 MPa [76]. Here, we are dealing with cavitation nanobubbles subjected to compression waves with high pressures ranging between 2-19 GPa (depending on the particle velocity), which is why nanobubbles collapse (see Figure 3.6 for more details). 45 Figure 3.6: Time evolution of the shape of the nanobubble of diameter D = 10 nm at u p = 3 km/s. At the onset of nanobubble collapse, we observe a sudden increase in the translational kinetic energy of molecules at the shock front and also an appreciable increase in the rotational energy of those molecules. In the final stages of nanobubble collapse, the molecular vibrational energy also increases to the same level as the rotational energy. For u p = 3 km/s and D = 10 nm, the maximum rotational and vibrational energies are about 0.2 eV, which is well below the self-dissociation energy [2H 2 O → OH - + H 3 O + ] in pure water. Figure 3.7 shows the effects of the shock amplitude (i.e., particle velocity) and nanobubble size on velocity vector fields of water molecules in and around the nanobubble just after it collapses. Panels (a), (b) and (c) display velocity vector fields for a fixed particle velocity, u p = 3 km/s, but different initial diameters D of the nanobubble. Evidently, the size of the jet (ellipsoidal red region) increases with an increase in the 46 nanobubble diameter. Panels (d), (e) and (f) show the effect of the particle velocity on the molecular velocity vector fields in systems containing the same size nanobubble (D = 10 nm). The high-speed region does not grow as much with an increase in the particle velocity as with the nanobubble size. Figure 3.7: Effects of particle velocity and nanobubble size on bubble collapse. Snapshots in panels (a)-(c) are taken for u p = 3 km/s and nanobubble diameter D = 6, 8 and 10 nm, respectively. Snapshots in panels (d)-(f) are for a fixed nanobubble size (10 nm) and u p = 2.5, 3.0 and 3.5 km/s, respectively. The size of the high-speed region (ellipsoidal red regions) increases more with an increase in D than with u p . All snapshots are taken immediately after the nanobubbles collapse. From the onset of nanojet formation and disintegration, we have determined the penetration length l jet and persistence time τ jet for the jet to examine how these quantities are affected by the particle velocity u p and initial bubble radius R. For all particle velocities we find l jet increases linearly with R. Experimentally, the length of jet tip in micron-to-millimeter size bubbles also scales linearly with R [55, 77]. This agreement is somewhat surprising because viscosity and surface tension may play larger roles in 47 nanobubbles than in micron-to-millimeter size bubbles. We find the nanojet persistence time is longer than the nanobubble collapse time by about 0.2 ps. Also the nanojet velocity, V jet = l jet /τ jet , is higher than the shock velocity u s (and hence also u p ), and V jet increases with an increase in the bubble size. Figure 3.8: [(a)-(c)] show the primary shock wave at t = 1.6 ps, 2.9 ps, and 3.3 ps, respectively. Here, the particle velocity, u p = 3 km/s and D = 10 nm. [(d)-(e)] show the secondary water hammer shock at t = 3.5 ps, 3.7 ps, and 3.9 ps, respectively. The color code corresponds to pressure and shock waves are identified by pressure discontinuities. The length scale bar for these panels is shown in (a). Figure 3.8 shows the local pressure from the Virial expression for a nanobubble of diameter D = 10 nm at u p = 3 km/s. The pressure in the compressed region remains constant and the wave front remains planar until the shock wave reaches the nanobubble; see Fig. 3.8(a). The pressure (18.8 GPa) is in reasonable accord with the estimate (19.5 GPa) from the jump condition, 𝑝−𝑝 ! = 𝜚 ! 𝑢 ! 𝑢 ! . Figure 3.8(b) shows the wave front becomes concave due to pressure gradient in the compressed region during the nanobubble collapse. As the nanobubble shrinks, the pressure at the shock front decreases. The shock front becomes planar again just before the complete collapse of the nanobubble; see Fig. 3.8(c). When the primary shock wave hits the proximal side of the 48 nanobubble, surrounding water molecules move into the nanobubble with high speeds. On reaching the distal side of the nanobubble, these high-speed molecules give rise to a secondary water hammer shock wave with a maximum pressure of 29 GPa [Fig. 3.8(d)] while the primary shock continues to move forward with a pressure of 15 GPa. The secondary water hammer shock propagates backward (opposite to the primary shock), spreading spherically with a velocity of 8 km/s as the pressure decreases; see Figs. 3.8(e) and 3.8(f) [54]. Snapshots in Fig. 3.8 show that the water hammer shock is not close to the lateral boundaries and so the interaction between mirror-image waves is negligible. Despite the differences in length scale and pulse duration, significant pressure amplification due to secondary water hammer shock wave and its rapid dissipation have also been observed in experiments [56] and continuum simulations [78]. 3.4 Conclusion In conclusion, the F-ReaxFF MD simulations of shock propagation in water are in good agreement with experimental results for the Hugoniot compression curve. We find molecular clusters with ice-VII like structure at particle velocity of 1 km/s. This is well supported by experimental observations of ice VII for particle velocities between 0.75 and 2 km/s. We have examined the effects of shock amplitude and the initial nanobubble size on the dynamics of nanobubble shrinkage and collapse. During shrinkage, we observe a focused nanojet whose length scales linearly with the nanobubble radius. This scaling relation has also been found experimentally in shock-induced collapse of micron- to-millimeter size bubbles. In the next chapter we will describe multimillion-atom MD simulations to examine the effect of nanojets from nanobubble collapse on lipid bilayers. The results 49 indicate that the nanojet impact creates a transient localized deformation of non-uniform width and poration in the lipid bilayer. In addition, we observe shear flow of water on the lipid bilayer which has also been reported in experiments on the interaction of cavitation bubbles with cells [31, 65]. Transient cell poration has potential applications in gene therapy and drug delivery. 50 4 Poration of Lipid Bilayers by Shock-Induced Nanobubble Collapse 4.1 Introduction In recent years, non-invasive drug- and gene-delivery approaches have garnered significant interest because of direct applications in cancer treatment and gene therapy. Much of the experimental effort is focused on designing a targeted approach that has both spatial and temporal specificities. Research in this area relies mainly on the use of electric fields or pressure waves to enhance the permeability of cell membranes. In one of the most commonly used techniques known as electroporation [79, 80], electric fields are applied across the cell to increase the cell-membrane permeability. Reversible electroporation, in which the cell permeability is enhanced temporarily, is used for drug delivery and gene therapy. Electric fields applied over a sufficiently long time can kill the cell because of temperature elevation resulting from Joule heating. This irreversible electroporation process is commonly used in the food industry to inactivate microbes and also in minimally invasive treatment of cancerous tissues. Sonoporation is another promising DNA-, protein- and drug-delivery approach [81-83]. To achieve high efficiency in sonoporation, in vivo gas bubbles are used in conjunction with diagnostic level ultrasound exposures [84-86]. Sonoporation experiments show that the collapse of bubbles by ultrasound generates water jets [55] whose impact on the cell membrane increases the permeability thereby allowing the intracellular delivery of drug/gene payload [31]. Shock waves in tandem with nanobubbles provide another promising approach to targeted delivery of drugs and genes. 51 Shock-wave phenomena [87], such as extracorporeal shock wave lithotripsy, have been used in living tissues. 4.2 Simulation Method Here we focus on poration of lipid bilayers by the interaction of shock waves with nanobubbles. We have performed molecular dynamics (MD) simulations to study the impact of shock waves on nanobubbles in the vicinity of a dipalmitoylphosphatidylcholine (DPPC) phospholipid bilayer embedded in water. We use GROMACS [44] to simulate the SPC model for water and DPPC bilayer. During equilibration the DPPC bilayer and water molecules are coupled to a thermostat and barostat, bond lengths are kept constant, and periodic boundary conditions are imposed in all directions. After equilibration, a nanobubble is created close to the lipid bilayer by removing water molecules within a sphere of diameter D. The distance between the right edge of the nanobubble and the left boundary of the lipid bilayer is 1.9 nm. The simulations were done for D = 10, 20, 40 nm and these systems contain about 1.9, 6.2, and 30.8 million atoms, respectively; the dimensions of the corresponding MD cells are 19×19×61, 34×34×82, and 64×64×92 nm 3 . The parameters for bonded and non-bonded interactions as well as partial charges on atoms in the DPPC bilayer system are taken from refs. [49, 88]. The force field for DPPC molecule is validated against experimental data for area per lipid and order parameter [88]. We choose the SPC model for water because it works better with simulations involving interfaces although other water models like SPC/E reproduce radial distribution function and diffusion more accurately. The main reason for using the SPC model is that it gives an area per lipid value and the width 52 of DPPC/water interface close to experimental values. Bending modulus and gel to liquid-crystalline phase transition temperature obtained from MD simulation are also in good agreement with experiments [25, 89]. The lateral diffusion coefficient for DPPC molecules in bilayers has been measured experimentally [90, 91]. In the fluid phase the value ranges between 0.6-2 × 10 -7 cm 2 /s, whereas in the gel phase it ranges between 0.04- 16 × 10 -9 cm 2 /s. In coarse-grained MD simulations based on all-atom MD simulation described here, the lateral diffusion coefficient is between 1-4 × 10 -7 cm 2 /s in the fluid phase and between 0.5-4 × 10 -9 cm 2 /s in the gel phase [92]. Figure 4.1: Shock velocity versus particle velocity. The simulation results for SPC water are in good agreement with experimental data. The inset shows the setup for shock simulation. The gray plate is the momentum mirror. 53 We checked the validity of the SPC model for water under shock. An MD simulation was performed to equilibrate the water system by itself for 1 ns using a time step of 2 fs. Periodic boundary conditions (PBC) were imposed during equilibration. To apply shock, we inserted a vacuum layer of thickness 2 nm at the end of the MD box in the x direction and moved the system with a constant particle velocity u p towards a momentum mirror [93], i.e. along –x in the inset of Fig. 4.1. The mirror reverses the x component of atomic velocity if an atom crosses the mirror plane located at x = 0.2 nm. This generates a planar shock in the +x direction. Pressure and temperature couplings as well as bond-length constraints are turned off in shock simulations. The shock velocity u s is determined by monitoring the shock front, i.e. discontinuity in pressure or mass density of water at two instants of time. Figure 4.1 shows that the MD results for u s as a function of u p are in good agreement with experimental data [72]. Next, we performed MD simulations to equilibrate initial configurations of lipid bilayers and water molecules at temperatures T i = 300 and 323 K and pressure P = 1 bar using a time step of 2 fs. After equilibration we created a bubble near the bilayer (see Fig. 4.2) and applied planar shock as discussed above. A number of simulations were performed for different nanobubble diameters (D = 10, 20 and 40 nm), particle velocities (u p = 0.4, 0.7 and 1 km/s), and initial temperatures (T i = 300 and 323 K). 54 Figure 4.2: Schematic of the MD box containing a nanobubble in water and a lipid bilayer. The gray plate is the momentum mirror. 4.3 Results When a planar shock front hits the proximal side of a nanobubble, water molecules from the bubble periphery accelerate towards the center of the bubble and form a nanojet. The size of the nanojet depends on the particle velocity and nanobubble diameter. As the particle velocity increases from 0.4 km/s to 1.0 km/s, the number of water molecules in the nanojet for a fixed value of D increases by an order of magnitude. Figure 4.3 displays instantaneous molecular velocities averaged in voxels of dimension 0.5 nm for a bubble of initial diameter D = 40 nm under the impact of a shock front moving with velocity u p = 0.7 km/s. Figure 4.3(a) shows that velocities of water molecules in the domain of the shrinking nanobubble are focused in the form of a nanojet. As the particle velocity increases from 0.4 km/s to 1.0 km/s, the average x- component of molecular velocities inside the nanojet increases from 2.6 km/s to 3.5 km/s 55 for all simulated nanobubble sizes. For D = 40 nm, we find the length of the nanojet l jet is 57 nm. Our results for other nanobubble sizes and particle velocities also indicate that the length of the jet scales as l jet ≈ 1.5D. Surprisingly, the same linear scaling has been observed in experimental studies of shock-induced collapse of micron-to-millimeter size bubbles [77, 94]. Figure 4.3: Snapshots of velocity profile for the system with D = 40 nm, T i = 300 K and u p = 0.7 km/s. Arrows show the direction of average molecular velocities and the velocity magnitudes are color-coded. Panel (a) shows a nanojet in the system at t = 20 ps. The white vertical region is the bilayer. Panel (b) shows a spreading flow at t = 24 ps resulting from the impact of the nanojet on the lipid bilayer. We have performed additional simulations of nanobubble collapse in water at a particle velocity of 3 km/s. We find the bubble collapse time to be 0.9, 1.2 and 1.5 ps for bubble diameter 6, 8 and 10 nm, respectively. Using the Rayleigh formula ( € τ =0.45D ρ ΔP , where ρ is the mass density and ΔP is the pressure difference across the bubble surface), we obtain τ to be 0.8, 1.1, and 1.3 ps for the three bubble sizes. The differences between our calculation and the estimates from the Rayleigh formula arise 56 from the facts that: (1) in Rayleigh collapse it is assumed that the bubble collapses within a fluid of uniform pressure and density, whereas in our simulations pressure and density become nonuniform due to the shock front; and (2) the Rayleigh equation does not include viscosity and surface tension effects which arise due to interatomic interactions. From the onset of nanojet formation and disintegration, we have determined that the persistence time τ jet for the jet exceeds the bubble collapse time by at least 0.2 ps. In Fig. 4.3(b) we show the interaction between water molecules in the nanojet and the DPPC molecules in the bilayer. Water molecules in the nanojet form a spreading flow after hitting the leaflet of the DPPC bilayer [31]. Figure 4.4 shows the averaged lateral velocity of water molecules in the vicinity of a lipid bilayer versus the distance from the center of the bilayer for a nanobubble of diameter 40 nm. Over half the nanobubble radius the lateral velocity is larger than the thermal velocity. The peak in the lateral-flow velocity appears when the nanojet hits the bilayer. The distance over which the lateral velocity is larger than the thermal velocity is half the nanobubble radius. Experiments [95] on millimeter size bubbles in the vicinity of a hard surface indicate that this distance is of the order of the bubble radius. The differences between experimental and our MD results are due to the fact that bubble sizes differ by several orders of magnitude and the surfaces are soft in MD simulation and hard in experiments. We also observe vortices in the collapsed bubble when water molecules bouncing back from the bilayer encounter other water molecules in the incoming shock wave. 57 Figure 4.4: Average lateral velocity of water molecules in the vicinity (1nm) of a lipid bilayer versus the distance from the center of the bilayer for a nanobubble of diameter 40 nm. The inset shows the results for a nanobubble of diameter 10 nm. In both systems u p = 0.7 km/s and T i =300 K. 58 Figure 4.5: (a) and (b) are snapshots of the density of water at t = 20 and 28 ps. Here D = 40 nm, T i = 300 K and u p = 0.7 km/s. The central blue region is the lipid bilayer. Panel (a) shows the nanojet traveling towards the distal side of the nanobubble. Panel (b) shows the deformed bilayer and water-hammer shock. Figure 4.5(a) shows the water density around the DPPC bilayer (blue region) just before the bilayer is hit by the water nanojet from a collapsing nanobubble of diameter 40 nm at u p = 0.7 km/s. The curved blue region to the left of the bilayer indicates that the bubble has not collapsed completely and the water density around the nanobubble is close to the normal density of water. The water density around the bilayer leaflet closer to the collapsed bubble increases to 1.5 g/cc. Figures 4.5(b) shows that the deformed bilayer is hemi-spherical. In addition, we observe water-hammer shock when water molecules in the nanojet hit the distal side of the nanobubble. This secondary water-hammer shock spreads spherically and its initial speed (until 4 ps after formation) is approximately 1.6 km/s. The amplitude of the secondary shock decreases but its velocity increases with time. Secondary water-hammer shocks have been observed in experiments [56] and continuum simulations [78]. 59 After the nanojet impact, the DPPC bilayer deforms and becomes significantly disordered. To characterize the evolution of disorder in the bilayer, we calculate the order parameter tensor € ˜ S defined as:  S ij = 1 2 3cosθ i cosθ j −δ ij , (4.1) where € θ i is the angle between the i th “molecular vector” and the bilayer normal. The “molecular vector” for the i th CH 2 unit has z component from C i-1 to C i+1 and the y component in the plane containing C i-1 , C i , and C i+1 .The brackets denote an ensemble average. From the diagonal elements of € ˜ S , we define the normalized order parameter at time t; S n (t)= 2  S xx (t)+  S yy (t) 2  S xx (0)+  S yy (0) , (4.2) where t = 0 is the time at which the nanobubble is inserted in the system. Figure 4.6 shows how rapidly the impact of the nanojet partially destroys order in the lipid bilayer. After impact, the order parameter within a radius of 5 nm from the center of the hemi- spherically deformed bilayer decreases by 60% in about 4 ps. As expected, the disorder in the bilayer decreases as the lateral distance from the deformation center increases. For example, the drop in the order parameter 20 nm from the center of the deformed bilayer region is 36%. 60 Figure 4.6: Temporal changes in the normalized order parameter for the lipid bilayer region directly impacted by the nanojet. At t = 16 ps, the density of water molecules around the bilayer starts to increase and hence the bilayer begins to show increasing disorder. The impact of the nanojet causes poration in the lipid bilayer. Figure 4.7 shows poration resulting from the impact of the collapsed nanobubble of initial diameter 40 nm at u p = 0.7 km/s. The poration was calculated by dividing the impacted region of the bilayer into pixels of size 0.1 nm and determining the area of empty pixels, i.e., those containing no lipid molecules. For the bilayer initially in the gel phase at T i = 300 K, the nanojet impact increases poration by a factor of 30 over its normal value before the nanojet impact; see Fig. 4.7(a). For the bilayer initially in the liquid phase at T i = 323 K the poration increases by another factor of 5 relative to the poration in the gel phase; see Fig. 4.7(b). In the liquid phase at 323K, the maximum nanopore size is 0.7 nm as compared to 0.4 nm in the gel phase. The poration varies with the particles velocity and 61 nanobubble diameter. At u p = 0.4 km/s, we do not observe any significant change in the porosity of the gel phase for the three bubble sizes we have considered. However, at u p = 1 km/s the maximum nanopore size increases to 0.3 nm for D = 10 nm and it increases linearly with the initial diameter of the nanobubble. Figure 4.7: Poration of lipid bilayers by collapsed nanobubbles. Here u p = 0.7 km/s and D = 40 nm. In (a), the bilayer was initially in the gel phase at T i = 300 K and in (b) it was in the liquid phase at T i = 323 K. In the deformed DPPC bilayer that was initially in the liquid phase, the pores are large enough (≥ 0.5 nm) to allow rapid translocation of water molecules. Translocation events are observed for u p = 1.0 km/s and D ≥ 10 nm and also for u p = 0.7 km/s and D = 40 nm; see fig. 4.8. Water molecules can diffuse through the lipid bilayer in the absence of shock, but the diffusion is almost four orders-of-magnitude slower than in bulk water. The poration by nanojet impact and the large pressure difference (~ 9 GPa) across the bilayer combine to shorten the average time of passage for water molecules by 6 orders of magnitude. 62 Figure 4.8: Passage of a water molecule (white dot) across the lipid bilayer for the system with D = 10 nm, T i = 323 K, and u p = 1.0 km/s. The water molecule, initially in the red hydrophilic region on the left side of the bilayer, crosses the hydrophobic region (light green) and ends up on the right side (red) of the bilayer. The bilayer poration is, however, temporary because the nanopores disappear and the bilayer heals after the passage of shock wave. Figure 4.9 shows the relaxation of a lipid bilayer after the impact of the nanojet. Here the MD simulation was performed on a configuration of a deformed bilayer (including surrounding water molecules in an 11 nm wide region) embedded in water under ambient conditions. The simulation was run for 5 ns in the microcanonical ensemble using periodic boundary conditions. After relaxation, the deformation in the bilayer disappears. 63 Figure 4.9: (a) Deformed and (b) relaxed lipid bilayer configurations. 4.4 Summary In summary, multimillion-atom MD simulations reveal the mechanism of transient poration in lipid bilayers by shock-induced collapse of nanobubbles. When a planar shock front strikes a nanobubble, water molecules from the bubble periphery accelerate towards the center of the bubble to form a nanojet. The length of the nanojet scales linearly with the initial nanobubble size which, surprisingly, is also observed in experimental studies of shock-induced collapse of micron-to-millimeter size bubbles. The 64 MD simulations reveal that the nanojet impact significantly deforms and thins the lipid bilayer and water molecules in the nanojet form a spreading flow pattern after impact. Deformation and thinning of bilayers combined with large pressure gradients across and spreading flow around the bilayers create transient nanochannels through which water molecules translocate across the bilayer. 65 5 Cholesterol Translocation in a Phospholipid Membrane 5.1 Background In the standard fluid mosaic model, cell membranes are dynamic structures in which molecular motions play a key biophysical role [96]. These motions include diffusion of proteins and lipids in the plane of the bilayer and flip-flop transport of lipids between the two leaflets composing the bilayer. While flip-flop motion of phospholipids is typically considered to be a slow process, with little relevance on physiological timescales [17], this is not the case for cholesterol (CHOL). While the transport rate of CHOL has been elusive from an experimental point of view, most studies suggest that it is relatively fast [18]. This rapid flip-flop transport has profound implications for the biological role of CHOL in cell membranes. The potential role of CHOL in controlling the mechanical behavior of membranes can be understood in the context of the bilayer couple hypothesis originally suggested by Sheetz and Singer, in which the overall curvature of a membrane is determined by the relative areas of the two leaflets [97]. The membrane curves away from the leaflet with the relative excess of surface area, much in the way a bimetallic couple curves away from the metal film with the greater thermal expansion modulus. Lange and coworkers first suggested that CHOL may act as a “surface area buffer” in cell membranes, stabilizing changes in membrane shape by rapidly flipping between the bilayer leaflets, decreasing the surface area of the relatively compressed leaflet while increasing the surface area of the relatively expanded leaflet [98]. This was consistent with their observations that 66 CHOL concentration determines the equilibrium curvature of membrane ghosts derived from red blood cells [99]. The concept of CHOL as a molecule to facilitate membrane curvature seems at odds with its observed effect on the mechanics of membrane bending, however. Multiple studies have demonstrated that as the CHOL content of a membrane increases, the membrane becomes more rigid, requiring more force to produce the same degree of curvature [100-103]. Recently, Bruckner and coworkers suggested that the effect of CHOL on membrane mechanics may be time-dependent, resulting in an increased instantaneous elastic modulus but allowing for relaxation of curved membranes at longer timescales by a flip-flop mechanism [20]. This is consistent with the results of Yanagisawa and coworkers, who observed that in membrane vesicles in which CHOL- rich regions are formed, deformation at long timescales occurred preferentially where CHOL concentration was the greatest [104]. Another potential role for CHOL flip-flop in cell physiology is related to the interleaflet coupling of lipid phase separation. The chemical composition of the eukaryotic plasma membrane is asymmetric [105, 106]; that is, if the bilayer membrane is considered as two monolayers (or “leaflets”) packed tail-to-tail, the leaflet facing the exterior of the cell contains different lipids at different concentrations than the leaflet facing the interior of the cell. In terms of lipid phase state, it is striking that the interior leaflet of a eukaryotic plasma membrane has a composition that does not form phase segregated liquid ordered (l o ) and liquid disordered (l d ) domains in vitro while the exterior leaflet will segregate into these two liquid phases [107]. However, evidence from biochemical assays [108] and recent in vitro experiments [109-112] has shown that the 67 phase separation in one leaflet of an asymmetric bilayer can induce domain formation in the opposing leaflet. Rapid CHOL transport between the leaflets may be the mechanism of this induction [113]. Understanding the role of CHOL in the membrane requires understanding the dynamics of CHOL flip-flop. Our hypothesis is that membrane bending stress can be relieved by cholesterol flip-flop across the membrane with a characteristic timescale in the microsecond range. Broad ranges have been reported for CHOL flip-flop rates, but the general trend has been for the upper limit on flip-flop half-time to decrease as analytical techniques have become more advanced. Early studies reported that this half time was definitively shorter than phospholipid flip-flop half time, but ranged broadly from 3 s [98] to tens of minutes [114]. More recent studies have clarified that the half time is certainly shorter than 1 s [115], and Bruckner et al. used an NMR-based technique to identify CHOL residence time on a single leaflet of the bilayer as less than 10 ms [20]. Wide discrepancies remain in experimental studies of CHOL flip-flop, however. In 2011, Garg and coworkers reported neutron scattering results consistent with a CHOL flip-flop half time in the range of hundreds of minutes [19]. Physical insights into the molecular details of flip-flop will be a useful tool in resolving these discrepancies. In fact, molecular dynamics (MD) simulations of CHOL flip-flop transport have generally identified fast flip-flop rates, with half times in the range of tens of ns to ms [116, 117]. There have also been MD simulations identifying relatively slow CHOL flip- flop, however, with recent work by Bennett and Tieleman supporting half times in the range of minutes for raft-like bilayers [118]. The MD simulations that have been 68 performed so far have used coarse-grain methods and/or steered approaches to capture flip-flop events. Here, we discuss an all-atom MD simulation study of spontaneous flip-flop events in a DPPC/CHOL bilayer at physiological CHOL concentrations. We have examined the effects of CHOL translocation on structural, dynamical and mechanical properties of the DPPC bilayer membrane. We performed two simulations: in one case, the initial configuration of the bilayer leaflets had a slight imbalance (31% and 29%) in CHOL concentrations; and in the other case, each leaflet had the same number of CHOL molecules in the initial configuration. We ran the first simulation for 15 µs and observed 24 CHOL flip-flop events. In the second case, we performed the simulation with the parallel replica method [119, 120] and observed 5 CHOL translocations over a period of 2.9 µs. On average, a CHOL molecule migrates across the lipid bilayer in about 73 ns after a flip-flop event is triggered. The potential-energy barrier for a CHOL translocation event is estimated to be 100 KJ/mole, whereas the free energy barrier is significantly lower (21.5 KJ/mol). We have also calculated the effects of CHOL translocation on mechanical stresses across the lipid bilayer and on the curvatures of bilayer leaflets. We find that the compressive lateral pressure and curvature of a leaflet decrease after a CHOL molecule leaves the leaflet. 5.2 Simulation Method Molecular dynamics simulations were performed with the GROMACS 4.5 package [44]. All-atom MD simulations were performed on a system consisting of 360 DPPC and 152 CHOL molecules, which corresponds to an average CHOL concentration 69 of 30%. The system also contained 14,620 simple point charge (SPC) water molecules. Altogether there were 66,268 atoms in the system. We used the DPPC force field of Tieleman et al. [88] and a slightly modified GROMOS force field for CHOL molecules [121, 122]. Periodic boundary conditions were applied in all directions. The cutoff for the non-bonded Lennard-Jones interaction was equal to 0.9 nm. The Coulomb interaction was calculated with the particle mesh Ewald (PME) method. In PME the real-space contribution was obtained with a cutoff of 0.9 nm and the reciprocal-space contribution was calculated using the fast Fourier transform technique with a grid spacing of 0.16 nm. We used the Nose-Hover scheme to maintain the temperature at 323K. DPPC, CHOL and water molecules were coupled independently with the thermostat using a time constant of 0.2 ps. We also used the Parrinello-Rahman scheme for semi-isotropic coupling in the x- y and z directions. The coupling constant was taken to be 2.5 ps and the pressure in the system was 1 bar. Equations of motions were integrated with the velocity-Verlet algorithm using a time step of 2 fs and we ran the simulation for 15 µs. We created the DPPC bilayer by placing DPPC molecules on a grid and then randomly replacing a fraction of DPPC molecules by CHOL molecules. The size of the MD box was chosen to give a reasonable initial value for the average area per lipid. Subsequently, we added water molecules randomly and equilibrated the entire system for one hundred nanoseconds. We also studied CHOL flip-flop dynamics with the parallel replica MD simulation scheme [119]. Parallel replica is an efficient approach to sample rare events, such as CHOL translocation across a membrane. In this scheme, one generates multiple replicas of an atomic configuration, i.e., atomic positions at a given instant of time. These 70 replica configurations are implemented on a parallel machine with identical processors. The number of available parallel processors p is partitioned into n groups, where n is equal to the number of replicas. Each replica configuration is implemented on p/n number of processors. Although the initial atomic configurations of the replicas are identical, the initial atomic velocities in the replicas are not the same. They are chosen randomly from the Maxwell distribution corresponding to the temperature of the system. The n replica simulations are run concurrently for a dephasing time on the order of a few ns, and after dephasing the parallel replica simulations are run until a rare event is observed in one of the replicas. At that instant, the simulations are stopped and the atomic configuration of the replica with the rare event is implemented on the other (p/n-1) groups of processors. Again, new Maxwellian random velocities are assigned to all the atoms in the n replicas and then replica simulations are run in parallel until the next rare event is observed. In our parallel replica simulation, we run five parallel replicas on five sets of identical processors. Whenever the difference between the maxima and minima of the z coordinate of the hydroxyl group of a CHOL molecule in a replica exceeds 3 nm (bilayer thickness is 4-5 nm), we record that as a “rare event” and advance the clock by five times the simulated time of the rare event. 71 5.3 Results To study CHOL flip-flop dynamics resulting from a small difference in CHOL concentrations in the bilayer leaflets, we set up an initial MD configuration in which the CHOL concentrations in the two leaflets were 31% and 29%. Hereafter the leaflets with 31% and 29% initial CHOL concentrations will be referred to as the upper leaflet (UL) and lower leaflet (LL), respectively. In this initial configuration, the CHOL molecules were placed randomly in the UL and LL. We ran the MD simulation for 15 µs and observed 22 CHOL translocation events from the UL to the LL and 2 of those 22 CHOL molecules translocated back to the UL. A translocation event is determined by monitoring the z-coordinate of the hydroxyl group of a CHOL molecule, z OH . If the change in z OH exceeds two-thirds the membrane thickness (~ 4-5 nm), we consider that a translocation or “flip-flop” event. For each flip-flop event we monitor residence times in the bilayer leaflets before the onset of the flip-flop event and after the trans-membrane migration of the CHOL molecule. We also monitor the CHOL residence time in the middle of the bilayer. Figure 5.1 (a)-(c) show snapshots of a CHOL molecule (red and cyan) migrating from the UL to the LL. The snapshots also show the surrounding DPPC molecules (blue and yellow) in the bilayer. For clarity, other CHOL or water molecules are not shown here. 72 Figure 5.1: Snapshots of a translocating CHOL molecule (head group is red and ring/tail structure is cyan) surrounded by DPPC molecules. For clarity, other CHOL and water molecules are not shown here. Panels (a), (b) and (c) show CHOL translocation from the UL to the LL. These snapshots were taken at time t = 166, 173 and 189 ns. Panel (d) shows how the tilt angle (angle between the bilayer normal and the line joining atoms C5 and C21 in a CHOL molecule; see Fig 5.2) changes with time. Panels (e), (f) and (g) are snapshots of the same CHOL molecule at t = 13407, 13412 and 13418 ns, respectively. These snapshots show translocation from the LL to the UL. Panel (h) shows changes in the CHOL tilt angle with time. 73 Figure 5.2: Structure of a cholesterol molecule with the C5 and C29 atoms highlighted. These atoms are used to define the tilt angle. Cholesterol molecules with their rigid hydrocarbon rings and flexible hydrocarbon tails are usually aligned parallel to the bilayer normal (z axis). When a CHOL molecule starts to migrate across the bilayer, the angle between its molecular axis (i.e., the line joining atoms C5 and C21 in a CHOL molecule; see Fig. 5.2) and the bilayer normal begins to increase. Near the middle of the bilayer, the translocating CHOL molecule lies almost entirely in the x-y plane; see Fig. 5.1(b). When the molecule leaves the middle of the bilayer and approaches the opposite leaflet, its molecular axis begins to tilt towards the bilayer normal with its hydroxyl group facing that leaflet; see Fig. 5.1(c). Figure 5.1(d) shows changes in the angle between the axis of the translocating CHOL molecule and the bilayer normal as the molecule migrates across the bilayer. We observe similar changes in the angle for every translocating CHOL molecule. Figure 5.1 (e)-(g) show the same CHOL molecule migrating back to the UL approximately 13 µs later. In this reverse translocation, the residence time for the molecule in the middle of the bilayer was only 7 ns. The molecule took 37 ns to diffuse from the LL to the middle of the bilayer and another 29 ns to diffuse to the UL. 74 Figure 5.3: (a) The number of CHOL molecules in the upper (red squares) and lower (blue circles) leaflets as the simulation progressed. (b) A histogram of total translocation times for the 24 flip-flop events we observed in 15 µs. 75 We find a wide range of time intervals between successive flop-flop events. As shown in Fig. 5.3(a), 13 CHOL translocation events occurred in the first 5 µs. The shortest time interval between consecutive flip-flop events was 20 ns and the longest 2 µs. The frequency of CHOL flip-flop events dropped after 5 µs. Between 5 and 10 µs, there were 7 UL → LL flip-flop events. After 10 µs, we observed 2 translocations from the UL to the LL and 2 in the opposite direction. The histogram in Fig. 5.3(b) shows the distribution of translocation event durations. The shortest translocation time is 20 ns and the longest 200 ns. The average translocation time for the 24 flip-flop events we observe is 73 ns. The first CHOL flip-flop event happened at 165 ns. In that event, the CHOL molecule moved from the UL to the middle of the bilayer in about 7 ns, stayed there for 10 ns, and then moved to the LL in about 7 ns. For translocating CHOL molecules, the residence time in the middle of the bilayer varied between 10 ns and 90 ns. Similarly, we observed a wide range of timescales for transitions from the UL to the middle of the bilayer (10 - 70 ns) and from there to the LL (5 - 45 ns). A couple of times during the simulation (15 µs), we observed a pair of CHOL molecules residing simultaneously in the middle of the bilayer while translocating from the UL to the LL. Snapshots in Fig. 5.4 show the transmembrane migration for a pair of molecules. The second flip-flop event was triggered about 20 ns after the first one. The two molecules resided in the middle of the bilayer for about 60 ns and 13 ns, and the second molecule reached the LL 11 ns before the first one. 76 Figure 5.4: A pair of CHOL molecules (head groups are red and ring/tail structures are cyan) translocating almost simultaneously from the UL to the LL. Snapshots (a), (b), (c) and (d) were taken at time t = 4340, 4385, 4465 and 4505 ns, respectively. We use the method suggested by Hosfaβ et al. [123] to calculate the area per DPPC molecule, A DPPC , and the area per CHOL molecule, A CHOL . We first calculate the bilayer height, h= V box −N W V W A box , (5.1) where V box and A box are respectively the volume and area of the simulation box, and V W is the volume of a water molecule. Following Hosfaβ et al., we set V W = 0.0312 nm 3 (N W = 14620) and calculate the volume of a DPPC molecule in the UL: V DPPC UL = 1 N DPPC UL V box −N W V W 2 −N CHOL UL V CHOL " # $ % & ' , (5.2) 77 where € N DPPC UL is half the total number of DPPC molecules and € N CHOL UL is the number of CHOL molecules in the upper leaflet, which varies with time. Hosfaβ et al. suggest a value of 0.593 nm 3 for € V CHOL , the crystalline volume of a CHOL molecule. The relation, A DPPC UL h 2 =V DPPC UL , (5.3) gives € A DPPC UL , and € A CHOL UL is obtained from A CHOL UL = 1 N CHOL UL A box −A DPPC UL N DPPC UL " # $ % . (5.4) We follow the same procedure to calculate the areas for DPPC and CHOL molecules in the lower leaflet. We find that A CHOL = 0.28 nm 2 at t = 100 ns and 15 µs. At 15 µs, the CHOL concentration in the UL and LL are 25% and 34% respectively. Hosfaβ et al. report A CHOL = 0.28 nm 2 at 25% CHOL concentration and 0.27 nm 2 at 40% CHOL concentration. At t = 15 µs the corresponding A DPPC are 0.57 nm 2 and 0.52 nm 2 . Hosfaβ et al. report A DPPC = 0.55 nm 2 at 25% CHOL concentration and 0.54 nm 2 at 40% CHOL concentration. We calculate the order parameter tensor € ˜ S defined as: € ˜ S ij = 1 2 < 3cosθ i cosθ j −δ ij > , (5.5) where € θ i is the angle between the i th “molecular vector” and the bilayer normal. The “molecular vector” for the i th CH 2 unit has z component from C i-1 to C i+1 and the y component in the plane containing C i-1 , C i , and C i+1 .The brackets denote an ensemble average. From the diagonal elements of € ˜ S , we obtain the deuterium order parameter S CD 78 € S CD = 2 3 S xx + 1 3 S yy . (5.6) Figure 5.5: Time evolution of order parameters for DPPC molecules in the upper (red squares) and lower (blue circles) leaflets. We investigated the effect of CHOL translocation on the deuterium order parameter, S CD , for the acyl chains of DPPC molecules. Figure 5.5 shows the effect of CHOL flip-flop on the order parameter. The results shown here were averaged over all the atoms of acyl chains in the UL and LL. At time t = 100 ns the value of S CD for the UL is higher than that for the LL because the UL has more CHOL molecules (80) than the LL (72). The order parameter for DPPC molecules in the LL increases with CHOL migration from the UL to the LL. This is consistent with experimental [124] and MD simulation [123, 125] studies, which show that an increase in CHOL concentration reduces the number of gauche defects in the hydrocarbon chains of DPPC molecules and 79 thus the value of the order parameter is increased. Figure 5.5 shows that the order parameter for the LL levels off after 10 µs, since there are very few flip-flop events in the last 5 µs of the simulation. We have also determined the lateral self-diffusion coefficients of DPPC and CHOL molecules in the two leaflets. The lateral self-diffusion coefficients for DPPC and CHOL molecules are obtained from the slopes of mean square displacements in the x-y plane as a function of time. The mean square displacements are averaged over 500 time frames over a time period of 100 ns. We find the diffusion coefficients for both CHOL and DPPC molecules are nearly the same when these molecules are in the same leaflet. However, the self-diffusion coefficients for molecules in opposite leaflets vary significantly because the two leaflets have different CHOL concentrations. For example, after 5 µs the diffusion coefficient was 35 % higher in the UL than in the LL. This was caused by the translocation of CHOL molecules, which made the CHOL concentration in the LL 5% higher than in the UL. The higher the CHOL concentration, the less mobile are the molecules in the leaflet [123]. Figure 5.6 shows changes in the potential energy of a CHOL molecule with time as the molecule translocates from the UL to the LL. The potential energy is calculated by summing up the bonded and non-bonded (Lennard-Jones and electrostatic) energies for the translocating CHOL molecule. The potential energy of the translocating molecule fluctuates around -350 kJ/mole when the molecule is in the UL or LL. In the middle of the bilayer (see inset), the potential energy of the molecule rises to -250 kJ/mole. Thus, the potential energy barrier to CHOL inter-leaflet migration is around 100 kJ/mole. In 80 every CHOL flip-flop event, we find almost the same value for the potential-energy barrier. Figure 5.6: Potential energy of a flip-flopping CHOL molecule as a function of time. The potential energy barrier is about 100 KJ/mol. The inset shows the cholesterol molecule (head group is red and ring/tail structure is cyan) in the middle of the bilayer. To correlate structural changes with the mechanical response of the membrane, we calculated the density profile and mechanical stresses across the membrane (i.e., along z, the bilayer normal) at various time intervals. The membrane thickness was calculated from the difference between the two highest peaks in the density profile and the lateral pressure, π(z), from the diagonal components of the stress tensor . ; . (5.7) Here m i is the mass of the i th atom,  v i is its velocity, and and are the position vector and force between atoms i and j, respectively. The Irving-Kirkwood contour method is € P(z) € P(z) = m i i ∑ v i ⊗v i + 1 V r ij i< j ∑ ⊗F ij € π(z) = (P xx + P yy )/2−P zz € r ij € F ij 81 used to calculate stress profiles normal to the DPPC bilayer (z direction). We divide the system into slabs of length 0.2 nm in the z direction and calculate stresses by averaging over 5000 frames spanning 50 ns. Figure 5.7: (a) and (b) are density and stress, Π(Z), profiles, respectively, across the bilayer at time t = 15 µs. Membrane midplane is at Z = 0. Figure 5.7(a) and 5.7(b) show the density and lateral pressure, respectively, at time t = 15 µs. The two leaflets have reached mechanical equilibrium because they have nearly the same stress profiles. Cholesterol molecules may flip-flop between the leaflets, 82 but the probability for the UL → LL translocation is balanced by the probability for the LL → UL translocation. Curvature plays an important role in the structure and mechanical behavior of biomembranes [126]. Thus, we calculated the curvature energies of UL and LL to determine if the CHOL flip-flop dynamics lead to membrane relaxation. The bending free energy is defined as, G bend = K b 2 dx 1 ∫ dx 2 (κ 1 (x 1 ,x 2 )+κ 2 (x 1 ,x 2 )) 2 , (5.8) where K b is the bending rigidity. κ 1 and κ 2 are principal curvatures, defined in terms of the eigenvalues of the curvature matrix: Κ= κ 11 κ 12 κ 21 κ 22 " # $ $ % & ' ' , (5.9) where € κ ij = ∂ 2 h ∂x i ∂x j and h is the height of the membrane. We calculated the curvature energy of the UL and LL at intervals of 2.5 µs. We divided the x-y plane into pixels of size 1.5 nm. This grid size was chosen to ensure that each pixel was occupied. The height of the membrane was determined from the z coordinate of either the N4 or P8 atom of DPPC molecules. The calculations involved averaging over 21 frames spanning 4 ns. Our results indicate that CHOL migration from the UL to the LL decreases the curvature energy of the UL. However, the curvature energy of the LL remains more or less the same over the entire duration of the simulation. At 15 µs, the curvature energies of the two leaflets are very close to each other. 83 Finally, we performed another all-atom MD simulation in which the two leaflets had the same number of CHOL molecules. In this case, we used the parallel replica method (PRM) to study CHOL flip-flop events. In PRM the initial configuration of the system was taken from the MD simulation after four CHOL molecules had migrated from the UL to the LL of the bilayer. So in the initial configuration of the parallel replica MD simulation, each leaflet had 76 CHOL molecules or 30% CHOL concentration. We created five replicas of this configuration and assigned each atom in the system a random velocity chosen from the Maxwell distribution corresponding to a temperature of 323K. Using a thermostat and keeping the pressure constant at 1 bar, the five replicas were run in parallel on 900 identical processors of a Linux cluster (180 processors for each replica) for a dephasing time of 5 ns. There were no flip-flop events during the dephasing time. After dephasing, the parallel replica simulations were continued until a flip-flop event occurred in one of the replicas. The total time for a flip-flop event is simply five times the elapsed MD simulation time after dephasing. Altogether we observed 5 translocation events from the UL to the LL in a span of 2.9 µs. 5.4 Discussion In the direct all-atom MD simulation we observed 24 CHOL flip-flop events over a span of 15 µs, and in the parallel replica MD simulation there were 5 translocation events over 2.9 µs. Based on these results, the average CHOL translocation frequency in the simulated system is about 1.6 ×10 6 s -1 and the rate constant for flip-flop events in the first 10 µs is found to be 30000 s -1 , corresponding to a flip-flop half time of 23 µs. The CHOL flip-flop rate drops off significantly with time as shown in Fig. 5.8. 84 Figure 5.8: CHOL flip-flop rates measured in inverse microseconds versus time. The rates are calculated using block (red circles) and running (blue squares) averages. The block size is 3 µs. Continuous lines are drawn through the data points to guide the eye. The flip-flop rates decrease by an order of magnitude (red curve) with time. We have estimated the free-energy barrier for flip-flop events from the calculation of diffusion maps of CHOL molecules [127]. We find the free-energy barrier (21.5 KJ/mol) to be significantly lower than the potential energy barrier shown in Fig. 5.6. Our estimate for the free-energy barrier is consistent with the calculated flip-flop rate. Bennett et al. report the free-energy barrier to be 25 KJ/mol for 0% and 41 KJ/mol for 40% CHOL concentration. 85 The average CHOL flip-flop rate in our MD simulations is faster than the rate reported for other MD simulations of DPPC bilayers. In their 2009 paper, Bennett et al. report the CHOL flip-flop event rate constant to be between 10 and 8100 s -1 for DPPC bilayers with CHOL concentrations between 20 and 40% at 323 K [117]. These rate constants are based on umbrella sampling free energy calculations with all-atom MD simulations. In the same work, Bennett et al. also observe direct flip-flop events in their 3 µs long coarse-grained MD simulations. They find that the flip-flop rate constant exceeds 14000 s -1 when the CHOL concentration is 40%. Jo and coworkers find a CHOL flip-flop rate constant between 1200 and 16000 s -1 based on an umbrella sampling MD simulation of DPPC bilayers with two cholesterol molecules at 323 K [116]. Perhaps the reason for the faster flip-flop rate in our all-atom MD simulation is that initially the UL and LL had slightly different initial CHOL concentrations, 31% and 29%, respectively. Fast flip-flop rates have been reported in recent experimental work; Steck et al. find that flip-flop has a half time of < 1 s [115] while Bruckner et al. put an upper limit of 10 ms on the half time. These results, as well as other recent MD results, suggest that these fast flip-flop measurements are certainly plausible. They also call for experimental approaches with finer temporal resolution, capable of probing the time scales simulated here. The fast flip-flop rates we see here have implications for the physiological mechanism that maintains CHOL concentration asymmetry in the plasma membrane. In the context of CHOL molecules that can migrate across the membrane so quickly, it makes no sense to postulate active mechanisms (mediated either by transport proteins or steady-state membrane remodeling) that maintain the asymmetry. Instead, the asymmetric CHOL composition must be built into the chemical physics of the membrane, 86 and is most likely a result of interactions between CHOL and other lipid species. In particular, interactions between CHOL and sphingomyelin (SM) seem a likely source. It is well known that these molecules associate strongly [128], and SM flip-flops slowly, with a long residence time in a single bilayer leaflet [129]. A recent report by Bennett and Tieleman shows MD simulations in which the presence of SM in a “raft-like” membrane slows CHOL flip-flop tremendously, yielding a rate constant of 4.4 × 10 -4 s -1 [118]. From the observed direct events in our simulations we find that on average CHOL translocation time is about 73 ns. Since CHOL can exchange between the leaflets rapidly, CHOL flip-flop is a likely mechanism for a range of physiological processes, including coupling of phase structure across the membrane and stress relaxation upon membrane bending. Our results provide evidence of an induced phase transition upon CHOL flip-flop. Figure 5.5 shows clearly that as CHOL molecules migrate from the UL to the LL, the order parameter of the LL increases. At the temperature at which this simulation was performed, binary membranes of CHOL and DPPC can take on either an l o or l d configuration, depending on the concentration of CHOL. Chiang et al. have identified l o /l d phase coexistence at ~20% CHOL with a transition to pure l o phase occurring at about 30% CHOL concentration [130]. Our simulation, then, was performed around the transition between the phase coexistence regime and the pure l o regime. Chiang et al. also saw that the order parameter for the l o phase was about 0.1 greater than that for the l d phase, roughly consistent with the difference between the order parameters for the two leaflets at the 15 µs time point, at which the LL is 34% CHOL while the UL is 25% CHOL. The constant increase of the order parameter of the LL in the first ~7.5 µs is 87 consistent with ordering induced by the 22 CHOL molecules that flip-flop into this leaflet, taking it from 29 to 34% CHOL concentration. Rapid flip-flop of CHOL from a CHOL-rich leaflet to a CHOL-poor leaflet can clearly induce lipid ordering. Our results for the order parameter are in good agreement with previous MD calculations of order parameter at various CHOL concentrations. At 15 µs we find the CHOL concentration in the UL is 25 % and the order parameter is 0.28. Hofsaβ et al. [123] reported a value of 0.3 for the order parameter at the same CHOL concentration. In the LL we calculate the order parameter to be 0.33 at t =15 µs and the corresponding CHOL concentration to be 34 %. Hofsaβ et al. reported a value of 0.31 for S CD at 40 % CHOL concentration. (Note in order to compare our results with those of Hofsaβ et al. we averaged over atoms 2-8 in the acyl chains). The slight discrepancy between their results and ours is likely due to the difference in hydration levels. Our results also provide direct evidence for the role of relaxation of membrane stress upon CHOL flip-flop. We found the total stress profile of the membrane changes as CHOL migrates from the upper leaflet to the lower leaflet. The total compressive stress in the upper leaflet (Z > 0 nm) decreases significantly, which suggests that membrane bending in the direction of this leaflet would require significantly less energy following the transport of CHOL to the lower leaflet; see Fig. 5.9. Cholesterol transport rate slows as the membrane relaxes, (see Fig. 5.8), indicating that the initially observed fast transport rates may be driven in part by membrane stress. This all suggests that rapid CHOL flip-flop may be an important mechanism by which a stressed membrane is relaxed. 88 Figure 5.9: Stress, Π(z), profiles across the bilayer at time t = 5 µs (red) and 15 µs (blue). Given the experimental challenges in directly addressing the time scales of flip- flop events shown here, the best approaches to further investigating the biophysics of CHOL transport in membranes likely revolve around examining the secondary effects of flip-flop, particularly the ordering and stress relaxation effects described above. One promising model system for experimental investigation is the asymmetric SM vesicle system described by London and coworkers [129, 131]. Such systems allow for the dynamics of CHOL to be modulated in each leaflet independently via the strong SM- CHOL interaction. Combined with detailed simulations of this interaction and the implications of this interaction for membrane properties, the examination of how SM 89 asymmetry alters membrane properties has the potential to provide key evidence for the role of CHOL flip-flop and its dynamics in key physiological processes. In summary, our all-atom MD simulations of DPPC-CHOL membrane reveal that a slight imbalance (2%) in CHOL concentration between the UL and LL leads to spontaneous CHOL flip-flop events, revealing a relatively large rate constant of ~ 30000 s -1 . The average time for a CHOL molecule to migrate across the DPPC bilayer is 73 ns. At the end of our 15 µs MD simulation, we find that the stress distributions in the UL and LL are nearly symmetric and that their curvature energies are approximately the same. Hence we conclude that CHOL flip-flop events tend to relax stressed membranes. 90 6 Small Interfering RNA Induces Liquid-to-Ripple Phase Transformation in a Phospholipid Membrane 6.1 Introduction Since the discovery of RNAi substantial effort has gone into developing siRNA as a therapeutic technology [26]. RNAi using siRNA drugs can directly destroy mRNA transcribed from oncogenes, potentially arresting pathological cell proliferation. RNAi also has the capacity to knock out other points in biochemical pathways, potentially crippling proliferative mechanisms at arbitrary “weak points”. Strategies based on siRNA treatment have been investigated for gastric, colon, prostate, breast, lung, bladder, and ovarian cancer [27]. Despite a great deal of experimental research, the nature of mechanical barriers to siRNA transfection through a cell membrane is not well understood. The common hypothesis is that a bare siRNA molecule is unable to cross a cell membrane because of its large molecular weight (~ 13 kDa) and anionic charge [27]. Here we examine molecular processes that give rise to siRNA-induced transfection barriers in a model biomembrane consisting of DPPC molecules. Using all-atom MD simulations, we probe the structure, dynamics and mechanical response of a DPPC bilayer in the presence of a siRNA molecule at a temperature of 323 K and atmospheric pressure. Under these conditions, the thermodynamic state of a fully hydrated DPPC bilayer without the siRNA is the liquid crystalline phase 𝐿 ! [8, 132, 133]. In this phase DPPC molecules diffuse laterally and their hydrocarbon chains have conformational disorder in the form of gauche defects. Upon cooling at atmospheric pressure, the bilayer transforms into a ripple 91 phase (labeled L R ) below a temperature of 315 K [8]. The structure of the ripple phase has been the subject of several experimental and theoretical investigations [5, 132, 133]. X- ray studies indicate that the ripple phase in a DPPC bilayer consists of major (M) and minor (m) regions separated by kinks [2-4]. Experimental values for thicknesses of major and minor domains are consistent with the thickness of gel phase and liquid crystalline phase, respectively [5]. From AFM experiments, Leonenko et al. [134] find that for fully hydrated DPPC bilayers, the thickness of the gel phase is 4.4 ± 0.2 nm and that of liquid crystalline phase is 3.6 ± 0.3 nm [134]. X-ray studies [2-4, 135, 136] have been used to characterize other structural features of the ripple phase, e.g., ripple length (λ r ), stacking repeat distance (d) and oblique angle (γ). Fourier transform infrared and 13 C-NMR measurements show a high degree of tail stretching of DPPC molecules in the ripple phase [6, 7]. The structure and dynamics of lipid molecules in kink regions are not well understood [5], although experiments and MD simulations indicate that hydrocarbon chains of lipid molecules may be disordered and lipids may be more mobile in kinks than in major and minor regions of the ripple phase. On lowering the temperature of the ripple phase, the DPPC bilayer first goes into a gel phase (labeled Lβ ' ) at T = 308 K and then crystallizes into a subgel structure at 280 K [8]. In the gel phase hydrocarbon chains are nearly parallel and the membrane is thicker and less permeable than the fluid phase. The characteristics of 𝐿 ! → L R and L R → 𝐿 ! ! phase transformations in lecithin have been studied both experimentally and by computer simulations [2-7, 134-147]. 92 6.2 Simulation Method All-atom MD simulations were performed on a system consisting of a DPPC bilayer and a bare siRNA molecule. The double-stranded (ds) siRNA usually has 21-25 nucleotides on each strand and 2 nucleotide overhangs at the 3'- ends. In our simulations, we chose the siRNA sequence to be ss(UU)-ds(GACAGCAUAUAUGCUGUC)-ss(UU). Figure 6.1 shows a snapshot of the simulated system comprising a siRNA, 648 DPPC and 216,100 water molecules. Forty sodium ions were distributed randomly around the siRNA to maintain charge neutrality in the system [148, 149]. All in all, the MD cell contained 733,917 atoms and the initial dimensions of the cell were 14.4×14.4×38.6 nm 3 . The pre-equilibrated DPPC membrane was taken from ref. [150] and replicated in the x and y directions to generate a membrane of cross-sectional area 14.4×14.4 nm 2 . The MD simulation was carried out with CHARMM 36 force field [151-153] using the GROMACS MD package [44]. Periodic boundary conditions were imposed and we used the particle mesh Ewald method to compute the long-range Coulomb interaction [154]. The non-bonded Lennard-Jones interaction was calculated with a cut-off of 1.2 nm. The real-space contributions to Coulomb potential and forces were calculated with a cutoff of 1.0 nm. The Fourier space contributions to Coulomb potential and forces were calculated on a mesh of 0.15 nm. The simulation was carried out in the isothermal-isobaric ensemble using a Nosé-Hoover thermostat [155] to maintain the temperature at 323 K and Parrinello-Rahman barostat to keep the pressure at 1 bar [53]. Interatomic bond lengths were constrained using the LINCS algorithm [52]. Equations of motion were integrated with the velocity-Verlet algorithm using a time step of 1 fs and the simulation was run for 500 ns. 93 Figure 6.1: A snapshot of the system at time t = 10 ns. The siRNA is shown in magenta, sodium ions are yellow, and carbon tails of proximal and distal leaflets of the DPPC bilayer are cyan and orange spheres, respectively. Blue and red spheres represent lipid head groups. Only a few water molecules (red and white dots) are shown here. 6.3 Results The siRNA induces significant changes in the structure of the DPPC bilayer. The snapshot in Fig. 6.2(a) gives an overall view of the membrane from the top of the distal leaflet (farther from the siRNA) at time t = 500 ns. The head groups and tails of lipid molecules in the distal leaflet are blue and cyan, respectively, and the lipid tails of the proximal leaflet are brown. The apparent patchiness in the x-y plane of the membrane indicates significant variation in the membrane thickness. To quantify this, we divided the membrane into 200×200 pixels in the x-y plane and calculated the bilayer thickness in each pixel [156]. (The bilayer thickness, D, in a pixel is defined as the maximum distance between the lipid head groups in the proximal and distal leaflets.) Figure 6.2(b) shows the 94 Figure 6.2: (a) Top view of the lipid bilayer from the distal side of the siRNA. The head groups and tails of lipid molecules in the distal leaflet are blue and cyan, respectively, and the lipid tails of the proximal leaflet are brown. (b) Bilayer-thickness distributions for the fluid phase at time t = 50 ns (without siRNA, in red) and the ripple phase at t = 500 ns (in blue). (c) Pixelated membrane thickness shows the thicker gel domain (orange and red) embedded in a thinner gel phase (purple). 95 distribution functions for membrane thickness at time t = 2 ns and t = 500 ns. At 2 ns the distribution function has a single peak, indicating the membrane is uniform and its thickness is about 3.7 nm. As time evolves, the membrane develops a bi-modal thickness distribution and at t = 500 ns the distribution function has a prominent peak around 4.7 nm and a smaller one around 3.0 nm; see Fig. 6.2(b). Figure 6.2(c) is a color-coded snapshot of the membrane thickness in 200×200 pixels at t = 500 ns. Here the membrane thicknesses in the major (orange) and minor (purple) regions are D M ≈ 4.7 nm and D m ≈ 3.0 nm, respectively. In Fig. 6.3, panels (a) and (b) present side views of the bilayer structure in minor and major regions, respectively, at time t = 500 ns. Here the lipid head groups and hydrocarbon tails in the proximal leaflet of the DPPC bilayer are blue and cyan, and in the distal leaflet they are red and brown, respectively. Figures 6.3(a) and (b) show that the minor region consists of interdigitated and the major region contains non-interdigitated lipid molecules. These snapshots also show that lipid molecules in major and minor regions are tilted relative to the bilayer normal, i.e., the z-axis. After t = 100 ns the average tilt angle 𝜃 ! is around 26° for the patch of noninterdigitated lipid molecules. Experimental value of 𝜃 ! in the gel phase of a fully hydrated DPPC bilayer is 32° [157]. The major and minor regions make up almost entirely the ripple phase of the DPPC bilayer. Figure 6.4 is a snapshot of the ripple phase taken at t = 500 ns. We have used such snapshots to estimate the average wavelength (λ r ), amplitude (A), tilt angle (γ), and widths of the major (d M ) and minor (d m ) regions. We find: λ r = 10.85 nm, A = 1.07 nm, γ = 104°, d M = 3.33 nm, and d m = 1.41 nm. 96 Figure 6.3: A cross-section of the DPPC bilayer at time t = 500 ns shows interdigitated (a) and non-interdigitated (b) lipid molecules in the siRNA-induced ripple phase. In the top leaflet the head groups and lipid tails are blue and cyan and in the bottom leaflet they are red and orange, respectively. Figure 6.4: The ripple phase where X (horizontal) is the ripple vector direction at t = 500 ns. The cross-section is obtained for Y between 6 and 8 nm. 97 A noticeable feature of the major and minor regions is that the hydrocarbon tails of lipid molecules are almost straight. In contrast, the hydrocarbon chains of lipid molecules in the intervening kink between the major and minor regions are disordered. To quantify this disorder, we have calculated gauche defects in the lipid tails and find peaks in the distribution function around ±68º. These peaks are similar to the peaks in the distribution function for gauche defects in the 𝐿 ! phase. Moreover, the planar (x-y) self- diffusion coefficient of lipid molecules has nearly the same value in the kink region as in the 𝐿 ! phase. Together, the gauche defect density and lipid diffusion indicate that the kink region may be akin to the liquid crystalline phase of the DPPC bilayer. The siRNA molecule also has a dramatic slowing-down effect on the dynamics of lipid molecules. The MD simulation shows that lipid molecules are much less mobile in the ripple phase than in the liquid crystalline phase at 323 K. Taking a cue from supercooled liquids; we analyzed the planar dynamics of lipid molecules in terms of a cage correlation function, C(t)≡ Θ c−n i out (0,t) ( ) , (6.1) where c = 1, n i out 0,t ( ) = l i (0) 2 −l i (0)•l i (t). (6.2) l i (t) is the number of molecules in the neighbor list of the i th molecule at time t. The calculations were performed by using the coordinates of N atom of the headgroups of DPPC molecule. We used 200 frames over a span of 40 ns to calculate C(t). The calculations where done independently for the two leaflets. β is reported by fitting C(t) to the stretched exponential function. 98 Figure 6.5 shows the temporal decay of the cage correlation function in the ripple phase. We find that the stretched exponential relaxation, exp− 𝑡 𝜏 ! , with β = 0.2 gives the best fit to C(t). Figure 6.5: Stretched exponential behavior for the lipids. Here N atoms of the lower leaflet (which is the proximal leaflet to siRNA) are chosen for calculations. The calculations are done using 200 frames from 455 ns to 505 ns. In the major and minor regions of the ripple phase, lipid tails are straighter and more closely packed than in the 𝐿 ! phase. As a result, we expect that stresses in the 99 hydrophobic tails of the lipid bilayer will be more compressive in the L R than in the 𝐿 ! phase. To quantify this, we calculated the stress tensor 𝑃 [54]:  P= m i i ∑  v i ⊗  v i + 1 V  r ij i<j ∑ ⊗  F ij , (6.3) where m i is the mass of the i th atom, 𝑣 ! is its velocity, and 𝑟 !" and 𝐹 !" are the position vector and force between atoms i and j, respectively. Stress profiles are calculated by dividing the lipid-bilayer region into voxels of length 0.2 nm. At each instant of time, we use Eq. (6.3) to calculate the stress tensor in each voxel. The components of the stress tensor are then averaged over 5000 time frames spanning 50 ns. The z-dependence (i.e., along the bilayer normal) of stress profiles are calculated by averaging over all the voxels in the x-y plane for a given value of z. 100 Figure 6.6: Lateral stress profile normal to the lipid membrane (z-direction). Figure 6.6 shows the variation of stress π(z) = (P xx + P yy )/2 - P zz in the ripple phase. Here the positive peaks indicate membrane expansion in the lipid head-group region. The negative peaks in π(z) imply membrane compression in lipid tails. Evidently, the siRNA reduces the positive lateral pressure in the head-group region while enhancing the compressive stress in the lipid tails. We observe that the lateral pressure, (P xx + P yy )/2, develops a large compressive peak in the middle of the bilayer. 101 Figure 6.7: Snapshots of a SMD simulation show how a siRNA molecule migrates across the DPPC bilayer. Here the head groups of DPPC molecules are shown in red and blue spheres, respectively, and the corresponding tails groups are orange and cyan. (a) At time t = 0 ns in the SMD simulation, the siRNA molecule (shown in yellow) along with a few sodium ions (shown in magenta) lies flat on the lipid membrane; (b) at t = 22 ns, the siRNA is in the middle of the lipid bilayer; and (c) at t = 32 ns the siRNA molecules is aligned parallel to the bilayer normal. We also performed steered molecular dynamics (SMD) simulations to estimate the force required to pull a siRNA molecule across the DPPC bilayer. In SMD the siRNA center-of-mass was coupled to a spring and pulled with a constant velocity v normal to the bilayer plane, i.e., along the z direction. The value of the spring constant k was chosen to be 10 kJ mol -1 Å -2 , and SMD simulations were performed with pulling velocities 2.5, 5.0, and 10.0 Å/ns. The results for v = 2.5 Å/ns are shown in Fig. 6.7. As the siRNA approaches the bilayer, the electrostatic interaction between the anionic siRNA and the zwitterionic head groups of lipid molecules tilts the siRNA parallel to the bilayer plane; see Fig. 6.7(a). When the siRNA is pulled further into the membrane, lipid molecules in the proximal leaflet move away to open up a pore around the siRNA. Figure 6.7(b) shows the siRNA in the middle of the bilayer. The hydrophilic head groups (red spheres) of lipid molecules displaced from the proximal leaflet line up the pore surface and exert a torque 102 that rotates the siRNA 90°; see Fig. 6.7(c). Subsequently, the siRNA remains normal to the bilayer plane as it penetrates further into the membrane and moves toward the distal leaflet. The pore in the proximal leaflet begins to close up after the siRNA rotates in the bilayer and moves toward the distal leaflet. The bilayer heals back to the ripple phase after the siRNA goes through the distal leaflet. Figure 6.8 shows the force, F = – k (Z spring – Z siRNA ), required to pull the siRNA across the bilayer as a function of time for v = 2.5 Å/ns. Here Z siRNA is the z-coordinate of the center-of-mass of the siRNA and Z spring is the z-coordinate of a point on the spring. The maximum force required to pull a siRNA through the DPPC bilayer is about 1.5 nN. It is at the maximum value of the applied force that the siRNA completes its 90º rotation in the bilayer. Subsequently, the applied force decreases because it is easier to pull the siRNA when it is aligned normal to the bilayer plane. For higher pulling velocities, v = 5.0 and 10.0 Å/ns, the results are similar to those shown in Figs. 6.7 and 6.8. 103 Figure 6.8: The force versus time curve for SMD at v = 2.5 Å/ns. 6.4 Discussion Our all-atom MD simulation reveals that the electrostatic interaction between the anionic siRNA and the zwitterionic head groups of lipid molecules induces a liquid crystalline-to-ripple phase transformation in the DPPC bilayer at T = 323 K. It is known that on reducing the temperature below the chain-melting transition the liquid crystalline phase transforms to a gel phase in which the hydrocarbon chains of lipid molecules are fully stretched and tilted with respect to the bilayer normal. Upon cooling the DPPC 104 bilayer, a ripple phase is found between the liquid crystalline and gel phases. Wittebort et al. [7] have shown that the ripple phase exhibits characteristics of both liquid crystalline and gel phases. Ripple phase has been well-studied using differential scanning calorimetry [158], X-ray [138] and neutron diffractions [4], scanning-tunneling microscopy [137], and freeze-fracture electron microscopy [4, 138]. These experiments indicate that the size of the hydrated head group, the tilt of acyl chains, the extent of the chain movement, and the strength of intraleaflet interactions are the key determinants of the ripple phase [138]. Sengupta et al. [3] have used electron density maps of ripple phases to establish that the average chain tilt along the direction of rippling is responsible for the ripples. Other models assume that ripples form to relieve packing frustrations when the relationship between the cross-sectional areas of the head group and apolar tails exceeds a certain threshold [139]. It is also suggested that periodic local spontaneous curvature in the lipid bilayer gives rise to ripples. The reasons for the origin of this local spontaneous curvature could be electrostatic coupling between water molecules and the polar lipid headgroups, coupling between membrane curvature and molecular tilt, or the generation of curvature by linear arrays of fluid state lipid molecules [139]. Several simulation studies have also reported the formation of ripple phases. Lenz et al. [143] used Monte Carlo simulations of a coarse-grained molecular model for lipid- solvent mixtures to study ripple phases. They showed that a mismatch between the head group and tail size and the possibility of interdigitation were important factors for the formation of ripple phases. Kranenburg et al. [146] performed mesoscopic simulations to study the ripple phase in lipid bilayers. They found that the frustration arising from the 105 surface area of head groups and the lateral density of the tails was a key factor in the formation of ripple phases. Our results for the average wavelength (λ r = 10.85 nm), amplitude (A = 1.07 nm), tilt angle (γ = 104), and widths of major (d M = 3.33 nm) and minor (d m = 1.41) regions of the ripple phase are consistent with experimental and all-atom MD simulation studies. Yao et al. [2] find λ r = 13.64 nm, and γ = 95 o for the ripple structure in DPPC. They also find that the hydrocarbon chains in the primary ripple structure are tilted relative to the normal to the bilayer plane. Freeze-fracture electron microscopy and X-ray diffraction also suggest that λ r ~ 13 nm in the ripple phase of DPPC [2, 4]. The non-interdigitated and interdigitated major and minor regions of the ripple phase in the presence of a siRNA are akin to those found by de Vries et al. [5] in their all- atom MD simulations of lecithin bilayers. They observe that the interdigitated patch connects the neighboring noninterdigitated region such that the upper leaflet of the bilayer on one side crosses over into the lower leaflet on the other side. They find that λ r = 12 – 16 nm, A = 2.4 nm, γ = 56 – 118 o , d M = 4.6 nm, and d m = 3.2 nm for the ripple phase in the DPPC membrane. In our simulation, we find the ripple phase consists of a 3.0-nm-thick minor region containing interdigitated lipid molecules and a 4.7-nm-thick major region containing noninterdigitated lipid molecules. There is experimental evidence for interdigitation of lipid molecules in the presence of anions: Lesslaue et al. [159] report that the anionic fluorescent probe 1-anilino-8-naphthalene sulfonate decreases the lipid bilayer thickness to 3 nm because of interdigitation of lipid chains; and X-ray diffraction data show that chlorpromazine also reduces the bilayer thickness from 5 nm to 3 nm 106 [160]. Wang et al.’s experiments indicate that anionic nanoparticles induce a transformation from a liquid crystalline to gel phase in pure phospholipid vesicles [161]. As far the dynamics of lipid molecules are concerned, the presence of siRNA significantly reduces the in-plane diffusivity of lipids in the ripple phase compared with diffusion in the liquid crystalline phase without the siRNA. The DPPC molecules in the presence of a siRNA display stretched exponential relaxation and their cage correlation function decays with time as exp− 𝑡 𝜏 ! , where β = 0.2. Brandt et al. [162] performed molecular dynamics simulations of 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) lipid bilayer over a broad range of length and time scales to show that the time correlation functions of structure fluctuations in DMPC lipid bilayer are stretched exponentials. They find that β = 0.45 for the DMPC lipids using the Berger united atom force field. Their simulations were performed at T = 300 K and atmospheric pressure. The barrier to siRNA transfection across the DPPC bilayer arises from the presence of large compressive stresses in the middle of the bilayer. Our steered MD simulations reveal that the maximum force required for the siRNA passage across the lipid ripple phase is about 1.5 nN. This should be compared with electroporation studies of siRNA using MD simulations. Tarek et al. [29, 30] applied high electric fields (several orders more than that used in experiments) to cause electroporation on the nanosecond timescale. They find that the electroporation is accompanied by the formation of water wires and channels across the membrane. These channels are stabilized by hydrophilic pores formed by lipid head groups and acyl chains. They also infer that electric field induces a significant lateral stress on the bilayer. Switching off the external transmembrane potential allows for complete resealing and reconstruction of the bilayer. 107 In an experimental electroporation study, Breton et al. [29] indicate that a single 10 ns high-voltage electric pulse can permeabilize lipid vesicles and allow delivery of siRNA to the cytoplasm. Also, in their MD simulation they apply an electric field of 140 kV/mm across the bilayer. With a -40e charge on the siRNA, this amounts to a constant force of ~ 0.9 nN on the siRNA and it takes about 10 ns for siRNA to go across the membrane. They show that low-magnitude nanopulses drag the siRNA towards the membrane where it interacts strongly with the polar head groups of lipid molecules. When the field strength is increased, small nanopores open up in the lipid membrane and the siRNA slides through these pores within nanoseconds. We have observed similar behavior in our constant force SMD simulation. 6.5 Conclusion In summary, our all-atom MD simulation shows that the siRNA induces a liquid crystalline-to-ripple phase transformation in the bilayer. The ripple phase consists of a major region of non-interdigitated and a minor region of interdigitated lipid molecules with an intervening kink. Slow dynamics of DPPC molecules in the ripple phase are analyzed in terms of a cage correlation function which exhibits stretched exponential relaxation. In the ripple phase hydrocarbon chains of lipid molecules have large compressive stresses, which present a considerable barrier to siRNA transfection across the bilayer. Steered MD simulations reveal that siRNA transfection through the DPPC bilayer requires a force of ~ 1.5 nN. SiRNA transfection can be facilitated by the fluidization of the cell membrane, which will lower the mechanical barrier arising from large compressive stresses in the ripple phase. Experiments by Wang et al. indicate that cationic nanoparticles can fluidize 108 the gel phase of lipid vesicles [161], and numerous experiments have shown that encapsulation in cationic liposomes enhances the transfection efficiency of siRNA molecules [27]. 109 7 Conclusions We have performed large scale and long time all-atom MD simulations to study: (1) poration of lipid bilayers by shock-induced nanobubble collapse, (2) cholesterol translocation in a phospholipid membrane and (3) small interfering RNA (siRNA) induced phase transformation in a phospholipid membrane. Our largest simulation involved 30 million atoms and our longest simulation reached 15 microseconds. We studied shock-induced structural changes in water using a reactive force field MD simulation approach that takes into account the formation and dissociation of chemical bonds. The observed Hugoniot curve is in good agreement with experimental data. We find molecular clusters with ice-VII-like structure at particle velocity of 1 km/s. This is well supported by experimental observations of ice VII for particle velocities between 0.75 and 2 km/s. The dynamics of shock-induced collapse of nanobubbles in water has been studied in detail. When a shock wave strikes a nanobubble, the nanobubble starts shrinking and we observe a focused nanojet whose length scales linearly with the nanobubble radius. This scaling relation has also been observed experimentally in shock-induced collapse of millimeter-to-micron sized bubbles. We find that the size of the nanojet increases with an increase in the nanobubble diameter and shock amplitude. When the nanojet hits the distal side of the nanobubble, a secondary water hammer shock is generated. This secondary shock propagates spherically backward with a speed higher than that of the primary shock. At the onset, the water hammer shock has a pressure larger than the primary shock, but the pressure drops rapidly as the secondary shock spreads. As a potential application, we studied poration in a lipid bilayer due to shock- 110 induced nanobubble collapse. Our multimillion-atom MD simulations reveal the mechanism of transient poration in lipid bilayers by shock-induced collapse of nanobubbles. The MD simulations indicate that the nanojet impact significantly deforms and thins the lipid bilayer and water molecules in the nanojet form a spreading flow pattern after impact. Deformation and thinning of bilayers combined with large pressure gradients across and spreading flow around the bilayers create transient nanochannels through which water molecules translocate across the bilayer. We have also examined the flip-flop dynamics of cholesterol (CHOL) molecules in a DPPC-CHOL bilayer. CHOL molecules play a key role in modulating the rigidity of cell membrane and in controlling intracellular transport and signal transduction. Using an all-atom MD approach, we study the process of CHOL interleaflet transport (flip-flop) in a DPPC–CHOL bilayer over a time period of 15 µs. Our simulations reveal that a slight imbalance (2%) in CHOL concentration between the leaflets leads to spontaneous CHOL flip-flop events, revealing a relatively large rate constant of ~ 30000 s -1 . The rapid flip- flop rates reported here have important implications for the role of CHOL in mechanical properties of cell membranes, formation of domains, and maintaining CHOL concentration asymmetry in plasma membrane. The average time for a CHOL molecule to migrate across the DPPC bilayer is 73 ns. At the end of our 15 µs MD simulation, we find that the stress distributions in the leaflets are nearly symmetric and that their curvature energies are approximately the same. Hence we conclude that CHOL flip-flop events tend to relax stressed membranes. These results have been verified by performing various independent parallel replica (PR) simulations. Our PR approach can reach submillisecond time scales and bridge the gap between MD simulations and Nuclear 111 Magnetic Resonance (NMR) experiments on CHOL flip-flop dynamics in membranes. Lastly, we examined siRNA transfection barriers in a DPPC bilayer. SiRNA molecule plays a pivotal role in gene silencing via the RNAi mechanism. A key limitation to the widespread implementation of siRNA therapeutics is the difficulty of delivering siRNA-based drugs to cells. We have performed all-atom MD simulation to examine changes in the structure, dynamics and mechanical response of a DPPC bilayer in the presence of a siRNA molecule. We find that the siRNA induces a liquid crystalline-to-ripple phase transformation in the bilayer. The ripple phase consists of a major region of non-interdigitated and a minor region of interdigitated lipid molecules with an intervening kink. Slow dynamics of DPPC molecules in the ripple phase are analyzed in terms of a cage correlation function which exhibits stretched exponential relaxation. In the ripple phase hydrocarbon chains of lipid molecules have large compressive stresses, which present a considerable barrier to siRNA transfection across the bilayer. Steered MD simulations reveal that siRNA transfection through the DPPC bilayer requires a force of ~ 1.5 nN. 112 References [1] V. A. Raghunathan and J. Katsaras, Phys Rev Lett 74, 4456 (1995). [2] H. Yao, S. Matuoka, B. Tenchov, and I. Hatta, Biophys J 59, 252 (1991). [3] K. Sengupta, V. A. Raghunathan, and J. Katsaras, Europhysics Letters 49 (2000). [4] Beth A. Cunningham, Ari-‐David Brown, David H. Wolfe, W. Patrick Williams, and A. Brain, Physical review. E, Statistical, nonlinear, and soft matter physics 58 (1998). [5] A. H. de Vries, S. Yefimov, A. E. Mark, and S. J. Marrink, P Natl Acad Sci USA 102, 5392 (2005). [6] D. G. Cameron, H. L. Casal, and H. H. Mantsch, Biochemistry 19, 3665 (1980). [7] R. J. Wittebort, C. F. Schmidt, and R. G. Griffin, Biochemistry 20, 4223 (1981). [8] S. Tristram-‐Nagle and J. F. Nagle, Chem Phys Lipids 127, 3 (2004). [9] M. B. Sankaram and T. E. Thompson, P Natl Acad Sci USA 88, 8686 (1991). [10] K. Simons and D. Toomre, Nat Rev Mol Cell Bio 1, 31 (2000). [11] C. Dietrich, Z. N. Volovyk, M. Levi, N. L. Thompson, and K. Jacobson, P Natl Acad Sci USA 98, 10642 (2001). [12] T. Baumgart, S. T. Hess, and W. W. Webb, Nature 425, 821 (2003). [13] N. Kahya, D. Scherfeld, K. Bacia, B. Poolman, and P. Schwille, J Biol Chem 278, 28109 (2003). [14] A. V. Samsonov, I. Mihalyov, and F. S. Cohen, Biophys J 81, 1486 (2001). [15] J. Zhao, J. Wu, F. A. Heberle, T. T. Mills, P. Klawitter, G. Huang, G. Costanza, and G. W. Feigenson, Bba-‐Biomembranes 1768, 2764 (2007). [16] H. J. Risselada and S. J. Marrink, P Natl Acad Sci USA 105, 17367 (2008). [17] M. Nakano, M. Fukuda, T. Kudo, N. Matsuzaki, T. Azuma, K. Sekine, H. Endo, and T. Handat, Journal of Physical Chemistry B 113, 6745 (2009). [18] J. A. Hamilton, Curr Opin Lipidol 14, 263 (2003). [19] S. Garg, L. Porcar, A. C. Woodka, P. D. Butler, and U. Perez-‐Salas, Biophys J 101, 370 (2011). [20] R. J. Bruckner, S. S. Mansy, A. Ricardo, L. Mahadevan, and J. W. Szostak, Biophys J 97, 3113 (2009). [21] A. J. Garcia-‐Saez and P. Schwille, Biochim Biophys Acta 1798, 766 (2010). [22] G. Lindblom and G. Oradd, Biochim Biophys Acta 1788, 234 (2009). [23] J. Perez-‐Gil, Biochim Biophys Acta 1798, 701 (2010). [24] J. P. Slotte, Biochim Biophys Acta 1788, 1 (2009). [25] E. Lindahl and O. Edholm, Biophys J 79, 426 (2000). [26] A. Fire, S. Q. Xu, M. K. Montgomery, S. A. Kostas, S. E. Driver, and C. C. Mello, Nature 391, 806 (1998). [27] K. A. Whitehead, R. Langer, and D. G. Anderson, Nat Rev Drug Discov 8, 516 (2009). [28] M. Parisien, X. Wang, G. Perdrizet, 2nd, C. Lamphear, C. A. Fierke, K. C. Maheshwari, M. J. Wilde, T. R. Sosnick, and T. Pan, Cell Rep 3, 1703 (2013). [29] M. Breton, L. Delemotte, A. Silve, L. M. Mir, and M. Tarek, J Am Chem Soc 134, 13938 (2012). 113 [30] M. Tarek, Biophys J 88, 4045 (2005). [31] C. D. Ohl, M. Arora, R. Ikink, N. de Jong, M. Versluis, M. Delius, and D. Lohse, Biophys J 91, 4285 (2006). [32] V. P. Torchilin, Nat Rev Drug Discov 4, 145 (2005). [33] K. K. Ewert, H. M. Evans, N. F. Bouxsein, and C. R. Safinya, Bioconjugate Chem 17, 877 (2006). [34] A. Ahmad, H. M. Evans, K. Ewert, C. X. George, C. E. Samuel, and C. R. Safinya, Journal of Gene Medicine 7, 739 (2005). [35] S. Bamrungsap, Z. L. Zhao, T. Chen, L. Wang, C. M. Li, T. Fu, and W. H. Tan, Nanomedicine-‐Uk 7, 1253 (2012). [36] W. Y. Qian, D. M. Sun, R. R. Zhu, X. L. Du, H. Liu, and S. L. Wang, Int J Nanomed 7, 5781 (2012). [37] S. J. T. Rezaei, M. R. Nabid, H. Niknejad, and A. A. Entezami, Polymer 53, 3485 (2012). [38] H. Wu, L. Zhu, and V. P. Torchilin, Biomaterials 34, 1213 (2013). [39] S. Cajot, K. Van Butsele, A. Paillard, C. Passirani, E. Garcion, J. P. Benoit, S. K. Varshney, and C. Jerome, Acta Biomater 8, 4215 (2012). [40] S. Donnini, F. Tegeler, G. Groenhof, and H. Grubmuller, J Chem Theory Comput 7, 1962 (2011). [41] M. Vedadi, A. Choubey, K. Nomura, R. K. Kalia, A. Nakano, P. Vashishta, and A. C. T. van Duin, Phys Rev Lett 105 (2010). [42] A. Choubey, M. Vedadi, K. Nomura, R. K. Kalia, A. Nakano, and P. Vashishta, Appl Phys Lett 98 (2011). [43] A. Choubey, R. K. Kalia, N. Malmstadt, A. Nakano, and P. Vashishta, Biophys J 104, 2429 (2013). [44] B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J Chem Theory Comput 4, 435 (2008). [45] T. Darden, D. York, and L. Pedersen, J Chem Phys 98, 10089 (1993). [46] D. P. Tieleman, J. L. Maccallum, W. L. Ash, C. Kandt, Z. Xu, and L. Monticelli, J Phys Condens Matter 18, S1221 (2006). [47] J. B. Klauda, R. M. Venable, J. A. Freites, J. W. O'Connor, D. J. Tobias, C. Mondragon-‐Ramirez, I. Vorobyov, A. D. MacKerell, Jr., and R. W. Pastor, J Phys Chem B 114, 7830 (2010). [48] L. Monticelli, E. Salonen, and D. P. Tieleman, Biophys J 98, 668A (2010). [49] O. Berger, O. Edholm, and F. Jahnig, Biophys J 72, 2002 (1997). [50] S. J. Marrink, A. H. de Vries, and A. E. Mark, Journal of Physical Chemistry B 108, 750 (2004). [51] D. P. Tieleman, J. L. MacCallum, W. L. Ash, C. Kandt, Z. T. Xu, and L. Monticelli, J Phys-‐Condens Mat 18, S1221 (2006). [52] B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije, J Comput Chem 18, 1463 (1997). [53] S. Nose and M. L. Klein, Mol Phys 50, 1055 (1983). [54] O. H. S. Ollila, H. J. Risselada, M. Louhivuori, E. Lindahl, I. Vattulainen, and S. J. Marrink, Phys Rev Lett 102 (2009). [55] C. D. Ohl and R. Ikink, Phys Rev Lett 90, 214502 (2003). 114 [56] E. A. Brujan, G. S. Keen, A. Vogel, and J. R. Blake, Physics of Fluids 14, 85 (2002). [57] Y. Tomita and T. Kodama, Journal of Applied Physics 94, 2809 (2003). [58] G. N. Sankin and P. Zhong, Physical Review E 74, 046304 (2006). [59] A. Tufaile and J. C. Sartorelli, Physical Review E 66, 056204 (2002). [60] E. Zwaan, S. Le Gac, K. Tsuji, and C. D. Ohl, Phys Rev Lett 98, 254051 (2007). [61] V. S. Ajaev and G. M. Homsy, Annual Review of Fluid Mechanics 38, 277 (2006). [62] C. E. Brennen, Cavitation and Bubble Dynamics (Oxford University Press, New York, 1995). [63] T. Kodama and K. Takayama, Ultrasound in Medicine and Biology 24, 723 (1998). [64] T. Kodama, M. R. Hamblin, and A. G. Doukas, Biophysical Journal 79, 1821 (2000). [65] C. D. Ohl, M. Arora, R. Dijkink, V. Janve, and D. Lohse, Applied Physics Letters 89, 074102 (2006). [66] M. Delius and G. Adams, Cancer Research 59, 5227 (1999). [67] T. Kodama, Y. Tomita, Y. Watanabe, K. Koshiyama, T. Yano, and S. Fujikawa, Journal of Biomechanical Science and Engineering 4, 124 (2009). [68] A. C. T. van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard, Journal of Physical Chemistry A 105, 9396 (2001). [69] A. K. Rappe and W. A. Goddard, Journal of Physical Chemistry 95, 3358 (1991). [70] A. Nakano et al., Computational Materials Science 38, 642 (2007). [71] K. I. Nomura, R. K. Kalia, A. Nakano, and P. Vashishta, Computer Physics Communications 178, 73 (2008). [72] A. P. Rybakov and I. A. Rybakov, Eur J Mech B-‐Fluid 14, 323 (1995). [73] J. D. Jorgensen and T. G. Worlton, J Chem Phys 83, 329 (1985). [74] S. Klotz, L. E. Bove, T. Strassle, T. C. Hansen, and A. M. Saitta, Nature Materials 8, 405 (2009). [75] M. Moseler and U. Landman, Science 289, 1165 (2000). [76] B. M. Borkent, S. M. Dammer, H. Schonherr, G. J. Vancso, and D. Lohse, Phys Rev Lett 98, 204502 (2007). [77] T. Kodama and Y. Tomita, Applied Physics B-‐Lasers and Optics 70, 139 (2000). [78] E. Johnsen and T. Colonius, Journal of the Acoustical Society of America 124, 2011 (2008). [79] A. Golberg and B. Rubinsky, Biomedical Engineering Online 9, 13 (2010). [80] P. T. Vernier, Y. H. Sun, and M. A. Gundersen, Bmc Cell Biology 7, 37 (2006). [81] S. Mitragotri, Nature Reviews Drug Discovery 4, 255 (2005). [82] I. Rosenthal, J. Z. Sostaric, and P. Riesz, Ultrasonics Sonochemistry 11, 349 (2004). [83] C. M. H. Newman and T. Bettinger, Gene Therapy 14, 465 (2007). [84] M. Tamagawa, I. Yamanoi, N. Ishimatsu, and S. Suetsugu, World Congress on Medical Physics and Biomedical Engineering 2006, Vol 14, Pts 1-‐6 14, 3236 (2007). [85] N. Kudo, K. Okada, and K. Yamamoto, Biophysical Journal 96, 4866 (2009). 115 [86] Y. Zhou, J. M. Cui, and C. X. Deng, Biophysical Journal 94, L51 (2008). [87] A. G. Doukas and N. Kollias, Advanced Drug Delivery Reviews 56, 559 (2004). [88] D. P. Tieleman. and H. J. C. Berendsen., J Chem Phys 105, 4871 (1996). [89] S. Leekumjorn and A. K. Sum, Biochim Biophys Acta 1768, 354 (2007). [90] A. L. Kuo and C. G. Wade, Biochemistry 18, 2300 (1979). [91] B. S. Lee, S. A. Mabry, A. Jonas, and J. Jonas, Chem Phys Lipids 78, 103 (1995). [92] S. J. Marrink, J. Risselada, and A. E. Mark, Chem Phys Lipids 135, 223 (2005). [93] K. I. Nomura, R. K. Kalia, A. Nakano, P. Vashishta, A. C. T. van Duin, and W. A. Goddard, Phys Rev Lett 99 (2007). [94] C. D. Ohl and R. Ikink, Phys Rev Lett 90 (2003). [95] C. D. Ohl, M. Arora, R. Dijkink, V. Janve, and D. Lohse, Appl Phys Lett 89 (2006). [96] S. J. Singer and G. L. Nicolson, Science 175, 720 (1972). [97] M. P. Sheetz and S. J. Singer, P Natl Acad Sci USA 71, 4457 (1974). [98] Y. Lange, J. Dolde, and T. L. Steck, J Biol Chem 256, 5321 (1981). [99] Y. Lange, H. B. Cutler, and T. L. Steck, J Biol Chem 255, 9331 (1980). [100] E. Evans and W. Rawicz, Phys Rev Lett 64, 2094 (1990). [101] J. B. Song and R. E. Waugh, Biophys J 64, 1967 (1993). [102] J. Henriksen, A. C. Rowat, E. Brief, Y. W. Hsueh, J. L. Thewalt, M. J. Zuckermann, and J. H. Ipsen, Biophys J 90, 1639 (2006). [103] J. Henriksen, A. C. Rowat, and J. H. Ipsen, Eur Biophys J Biophy 33, 732 (2004). [104] M. Yanagisawa, M. Imai, and T. Taniguchi, Phys Rev Lett 100, 148102 (2008). [105] M. Bretsche, Nature-‐New Biology 236, 11 (1972). [106] P. F. DEVAUX, Biochemistry 30, 1163 (1991). [107] T. Y. Wang and J. R. Silvius, Biophys J 81, 2762 (2001). [108] B. Baird, E. D. Sheets, and D. Holowka, Biophysical Chemistry 82, 109 (1999). [109] V. Kiessling, J. M. Crane, and L. K. Tamm, Biophys J 91, 3313 (2006). [110] M. D. Collins and S. L. Keller, Proc Natl Acad Sci USA 105, 124 (2008). [111] V. Kiessling, C. Wan, and L. K. Tamm, Bba-‐Biomembranes 1788, 64 (2009). [112] C. Wan, V. Kiessling, and L. K. Tamm, Biochemistry 47, 2190 (2008). [113] M. D. Collins, Biophys J 94, L32 (2008). [114] Y. Lange, C. M. Cohen, and M. J. Poznansky, P Natl Acad Sci USA 74, 1538 (1977). [115] T. L. Steck, J. Ye, and Y. Lange, Biophys J 83, 2118 (2002). [116] S. Jo, H. A. Rui, J. B. Lim, J. B. Klauda, and W. Im, Journal of Physical Chemistry B 114, 13342 (2010). [117] W. F. D. Bennett, J. L. MacCallum, M. J. Hinner, S. J. Marrink, and D. P. Tieleman, J Am Chem Soc 131, 12714 (2009). [118] W. F. D. Bennett and D. P. Tieleman, J Lipid Res 53, 421 (2012). [119] A. F. Voter, Phys Rev B 57, 13985 (1998). [120] E. Cubuk, A. Waterland, and E. Kaxiras, American Physical Society, APS March Meeting, S1.189 (2012). [121] M. Höltje., T. Förster., B. Brandt., T. Engels., W. v. Rybinski., and H.-‐D. Höltje., Biochimica et Biophysica Acta (BBA) -‐ Biomembranes 1511, 156 (2001). 116 [122] S. A. Pandit, S. W. Chiu, E. Jakobsson, A. Grama, and H. L. Scott, Langmuir 24, 6858 (2008). [123] C. Hofsass, E. Lindahl, and O. Edholm, Biophys J 84, 2192 (2003). [124] M. Lafleur, P. R. Cullis, and M. Bloom, Eur Biophys J 19, 55 (1990). [125] S. W. Chiu, E. Jakobsson, R. J. Mashl, and H. L. Scott, Biophys J 83, 1842 (2002). [126] R. Parthasarathy and J. T. Groves, Soft Matter 3, 24 (2007). [127] A. L. Ferguson, A. Z. Panagiotopoulos, I. G. Kevrekidis, and P. G. Debenedetti, Chem Phys Lett 509, 1 (2011). [128] M. P. Veiga, J. L. R. Arrondo, F. M. Goni, A. Alonso, and D. Marsh, Biochemistry 40, 2614 (2001). [129] H. T. Cheng, Megha, and E. London, J Biol Chem 284, 6079 (2009). [130] Y. W. Chiang, A. J. Costa-‐Filho, and J. H. Freed, Journal of Physical Chemistry B 111, 11260 (2007). [131] S. Chiantia, P. Schwille, A. S. Klymchenko, and E. London, Biophys J 100, L1 (2011). [132] W. J. Sun, S. TristramNagle, R. M. Suter, and J. F. Nagle, P Natl Acad Sci USA 93, 7008 (1996). [133] M. Rappolt, G. Pabst, G. Rapp, M. Kriechbaum, H. Amenitsch, C. Krenn, S. Bernstorff, and P. Laggner, Eur Biophys J Biophy 29, 125 (2000). [134] Z. V. Leonenko, E. Finot, H. Ma, T. E. Dahms, and D. T. Cramb, Biophys J 86, 3783 (2004). [135] D. V. Soloviov, Y. E. Gorshkova, O. I. Ivankov, A. N. Zhigunov, L. A. Bulavin, V. I. Gordeliy, and A. I. Kuklin, J Phys Conf Ser 351 (2012). [136] X. Y. Wang, K. Semmler, W. Richter, and P. J. Quinn, Arch Biochem Biophys 377, 304 (2000). [137] J. T. Woodward and J. A. Zasadzinski, Biophys J 72, 964 (1997). [138] D. M. Czajkowsky, C. Huang, and Z. F. Shao, Biochemistry 34, 12501 (1995). [139] T. Kaasgaard, C. Leidy, J. H. Crowe, O. G. Mouritsen, and K. Jorgensen, Biophys J 85, 350 (2003). [140] C. Cametti, F. Deluca, A. Dilario, M. A. Macri, G. Briganti, and B. Maraviglia, Prog Coll Pol Sci S 84, 465 (1991). [141] K. Honda and H. Kimura, J Phys Soc Jpn 60, 1212 (1991). [142] M. L. Belaya, M. V. Feigelman, and V. G. Levadny, J Phys Ii 1, 375 (1991). [143] O. Lenz and F. Schmid, Phys Rev Lett 98 (2007). [144] W. S. Mccullough and H. L. Scott, Biophys J 57, A473 (1990). [145] P. A. Pearce and H. L. Scott, J Chem Phys 77, 951 (1982). [146] M. Kranenburg, C. Laforge, and B. Smit, Phys Chem Chem Phys 6, 4531 (2004). [147] M. A. Kamal, A. Pal, V. A. Raghunathan, and M. Rao, Epl-‐Europhys Lett 95 (2011). [148] P. K. Maiti, T. A. Pascal, N. Vaidehi, J. Heo, and W. A. Goddard, Biophys J 90, 1463 (2006). [149] S. A. Pandit and M. L. Berkowitz, Biophys J 82, 1818 (2002). [150] http://terpconnect.umd.edu/~jbklauda/research/download.html 117 [151] P. Bjelkmar, P. Larsson, M. A. Cuendet, B. Hess, and E. Lindahl, J Chem Theory Comput 6, 459 (2010). [152] S. E. Feller and A. D. MacKerell, Journal of Physical Chemistry B 104, 7510 (2000). [153] A. D. MacKerell and N. K. Banavali, J Comput Chem 21, 105 (2000). [154] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J Chem Phys 103, 8577 (1995). [155] W. G. Hoover, Phys Rev A 31, 1695 (1985). [156] W. J. Allen, J. A. Lemkul, and D. R. Bevan, J Comput Chem 30, 1952 (2009). [157] J. F. Nagle and S. Tristram-‐Nagle, Bba-‐Rev Biomembranes 1469, 159 (2000). [158] R. N. A. H. Lewis, N. Mak, and R. N. Mcelhaney, Biochemistry-‐Us 26, 6118 (1987). [159] Lesslaue.W, J. E. Cain, and J. K. Blasie, P Natl Acad Sci USA 69, 1499 (1972). [160] T. J. Mcintosh, R. V. Mcdaniel, and S. A. Simon, Biochim Biophys Acta 731, 109 (1983). [161] B. Wang, L. F. Zhang, S. C. Bae, and S. Granick, P Natl Acad Sci USA 105, 18171 (2008). [162] E. G. Brandt and O. Edholm, J Chem Phys 133, 115101 (2010).

Abstract (if available)

Abstract Biological cell membranes provide mechanical stability to cells and understanding their structure, dynamics and mechanics are important biophysics problems. Experiments coupled with computational methods such as molecular dynamics (MD) have provided insight into the physics of membranes. We use long‐time and large‐scale MD simulations to study the structure, dynamics and mechanical behavior of membranes. ❧ We investigate shock‐induced collapse of nanobubbles in water using MD simulations based on a reactive force field. We observe a focused jet at the onset of bubble shrinkage and a secondary shock wave upon bubble collapse. The jet length scales linearly with the nanobubble radius, as observed in experiments on micron‐to‐millimeter size bubbles. Shock induces dramatic structural changes, including an ice‐VII‐like structural motif at a particle velocity of 1 km/s. The incipient ice VII formation and the calculated Hugoniot curve are in good agreement with experimental results. We also investigate molecular mechanisms of poration in lipid bilayers due to shock‐induced collapse of nanobubbles. Our multimillion‐atom MD simulations reveal that the jet impact generates shear flow of water on bilayer leaflets and pressure gradients across them. This transiently enhances the bilayer permeability by creating nanopores through which water molecules translocate rapidly across the bilayer. Effects of nanobubble size and temperature on the porosity of lipid bilayers are examined. ❧ The second research project focuses on cholesterol (CHOL) dynamics in phospholipid bilayers. Several experimental and computational studies have been performed on lipid bilayers consisting of dipalmitoylphosphatidylcholine (DPPC) and CHOL molecules. CHOL interleaflet transport (flip-flop) plays an important role in interleaflet coupling and determining CHOL flip‐flop rate has been elusive. Various studies report that the rate ranges between milliseconds to seconds. We calculate CHOL flip‐flop rates by performing a 15 μs all‐atom MD simulation of a DPPC–CHOL bilayer. We find that the CHOL flip‐flop rates are on the sub microsecond timescale. These results are verified by performing various independent parallel replica (PR) simulations. Our PR simulations provide significant boost in sampling of the flip‐flop events. We observe that the CHOL flip‐flop can induce membrane order, regulate membrane‐bending energy, and facilitate membrane relaxation. The rapid flip‐flop rates reported here have important implications for the role of CHOL in mechanical properties of cell membranes, formation of domains, and maintaining CHOL concentration asymmetry in plasma membrane. Our PR approach can reach submillisecond time scales and bridge the gap between MD simulations and Nuclear Magnetic Resonance (NMR) experiments on CHOL flip‐flop dynamics in membranes. ❧ The last project deals with transfection barriers encountered by a bare small interfering RNA (siRNA) in a phospholipid bilayer. SiRNA molecules play a pivotal role in therapeutic applications. A key limitation to the widespread implementation of siRNA‐based therapeutics is the difficulty of delivering siRNA‐based drugs to cells. We have examined structural and mechanical barriers to siRNA passage across a phospholipid bilayer using all‐atom MD simulations. We find that the electrostatic interaction between the anionic siRNA and head groups of phospholipid molecules induces a phase transformation from the liquid crystalline to ripple phase. Steered MD simulations reveal that the siRNA transfection through the ripple phase requires a force of ~ 1.5 nN.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Shock-induced nanobubble collapse and its applications

PDF

Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations

PDF

Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations

$Multimillion-to-billion atom molecular dynamics simulations of deformation, damage, nanoindentation, and fracture in silica glass and energetic materials$

PDF

Multimillion-to-billion atom molecular dynamics simulations of deformation, damage, nanoindentation, and fracture in silica glass and energetic materials

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras

PDF

Large-scale molecular dynamics simulations of nano-structured materials

PDF

High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks

PDF

Molecular dynamics study of energetic materials

Asset Metadata

Creator Choubey, Amit (author)

Core Title Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 02/18/2014

Defense Date 02/05/2014

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag cholesterol flip-flop,DPPC bilayer,molecular dynamics,nanobubble collapse,OAI-PMH Harvest,poration,shock,siRNA

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Kalia, Rajiv K. (committee chair), El-Naggar, Mohamed Y. (committee member), Malmstadt, Noah (committee member), Nakano, Aiichiro (committee member), Vashishta, Priya (committee member)

Creator Email choubey@usc.edu,kgp.amit@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-363504

Unique identifier UC11296412

Identifier etd-ChoubeyAmi-2248.pdf (filename),usctheses-c3-363504 (legacy record id)

Legacy Identifier etd-ChoubeyAmi-2248.pdf

Dmrecord 363504

Document Type Dissertation

Format application/pdf (imt)

Rights Choubey, Amit

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA