Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Molecular dynamics simulations of lipid bilayers in megavolt per meter electric fields
(USC Thesis Other)
Molecular dynamics simulations of lipid bilayers in megavolt per meter electric fields
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
MOLECULAR DYNAMICS SIMULATIONS OF LIPID BILAYERS IN MEGAVOLT
PER METER ELECTRIC FIELDS
by
Matthew James Ziegler
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(CHEMICAL ENGINEERING)
December 2008
Copyright 2008 Matthew James Ziegler
ii
DEDICATION
To Leslie.
iii
ACKNOWLEDGEMENTS
It is a pleasure to thank all the people who made this thesis possible.
First and foremost I thank my adviser Dr. Vernier. Starting with only an idea and a
computer, we were able to make several significant contributions to the research
community and get the attention of some of the top people in the field of molecular
dynamics, electroporation, and lipid research.
I thank my committee members — Dr. Gundersen for serving as chair of my non-
traditional committee, his encouragement, and his probing questions; Dr. Shing for her
honesty and guidance; and Dr. Warshel for his expertise and positive feedback.
I thank MOSIS and Cesar Piña for taking a chance on a Buckeye and providing financial
support, and Dr. Adebe for help with obscure math questions.
I thank Garrick Staples and the staff at the HPCC, Brad Cleaver, and Tom Wisniewski,
for keeping the computers humming along.
I thank our collaborators Dr. Tieleman and Dr. Dimova for their time and efforts and
contributions to my research.
iv
TABLE OF CONTENTS
DEDICATION................................................................................................................... ii
ACKNOWLEDGEMENTS ............................................................................................ iii
LIST OF TABLES ......................................................................................................... viii
LIST OF FIGURES ......................................................................................................... ix
ABSTRACT ..................................................................................................................... xii
CHAPTER 1: INTRODUCTION .................................................................................... 1
1.1 Overview ................................................................................................................... 1
1.2 Molecular Dynamics ................................................................................................. 5
1.2.1 Are molecular dynamics valid? ......................................................................... 5
1.2.2 Alternatives to Molecular Dynamics ................................................................. 8
1.2.3 Formalism of Molecular Dynamics ................................................................... 9
1.2.4 GROMACS Input Files.................................................................................... 12
1.2.5 Energy Minimization ....................................................................................... 13
1.2.6 Water Models ................................................................................................... 15
1.2.7 Phospholipid Parameters .................................................................................. 16
1.2.8 Electrostatics .................................................................................................... 17
1.2.9 Temperature and Pressure coupling ................................................................. 19
1.2.10 Bond and Angle restraints (LINCS, SHAKE, & SETTLE) ........................... 19
1.2.11 Analysis and Analysis Tools .......................................................................... 20
CHAPTER 2: PHOSPHATIDYLSERINE EXTERNALIZATION .......................... 21
2.1 Abstract ................................................................................................................... 21
2.2 Results and Discussion ........................................................................................... 22
2.3 Methods................................................................................................................... 30
2.3.1 Cell lines and culture conditions. ..................................................................... 30
2.3.3 Fluorescence microscopy. ................................................................................ 31
2.3.4 Molecular dynamics simulations. .................................................................... 31
2.4 Acknowledgment .................................................................................................... 33
CHAPTER 3: PHOSPHATIDYLSERINE EXTERNALIZATION PART 2 ........... 34
3.1 Abstract ................................................................................................................... 34
3.2 Introduction ............................................................................................................. 35
3.3 Materials and methods ............................................................................................ 40
3.3.1 Cell lines and culture conditions ...................................................................... 40
3.3.2 Pulse generator and pulse exposures ................................................................ 40
v
3.3.3 Fluorescence microscopy and microphotometry ............................................. 41
3.3.4 Molecular dynamics simulations ..................................................................... 41
3.3.5 Single bilayer ................................................................................................... 42
3.3.6 Double bilayer .................................................................................................. 43
3.4 Results and Discussion ........................................................................................... 44
3.4.1 Imaging observations of nanoelectropulse-induced PS externalization in
living cells ................................................................................................................. 45
3.4.2 Imaging observations of nanoelectropulse-induced membrane
permeabilization in living cells ................................................................................. 48
3.4.3 Molecular dynamics simulations of an asymmetric phospholipid bilayer in
a high electric field .................................................................................................... 50
3.4.4 Molecular dynamics simulations of an asymmetric phospholipid bilayer
with a porating transmembrane potential generated by charge separation ............... 57
3.4.5 Equivalence of MD simulations of nanoelectroporation with
transmembrane potential generated by global electric field and charge
imbalance .................................................................................................................. 66
3.4.6 Water dipole moment at the membrane interface — cause or effect? ............. 70
3.4.7 Are the electric fields and transmembrane potentials in these MD
simulations realistic? ................................................................................................. 75
3.5 Conclusions and Outlook ........................................................................................ 77
3.6 Acknowledgment .................................................................................................... 80
CHAPTER 4: THE MINIMUM PORATING ELECTRIC FIELD .......................... 81
4.1 Abstract ................................................................................................................... 81
4.2 Introduction ............................................................................................................. 82
4.3 Methods................................................................................................................... 88
4.3.1 Determination of the minimum porating field ................................................. 90
4.3.2 Dipole-field angle ............................................................................................ 91
4.3.3 Images .............................................................................................................. 91
4.4 Results and Discussion ........................................................................................... 91
4.4.1 Minimum porating electric field varies with phospholipid properties. ............ 94
4.4.2 Common poration sequence ............................................................................. 97
4.4.3 Location of pore-initiating defects ................................................................... 98
4.4.4 Head group tilt not a factor in electroporation ............................................... 106
4.5 Conclusions ........................................................................................................... 110
4.6 Acknowlegement .................................................................................................. 112
CHAPTER 5: STRUCTURE OF AN ELECTROPORE .......................................... 113
5.1 Abstract ................................................................................................................. 113
5.2 Introduction ........................................................................................................... 114
5.3 Simulation Methods .............................................................................................. 115
5.3.1 Molecular dynamics simulations ................................................................... 115
vi
5.4 Results and Discussion ......................................................................................... 117
5.5 Conclusions ........................................................................................................... 124
5.6 Acknowledgement ................................................................................................ 126
CHAPTER 6: EQUILIBRIUM PROPERTIES OF CALCIUM IN LIPID
BILAYERS .................................................................................................................... 127
6.1 Abstract ................................................................................................................. 127
6.3 Materials and Methods .......................................................................................... 131
6.3.1 Molecular dynamics simulations ................................................................... 131
6.3.2 Structures ....................................................................................................... 132
6.3.3 Calcium binding criteria ................................................................................ 133
6.4 Results and Discussion ......................................................................................... 134
6.4.1 System equilibration time, area per lipid, bilayer thickness .......................... 134
6.4.2 Calcium distribution in DOPC bilayers ......................................................... 136
6.4.3 Calcium distribution in DOPC:DOPS bilayers .............................................. 139
6.4.4 Calcium and chloride interfacial distribution in DOPC and DOPC:DOPS ... 139
6.4.5 Calcium density in DOPC and DOPC:DOPS systems .................................. 142
6.4.6 Water density in DOPC and DOPC:DOPS systems ...................................... 145
6.4.7 Bound, interface, and bulk calcium ............................................................... 146
6.4.8 Calcium-PS complex ..................................................................................... 152
6.4.9 Calcium and sodium in the bilayer interface ................................................. 155
6.4.10 Calcium and PC and PS head group dipole angle........................................ 159
6.5 Conclusions ........................................................................................................... 164
6.6 Acknowledgement ................................................................................................ 167
CHAPTER 7: CALCIUM BINDING KINETICS ..................................................... 168
7.1 Abstract ................................................................................................................. 168
7.2 Introduction ........................................................................................................... 169
7.3 Methods................................................................................................................. 171
7.3.1 Simulation Parameters ................................................................................... 171
7.3.2 Structures ....................................................................................................... 172
7.3.4 Kinetic Models ............................................................................................... 177
7.3.5 Equilibrium Properties ................................................................................... 178
7.4 Results ................................................................................................................... 180
7.4.1 Qualitative observations from the trajectory .................................................. 180
7.4.2 Time evolution of the area per lipid ............................................................... 185
7.4.3 Calcium binding to carbonyl oxygen and associated kinetics ....................... 186
7.4.4 Generalized calcium binding ......................................................................... 190
7.4.5 Calcium-carboxyl complex Ca·COO
-
............................................................ 191
7.4.6 Calcium-carbonyl complex Ca·CO ................................................................ 192
7.4.7 Calcium-phosphatidyl oxygen complex Ca·PO
4
2-
......................................... 192
7.4.8 Combination complexes ................................................................................. 192
vii
7.4.9 Calcium and the head group dipole angle ...................................................... 194
7.4.10 Updated kinetic model ................................................................................. 195
7.4.11 Equilibrium binding constants ..................................................................... 200
7.4.12 Water hydration shell ................................................................................... 203
7.4.13 Water organization at the lipid interface ...................................................... 205
7.4.14 Diffusion rates of bound and free calcium ................................................... 206
7.4.15 Diffusion rates of lipids ............................................................................... 208
7.5 Conclusions ........................................................................................................... 210
7.6 Acknowledgement ................................................................................................ 211
REFERENCES .............................................................................................................. 212
viii
LIST OF TABLES
Table 4.1. Relationship between threshold field and bilayer properties .......................... 94
Table 4.2. Phosphorus-nitrogen (phosphatidylcholine head group dipole) vector
without an external electric field and at the minimum porating field ............................. 106
Table 6.1. System configurations .................................................................................... 133
Table 6.2. Variation in bilayer area per lipid and thickness with number of calcium
ions and PS molecules ................................................................................................... 136
Table 6.3. Calcium density in phosphate and carbonyl regions of the phospholipid
interface in the DOPC and DOPC:DOPS leaflets. .......................................................... 143
Table 6.4. Water density in phosphate and carbonyl regions of the phospholipid
interface ........................................................................................................................... 146
Table 6.5. Distribution of calcium in DOPC and DOPC:DOPS bilayers ....................... 147
Table 6.6. Equilibrium calcium binding in DOPC and DOPC:DOPS systems .............. 148
Table 6.7. Calcium coordination with DOPS carboxyl oxygens in DOPC:DOPS. ........ 153
Table 6.8. Sodium density in phosphate and carbonyl regions of the phospholipid
interface in the DOPC and DOPC:DOPS leaflets. .......................................................... 158
Table 7.1. System configurations .................................................................................... 173
Table 7.2. Calcium·Carbonyl Binding Rate Constants ................................................... 188
Table 7.3. Generalized Binding Rate Constants ............................................................. 200
Table 7.4. Equilibrium binding constants ....................................................................... 201
Table 7.5. Distribution of the coordination number of bulk calcium with wate ............. 204
Table 7.6. Equilibrium binding constants (Loss of 50% Hydration Shell) ..................... 205
Table 7.7. Diffusion Coefficients for Bulk and Bound Calcium .................................... 207
Table 7.8. Diffusion Coefficients for DOPS and DOPS ................................................. 209
ix
LIST OF FIGURES
Figure 2.1. PS externalization in RPMI 8226 human multiple myeloma cells after
exposure to 4 ns, 8 MV/m pulses ...................................................................................... 24
Figure 2.2. Nanoporation and PS translocationcation in a DOPC-DOPS bilayer in a
450 mV/nm electric field .................................................................................................. 26
Figure 2.3. Water dipoles aligned in the low dielectric permittivity of the membrane
interior during the formation of a hydrophilic pore .......................................................... 28
Figure 3.1. Membrane phospholipid bilayer simulation — double bilayer with ion
imbalance .......................................................................................................................... 44
Figure 3.2. Nanoelectropulse-induced membrane perturbation and
phosphatidylserine externalization.................................................................................... 46
Figure 3.3. Nanoelectropulse-induced membrane permeabilization ................................ 49
Figure 3.4. Membrane electric field versus external electric field vector ........................ 52
Figure 3.5. Resealing of hydrophobic pore ....................................................................... 55
Figure 3.6. Mass density, electric field, electric potential, and mean water dipole
moment for one bilayer in a double bilayer simulation with charge separation ............... 60
Figure 3.7. Pore formation and PS translocation in one bilayer in a double bilayer
simulation with a 5-ion charge separation ........................................................................ 64
Figure 3.8. Equivalence of global field and charge imbalance methods for generating
transmembrane potential ................................................................................................... 68
Figure 3.9. Effect of transmembrane charge imbalance on the mean water dipole
moment profile across a DOPC-DOPS bilayer before poration ....................................... 73
Figure 4.1. The average time to the formation of a hydrophilic pore based on three
independent simulations is plotted against the applied external field ............................... 93
Figure 4.2. Pore formation sequence for DOPC and POPC ............................................. 96
Figure 4.3. Membrane properties for 0.1 nm slices normal to the membrane surface ... 101
Figure 4.4. The field profile for POPC without a field and at the minimum field for
poration (320 mV/nm) .................................................................................................... 105
x
Figure 4.5. The distribution of the phosphorus-nitrogen vector ..................................... 107
Figure 5.1. Dipole-electric field alignment during pore formation in a POPC bilayer .. 118
Figure 5.2. Snapshots from the POPC bilayer electroporation simulation ..................... 120
Figure 5.3. Distribution of head group P →N dipole ...................................................... 121
Figure 5.4. P →N dipole angles with bilayer normal versus time for representative
lipid molecules ................................................................................................................ 122
Figure 5.5. Field-driven head group dipole alignment. .................................................. 123
Figure 6.1. Time Evolution of Area per Lipid and Calcium Pairs .................................. 135
Figure 6.2 Number density (nm-3) profiles for selected components of DOPC and
DOPC:DOPS................................................................................................................... 137
Figure 6.3. Number density (nm-3) profiles for the interface regions ............................ 141
Figure 6.4 Radial distribution functions ......................................................................... 150
Figure 6.5. Calcium-DOPS-water complex .................................................................... 154
Figure 6.6. Number density (nm-3) profiles ................................................................... 156
Figure 6.7. Calcium-induced changes in PC and PS head group orientation in DOPC
and DOPC:DOPS bilayers .............................................................................................. 161
Figure 7.1. Pair correlation functions ............................................................................. 174
Figure 7.2. Potential binding atoms for DOPC and DOPS ............................................. 176
Figure 7.3. Snapshots from the 0 PS 100 Ca simulation ................................................ 179
Figure 7.4. Progression of free to bound calcium ........................................................... 181
Figure 7.5. Time course of the equilibration of area per lipid in asymmetric
DOPC:DOPS bilayers with calcium ions ....................................................................... 183
Figure 7.6. Snapshots showing the migration of calcium from bulk water to the
membrane for the 0 PS 10 Ca system ............................................................................. 184
Figure 7.7. Fit of carbonyl-only binding model for the 0 PS 10 Ca, system .................. 190
xi
Figure 7.8. Calcium complexes ...................................................................................... 191
Figure 7.9. Time evolution of the P →N dipole .............................................................. 195
Figure 7.10. Number of calcium complexes based on coordination number ................. 196
Figure 7.11. Number of calcium complexes ................................................................... 198
Figure 7.12. Snapshot of bulk calcium coordinated with 8 water molecules ................. 203
Figure 7.13. The bilayer normal component of the water dipole moment ..................... 206
xii
ABSTRACT
Recent advances in computing technology have facilitated the application of simulations
to studying biological systems at the atomic level. In particular atomistic molecular
dynamics provide an opportunity to model systems that are unobservable through
conventional experimental methods as well as supplement understanding of observations.
In this thesis molecular dynamics were applied to study biological cell membranes,
specifically lipid bilayers, the primary constituent of the cell membrane, in electric fields,
and to understand the mechanism and events associated with electroporation.
Electroporation is a widely used experimental and commercial technique for introducing
normally excluded compounds such as DNA, RNA, ions, drugs, and other chemicals into
cells. Traditional electroporation utilizes kilovolt-per-meter electric fields applied on the
order of microseconds that disrupt and scramble the cell membrane and allow diffusive
entry, however, ultra-short nanosecond pulses at megavolt-per-meter fields produce
different effects in cells which remain largely uncharacterized. One such effect, the
migration of the negatively charged lipid phosphatidylserine from the inner to outer
leaflet of the cell, is of particular biological interest because of its association with cell
apoptosis, or programmable cell death. Control of such an event could be useful in
developing a targeted treatment for removing unwanted cells such as tumors or
melanoma.
xiii
In this thesis I begin by introducing electroporation and its history and explain how
molecular dynamics and its techniques can help advance our understanding of the field.
In the following chapters, consisting of peer reviewed and submitted journal articles
loosely tied together, I explore the mechanism of phosphatidylserine translocation
induced by nanosecond pulses in megavolt-per-meter electric fields, and correlate
experimental data and anecdotal evidence of phosphatidylserine translocation in vitro
with a detailed molecular mechanism provided by simulations. Upon developing this
relationship, I further explore the more generic mechanism of electroporation by
simulating the behavior of lipids of differing composition and introduce the concept of a
minimum porating electric field to aid comparison. Next I probe into the detailed
structure of an electropore, focusing specifically on the alignment of both the lipid
headgroups and water dipoles along the pore walls and unperturbed sections of the
bilayer. I conclude by studying the lipid structural changes caused by the introduction of
calcium ions to the system, as well as the kinetics of calcium binding to lipid bilayers.
1
CHAPTER 1: INTRODUCTION
1.1 Overview
Biological membranes, universal to all life, provide structural integrity to cells and
compartmentalize and isolate different functional units. They form the barrier between
the cytoplasm and extracellular fluid and separate organelles. Comprised primarily of
various phospholipids of amphipathic character aligned tail to tail, the cell membrane is
more than just a physical barrier. The membrane is an active and functional part of a cell
and responds dynamically to changes in physiological conditions by utilizing its
embedded proteins, modifying its structural arrangement, or changing its composition. It
helps regulate ion transport and osmotic pressure and is involved in cell signaling and
communications. For example, phosphatidylserine, a lipid relegated to the membrane
interior of a healthy cell, when externalized is a symptom of apoptosis or programmable
cell death (202). The exposure of negative charge on the cell exterior acts as a signaling
mechanism which is used to trigger macrophages in the blood stream to dispose of the
distressed cell (200).
Due to the large number of functions that the cell membrane helps manage,
understanding the delicate balance of the cell membrane is important. Furthermore, a
rigorous understanding of the behavior of the membrane under a variety of conditions can
give insight into ways to engineer the cell membrane to do useful things. A central
2
problem of medicine is the delivery of therapeutic compounds to a specific cell. Many
potentially curative compounds go unutilized due to the inability to deliver the drug to the
proper location. A drug might not be soluble in the correct tissues or it may be harmful if
absorbed in the wrong part of the body. In nature, functional peptides, ion channels, and
viruses such as HIV have evolved to overcome the barrier function of the membrane and
allow specific entry and exit of compounds. Man, being much less sophisticated than
Mother Nature, has developed other techniques for circumventing the barrier function of
membranes. For instance, mechanical and electrical stresses can be applied to the
membrane to open pores, and to introduce excluded substances inside a cell.
Electroporation, a technique in which a large electrical potential is applied across a cell,
results in an increase in the conductivity and permeability of the membrane. This
process, which can be both reversible and irreversible, causes pores and is used to
introduce probes, DNA, plasmids, drugs, or other species inside the cell. The precise
effect of an electric field on a cell membrane or lipid bilayer, however, is dependent on
many factors including the strength and duration of the electric field, the ionic strength of
the surrounding solution, and the type and composition of the cells or lipid bilayers.
Understanding the complex interplay of these factors will continue to improve the utility
of electroporation and its applications in helping cure human disease.
Much of the early work characterizing the electrical properties of cells has been attributed
to Schwan (162) and Cole (31). They observed the frequency dependence of the
3
membrane dielectric and conductivity. In parallel, Goldman measured the voltage-current
relationship of plant viruses (55) and observed an increase in membrane conductance
with an applied field. Later, Neumann and Rosenheck (124) would observe temporary
increases in the permeability of chromaffin granules using kilovolt-per-centimeter fields,
which decayed after 150 microseconds and resulted in limited cell destruction. This
reversible permeabilization or electrical breakdown coincided with a large increase in the
conductance of the cells (1, 10, 27, 92). Non reversible electropermeabilization was first
described in detail by Sale and Hamilton (63, 155, 156) in a series of three papers that
described the non-thermal killing of bacteria, yeast, and erythrocytes when exposed to
applied fields of 30 kilovolts-per-centimeter. Lysis was observed when the potential
difference across the membrane reached between 0.3 V and 1.15 V.
Direct observation of the processes described above were not experimentally feasible,
consequently a combinational of mathematical models (35, 36, 56, 176, 218) and indirect
observations (51, 70, 119, 125, 149) were correlated in order to explain this process
coined electroporation. Experimental systems consist primarily of cell suspensions,
individual cells, and planar lipid bilayers and the experiments are divided into two main
groups — charge pulse and voltage clamp. Charge pulse experiments apply a constant
electric field in doses and measure the resulting electrical properties of the system.
Voltage clamp experiments use a feedback system to maintain a constant voltage over the
membrane. In charge pulse experiments four outcomes are possible for the membrane
(6): 1) a small increase in electrical conductance, 2) mechanical rupture, 3) incomplete
4
reversible electrical breakdown, or 4) complete reversible electrical breakdown, the latter
two differing by whether or not the membrane is completely discharged. Charge pulse
experiments take place on a shorter time scale than voltage clamp since the discharge
takes place rapidly in less than one microsecond.
As explained previously much of the interest in electroporation came from its transport
applications, however, recent observations indicate that electrical pulses at the megavolt-
per-meter scale induce different biological responses that are different from the four
mentioned above. Nanosecond megavolt-per-meter pulses do not cause mechanical or
electrical breakdown of the cell, nor significantly increase the electrical conductance or
temperature. However they do allow entry of small molecules such as YO-PRO-1 (206)
and cause calcium bursts (202, 219) and migration of some cells to apoptosis (8, 9, 201).
These observations indicate fundamental differences in microsecond-kilovolt-per-meter
versus nanosecond-megavolt-per-meter electroporation. Direct observation of the events
causing these results is difficult, however, the shortened time scale of nanosecond
megavolt-per-meter pulses is within the range of real time atomistic molecular dynamics
simulations.
5
1.2 Molecular Dynamics
1.2.1 Are molecular dynamics valid?
Studies of membranes and lipid bilayers date back to the 18
th
century and first
characterizations of the electrical properties of the cell date back to 1925 when Fricke
measured the capacitance of red cell membrane (50). He measured it to be 1 microfarad
per square centimeter, from which he estimated the thickness of the membrane, 33
angstroms. Advances in electron microscopy technique later gave more accurate
measurements of the membrane thickness and atomic force microscopy can now probe
individual atoms and provide high resolution images of the membrane surface.
Diffraction methods give the areas per lipid down to square angstrom level precision, and
NMR can describe the ordering of lipid tails and verify lipid area calculations determined
by X-ray diffraction. These experiments would provide benchmarks for lipid models.
In 1981, Ploeg and Berendsen published the first molecular dynamics simulation of a
bilayer system (197). The membrane was represented by two layers of 16 decane
molecules restrained by a dihedral potential function and was used to estimate the order
parameter S
CD
, which measures the carbon deuterium bond direction with respect to the
bilayer normal. The order parameter, a second order Legendre polynomial, was first used
by Seelig to model the rigidity and flexibility of the lipid bilayer and determine the
structure of the bilayer tails (164, 165). The 32 molecule simulation lasted 80
picoseconds and took 10 hours. Order parameters were successfully reproduced,
6
verifying that lipid rigidity decreases further down the hydrocarbon chain, and showed
that molecular dynamics of lipids could provide results consistent with experimental data.
Time scales and system sizes, however, were limited to a few square nanometers and
under a nanosecond, which meant studies were limited to order parameters, density and
electron density profiles, and lipid tilt angles.
As time progressed and computers became faster the complexity and duration of
simulations increased. An explosion in computational power meant more atoms could be
simulated for longer times. Solvents could be treated explicitly instead of using a crude
potential function, better methods for calculating long range electrostatic interactions
were possible, and more extensive parameterization and force fields could be used. This
combination of factors now allows simulations to make predictions about physical
systems that are beyond the reach of current experimental techniques. For example ion
transport could be observed, lipids could be completely solvated (fully hydrated),
dynamic processes like electroporation could be simulated in near all atom detail (193,
194), and even free energy profiles could be obtained (195). Simulating non-equilibrium
events is an advancement over simulating basic lipid properties. Still, the reliability of
any predictions must be qualified, and the selection of parameters must be verified and
grounded by comparison with available experimental data.
All molecular dynamics simulations are dependent on the choice of force field and its
parameters which define the bonds, bond angles, electrostatic interactions, and van der
7
Waals interactions with a combination of potential functions. The most widely used
parameters for lipid simulations using GROMACS use a combination of parameters from
Berger et. al (14), the OPLS force field (86), the SPC water model (11), and specialized
parameters for headgroups which are obtained from ab initio calculations. Tieleman et. al
performed a detailed study of the effect of different ensembles, water models, boundary
conditions and electrostatics conditions on lipid simulations (189). The authors concluded
that constant pressure simulations are preferred to their constant volume counterparts,
SPC, although one of the simplest water models, better reproduces the lipid water
interface, and that full atomic charges could be used. Further improvements would
require a force field capable of handling polarizability rather than extensive
reparameterization of the current ensemble. A more recent study by Anezo et. al.
confirmed an earlier conclusion based on longer 150 ns simulations (4). A limited
number of properties are available both experimentally and through simulation including
electron and mass density, order parameters, dipole potential, and lipid areas. The de
facto standard for determining whether or not to trust lipid simulations has become the
area per lipid, which is easily computed by dividing the surface area of the bilayer (XY)
plane with half the number of lipids. Areas for fully hydrated lipids are available from
work done by Nagle and Tristram-Nagle including but not limited to DPPC (98, 121),
DMPC (97), DLPC (97), DOPC (99), DOPS (141), and POPC (99). Systems are
considered production ready when the fluctuation in equilibrium system area is less than
1% for an extended period of time typically taken to be at least 10 ns. With the membrane
8
structure and properties well defined and consistent between experiment and simulations,
new efforts focused on transport properties of water and ions, as well as small molecule
permeability. Other simulations also studied transmembrane species embedded in a lipid
bilayer including peptides (88, 89), porins (157, 158, 191), and cholesterol (44, 73, 120,
137) and their effect on the bilayer. Simulations could also help answer questions
regarding permeability, lipid flip flop, diffusion, ion binding, and electroporation.
1.2.2 Alternatives to Molecular Dynamics
Prior to molecular dynamics many mathematical models were used to model membrane
systems and continue to be refined. The simplest models treat lipid bilayers as a smooth
hydrophobic surface and apply Poisson's equation or Nerst-Planck relations to determine
the potential at the surface. Double layer models such as the Gouy-Chapman-Stern model
are better yet, but the lipid headgroups aren't a smooth surface and have their own
charges which complicate the electrostatic environment.
These base models for calculating the potential at the bilayer surface were extended into
other models of electroporation. Electroporation models fall into two categories 1)
surface potential models, which calculate the potential profile at the bilayer and combine
it with a model of bilayer stability and 2) stochastic pore models which assume the
existence of pores and that the membrane oscillates between stable and porated phases.
9
Abidor, Chizmadzhev, and Chernomordik pioneered models of the first type (1, 28) while
Sugar and Neumann developed the second type (177, 178).
Modern models use aspects of both techniques. DeBruin and Krassowska developed
models for electroporation that combined electrical potential, ion current, and pore
density (34-36). The most recent model also includes terms for concentrations of
different ionic species. The terms are derived from a Poisson’s equation application of an
electrical cell, Nerst-Planck relations, a first order differential equation for the pore
density, and physical constants based largely on the potential induced on sea urchin eggs
(71).
Transport lattice models (56, 176) pioneered by Weaver and coworkers have also been
applied to the electroporation problem and present an alternative to the differential
equation based models. Each unit cell is treated as a circuit element connected to 4 (2D)
or 6 (3D) other nodes. Each unit cell can contain capacitors and resistors and have
different values depending on whether the nodes are membrane or electrolyte.
1.2.3 Formalism of Molecular Dynamics
Molecular dynamics solve an adjusted set of Newton’s equations of motion for all atoms
in a system
10
1 …
(1.1)
where forces in this equation are the negative derivatives of a potential function with
respect to each atom.
(1.2)
The underlying concept is simple yet powerful. With an understanding of the position of
all atoms in the system for the entire duration of a simulation, nearly any physical
property can be calculated using statistical mechanics or other mathematical
relationships. In order to improve the realism of the system, several adjustments are made
in order to account for physical constraints. The internal energy and pressure of the
system are kept constant by the use of a thermostat and barostat respectively and bond
lengths and angles are treated as constraints on Equations 1.1 and 1.2.
Some interactions cannot be described by classical mechanics. For example when the
resonance frequency ν of a bond is near to or greater than K
b
T/h (where K
b
is the
Boltzmann constant, T is the temperature, and h is Planck's constant), as is the case with
hydrogen and helium bonds, the quantum mechanical oscillator energy deviates from the
classical version. Either a correction term must be applied to the energy or united atoms
can be implemented. Carbon bonded to one or more hydrogen atoms can be modeled as a
11
single united atom. This circumvents the discrepancy between classical and quantum
energies and allows longer time steps (typically 2 or 4 picoseconds) to be used.
Molecular dynamics have additional limitations. All atoms are assumed to be in the
ground state and simulating chemical reactions, which are best described by quantum
mechanics, is difficult. Computational power limits the number of atoms in the system
and forces the use of cutoff terms for interactions in order to reduce calculations. One
must then choose between surrounding the system with a vacuum or utilizing periodic
boundary conditions. Both choices can results in simulation artifacts including excess
ordering and net system dipoles.
The general simulation algorithm conceptually isn’t very complex. The inputs include the
positions and current velocities of all atoms in the system and a corresponding potential
function that describes their interactions. For every atom in the system the force is
calculated and then Newton’s equations are solved numerically. Finally the new
positions and velocities are written and the process is repeated. Many packages are
available that perform these calculations including CHARMM (22), AMBER (24),
GROMACS (13, 105, 198), MOLARIS (169), and NAMD (143). They vary primarily in
the way force fields are constructed — all atom versus atomic groups, their distribution
method, and the set of built in analysis tools. While there has been some effort for
interoperability between programs, certain applications are better suited for one package
than another. For the purpose of studying lipid bilayers a robust set of force fields are
12
available for the GROMACS package that have been tested and validated over many
years of work. Furthermore, GROMACS is distributed under the GNU General Public
License (GPL) which makes it free to operate and modify. Additionally, a set of analysis
tools relevant to studying bilayers is included along with source code. Due to these
advantages GROMACS was chosen as the platform for all work and all references to
simulations should be assumed to apply to GROMACS unless otherwise noted.
1.2.4 GROMACS Input Files
The GROnigen MAchine for Chemical Simulation, or GROMACS for short, is a
complete molecular dynamics and analysis suite. In order to start a simulation with
GROMACS very few items are needed. Molecular coordinates, in either native .gro or
.pdb (protein databank) format, a force field and molecular topology, and a set of
simulation parameters including simulation time, pressure, temperature, time step, cut off
radii, output options, etc. are all that are required for simulations. Many of the files can be
automatically generated by the built-in GROMACS tools if one uses one of the supported
force fields. The output of the primary program mdrun is a .xtc or .trr trajectory file,
which contains the positions and velocities of all atoms at all times. Simulation
structures, however, are difficult to create and very few options are available for building
them besides manual coding. Software packages such as Molden or Insight II lack
functionality and can create structures which are inconsistent. This leads to unnatural
13
angles, bond lengths, overlapping atoms, and naming inconsistencies. Other alternatives
include the NCBI Protein Database which has many structures, but like the software
packages mentioned above, give structures that are only consistent with a limited number
of force fields. Lipid structures in particular are difficult to find. Most structures are
therefore made by manually adjusting coordinates and working from existing models to
create more complex structures. Incremental changes are made and then fed into the
energy minimization routine. The goal is to create a structure in which the energy
minimization algorithm can fix the bad contacts and high energy regions in order to
produce a workable structure.
1.2.5 Energy Minimization
Systems are energy minimized using a steepest descents method. New positions are
calculated by Equation 1.3 —
max| |
(1.3)
where F is the force, r is the position of the atom, and h is the maximum displacement, a
parameter of the algorithm. F is the negative gradient of the potential V defined by
Equation 1.4.
14
4
1
2
1
2
(1.4)
1
The terms in order above represent electrostatic interactions, van der Waals interactions,
direct bonded interactions, angles, dihedral angles, and improper dihedral angles. The
target |Fn| and number of iterations are specified and the process continues until one of
these conditions is met. The component forces (Equations 1.5-1.11) are calculated from
the derivative of the potential function.
4
(1.5)
12
6
(1.6)
(1.7)
(1.8)
15
(1.9)
(1.10)
1.2.6 Water Models
Many different water models exist and vary in complexity. The most commonly used
water model is the simple point charge model or SPC model (11). The SPC water model
treats each atom as a center of charge. Each hydrogen has a charge of 0.41 e and the
oxygen has a charge of -0.82 e. The assumption of point charges result in an incorrect
dipole moment, so the angle H-O-H is changed to 109.42° versus the experimentally
determined value of 104.45° to reproduce the correct permanent dipole. Because SPC
water is missing two lone electron pairs, it diffuses faster than real water.
This effect decreases as temperature increases. The SPC/E and TIP3P models adjust the
charges relative to SPC, while the TIP4P creates a tetrahedral with a dummy location. A
comprehensive study of water models can be found on Chaplin's website (26). For the
purpose of studying water interfaces including lipid water interfaces, the SPC model
provides results closest to experimental observations, and has the added advantage of
being computationally expedient. SETTLE, an analytical version of the SHAKE
16
algorithm for restraining bond lengths, exists for SPC water provides a 3 to 4 fold
increase in performance over iterative methods.
1.2.7 Phospholipid Parameters
The phospholipid parameters used in most present day lipid simulations have roots in the
GROMOS and OPLS (optimized parameters for liquid systems) based force fields. The
primary feature of the force field was that it was designed to replicate the density of the
lipids, known experimentally within 1% accuracy. This was done by adjusting the
Lennard-Jones parameters of a pentadecane system until the density and heat of
vaporization agreed with experiment. The parameters were used in order to estimate the
area per lipid of dipalmitoylphosphatidylcholine (DPPC) bilayers. Area per lipid is
calculated by simulating the system for a period of time until the box dimensions normal
to the bilayer are constant for a period of many nanoseconds. The box dimensions are
multiplied together and divided by the number of lipids in a single leaflet. It is
recommended that the system be allowed to equilibrate for more than 25 ns for the best
reproduction of properties.
Much work has been done optimizing parameters for applications to specific simulation
systems; however, most current papers no longer validate simulation parameters beyond a
few references. The de facto standard for testing simulation parameters is the area per
lipid test. Most lipids fall between 0.56 and 0.72 nm
2
. Better estimates are available for
17
pure lipids, however, mixture data is often not available. Deviations in areas per lipid can
vary as much as 25% from the pure component for mixtures of DMPC with DMTAP and
for DPPC with cholesterol.
Many choices exist for simulation conditions however some are more optimal than
others. One of which is the choice of ensemble. Common choices include constant
volume (NVT) and constant isotropic pressure (NPT) and constant anisotropic pressure.
For the purpose of determining lipid properties the NVT ensemble is a poor choice since
the area per lipid should be determined by the simulation rather than set initially. Berger
concludes for lipid simulations NPT is the best choice (14).
1.2.8 Electrostatics
Electrostatic interactions comprise the majority of the computational resources involved
in the simulation, and are also the most important factor in determining the simulation
results. Electronic interactions scale with the inverse of distance, consequently, they
cannot be cut-off arbitrarily. A twin range cutoff technique is used. Short range
electrostatic interactions are calculated every time step and are cut off at 1 nm. Long
range electrostatic interactions are calculated every 10 time steps using the Particle Mesh
Ewald (PME) method. This combination provides exact results when no molecules pass
within the 1 nm neighbor-list cutoff.
18
Simulations of lipid bilayers are inherently periodic due to the structure of the membrane
which requires special treatment of the electrostatics. In order for the PME method to
converge when using periodic boundary conditions the system must be electrically
neutral. If this condition is not made an opposing charge is evenly distributed around the
system. Additionally for periodic systems with electric fields, PME calculations for the
long range electrostatics introduce error in the total energy if the system has a net dipole.
In periodic systems with an externally applied electric field (a simulation parameter
resulting in an additional force term in the potential) this has the consequence of
introducing a system size dependence on the field experienced by atoms. The exact
relationship between the applied and experienced fields is non-linear and difficult to
predict since it depends on the distribution and properties of atoms in the system.
Spherical boundary conditions produce the most physically consistent results. In 3-d
periodic boundary conditions with either electric fields or dipoles, non physical image
charges in the periodic cells are introduced in order to allow a convergence path for the
PME treatment of electrostatics (213) which can cause artifacts in the system. However,
simulations of the free energy required to rotate a water dipole in bulk water, the lipid
headgroup region, and the low dielectric membrane interior with proper spherical
treatment of the electrostatic interactions did yield any significant differences to
calculations made using PME with periodic boundary conditions.
19
1.2.9 Temperature and Pressure coupling
Temperature and pressure coupling are each handled by a Berendsen thermostat or
barostat of the following form
(1.11)
where T is the temperature in Kelvin at a simulation time t. T
0
is target simulation
temperature and τ is the relaxation time in picoseconds. Pressure coupling is handled in
the same manner by substituting P for T. Since pressure is a tensor rather than a scalar
the different dimensions can be treated either anisotropically, isotropically, or
semiisotropically.
1.2.10 Bond and Angle restraints (LINCS, SHAKE, & SETTLE)
The time step in classical MD simulations is limited by the bond oscillations. In order to
accelerate the solution all bonds are assumed to be in the ground state and modeled using
constraints. After each time step the length of the bonds and angles need to be returned
to acceptable values. Since the bonds and angles are all intertwined the solution is
nonlinear. Analytical solutions exist for small molecules such as water, and the SETTLE
algorithm is the fastest. For larger molecules an iterative method is used, typically either
20
LINCS or SHAKE. LINCS is preferred since it is 3-4 times faster and has a wider
convergence range.
1.2.11 Analysis and Analysis Tools
Analysis tools read the binary output files which contain information regarding the
positions, forces, velocities, and energies of all atoms in the system at specified intervals.
With such information it is possible to calculate basic properties such as density (both
electron and mass), charge distributions, average orientations, and others. It is also
possible to derive thermodynamic properties from the radial distribution function such as
the diffusion coefficient. These tools are provided within the GROMACs set of
programs. It is also possible to write specialized applications for specific systems
including DNA and polymers. Some of these applications are available at
http://www.gromacs.org/contributed_by_users/ while others remain proprietary to the
author.
21
CHAPTER 2: PHOSPHATIDYLSERINE EXTERNALIZATION
Portions of this chapter are adapted from:
Vernier, P. T., M. J. Ziegler, Y. Sun, W. V. Chang, M. A. Gundersen, and D. P.
Tieleman. 2006. Nanopore formation and phosphatidylserine externalization in a
phospholipid bilayer at high transmembrane potential. Journal of the American Chemical
Society 128:6288-6289.
2.1 Abstract
Atomic-resolution molecular dynamics simulations of lipid bilayers containing 7%
phosphatidylserine (PS) on one leaflet are consistent with experimental observations of
membrane poration and PS externalization in living cells exposed to nanosecond,
megavolt-per-meter electric pulses. Nanometer-diameter aqueous pores develop within
nanoseconds after applying an electric field of 450 mV/nm, and electrophoretic transport
of the anionic PS headgroup along the newly constructed hydrophilic pore surface
commences even while pore formation is still in progress.
22
2.2 Results and Discussion
Externalization of phosphatidylserine (PS), a phospholipid normally confined to the inner
leaflet of the plasma membrane, follows exposure of biological cells to high-field electric
pulses of less than 100 ns duration (203). Nanosecond, megavolt-per-meter pulses cause
intracellular calcium release (202, 219) and induce apoptosis(9), but membrane
permeability to propidium iodide, a standard electroporation indicator, is not observed.
Nanoelectropulse-driven PS translocation is immediate and distinct from the PS
externalization that is symptomatic of apoptosis (201). Although nanosecond pulses do
not cause conventional electroporation, we have recently observed YO-PRO-1 influx, a
sensitive fluorescent reporter of membrane integrity, after exposure of cells to multiple 4
ns, 8 MV/m pulses, indicating a perturbative but not disruptive breaching of the
membrane hydrophobic barrier(205).
We present here observations of localized PS externalization in response to nanosecond
electric pulses (imaged for the first time with annexin V binding), and we correlate this
experimental data with molecular dynamics (MD) simulations of
dioleoylphosphatidylcholine-dioleoylphosphatidylserine (DOPC-DOPS) asymmetric
membranes in high electric fields. We demonstrate the association between field-driven
poration and PS translocation, and we draw a line of continuity from experiments on
living cells to mathematical models of pore formation in lipid bilayers, leading to the
23
conclusion that nanoelectropulse-driven PS externalization: 1. occurs on a nanosecond
time scale; 2. is restricted to the anode-facing pole of the cell; 3. is intimately connected
with and dependent on nanopore formation.
Fluorescence microscopy with the fluorescent-tagged PS-binding protein annexin V and
FM1-43, a cationic fluorochrome that partitions between the aqueous medium and the
outer leaflet of the cytoplasmic membrane, where its fluorescence greatly increases (228),
reveals nanoelectropulse-induced perturbations associated with PS externalization
localized exclusively at the anode-facing pole of the cell immediately after pulse
exposure (Figure 2.1). Because the kinetics and sensitivity of these imaging methods
permit only micrometer- and millisecond-scale resolution of post-pulse membrane
restructuring (194) and cannot be used to visualize molecular events with nanosecond
resolution, we turned to MD to provide insight into the mechanisms of nanosecond pulse-
driven poration and PS translocation.
24
Figure 2.1. PS externalization in RPMI 8226 human multiple myeloma
cells after exposure to 4 ns, 8 MV/m pulses. Separate cells stained with
annexin V -FITC (a) and FM1-43 (b) are shown before and 5 s after
pulse delivery (c and d). The anode (positive terminal) is at the top of
each image; the cathode is at the bottom. Fluorescence indicating PS
translocation appears at the anode-facing pole of the cell. Note
correspondence of the regions of pulse-induced annexin V binding and
FM1-43 fluorescence intensification.
Anode
Cathode
0 s
5 s
0 s
5 s
ab
d c
Previous MD electroporation studies show that with electric fields of 400 mV/nm water
intrusions into the bilayer extend across the membrane within nanoseconds, with
phospholipid headgroups following the initial water chain to line the resulting
hydrophilic pore (76, 182). In coarse-grained simulations, poration and PS externalization
occur during a single 10 ns pulse (193).
Using Gromacs-generated atomic resolution representations of a DOPC-DOPS bilayer in
water at 300 K, 1 bar, we ran simulations with varying field strengths and polarities. One
Figure 2.1. PS externalization in RPMI 8226 human multiple myeloma
cells after exposure to 4 ns, 8 MV/m pulses. Separate cells stained with
annexin V-FITC (a) and FM1-43 (b) are shown before and 5 s after
pulse delivery (c and d). The anode (positive terminal) is at the top of
each image; the cathode is at the bottom. Fluorescence indicating PS
translocation appears at the anode-facing pole of the cell. Note
correspondence of the regions of pulse-induced annexin V binding and
FM1-43 fluorescence intensification.
25
bilayer face contains 64 DOPC molecules; the other, representing a cell membrane
internal leaflet, contains 4 DOPS and 60 DOPC molecules, with 4 Na
+
ions to balance the
DOPS charge. Parameters and simulation methods are as previously described for
homogeneous bilayers (182, 204). Snapshots from a typical simulation are shown in
Figure 2.2.
26
Figure 2.2. Nanoporation and PS translocation in a DOPC-DOPS
bilayer in a 450 mV/nm electric field, directed from top to bottom, with
the positive electrode at the top. Water atoms are red (O) and white
(H); phospholipid tails are gray; glycerophosphate atoms are yellow ,
choline atoms green, serine atoms blue. Serine atoms are enlarged to
enhance visibility.
5.0 ns 5.5 ns
ab
c
6.0 ns 6.5 ns
d
Figure 2.2. Nanoporation and PS translocationcation in a DOPC-DOPS
bilayer in a 450 mV/nm electric field, directed from top to bottom, with
the positive electrode at the top. Water atoms are red (O) and white
(H); phospholipid tails are gray; glycerophosphate atoms are yellow,
choline atoms green, serine atoms blue. Serine atoms are enlarged to
enhance visibility.
27
With an electric field of 450 mV/nm, water intrusion occurs in about 5 ns (Figure 2.2a),
from the cathode side of the membrane. Within 0.5 ns a column of dipole-aligned water
molecules (Figure 2.3) spans the membrane interior, and phospholipid headgroups begin
to realign along the new aqueous surface (Figure 2.2b). By 6 ns PS headgroups are
migrating along the walls of the incipient pore (Figure 2.2c) toward the anode (Figure 2-
2d). In less than 7 ns the first PS molecule crosses the 5 nm span of the lipid membrane
interior. In 30 simulations with fields from 200 to 800 mV/nm, water defects appear
predominantly at the cathode-facing leaflet of the bilayer, regardless of the location of
PS, and PS translocates only after pore formation and always to the side of the membrane
facing the positive electrode, as in live cell observations (194) (Figure 2.1). Pore
formation occurs more rapidly as the field increases, and in some higher field cases water
intrusion occurs also at the anode side of the membrane. Poration is observed in less than
5 ns in 21 out of 24 simulations when the field is 450 mV/nm or greater. PS translocation
occurs only along pore walls and only when the field direction is as indicated in Figure
2.2. PS does not facilitate water pore formation, consistent with the suggestion that
interactions in the hydrophobic interior of the membrane dominate the relative ease of
propagation of a water column into the bilayer (182).
Pulse-driven PS externalization appears to be an electrophoretic transport of PS
headgroups through hydrophilic pores. Switching the field polarity during a simulation
reverses the direction of this interleaflet electrophoresis.
28
Figure 3. Water dipoles aligned in the low dielectric permittivity of the
membrane interior during the formation of a hydrophilic pore. This is
an enlarged view of the center of Fig. 2b (simulation time 5.5 ns). The
field direction is from top to bottom, with the positive electrode at the
top.
Why such a large electric field? At a porating potential of 1 V the area occupied by pores
is less than 1% of the cell surface area(27). To insure that we observe pore formation
during a simulation we must increase either the simulated area (22 nm
2
) or the simulation
time — both computationally expensive — or, since pore density in the stochastic model
(178) is dependent on transmembrane voltage (Δψ
m
) (34), increase Δψ
m
.
Figure 2.3. Water dipoles aligned in the low dielectric permittivity of
the membrane interior during the formation of a hydrophilic pore. This
is an enlarged view of the center of Figure 2.2b (simulation time 5.5
ns). The field direction is from top to bottom, with the positive
electrode at the top.
29
On a microsecond scale the conductive pores produced by supraphysiological
transmembrane potentials clamp the voltage at approximately 1 V (185), but current
models allow higher potentials for at least a few nanoseconds (36, 176). A simulated ψ
m
of 3.4 V (450 mV/nm)
8
produces poration and PS translocation in less than 10 ns.
Simulations that represent the poration-driven drop in ψ
m
will be required to establish
how well this reflects the experimental situation. At 1.5 V (200 mV/nm) we see transitory
water defects but no pores, even in multiple repetitions, consistent with the exponential
dependence on Δψ
m
in the stochastic pore model.
These simulations of phospholipid bilayers, combined with evidence from fluorescence
microscopic imaging observations of living cells, demonstrate that a transmembrane
electric field can produce nanometer-diameter hydrophilic pores within nanoseconds and
that the anionic phospholipid PS is driven electrophoretically through the pores. Further
explication of the energetic and kinetic details of these processes, with systems that more
closely represent the complexity of a living cell membrane, will lead to improved
methods for electroporation in research and clinical applications, and to a better
understanding of PS translocation mechanisms in apoptosis and signal transduction.
30
2.3 Methods
2.3.1 Cell lines and culture conditions.
Human RPMI 8226 multiple myeloma cells (ATCC CCL-155) were grown in RPMI
1640 medium (Irvine Scientific, Irvine, CA) containing 10% heat-inactivated fetal bovine
serum (FBS; Gibco, Carlsbad, CA), 2 mM L-glutamine (Gibco, Carlsbad, CA), 50
units/mL penicillin (Gibco, Carlsbad, CA), and 50 μg/mL streptomycin (Gibco, Carlsbad,
CA) at 37 C in a humidified, 5% carbon dioxide atmosphere.
2.3.2 Pulse generators and pulse exposures.
For microscopic observations, cells were placed in a rectangular channel 100 μm wide,
30 μm deep, and 12 mm long, with gold-plated electrode walls, microfabricated with
photolithographic methods on a glass microscope slide. A USC-designed and fabricated
NanoPulser with fast recovery diode switching was mounted on the microscope stage for
delivery of 4 ns pulses directly to the microchamber electrodes in ambient atmosphere at
room temperature (100).
31
2.3.3 Fluorescence microscopy.
Observations of live cells during pulse exposures were made with a Zeiss Axiovert 200
epifluorescence microscope. For visualization of PS translocation with FM1-43
(Molecular Probes, Eugene OR; λ
ex
= 480 nm, λ
em
= 535 nm), cells were washed twice,
resuspended in growth medium containing 20 μM FM1-43, incubated for 20 minutes at
37 °C, and observed without additional washing. Real-time binding of annexin V was
observed using reagents from an apoptosis detection kit from Oncogene Research
Products (San Diego). Cells were suspended in RPMI 1640 medium containing 2.5 mM
CaCl
2
and incubated before pulse exposure with media binding reagent and annexin V-
FITC for 15 min in the dark according to the manufacturer's instructions, except that in
order to increase the sensitivity of the fluorescence microscopic imaging the annexin V
concentration was 16 times the value recommended for flow cytometric analysis. Images
were captured and analyzed with a LaVision Imager QE camera and software (LaVision,
Goettingen, Germany).
2.3.4 Molecular dynamics simulations.
The bilayer test system consists of 128 lipids, 2873 Simple Point Charge water molecules
(11), and 4 sodium ions. One of the two leaflets of the lipid bilayer includes 4
dioleoylphosphatidylserine (DOPS) molecules and 60 dioleoylphosphatidylcholine
32
molecules (DOPC). The other bilayer leaflet contains 64 DOPC molecules. CH
2
and CH
3
groups are treated as a united atom. Test systems were energy minimized using a steepest
descent algorithm and then equilibrated for 25 ns. Production simulations with applied
electric fields were run for times from 1 to 25 ns with a time step of 2 fs using periodic
boundary conditions. Bond lengths were constrained using LINCS (68) for lipids and
SETTLE (117) for water. Electrostatic interactions were calculated with a Particle Mesh
Ewald algorithm (42) using fast Fourier transforms with a cutoff radius of 1 nm for
coulombic and Lennard-Jones interactions. Solvent, DOPC, DOPS, and ions were
coupled separately to temperature baths of 300 K with a coupling constant of 0.1 ps
(12). Pressure was maintained at 1 bar also using a Berendsen algorithm with a relaxation
time of 1 ps. The system size was allowed to fluctuate anisotropically with a
compressibility of 4.5 x 10
-5
/bar in the z direction and 0 in the x and y directions to
maintain the bilayer shape throughout the simulation. Systems were subjected to an
external electric field from 0.2 V/nm and 0.8 V/nm. The phospholipid parameters were
derived from the OPLS-united atom based parameters of Berger et al. (14) with
additional parameters for the PS headgroup from Mukhopadhyay et al. (118) and newly
derived charges (226). The behavior of DOPC and POPS has been well-tested, and minor
differences in the PS head group for this study are not important. All simulations were
run using GROMACS version 3.2.1 (13, 105) on the University of Southern California
High Performance Computing and Communications Linux cluster
(http://www.usc.edu/hpcc/). Molecular graphics images were generated with VMD (80).
33
2.4 Acknowledgment
This work was supported by the Air Force Office of Scientific Research. Computing
resources were provided by the USC Center for High Performance Computing and
Communications (http://www.usc.edu/hpcc/). PTV and MJZ are supported by MOSIS,
Information Sciences Institute, USC. Work in DPT’s group is supported by NSERC. DPT
is an AHFMR Senior Scholar, CIHR New Investigator, and Sloan Foundation Fellow.
34
CHAPTER 3: PHOSPHATIDYLSERINE EXTERNALIZATION PART 2
Portions of this chapter are adapted from:
Vernier, P. T., M. J. Ziegler, Y. Sun, M. A. Gundersen, and D. P. Tieleman. 2006.
Nanopore-facilitated, voltage-driven phosphatidylserine translocation in lipid bilayers--in
cells and in silico. Physical biology 3:233-247.
3.1 Abstract
Nanosecond, megavolt-per-meter pulses — higher power but lower total energy than the
electroporative pulses used to introduce normally excluded material into biological cells
— produce large intracellular electric fields without destructively charging the plasma
membrane. Nanoelectropulse perturbation of mammalian cells causes translocation of
phosphatidylserine (PS) to the outer face of the cell, intracellular calcium release, and in
some cell types a subsequent progression to apoptosis. Experimental observations and
molecular dynamics (MD) simulations of membranes in pulsed electric fields presented
here support the hypothesis that nanoelectropulse-induced PS externalization is driven by
the electric potential that appears across the lipid bilayer during a pulse and is facilitated
by the poration of the membrane that occurs even during pulses as brief as 3 ns. MD
simulations of phospholipid bilayers in supraphysiological electric fields show a tight
association between PS externalization and membrane pore formation on a nanosecond
35
time scale that is consistent with experimental evidence for electropermeabilization and
anode-directed PS translocation after nanosecond electric pulse exposure, suggesting a
molecular mechanism for nanoelectroporation and nanosecond PS externalization:
electrophoretic migration of the negatively charged PS head group along the surface of
nanometer-diameter electropores initiated by field-driven alignment of water dipoles at
the membrane interface.
3.2 Introduction
A biological cell in a time-varying external electric field develops a time-dependent
potential across the cell membrane, proportional to the magnitude of the field and the
radius of the cell, and a function also of the specific dielectric properties of the various
components of the system (144). When a cell is exposed to a long electrical pulse ( μs to
ms) or a low-frequency AC signal (Hz to kHz), the energy is deposited primarily in the
external membrane, which has a capacitive charging time constant in physiological media
of tens to hundreds of nanoseconds, depending on the cell type (31, 162). Microsecond,
kilovolt-per-meter pulses produce conductive pores in the cytoplasmic membrane,
permitting migration of charged species and large molecules across the mechanical and
dielectric barrier presented by the lipid bilayer. Electroporation technology, which
operates in the low-field, long-pulse regime, is widely used to facilitate transport of
nucleic acids, pharmaceutical compounds, and other materials into the cytoplasm of
36
living cells (65, 116, 125, 149).
Fast-rising electric pulses shorter than the charging time of the external membrane, and
electric fields varying at high frequencies, are felt in every part of the cell (38, 170).
Nanosecond pulses and gigahertz frequencies bypass the plasma membrane capacitance
and reach the membrane-bound structures and biomolecular assemblies of the cytoplasm
and nucleus (159, 201), opening promising new possibilities for biophysical
investigations and therapeutic interventions (160). Intracellular responses of biological
systems to ultra-short (ns), high-field (MV/m) electric pulses, which can be of megawatt
intensity but do not produce significant thermal effects because the total energy delivered
is small (joules per gram), include calcium bursts (202, 219), vacuole
permeabilization(159, 188), and the induction of apoptosis (8, 9, 201). Cancer cell killing
in vitro and tumor growth inhibition and regression in vivo have also been demonstrated
(8, 127) pointing to potential clinical applications.
Because nanoelectropulses are by definition very brief events, identifying the mechanism
or mechanisms through which they produce effects of biological consequence requires
investigation of the earliest observable pulse-induced changes in cells: increases in
transmembrane potential(48), elevation of intracellular calcium levels (202, 219),
phosphatidylserine (PS) externalization (204), and small molecule permeabilization(205).
The latter two phenomena, which are closely related, are the subject of this report.
37
Although the primary target of nanoelectropulse perturbation in most studies has been the
cell interior, PS externalization — translocation of PS from its normal location on the
inner leaflet of the plasma membrane to the cell exterior — always accompanies high-
field, nanosecond electric pulse exposure (203). Confinement of PS to the cytoplasmic
face of the lipid bilayer is regulated by transport proteins (7, 32, 168, 200) and enforced
by the energy barrier of approximately 100 kJ/mol (75) that impedes the transport of
charged species, including phospholipid head groups, across the low dielectric
(hydrocarbon) interior of the membrane (146, 184). Scrambling of the normal asymmetry
of the phospholipid bilayer and the resultant exposure of PS to the external medium
serves several recognition, binding, and signaling functions in healthy organisms:
promotion of platelet aggregation and catalysis of clot formation during blood
coagulation (5), marking of apoptotic cells, which flags them for phagocytosis (43), and
mediation of intramembrane signal transduction in lymphocytes (41). Because PS
externalization is a discrete, well-defined event that is directly and immediately
associated with nanoelectropulse treatment and is physiologically significant, it provides
a useful focal point for efforts to characterize the interactions of pulsed electric fields
with biological systems.
In earlier studies we presented evidence from flow cytometry and fluorescence
microscopy that the potential that develops across the lipid bilayer during an electric
pulse drives PS externalization, and we concluded on the basis of these and more recent
38
observations that the most likely mechanism for nanosecond pulse-induced PS
translocation involves facilitated diffusion or electrophoresis of PS along the walls of
short-lived, nanometer-diameter membrane pores, which are themselves caused by the
increased electric field in the membrane (203). Molecular dynamics (MD) simulations are
consistent with this hypothesis (76, 208). In addition, we have recently reported
observations of the immediate permeabilization of the plasma membrane by nanosecond,
megavolt-per-meter pulses (205), experimental confirmation of pore formation caused by
nanosecond electrical pulses. This "nanoporation" can be distinguished from
conventional, long-pulse electroporation because it permits influx of the fluorochrome
YO-PRO-1 while still excluding observable quantities of propidium iodide (PI), the dye
used in standard tests for membrane poration.
We report here experimental observations and MD simulations of phospholipid
membrane permeabilization and PS externalization in nanosecond pulsed electric fields,
and we propose a mechanism for nanoelectropulse-induced PS translocation that is
consistent with these findings and those reported by others. This work has been guided in
equal measure by experimental data and the predictions of models of cells in electric
fields; our understanding of each has been improved and extended by the other. The
results we present here are compatible with first and second order analyses of the
electrophysical dielectric shell model of the cell (49, 95) and with recent implementations
of the stochastic pore model for electroporation (145, 178) in a transport lattice (56, 174,
39
176, 199) , and we show an intersection as well with MD simulations of phospholipid
bilayers in electric fields (77, 182, 193, 194, 208).
In previous MD simulations of PS translocation in an electric field (208) we represented
the transmembrane potential by applying a global electric field vector normal to the lipid
bilayer to each atom in the simulation volume, as has been done also for MD studies of
electroporation (182, 193), and we report here additional results using that method. A
potentially more flexible and realistic approach, especially for simulations extending for
several nanoseconds or more, generates a transmembrane voltage by establishing an
asymmetric ion distribution across the membrane (39, 60, 152) . This method, which
requires two separate bilayers in order to maintain charge separation, has been employed
successfully for the simulation of electroporation (61), and we describe here the results of
double bilayer ion imbalance simulations demonstrating the close linkage of electric
field-driven water dipole alignment, electropore formation, and PS translocation, all on a
nanosecond time scale.
40
3.3 Materials and methods
3.3.1 Cell lines and culture conditions
RPMI 8226 human multiple myeloma cells (ATCC CCL-155) and Jurkat T lymphoblasts
(ATCC TIB-152) were grown in RPMI 1640 medium (Irvine Scientific, Irvine, CA)
containing 10% heat-inactivated fetal bovine serum (FBS; Gibco, Carlsbad, CA), 2 mM
L-glutamine (Gibco, Carlsbad, CA), 50 units/mL penicillin (Gibco, Carlsbad, CA), and
50 μg/mL streptomycin (Gibco, Carlsbad, CA) at 37 C in a humidified, 5% carbon
dioxide atmosphere.
3.3.2 Pulse generator and pulse exposures
For live-cell, fluorescence microscopy, cell suspensions were placed in a rectangular
channel 100 μm wide, 30 μm deep, and 12 mm long, with gold-plated electrode walls,
microfabricated with photolithographic methods on a glass microscope slide(179). A
USC-designed and fabricated fast recovery diode switching NanoPulser (100) was
mounted on the microscope stage for delivery of 4 ns pulses directly to the microchamber
electrodes in ambient atmosphere at room temperature as previously described (204).
41
3.3.3 Fluorescence microscopy and microphotometry
Observations of live cells during pulse exposures were made with a Zeiss Axiovert 200M
epifluorescence microscope. For visualization of PS translocation with FM1-43
(Molecular Probes–Invitrogen, Eugene OR; λ
ex
= 480 nm, λ
em
= 535 nm), cells were
washed twice, suspended in growth medium containing 20 μM FM1-43, incubated for 20
minutes at 37 °C, and observed without additional washing. Real-time binding of annexin
V was observed using reagents from an apoptosis detection kit from Oncogene Research
Products (San Diego). Cells were suspended in RPMI 1640 medium containing 2.5 mM
CaCl
2
and incubated before pulse exposure with media binding reagent and annexin V-
FITC for 15 min in the dark according to the manufacturer's instructions, except that in
order to increase the sensitivity of the fluorescence microscope imaging, the annexin V
concentration was 16 times the value recommended for flow cytometric analysis. YO-
PRO-1 ( λ
ex
= 480 nm, λ
em
= 535 nm) and PI ( λ
ex
= 540 nm, λ
em
= 580 nm) were from
Molecular Probes (Invitrogen, Eugene, OR). Images were captured and analyzed with a
LaVision Imager QE camera and DaVis software (LaVision, Goettingen, Germany).
3.3.4 Molecular dynamics simulations
Single and double bilayer systems were assembled and simulated using GROMACS
version 3.2.1 (13, 105) on the University of Southern California High Performance
42
Computing and Communications Linux cluster (http://www.usc.edu/hpcc/). Test systems
were energy minimized using a steepest descent algorithm and then equilibrated for 25
ns. Production simulations were 1 to 25 ns with a time step of 2 fs using periodic
boundary conditions. Bond lengths were constrained using LINCS (68) for lipids and
SETTLE (117) for water. Electrostatic interactions were calculated with a Particle Mesh
Ewald algorithm (42) using fast Fourier transforms with conductive boundary conditions
and a real space cutoff radius of 1 nm for coulombic interactions and Lennard-Jones
interactions. Solvent, lipids, and ions were coupled separately to temperature baths of 300
K with a coupling constant of 0.1 ps (12). Pressure was maintained at 100 kPa also using
a Berendsen algorithm with a relaxation time of 1 ps. The system size was allowed to
fluctuate anisotropically with a compressibility of 4.5 x 10
-5
/bar in the z direction for
single bilayer simulations and in x, y, and z for double bilayer simulations ( ΔV <
1%). Phospholipid parameters were derived from the OPLS, united-atom parameters of
Berger et al. (14) with additional parameters for the PS headgroup from Mukhopadhyay
et al. (118) and newly derived charges (226), and we used the Simple Point Charge (SPC)
model for water (11). Molecular graphics images were generated with VMD (80).
3.3.5 Single bilayer
The single bilayer system consists of 72 lipid molecules, 4 sodium ions (counterions for
the anionic PS head groups), and 2873 SPC water molecules. One of the two leaflets of
43
the lipid bilayer includes 4 dioleoylphosphatidylserine (DOPS) molecules and 32
dioleoylphosphatidylcholine (DOPC) molecules. The other bilayer leaflet contains 36
DOPC molecules. The application of an external electric field from 200 mV/nm to 800
mV/nm was represented by applying an electric field vector to all atoms in the system,
resulting in an additional force q
i
E on each atom.
3.3.6 Double bilayer
The double bilayer system consists of 144 lipid molecules, 8 sodium ions, and 5176 SPC
water molecules. One of the two leaflets of each bilayer contains 4 DOPS molecules and
32 DOPC molecules. The other leaflet in each case contains 36 DOPC molecules. The
two bilayers separate two volumes of water — one between the two bilayers and one
external to the two bilayers (but continuous in the z direction because of periodic
boundary conditions, see Figure 3.1). A transmembrane potential is established by adding
sodium and chloride ions to the inner and outer volumes of water to create a charge
imbalance while maintaining overall charge neutrality. For example, adding 4 sodium
ions to the outer volume and 4 chloride ions to the inner volume creates a 4-charge
imbalance, distributed over the total area of the two bilayers. Ions were introduced into
fully equilibrated double bilayer systems from a table of randomly generated positions
and then an additional energy minimization was performed to properly associate the ions
with surrounding water molecules.
44
3.4 Results and Discussion
Figure 3.1. Membrane phospholipid bilayer simulation — double bilayer
with ion imbalance. The simulation contains 136 DOPC, 8 DOPS, 12 Na
+
,
4 Cl
-
, and 5176 H
2
O in a rectangular volume approximately 4.4 nm x 5.4
nm x 14.7 nm. The z-axis (normal to the membrane-water interface) is the
horizontal axis in this view. Atoms and atom groups: red – water oxygen;
white – water hydrogen; yellow – choline; teal – serine; orange –
phosphate; blue – glycerol and carbonyls; gray – lipid HC tails; black –
sodium; pink – chloride. Not all atoms are visible in a single two-
dimensional snapshot. The lipid bilayer leaflets facing the inner water
compartment each contain 4 DOPS molecules and 32 DOPC molecules.
The outer leaflets each contain 36 DOPC molecules. The inner water
compartment contains eight sodium counter-ions for the eight negatively
charged PS head groups and four chloride ions. The outer water
compartments, which are actually a single contiguous compartment
because of periodic boundary conditions, contain four sodium ions. The
membranes separate the chloride ions in the inner compartment from the
sodium ions in the outer compartment, creating a transmembrane electrical
potential. Each bilayer, from water interface to water interface, is about 4
nm in the z-direction.
45
3.4.1 Imaging observations of nanoelectropulse-induced PS externalization in living cells
For flow cytometry and assays that are not time-critical, PS externalization is
conventionally detected with fluorescein-tagged annexin V (annexin V-FITC), a protein
that binds specifically to phospholipid membrane surfaces containing PS (94), and we
have previously characterized PS translocation at early times (minutes) after
nanoelectropulse exposure using standard annexin V-FITC flow cytometry methods
(201-203). To avoid the delays associated with processing samples for flow cytometry
and observe directly the immediate effects of nanosecond electric fields on living cells,
we integrated a pulse generator and exposure chamber with a fluorescence microscope
(100, 179).
Although real-time fluorescence microscopy imaging of annexin V-FITC binding to cells
during and immediately following pulse exposure is compromised by the relatively slow
kinetics of annexin V binding to externalized PS and by the background fluorescence of
the tagged annexin V reagent in the medium, we are able to demonstrate with annexin V-
FITC the appearance of externalized PS at the anode-facing pole of cells within 1 s after
exposure to 4 ns, 8 MV/m electric pulses (Figure 3.2a–d).
46
This annexin V binding pattern is similar to the anodic fluorescence intensification that
occurs when cells are pulsed in the presence of FM1-43, a cationic styryl fluorochrome
that equilibrates between a weakly fluorescent form in the aqueous medium and a
strongly fluorescent form in the external leaflet of the cytoplasmic membrane (161).
Figure 3.2. Nanoelectropulse-induced membrane perturbation and
phosphatidylserine externalization. A fluorescent cap appears at the
anode-facing pole of RPMI 8226 human multiple myeloma cells with
both annexin V-FITC (a–d) and FM1-43 (e–h) within 1 s after
exposure to 4 ns, 8 MV/m electric pulses delivered at a 1 kHz
repetition rate, indicating phosphatidylserine translocation and possibly
other associated membrane structural rearrangements. To overcome the
more limited sensitivity of the annexin V-FITC observations imposed
by binding kinetics and background fluorescence, a 100-pulse dose was
used for the annexin V-FITC images; the FM1-43 response followed a
30-pulse exposure. Images are representative of the responses of
hundreds of cells individually observed and recorded.
47
Externalization of PS shifts the aqueous-lipid FM1-43 distribution toward the lipid
membrane, resulting in a fluorescence increase (161, 227) and providing a means for
visualizing PS translocation in living cells that is faster and more sensitive than annexin
V binding (204, 220, 222). Figure 3.2e–h shows the anodic cap of increased FM1-43
fluorescence that appears immediately after nanosecond pulse exposure.
With both annexin V and FM1-43, the fluorescence increase that indicates
nanoelectropulse-induced membrane perturbation and PS externalization occurs only in
the membrane at the anode pole of the cell, never at the cathode-facing portion of the cell,
and never in the cell interior.
Each of these two methods for imaging PS externalization has limitations. Because
annexin V binding requires calcium, it cannot be used to investigate the role of calcium in
nanosecond pulse-induced PS translocation (204). The background fluorescence of
annexin V-FITC substantially reduces the sensitivity of imaging detection. Annexin V
binding also appears to inhibit additional, subsequent PS externalization, and it interferes
with the normal lateral diffusion of PS in the membrane (106, 107). FM1-43 is faster and
more sensitive, but it lacks the molecular binding specificity of annexin V for
externalized PS.
48
Taken together, the imaging evidence from experiments with these two different reagents
is consistent with a pulse count- and field-dependent (203, 204) restructuring of the
phospholipid bilayer of the cell membrane, including externalization of PS at the anode-
facing pole of the cell, during exposure to nanosecond, megavolt-per-meter electric
pulses.
3.4.2 Imaging observations of nanoelectropulse-induced membrane permeabilization in
living cells
Influx of normally excluded fluorescent dyes into the cytoplasm follows long-pulse
electroporation, and mammalian cell permeabilization in low conductivity buffer has
been reported with pulses as short as 10 ns using PI as the indicator (119). Previous
efforts to detect membrane poration with intracellular PI fluorescence after nanosecond
electric pulse exposure of cells in growth medium produced negative results (201, 202),
but in recent experiments with higher fields (> 6 MV/m) and faster repetition rates (1
kHz) we have observed pulse-induced membrane permeability to YO-PRO-1, an
indicator of early membrane damage associated with apoptosis (82) and also of P2X
7
(purinergic) receptor channel activation (72). With even higher fields (8 MV/m) and
repetition rates (10 kHz), detectable levels of PI enter the cytoplasm (205). Fluorescence
microscope images demonstrating electropermeabilization of Jurkat T cells by 4 ns pulses
and the differential influx of YO-PRO-1 and PI are shown in Figures 3.2 and 3.3.
49
YO-PRO-1 and PI are both small molecules, roughly similar in size, doubly positively
charged at physiological pH. Although YO-PRO-1 and PI are two of the three primary
components of a commercially available apoptosis assay kit (along with Hoechst 33342),
the uptake mechanism for YO-PRO-1 in apoptotic cells is not well characterized. Our
observation that nanoelectropulse-induced cell membrane permeabilization to PI requires
higher fields and repetition rates than are necessary for YO-PRO-1 suggests that more
than the simple opening of hydrophilic pores may be involved. A field-driven activation
of P2X
7
-like channels could account for the YO-PRO-1 influx in our experiments, for
example, or it may be that the slightly larger and bulkier propidium moiety (f.w. 415)
Figure 3.3. Nanoelectropulse-induced membrane permeabilization.
Fluorescence microscopic images of individual Jurkat T cells in growth
medium containing 5.0 μM YO-PRO-1 (a–d) or 7.5 μM PI (e–h)
immediately before and 5 min after delivery of multiple 4 ns, 8 MV/m
pulses at 1kHz show YO-PRO-1 (a–b) but not PI (e-f) influx after a 30
pulse exposure. 100 pulses produce a much stronger YO-PRO-1
response (c–d) and a detectable influx of PI (g–h).
50
does not pass as readily through a nanopulse-permeabilized membrane as the slimmer
YO-PRO-1 (f.w. 376). On the other hand, the qualitative nature of our fluorescence
imaging analysis does not permit us to exclude the possibility that our system is simply
more sensitive to YO-PRO-1, so that our baseline level of detection for intracellular YO-
PRO-1 is lower than that for PI. The development of methods for calibrating
electropermeabilization and electroporation — quantifying the molecular influx of
substances into an electrically pulsed cell — would help to answer this question and
contribute to the largely uncharacterized dosimetry of therapeutic delivery of
pharmaceutical agents and genetic material to living cells.
3.4.3 Molecular dynamics simulations of an asymmetric phospholipid bilayer in a high
electric field
In previous MD simulations of DOPC-DOPS bilayers we showed formation of
nanometer-diameter hydrophilic pores within nanoseconds when a large electric field
vector is applied to every atom in the system, and we demonstrated the rapid,
electrophoretically driven translocation of PS along the walls of the nanopores (208). The
net electric field at any point in the bulk water in this system is relatively low — the sum
of the globally applied external field vector and the large induced polarization vector of
the water dipoles. The net field is much higher in the low permittivity lipid bilayer. At the
membrane interface the charged atoms of the phospholipid head groups, their counter-
ions, and the hydrating water molecules all contribute to the net field.
51
Profiles for the electric field — obtained by integrating the charge density across the
simulation volume — are shown in Figure 3.2-4a for external electric field vectors of
250, 450, and 650 mV/nm. Figure 3.2-4b shows the field in the membrane interior as a
function of the applied field, with additional points taken from simulations with different
values for the external field. The net field in the bulk medium, as mentioned above, is
much smaller than the field in the low dielectric interior of the bilayer. Note from the
contours of the electric field profiles that one cannot directly derive the transmembrane
potential, a key parameter in electroporation studies, from any single value for the field
anywhere in the system.
52
Figure 3.4. Membrane electric field versus external electric field
vector. (a) Electric field profiles across a DOPC-DOPS lipid bilayer
with three external electric fields — 250, 450, 650 mV/nm. The
midpoint of the bilayer is at about 4.1 nm. (b) Electric field in the
membrane interior (obtained by integrating the simulated charge
density across the simulation volume) as a function of the electric field
vector applied to every atom in the simulation.
53
In simulations of 25 ns duration with global electric field vectors applied to systems of
single bilayers, we observed pore formation 21 out of 24 times when the applied field
was 450 mV/nm or greater, and in no case did we observe pore formation when the
applied field was 400 mV/nm or less (208). Although the periodic boundary conditions of
these simulations constrain the extraction of the membrane voltage, we estimate that this
corresponds to a porating transmembrane potential substantially greater than 2 V, larger
than expected from experimental studies and from other models. We discuss this below.
To characterize the relation between the time required for pore formation and the applied
electric field we define hydrophobic and hydrophilic pores on the basis of the presence of
water molecules in the hydrocarbon tail region of the lipid bilayer and on the mean dipole
moment of the water that appears in that region. During a simulation, if the
transmembrane potential is high enough to lead to poration, one or more water molecules
sporadically enter the membrane interior and then exit within a few picoseconds.
Eventually a quasi-stable column of dipole-aligned water forms, always beginning in our
simulations at the side of the membrane that has the more negative potential and
extending within a few hundred picoseconds to reach the other side of the membrane (0
ps panel of Figure 3.2-5). The structure of the pore at this point corresponds to the
hydrophobic pore of some electroporation models (1, 54). If the driving electric field is
maintained, the hydrophilic phospholipid head groups adjacent to the incipient pore
54
rearrange and follow the water, surrounding the aqueous column to form a hydrophilic
pore (178, 216) within a few nanoseconds.
55
Figure 3.5. Resealing of hydrophobic pore. At the beginning (0 ps) of
this representative sequence from a simulation of a single DOPC-DOPS
bilayer, a dipole-aligned column of water molecules has spanned the
membrane in the presence of an electric field pointing from left
(positive) to right. Within 50 ps of the removal of the electric field the
water column has disappeared. By 500 ps the phospholipid head groups
that had begun to follow the water into the membrane interior have
begun to settle back into place. By 2500 ps the membrane surface has
returned to normal. Color scheme is the same as figure 1. Hydrocarbon
tails are omitted for clarity.
56
The mean dipole moment of the water in the pore, initially oriented in the electric field,
reaches a maximum value and then declines as the water associates with the phospholipid
head groups that form the wall of the new hydrophilic pore, and as additional water
molecules entering the enlarging volume of the pore become more "bulk-like". This mean
water dipole moment signature, a quick rise, then a decline to a non-zero value, marks the
hydrophobic pore stage of nanosecond electropore formation. We use this signature to
identify the first persistent presence of water in the low dielectric region of the membrane
as the time of pore formation. This can only be determined by looking backward in a
simulation. Sometimes, for reasons unknown at this time but obviously of interest,
incipient water columns form and then collapse, with a signature of a single spike of non-
zero water dipole moment in the membrane interior. We have not yet identified any
predictive elements associated with intrusive water defects, or any feature of the porating
membrane, for example head group dipole orientation, that represents a point of
commitment to the establishment of a membrane-spanning bridge of water. The process
appears to proceed as a statistical ensemble of events, with one or more energy barriers
lowered in the presence of an electric field of sufficient strength.
The hydrophobic pore is not stable. If the electric field is returned to zero before the
hydrophilic walls are assembled, the water quickly exits the hydrocarbon tail region, the
head groups which have begun to follow the water reassemble into a smooth layer, and
within a few nanoseconds the defect is repaired. The pore sealing sequence shown in
57
Figure 3.5 begins just after a single chain of dipole-aligned water molecules has spanned
the membrane. Within 50 ps after the field is switched off the water has completely left
the membrane interior, and in less than 2.5 ns the normal membrane configuration has
been restored.
3.4.4 Molecular dynamics simulations of an asymmetric phospholipid bilayer with a
porating transmembrane potential generated by charge separation
Although the global electric field vector approach provides a useful representation of a
cell membrane exposed to a pulsed external electric field for perhaps a few hundred
picoseconds, membrane electrostatics are dominated over longer time intervals by the
migration of charged species in the growth medium and the cytoplasm and the
rearrangement of phospholipid and water dipoles at the membrane-water interface(30,
31). In addition, the normal resting transmembrane potential of 50–100 mV in a living
cell is generated by charge imbalances maintained at the expense of metabolic energy by
ion pumps in the membrane. In order to carry out MD simulations with transmembrane
potentials generated by charge distribution rather than by a uniform electric field vector,
we employ a double bilayer with two water compartments, one containing an excess of
positive ions, the other containing a corresponding excess of negative ions, maintaining
overall charge neutrality (Figure 3.1). Keeping the periodic boundary conditions of the
58
simulation in mind, the situation is a direct analog of a biological cell, with the membrane
separating the medium from the cytoplasm.
Figure 3.6 shows mass density profiles across the left bilayer in Figure 3.1, before
poration, for different components of the system, and profiles of the electric field, the
electric potential, and the mean water dipole moment. The solid lines represent the zero-
field case (no ion imbalance), and the dashed lines represent a simulation with a 4-ion
imbalance — 4 excess Na
+
in the left (outer) compartment, 4 excess Cl
-
in the right
(inner) compartment. Profiles from the right bilayer (not shown) are essentially a mirror
image of Figure 3.6, with minor differences arising from random (thermal) variations in
the positions of the sodium and chloride ions which are the source of the electric field.
Electric field values are determined by integration of the charge density along the z-axis
of the system (the axis normal to the membrane surface); the potential is obtained by
solving Poisson's equation across the bilayer, applying zero Dirichlet boundary
conditions and using the bvp4c algorithm in Matlab, a finite difference method that
implements the three-stage Lobatto IIIA implicit Runge-Kutta formula (91). The mean
dipole moment is calculated with the GROMACS program g_h2order. A dipole moment
that points to the right in our diagrams is positive. Dipole moment data points for the
isolated water molecules that occasionally enter the membrane interior are excluded to
improve the clarity of the plot; the water dipole moment is set to 0 where the time-
averaged water density is less than 25 g/L. Distances along the z-axis are normalized to
59
the zero-field dimensions to facilitate comparison of the bilayer profiles with different
transmembrane potentials. With an eight-ion imbalance the membrane compresses
approximately 15% in the direction normal to the electric field.
60
Figure 3.6. Mass density, electric field, electric potential, and mean
water dipole moment for one bilayer in a double bilayer simulation
with charge separation. Profiles along the z-axis (normal to the
membrane-medium interface) are plotted for simulations with a
balanced ion distribution (solid lines) and with a 4-ion imbalance
(dashed lines); to create the imbalance, 4 Na
+
ions in the outer
compartment are separated from 4 Cl
-
in the inner compartment by the
DOPC-DOPS membrane. For simplicity only the left bilayer is shown
here. Left and right bilayers have minor, random differences but are
functionally equivalent in all respects.
61
The electric field profile at the membrane-water interface, extending inward to about 1.5
nm from the bilayer center, provides one view of the complexity arising from the
combined phosphate-nitrogen, carbonyl, and water dipoles(53, 135, 225), and from the
anionic carboxyl of PS, that must be addressed by ongoing efforts to map this landscape
experimentally (30, 129). Potential-sensitive (Stark effect) fluorescent dyes, which are
one of the primary tools for measuring membrane potential, may be expected, to the
extent that these plots are representative, to respond dramatically to atomic-scale shifts in
position in the lipid bilayer.
From Figure 3.6 we can estimate that the simulated transmembrane potential generated
by a 4-ion imbalance across a total membrane area of 48 nm
2
is about 3.1 V. This is close
to the potential we would expect across a simple planar capacitor of these dimensions
with a dielectric thickness of 5 nm and a relative permittivity 2.5 (V = q·d / ε·A), a
verification of the physical consistency of the simulation and the validity of the
integration method. (Note that the "average" membrane electric field obtained by
dividing the potential by the membrane thickness is a quantity of dubious utility, given
the contours of the electric field profile.)
The pre-poration mean water dipole moment profile in Figure 3.6 shows the alignment of
water molecules at the membrane-water interface (about 2.2 nm from the center of the
bilayer), with the hydrogen (positive) end of the water dipoles oriented toward the
62
membrane and the oxygen (negative) end facing the bulk solvent. Although the peaks
closer to the center of the membrane represent a much smaller number of water
molecules — note the water density profile at the top of Figure 3.6 — they contribute
significantly to the field in the membrane interior (154). The zero-field asymmetry in the
profiles in the two halves of the bilayer is associated with the presence of PS in the leaflet
on the right side of the diagram in Figure 3.6 and is not observed in DOPC bilayers.
The effect of an electric field normal to the membrane surface on this thermodynamically
favorable configuration is apparent from the 4-ion imbalance data (dashed line) in the
water dipole moment panel of Figure 3.6. At the positive (anode-facing) side of the
membrane (the left side of the plots in Figure 3.6), the water dipole alignment is
reinforced by the electric field. At the negative (cathode-facing) side of the membrane,
the magnitude of the mean water dipole alignment is reduced. This effect of the electric
field on the arrangement of water at the membrane interface — stabilizing on one side
and destabilizing on the other — is consistent with our observation that in MD
simulations of nanoelectroporation the pore is initiated almost exclusively from the
negative side, the side on which the electric field counteracts the water dipole alignment.
The pre-poration profiles in Figure 3.6 are qualitatively similar in all respects to those
produced in simulations with global electric field vectors instead of charge separation.
63
Poration and PS externalization in an asymmetric phospholipid bilayer with an ion
imbalance-generated transmembrane potential. MD simulations of pore formation and PS
translocation in double bilayers with membrane potential generated by charge
separation are similar to those with single bilayers where the transmembrane potential is
generated by applying a global electric field vector. The sequence in Figure 3.7 is taken
from the right bilayer in a double bilayer simulation with a 5-ion charge separation —
Na
+
in the right (outer) compartment, Cl
-
in the left (inner) compartment. Water intrusion
into the membrane is always initiated in our systems on the negative (cathode-facing)
leaflet, as shown in the 1.0 ns panel of Figure 3.7. The phospholipid head groups
immediately begin rearranging behind the water, which forms a dipole-aligned column
across the membrane in this simulation after 1.3 ns. By 2.0 ns we see a roughly
cylindrical structure that may be considered an incipient hydrophilic pore, lined by
phospholipid head groups and filled with water. The amphiphilic PC head groups move
randomly in and out of the pore area and along the pore wall and in general do not cross
the bilayer, but the anionic PS head groups migrate consistently in the electric field
across the pore toward the positive (anode-facing) leaflet. In the simulation shown in
Figure 3.7 a PS head group is visible in the pore by 5.0 ns and has translocated by 8.0
ns.
64
Poration was not observed during double bilayer simulations (18 independent 10 ns
simulations) when the system charge imbalance was 3 or less ( ≤ 3 separated ion pairs/48
nm
2
). Poration occurred in less than 10 ns in either the left or the right bilayer in 29 of 30
simulations when the ion imbalance was 4 or more. The greater the ion imbalance, the
Figure 3.7. Pore formation and PS translocation in one bilayer in a
double bilayer simulation with a 5-ion charge separation. The right
(outer) compartment contains excess Na
+
. The left (inner) compartment
contains excess Cl
-
. Water intrusion from the negative (excess Cl
-
) face
of the bilayer into the hydrophobic interior of the membrane is evident
after 1.0 ns. Phospholipid head groups follow the water, enlarging the
defect, and a water column (arrow) spans the membrane by 1.3 ns. At
2.0 ns a cylindrical array of head groups connects the two leaflets of the
bilayer, forming a nanometer diameter pore, and by 5 ns a PS head
group (circled) can be seen migrating along the nanopore wall, reaching
the other (positive, excess Na
+
) side after 8 ns. Color scheme is the
same as figure 1.
65
sooner pores are formed. At the highest transmembrane potentials (largest charge
imbalances) pores form in both bilayers. Pore formation appears to be driven entirely by
the electric field-induced rearrangements and realignments of dipoles and charges at the
membrane-water interface, and pore development is a molecular process, a step-by-step,
atom group by atom group, extension of a defect. Mechanical forces on the membrane
(flexure, tension, compression) are not required for poration in these simulations, beyond
the electrostatic torques, repulsions, and attractions expressed at the molecular level, nor
are large statistical fluctuations involving the entire transmembrane region (221) required
for the nanometer scale electroporation observed in our systems.
Once a pore has formed, and usually beginning while the pore is still developing, PS
translocates electrophoretically along the hydrophilic pore walls toward the leaflet at the
more positive potential, never in the reverse direction. Nanopore-facilitated
phosphatidylcholine (PC) translocation is a much slower, diffusive (non-electrophoretic)
process observed only rarely in our 10 ns double bilayer simulations.
In all simulations where poration occurs, the water intrusion and subsequent pore
development begin on the negative side of the membrane. This is true whether or not the
leaflet at negative potential contains PS and regardless of whether the positive ion
imbalance is in the outer compartment or the inner compartment of the double bilayer.
66
This observation differs from reports based on coarse-grained MD simulations of DPPC
(dipalmitoyl PC) and DPPC-DPPS bilayers (76, 77), where pore initiation was seen on
the anode (more positive) face of the bilayer and attributed to the realignment of the
phospholipid head group dipoles in the electric field. Although we cannot directly
compare these results because of the many differences in the simulations, the discrepancy
points to the need to calibrate simulations with experimental observations, a continuing
challenge.
After a hydrophilic pore has formed we observe the migration of sodium ions, and less
frequently chloride ions, through the pore, in the direction that decreases the ion
imbalance (which is of course also the direction favored by the electric field).
3.4.5 Equivalence of MD simulations of nanoelectroporation with transmembrane
potential generated by global electric field and charge imbalance
Although the two methods for generating the potential gradients that lead to poration are
different, the simulation results obtained with the application of an explicit electric field
vector to each atom are qualitatively and quantitatively equivalent to those obtained by
establishing a charge imbalance across the membrane. A comparison that is both
illustrative and functional is the relation between the electric field in the membrane
interior and the time it takes to initiate formation of a pore. We calculate the field in the
67
central portion of the bilayer from the charge density data collected before the formation
of a hydrophobic pore (as described above), using an adaptive quadrature method for the
integration. For charge imbalance simulations the data are the average of several runs; for
explicit field vector simulations each data point is the result of one run. Figure 3.8a
shows the inverse relation between the electric field and the poration time and at the same
time demonstrates in a summary manner how the two methods produce results that are
not distinguishable within the scatter of the data.
68
Figure 3.8. Equivalence of global field and charge imbalance methods
for generating transmembrane potential. (a) Pore formation time versus
electric field. The time elapsed from the beginning of a simulation until
the formation of a hydrophobic pore (dipole-aligned, membrane-
spanning water column) is plotted against the electric field at the center
of the lipid bilayer, calculated by integration of the charge density in
the 50 ps immediately preceding the appearance of the pore. Data from
global field simulations is indicated with a square symbol; data from
charge imbalance simulations is indicated with a circle symbol. (b)
Membrane interior electric field produced by global field and charge
imbalance methods. The charge imbalance is calculated by dividing the
number of separated ion pairs (Na
+
, Cl
-
) by the area of the double
bilayer (47.5 nm
2
).
69
70
In Figure 3.8b, to further facilitate the comparison, the charge imbalance and applied
field are plotted together against the resulting internal membrane electric field. Since we
never observe poration with an applied field less than 450 mV/nm or with an ion
imbalance less than 4 (0.084 elementary charge units per nm
2
), Figure 3.8b shows that
the threshold for poration in our simulations, gauged by the magnitude of the field in the
membrane interior, is approximately the same for the two methods — between 1.0 and
1.5 V/nm.
3.4.6 Water dipole moment at the membrane interface — cause or effect?
A molecular mechanism for electropermeabilization or electroporation has not been
clearly established, and bridging the gaps between theory, models, and experiments
remains a major challenge, but the growing body of evidence for membrane restructuring
and poration on a nanosecond time scale narrows the range of possibilities for candidate
mechanisms and places constraints on the physical boundaries of the process. In all of our
simulations of high-field, nanosecond poration and PS translocation we observe two
polarized responses. First, PS externalization occurs preferentially at the anode-facing
pole of the cell. We have demonstrated both experimentally and with the methods of MD
that this is a nanosecond-scale event that is consistent with interleaflet electrophoresis of
the PS anionic head group along the walls of hydrophilic membrane pores. The role of
the electric field in the translocation is clear. Second, and for this we have now only the
71
conclusions we can draw from MD simulations, the initiation and extension of water
defects into hydrophobic and then hydrophilic, nanometer-diameter pores occurs
preferentially on the cathode-facing (more negative) side of the bilayer. This process of
pore development also takes place within nanoseconds (it precedes and is required for PS
translocation) and involves apparently only a limited number of players: water,
phospholipid head groups, and the hydrocarbon phospholipid tails. Furthermore,
electroporation is observed even in octane bilayers (194), which suggests that the lipid
headgroups play a secondary role.
What is different about the water-membrane interface on the two sides of the bilayer in
an electric field, and how could this account for the porating penetration of water into the
membrane interior exclusively from the negative side? Water and phospholipid head
group dipoles, which are oriented by the structure of the bilayer and by entropic and
enthalpic considerations at the membrane interface (109, 154), are the most likely
candidates. We could not extract convincing evidence for a significant pre-poration
reorientation of the head group dipoles from our simulation data, but the effect of the
electric field on the mean water dipole moment at the interface before poration is striking.
The view presented in Figure 3.6 is expanded in Figure 3.9 to show more clearly the
effect of varying the transmembrane charge separation. The field-dependence of the
dipole reorientation and the very different situations at the positive- and negative-facing
72
leaflets are evident. The water dipole orientation at the membrane surfaces (±2.4 nm from
the bilayer center) is reinforced by the electric field at the more positive side of the
bilayer and destabilized at the more negative side. For the water molecules at the
boundary of the hydrophobic membrane interior (±1.3 nm from the bilayer center), the
effect is the opposite. The orientation of the interior waters is disturbed by the field at the
positive side of the membrane (where the dipole moment in zero-field conditions points
outward, toward the bulk water, the opposite of the water layer a few atoms away at the
aqueous boundary) and reinforced at the negative side.
73
Figure 3.9. Effect of transmembrane charge imbalance on the mean
water dipole moment profile across a DOPC-DOPS bilayer before
poration. The dipole moment is the time-averaged value from the
beginning of a simulation until 50 ps before the initiation of pore
formation. For reference the 10% and 90% water density profile points
(relative to bulk water) are marked with arrows near the z-axis.
74
Although we cannot propose at this time a straightforward analytical description of the
sum of these factors, we can state an intuitive and physically sensible interpretation of the
MD results represented by the profiles in Figure 3.9, extending earlier results and
conclusions (194). In an electric field, the water dipole orientation associated with a
stable phospholipid membrane surface is favored at the positive-facing leaflet and
disrupted at the negative-facing leaflet. In conjunction with this surface effect, the entry
of water into the phosphoglycerate region of the membrane interior, which is stabilized
by a dipole orientation opposite from that at the bulk water interface, is favored at the
negative-facing leaflet and disfavored at the positive-facing leaflet.
This water dipole moment hypothesis is consistent with our observation that the presence
of PS in the membrane does not appear to affect pore initiation or formation, and with a
reported lack of correlation of the potential required for electroporation with the charge
on the phospholipid head groups (37). We should also expect, to the extent that
molecular dipole orientation at the interface is a key factor in the pore initiation process,
that water might be the leading partner, since water dipole rotation is considerably faster
than phospholipid head group rotation. Missing from this analysis is any consideration of
the critical interaction of the intruding water molecules with the lipid hydrocarbon tails,
which must be left for future work.
75
3.4.7 Are the electric fields and transmembrane potentials in these MD simulations
realistic?
Experimental data (48, 111, 185), continuum electrical models (57, 145, 178) and coarse-
grained MD simulations (77) indicate that the transmembrane potential in living cells in
external electric fields is limited by the formation of conductive pores to a maximum of
about 1.5 V, even for pulsed fields lasting less than 10 ns, yet all-atom MD simulations
seem to require membrane voltages higher than this before significant poration occurs, as
shown here and elsewhere (61, 182, 194, 208). We suggest two likely sources for this
apparent discrepancy.
First, the number of atoms in our simulations corresponds to a relatively small membrane
area, < 25 nm
2
per bilayer, barely large enough to permit the molecular rearrangements
required for pore formation. To increase the probability of a poration event during a 10 ns
simulation we can either increase the area — computationally expensive — or lower the
energy barrier for pore formation by increasing the transmembrane potential (27, 177,
215). For these initial studies, resource limits led us to choose the latter.
Second, the membrane in our simulations consists only of the phospholipids DOPC and
DOPS. A living cell membrane contains a variety of lipids and numerous proteins and
other constituents, and these components may affect the energy threshold of pore
76
formation and the transmembrane potential of a porated cell in ways that our simple
simulations cannot yet represent.
In addition, we concentrate on the membrane before poration. Limited computational
resources have prevented us from modeling larger areas and volumes for longer times,
which would permit us to take into account, as continuum models and experimental
observations do, the effect on transmembrane field and potential profiles of opening
conductive pores in the membrane during the application of the electrical pulse.
Finally, accurate simulations of biological cells and membranes in very short pulsed
electric fields must incorporate the relative "amplification" of the field in the membrane
that results from the 6–8 ps relaxation of water in the water-membrane-water dielectric
stack (23, 163), producing an electric field in the low dielectric membrane as much as 40
times greater than the field in the aqueous medium. For electric pulse durations greater
than a few nanoseconds, the migration of charged species in the cytoplasm and external
medium to the dielectric barrier of the cell membrane results in an amplification of the
external field by a factor of a thousand or more in the membrane, which permits us to
ignore the smaller dielectric effect, but for fast-rising, nanosecond pulses it should not be
neglected (59, 95, 196). For a 10 MV/m pulse for example, assuming relative
permittivities of 80 and 2 for water and membrane, respectively, the field in the
membrane will follow the rising edge of the external field to 400 MV/m.
77
3.5 Conclusions and Outlook
We have provided evidence from experimental observations with living cells and from
MD simulations to support the hypothesis that the PS externalization that immediately
follows exposure to nanosecond, megavolt-per-meter electric fields results from
nanopore-facilitated, electrophoretically driven translocation of the anionic PS head
group from the inner leaflet of the plasma membrane to the external face of the cell. We
have also shown that comparable results can be obtained from two different methods for
generating porating transmembrane voltages in MD simulations: 1. application of a
global electric force term to all atoms in the simulation; 2. separation of charges across
the membrane. Applying both of these methods to simulations of asymmetric DOPC-
DOPS bilayers in electric fields, we find that with sufficiently high electric fields,
nanometer-diameter hydrophilic pores form within nanoseconds, always beginning on the
face of the membrane at the more negative potential, providing a path for the
electrophoretic migration of PS across the dielectric barrier of the membrane core. PS
translocation is not observed in the absence of poration. Although the precise molecular
mechanism for nanosecond electroporation has not yet been elucidated, our simulations
suggest that the realignment of water dipoles at the membrane surface and at the level of
the phospholipid glycerol backbones is a primary event in this process.
78
The data and conclusions presented here are consistent with previous experimental
investigations of nanoelectropulse-induced PS externalization and with both continuum
models and MD simulations of membrane restructuring in high electric fields, and it is
clear from all of this that megavolt-per-meter fields, even when delivered in very short
electric pulses, permeabilize and porate the cell membrane, and perhaps affect it in other
ways.
Among many unanswered questions relating to mechanisms for interactions of
nanosecond, megavolt-per-meter electric pulses with biological systems, one that may be
more amenable to experimental investigation is the reason for the different
electropermeabilization responses observed with YO-PRO-1 and PI. If nanoelectropulse-
induced YO-PRO-1 influx is simply the result of the formation of pores through which
YO-PRO-1 passes more easily than PI, then understanding how two similar molecules
are treated differently by what we believe to be a relatively simple structure is likely to be
informative regarding the size of the nanosecond pulse-induced pores and their ability to
permit the passage of charged species. It is also possible, and more intriguing, that the
YO-PRO-1 entry we observe after pulse exposure is the result of the activation of
purinergic receptor channels like the P2X
7
receptor channels (72, 212), in which case we
may be flipping a membrane protein switch with potentially nanosecond precision, and
this may lead us to the still unknown mechanism(s) for nanoelectropulse-induced
apoptosis.
79
The ability to remotely manipulate the distribution of phospholipids in the membrane, the
primary focus of this work, suggests potentially fruitful explorations at several levels. As
a relatively simple yet physiologically relevant biophysical system, the asymmetric
phospholipid bilayer can be both modeled and experimentally tested, and it can be
extended to more complex assemblies of lipids, proteins, and polysaccharides. A better
understanding of nanoelectropulse-induced PS translocation may have ramifications in
the clinic as well. Systems for the nondestructive marking of cells with the semaphore of
externalized PS are under development for targeting specific cells and even subcellular
regions with a combination of pulse amplitudes, durations, and patterns, and optimization
of these therapies will depend not only on feedback from the clinic but also on the
strengthening of the fundamental knowledge on which this technology is constructed.
80
3.6 Acknowledgment
We gratefully recognize Andras Kuthi and Tao Tang for pulse generator engineering,
Christin Chong and Meng-Tse Chen for cell culture and technical assistance with
fluorescence microscopic imaging, and Nathan Baker for Matlab charge density
integration routines, and we thank James C. Weaver and Arieh Warshel for stimulating
and valuable discussions.
This work was made possible in part by support from the Air Force Office of Scientific
Research. Computation resources were provided by the University of Southern
California Center for High Performance Computing and Communications
(www.usc.edu/hpcc/). PTV and MJZ are supported by MOSIS, Information Sciences
Institute, Viterbi School of Engineering, University of Southern California. DPT is an
AHFMR Senior Scholar, CIHR New Investigator, and Sloan Foundation Fellow. Work in
his group is supported by NSERC.
81
CHAPTER 4: THE MINIMUM PORATING ELECTRIC FIELD
Portions of this chapter are adapted from:
Ziegler, M. J., and Vernier, P. T. 2008. Interface water dynamics and porating electric
fields for phospholipid bilayers. Journal of Physical Chemistry Submitted
4.1 Abstract
Lipid bilayers, normally a barrier to charged species and large molecules, are
permeabilized by electric fields, a phenomenon exploited by cell biologists and
geneticists for porating and transfecting cells and tissues. Recent molecular simulation
studies have advanced our understanding of electroporation, but the relative contributions
of atomically local details (interface water and head group dipole and counter-ion
configurations) and medium- and long-range electrostatic gradients and changes in
membrane structural shifts to the initiating conditions and mechanisms of pore formation
remain unclear. Molecular dynamics simulations of electroporation in several lipid
systems presented here reveal the effects of lipid hydrocarbon tail length and composition
on the magnitude of the field required for poration and on the location of the initial sites
of field-driven water intrusion into the bilayer. Minimum porating external fields of 260
mV nm
-1
, 280 mV nm
-1
, 320 mV nm
-1
, and 380 mV nm
-1
were found for 1,2-dilauroyl-sn-
glycero-3-phosphatidylcholine (DLPC), 1,2-dipalmitoyl-sn-glycero-3-
phosphatidylcholine (DPPC), 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine
82
(POPC), and 1,2-dioleoyl-sn-glycero-3-phosphatidylcholine (DOPC) bilayers,
respectively, and correlated most strongly with the bilayer thickness. These phospholipid
systems share several common features including a wide, dynamic distribution of the
head group dipole angle with the bilayer normal ranging from 0 to 155 degrees that is
only slightly shifted in a porating electric field, and similar electric field-induced shifts in
water dipole orientation, although the mean water dipole moment profile at the aqueous-
membrane interface is more sensitive to the electric field for DOPC than for the other
phospholipids. The location of pore initiation — at the anode- or cathode-facing leaflet
— varies with the composition of the bilayer and correlates with a change in the polarity
of the localized electric field at the interface.
4.2 Introduction
Biological cell membranes maintain the osmotic integrity of cells by selectively blocking
the passage of charged species and medium and large molecular weight compounds into
or out of the cell. Lipids, the defining constituents of the bilayer membrane, contain a
hydrophilic, polar head group and a hydrophobic, nonpolar tail. Assembled in a double
layer, with head groups facing outward, they form a barrier that separates the cytoplasm
from the extracellular medium and isolates internal membrane-bound compartments of
the cell. Many techniques have been developed for penetrating this barrier.
Electroporation utilizes electric pulses of varying amplitude, duration, and frequency to
83
permeabilize the cell membrane, temporarily allowing transport of normally excluded
substances into and out of the cell. Common laboratory and commercial electroporators
deliver kilovolt-per-meter pulses that last from microseconds to milliseconds, depending
on the application. The relatively large amount of energy (joules) associated with these
pulses can damage the cell membrane and cause necrosis and cell death if not carefully
controlled. In contrast, megavolt-per-meter, nanosecond pulses deliver several orders of
magnitude less energy and have profoundly different effects on cells, not only because
they are higher power but lower energy, but also because the high frequency components
of nanosecond pulses reach the interior of the cell rather than affecting primarily the
external membrane. Nanosecond, megavolt-per-meter pulses cause phosphatidylserine
externalization (201, 203), intracellular calcium release (202, 219), and the induction of
apoptosis(8, 9, 201), and may become a new tool in the cancer therapeutic arsenal (52,
127). Progress in investigations of this new technology, however, requires a better
understanding of the electric field-driven interactions among the constituents of the cell
membrane and solvent and solute species in the cytoplasm and external medium.
Molecular dynamics (MD) simulations, when they are anchored in a verifiably accurate
representation of a bilayer system (gauged by stability and area per phospholipid),
provide an atomic level view of the electroporation process that is beyond the reach of
current experimental methods and which can both guide and supplement experiment.
Simulations of phospholipid bilayers in porating electric fields have shown the
importance of interface water dipole orientation and the resultant charge density gradients
84
(193, 194, 207) during pore initiation and have provided a sequence of pore formation
events that is common to several different bilayer systems (182, 193, 194, 207). Other
studies have correlated simulation results with experimental observations of
phosphatidylserine translocation (77, 78, 207, 208) and ion transport (61, 90).
Pore initiation, the focus of this report, begins when a sufficient electric field is applied to
a lipid bilayer, disordering the aqueous-membrane interface and facilitating the intrusion
of water into the membrane interior. With increasing probability at higher fields, a water
defect can develop into a water column spanning the membrane — a hydrophobic pore.
MD simulations indicate that the initial steps required to form an electropore are similar
for lipid bilayers of varying compositions and even for octane slabs (no polar head
groups) (194). The force that drives water molecules into the membrane interior seems to
originate from localized field gradients that are dependent on the dynamic orientation of
water dipoles at the interface and to an unknown extent on the correlated interactions of
the phospholipid head groups with interfacial water (194). In the absence of an external
field interfacial water dipoles generally point toward the hydrophobic barrier presented
by a lipid bilayer or a non-miscible hydrocarbon like octane (25). That is, the water
molecules tend to orient with the hydrogens facing the barrier and the oxygens facing the
bulk water. This preferred orientation is due almost entirely to entropic and enthalpic
considerations arising from hydrogen bonding, since hydrogen bond associations (>10
kJ/mol) have a greater energy than van der Waals interactions and the energy associated
85
with aligning water dipoles with the local fields of 1 V nm
-1
(< 4 kJ/mol). When an
electric field is applied across a membrane, it reinforces the water orientation at the
anode-facing interface; the field and dipole moment directions are the same. But it
disrupts the interfacial water orientation on the cathode side as a result of the field-
generated torque on the water dipoles which tends to rotate them 180 degrees from their
zero-field alignment. This creates different local field environments on the two sides of
the membrane, a key to understanding the symmetries and asymmetries in pore formation
observed with systems constructed from various phospholipids (and octane).
Although molecular simulations of electroporation to this point have been instructive,
there are still discrepancies that need to be understood. In previous studies the simulated
transmembrane potentials associated with pore formation were much higher — as high as
3 V (207, 208) — than those typically encountered experimentally (185) and in other
models (36, 218). This discrepancy was explained as the result of a compromise between
limited computational power and the statistical nature of pore formation. Because
accurate all-atom MD simulations of fully hydrated lipid bilayers in the liquid phase
require at least 30-35 water molecules per lipid (122), practical computing resource limits
constrain bilayer system sizes to roughly 1000 lipids (~650 nm
2
) and simulation times to
fractions of a microsecond. This is not much larger than the area occupied by a
membrane pore and not much longer than the minimum time required for pore formation
(10), and it was suggested that higher fields and consequently higher simulated
86
transmembrane potentials were needed in order to observe pore formation events within
reasonable simulation times.
Direct comparisons of the values of MD-derived transmembrane potentials with those
obtained from other models and experimental observations, however, may not always be
valid. MD simulations of lipid bilayers are usually performed using periodic boundary
conditions to mitigate size effects, and this requires methods for calculating van der
Waals and electrostatic interactions that are compatible with this approach. Van der
Waals interactions decay rapidly with distance and can safely be ignored beyond a few
atomic diameters. Long-range electrostatic interactions are significant, however and must
be calculated. The three-dimensional Particle Mesh Ewald (PME) method is most often
employed because it is orders of magnitude faster than direct calculation and two-
dimensional PME methods, but it can introduce artifacts in the charge profile and has
been shown to be sensitive to the length of the system perpendicular to the plane of the
bilayer (175). The calculated value of the transmembrane potential is dependent on the
accuracy of the charge profile, and small errors in calculating long-range interactions can
result in non-trivial differences in the calculated potential that are magnified after
integration of Poisson’s equation. This has resulted in the use of Gaussian filters to
smooth the charge profile prior to integration (93). Another consequence of using
periodic boundary conditions, specifically the “tin-foil” boundary conditions commonly
employed in lipid bilayer simulations, is that the potential drop across the simulation box,
87
which is directly proportional to the length of the box in the field direction, cannot be
directly compared to the potential on the macroscopic charged electrodes employed in
experiments (20).
The results presented here circumvent some of these difficulties by comparing systems
based on the value of the external (applied) electric field rather than the transmembrane
potential. The external electric field, which represents the influence of free charge in a
real system, adds a force term (qE) to any particle in the system with a nonzero charge.
The minimum porating field, discussed in greater detail below, is used as a metric for
comparison of the different bilayer systems along with other structural properties of
several lipids.
The external fields in our simulations should not be confused with the electric fields
calculated for laboratory experimental apparatus by dividing the potential difference on a
pair of electrodes by the distance between them. The latter value, which is not directly
accessible in our simulations, represents in homogeneous systems the net electric field
that results from summing the external field arising from free charge (at the electrodes)
with the counteracting field arising from the dielectric polarization of the medium. The
magnitude of the external field in our simulations corresponds to the value of the field
that would be generated from the charge on the electrodes in the absence of an
intervening dielectric. The actual effective (net) electric field in the medium can be
88
calculated from the local water dipole moment (93, 190). A different approximation can
be obtained by considering the system from the viewpoint of continuum electrostatics and
reducing the applied field by a factor equal to the relative permittivity of the medium. For
example, given that the relative permittivity of SPC water is 61 (66), applying an external
field of 320 mV nm
-1
(MV m
-1
) in our simulations (the minimum porating electric field
for POPC (209)) produces an effective (net) field in the aqueous medium of 4.9 mV nm
-1
(MV m
-1
), comparable to the fields applied in experimental studies of nanosecond
electropermeabilization (206). This same “method of dielectrics”, assuming a
homogeneous membrane slab with a relative dielectric permittivity of 4 and a thickness
of 5 nm, yields a porating transmembrane potential of 0.4 V, again comparable to
experimental values (186). Although analytical limitations imposed by the “tin-foil”
boundary conditions mentioned above preclude a direct correspondence between these
simplified electromagnetics-based approximations (which also ignore the electrostatic
complexity of the membrane interface region (128)) and the values generated in the
simulations, the reasonable ranges for the fields and potentials in the comparisons
contribute to our confidence in our interpretation of the simulation results.
4.3 Methods
All simulations were performed using the GROMACS set of programs version 3.3.0 or
3.3.1 (13, 105, 198) on the University of Southern California High Performance
89
Computing and Communications Linux cluster (http://www.usc.edu/hpcc/). Lipid
parameters were derived from the OPLS, united-atom parameters of Berger et al. (14).
The Simple Point Charge (SPC) model (11) for water was chosen for compatibility with
the lipid models, computational advantages, and the ability to produce the correct area
per lipid. Systems were coupled to a temperature bath at 310 K with a relaxation time of
0.1 ps and a pressure bath at 1 bar with a relaxation time of 1 ps, each using a weak
coupling algorithm (12). Pressure was scaled semi-isotropically for lipid systems with a
compressibility of 4.5 x 10
-5
bar
-1
in the plane of the membrane and 4.5 x 10
-5
bar
-1
perpendicular to the membrane. For octane systems the compressibility perpendicular to
the membrane was set to 0.0 bar
-1
. Bond lengths were constrained using LINCS (68) for
lipids and SETTLE (117) for water. Short-range electrostatics and Lennard-Jones
interactions were cut off at 1 nm. Long-range electrostatics were calculated by the PME
algorithm (42) using fast Fourier transforms and tin-foil boundary conditions. Real-space
interactions were cut off at 1.0 nm, and reciprocal-space interactions were evaluated on a
0.12 nm grid with fourth order B-spline interpolation. The parameter ewald_rtol which
controls the relative error for the Ewald sum in the direct and reciprocal space was set to
10
-5
. Periodic boundary conditions were employed to mitigate system size effects.
Each lipid system contains 128 lipids and 4480 water molecules equilibrated for 25 ns.
Sample structures for DPPC and POPC were obtained from Peter Tieleman’s web site
(http://moose.bio.ucalgary.ca/) and modified for lipid and water content. DMPC and
DLPC structures were derived from the finished DPPC structure by removing carbons
90
from DPPC and re-equilibrating. DOPC structures were obtained from the Protein
Databank (http://www.rcsb.org/) and modified for lipid and water content. Octane
systems contained 313 octane molecules arranged in a slab formation with 3774 water
molecules. Production simulations were done in sets of 3 for each value of the electric
field. Each simulation was 25 ns with a 2 fs integration time step. The simulations were
extended up to 2 ns when hydrophilic pores were not fully developed after 25 ns.
4.3.1 Determination of the minimum porating field
The minimum porating electric field was determined from sets of triplicate simulations in
which the value of the external electric field was increased in steps of 20 mV nm
-1
until at
least 1 of 3 simulations formed a pore within 25 ns. This field was designated the
minimum porating (threshold) field. In cases where 0 of 3 simulations formed pores at a
given step and 3 of 3 simulations formed pores at the next higher step, the lower field was
marked as the minimum porating field. A minimum porating field for octane was not
determined. The highest field simulated for octane was 500 mV/nm.
In simulations with an applied field (or net dipole moment) that utilizes both periodic
boundary conditions and Ewald summation for long range electrostatics the potential
drop over the simulation box is proportional to the box length. Therefore, for comparative
purposes, simulations should maintain a similar box length. Additional simulations of
91
POPC with 5100 instead of 4480 water molecules, which increased the system size, did
not result in any appreciable effect on the results presented here.
4.3.2 Dipole-field angle
The dipole-field angle is the angle between the phosphorus-nitrogen (P →N) vector of the
PC lipid head group and the plane normal to the electric field.
4.3.3 Images
Molecular graphics images were generated with Visual Molecular Dynamics (VMD)
(80).
4.4 Results and Discussion
In the transient aqueous pore model, on which productive macroscopic and continuum
representations of electroporation have been constructed, pore formation is a stochastic
event (6, 178), and the relationship between the applied field and the amount of time
required to form pores is exponential (6, 10, 36, 176). Pores are formed, albeit at a very
low frequency and density, even in the absence of an external electric field, and are
92
thought to be at least partially responsible for the low basal rates of ion permeation
through cell membranes (62, 195, 217). These pores are reversible and do not result in a
sustained increase in the conductivity of the membrane. They must be distinguished from
pores formed via electroporation which can lead to either reversible or irreversible
electrical breakdown depending on the poration conditions (6, 178).
As a reference point for characterizing the electroporation of lipid bilayers of different
compositions and under different conditions, we have found it useful to find for each
system the minimum external electric field which can be expected to produce poration in
a simulation within a reasonable amount of time, where the number of lipids in the
simulation and the length of time that is “reasonable” are constrained by the computing
resources available to us. The threshold porating electric field is determined over a set of
simulations by increasing the value of the external electric field by 20 mV nm
-1
until a
pore forms within 25 ns in at least 1 of 3 independent simulations. In cases where no
pores were formed at one field step and pores formed in all three simulations at the next
step, we call the lower step the minimum porating (threshold) field.
93
Figure 4.1 - The average time to the formation of a hydrophilic pore
based on three independent simulations is plotted against the applied
external field. A hydrophilic pore is a quasi-stable structure defined as
a continuous water column surrounded by solvated lipid headgroups.
As noted above, the probability of pore formation is not a linear function of the electric
field. We observe that fields below the threshold do not result in pore formation even
after 100 ns. Increasing the field as little as 50 mV nm
-1
above the minimum porating
value can result in pore formation and disruptive reorganization of the bilayer in 10 ns or
less while fields 200 mV nm
-1
above the threshold result in near instantaneous collapse of
the bilayer (Figure 4.1).
Figure 4.1. The average time to the formation of a hydrophilic pore
based on three independent simulations is plotted against the applied
external field. A hydrophilic pore is a quasi-stable structure defined as
a continuous water column surrounded by solvated lipid headgroups.
94
Table 4.1. Relationship between threshold field and
bilayer properties
DLPC DPPC POPC DOPC
Threshold Field (V/nm)
0.26 0.28 0.32 0.38
Area per Lipid (nm
2
)
Experimental
Simulated
0.63
1
0.66
0.64
2
0.64
0.68
3
0.65
0.72
3
0.66
Bilayer Width (nm)
Experimental
Simulated
3.08
1
2.93
3.78
2
3.75
3.70
3
3.81
3.67
3
3.92
Hydrocarbon Thickness (nm)
Experimental
Simulated (sn-1 tail)
Simulated (sn-2 tail)
2.09
1
2.06
1.78
2.79
2
2.84
2.69
2.71
3
2.85
2.74
2.68
3
2.88
2.91
Carbons in Tail
12 16 18 18
Unsaturated Bonds
0 0 1 2
a
Kucerka, N.; Liu, Y. F.; Chu, N. J.; Petrache, H. I.; Tristram-Nagle, S.; Nagle, J.
F. Biophysical Journal 2005, 88, 245a.
b
Kucerka, N.; Tristram-Nagle, S.; Nagle, J. F. Biophysical Journal 2006, 90, L83.
c
Kucerka, N.; Tristram-Nagle, S.; Nagle, J. F. Journal of Membrane Biology
2006, 208, 193.
4.4.1 Minimum porating electric field varies with phospholipid properties.
Minimum porating fields for DLPC, DPPC, POPC, and DOPC are given in Table 4.1. An
increase in the threshold field correlates with an increase in the thickness of the
hydrocarbon interior of the bilayer (the mean distance between the C-1 atoms of the
anode- and cathode-side fatty acid tails, reported for both sn-1 and sn-2 tails) and the total
membrane thickness. No correlation is observed between simulated values for the
threshold field and area per lipid or temperatures between 280 K and 340 K (data not
shown). Discrepancies in simulated and experimental values for area per lipid result from
several factors. Simulated areas are obtained from the product of the simulation box
dimensions in the bilayer plane divided by the number of lipids in a leaflet; experimental
values take into account the volume of the tails as well as the headgroup volume. (122)
95
Area also varies considerably with temperature. A 15 K change in the temperature results
in a change in the area of 0.031 nm
2
to 0.033 nm
2
which is on the order of the difference
in area between lipids (133). Both all atom and coarse-grained models of lipids
underestimate the phase transition temperature of DPPC by 15–20 K (103, 110), which
means that the simulated areas presented here may be representative of experiments of
membranes at cooler temperatures.
For small pores the amount of work required to form a pore of radius r is approximated
by the equation below (1)
∆ 2 1 2
(4.1)
where F is the free energy, γ is the linear tension in the pore, σ is the tension of the
membrane, U is the voltage across the membrane, and C
m
is the membrane capacitance
defined by
(4.2)
where A is the membrane area and d is the hydrophobic width. The simple physical
model summarized by Equations 4.1 and 4.2 is valid for pores with a radius less than 1
nm (123) and is consistent with our observation that a thicker bilayer in a simulation
requires a larger field for poration. In order to form a pore of similar radius in a thicker
membrane, with its associated lower capacitance, a larger transmembrane potential and
96
therefore larger external field is required. The results are also consistent with the
expectation from basic kinetic considerations that water will have a lower probability of
penetrating a thicker bilayer.
Figure 4.2 – Pore formation sequence for DOPC and POPC T op) At
0.0 ns the lipid bilayer is in an equilibrium state. At 23.7 ns in this
simulation an aggregate of lipid head groups forms a cavity in the
water phase and a bump in the hydrocarbon phase. After 23.9 ns a
hydrophobic pore is formed. Water molecules further penetrate into the
hydrocarbon phase with their dipoles aligned with the electric field. At
24.2 ns solvated lipid head groups migrate into the hydrocarbon
interior. At 25.0 ns a hydrophilic pore with a stable interface has
formed. Bottom) Pore initiation for the DOPC bilayer is driven by
cathode-side water intrusion with little activity on the anode leaflet.
Figure 4.2. Pore formation sequence for DOPC and POPC Top) At 0.0
ns the lipid bilayer is in an equilibrium state. At 23.7 ns in this
simulation an aggregate of lipid head groups forms a cavity in the water
phase and a bump in the hydrocarbon phase. After 23.9 ns a
hydrophobic pore is formed. Water molecules further penetrate into the
hydrocarbon phase with their dipoles aligned with the electric field. At
24.2 ns solvated lipid head groups migrate into the hydrocarbon
interior. At 25.0 ns a hydrophilic pore with a stable interface has
formed. Bottom) Pore initiation for the DOPC bilayer is driven by
cathode-side water intrusion with little activity on the anode leaflet.
97
4.4.2 Common poration sequence
The sequence of pore formation events is similar for all lipids in this study. Figure 4.2
shows typical pore formation sequences near the threshold field for POPC and DOPC. At
the onset of pore formation a cluster of 6-10 lipids with inward-facing head group dipoles
appears, with head group positions closer to the bilayer center than the zero-field mean.
Water molecules associated with this defect (bump or dimple, depending on one’s
perspective) in the membrane surface gain access to the hydrocarbon interior. Water
dipoles in this low dielectric region are strongly oriented in the electric field, and the head
group dipoles along the curved surface of the bump also align with the electric field
(209). The stochastic motions of lipids and waters in these defects lead in some cases to
the disappearance of the bump and in others to a membrane-spanning column of water —
a hydrophobic pore. We have been unable to identify a definitive pre-poration event that
would allow us to predict which defects will lead to hydrophobic pore formation, but
once this stage is reached the probability of progression to full pore formation is very
high.
The time between the appearance of a pore-forming bump and the formation of a
hydrophobic pore is typically less than 200 ps. In the next few hundred picoseconds lipid
molecules and additional water molecules migrate into the hydrocarbon interior, forming
a cone of lipid head groups surrounding the water molecules. As the pore expands it
becomes more cylindrical and creates a new, highly curved, quasi-stable, lipid-water
98
interface — a hydrophilic pore. Hydrophilic pores continue to expand when the
threshold field is maintained. If the applied field is reduced to half of the threshold field
at the hydrophilic pore stage, the pore radius remains constant, and the head group
dipoles remain aligned with the field (which is tangential to the cylindrical pore wall)
(209). Reducing the field to zero results in head group dipole angle randomization and a
decreasing hydrophilic pore radius, with eventual collapse of the pore and repair of the
bilayer. Pore sealing takes much longer than pore formation. Hydrophilic pores in most
simulations persist for more than 100 ns in the absence of an external electric field and
allow for diffusive lipid flip flop.
4.4.3 Location of pore-initiating defects
Pore-initiating defects appear almost exclusively on the cathode-facing leaflet of DOPC
bilayers and on the anode side of octane slabs. In DLPC, DPPC, and POPC the head
group bumps develop on both the cathode and anode leaflets. When a pore forms in a
DOPC bilayer the anode leaflet remains near its equilibrium arrangement during the first
few nanoseconds of a pore formation event while both water and phospholipid head
groups enter the membrane interior from the cathode side (Figure 4.2, bottom panels).
Water molecules crossing the DOPC bilayer seek out waters near carbonyl groups on the
opposite leaflet (where their energy in the low-permittivity environment is minimized by
hydrogen bonding), which can lead to pores that are initially not perpendicular to the
bilayer surface, an indicator of the importance of hydrogen bonding for pore stabilization.
99
The pore initiation sequence for octane slabs is similar to the DOPC sequence except that
pore initiation is always on the anode side (194).
To investigate this bias further we simulated a chargeless DOPC system. In this system
the DOPC bilayer was simulated with the same parameters described in the methods
except that all the DOPC partial charges in the system were set to zero, so that the
molecules are entirely non-polar. Pore initiation in this chargeless DOPC bilayer occurred
in the same manner as in octane slabs — from the anode side. The hypothesis that
emerges from these simulations is that the more hydrocarbon-like a system is, the more
likely it is to produce pores from the anode side. We show here for our systems, as has
been reported previously for octane and DPPC (194), that this is associated with a force
asymmetry on the water molecules in the system.
Electric field gradient forces on interface water dipoles. The force profile on octane, a
point of reference for more complex lipid bilayers, is illustrated in Figure 4.3d (bottom
panel). The force F on any water molecule in an electric field gradient is given by the
equation
· (4.3)
where µ is the dipole moment of water, and E is the local electric field. Non-uniformities
in the local electric field arise from the distribution of charges in the system, particularly
from oriented dipoles such as interfacial water (and lipid head group dipoles in lipid
systems). In the absence of an external field the force on water molecules on each side of
the membrane is symmetrical. An external field places a torque on the water molecules at
100
the interface, which seek to establish an equilibrium alignment with the hydrogen atoms
adjacent to the interface in order to maximize hydrogen bonding, a phenomenon known
as the hydrophobic effect (25). The torque reinforces the water alignment on the anode
side but randomizes water molecules at the cathode side, which reduces the local field
gradient there and consequently the force on the dipoles. A significant driving force on
the water molecules remains at the anode side, which explains the anodic pore initiation
preference for octane slabs and mock hydrocarbon bilayers such as chargeless DOPC.
101
Figure 4.3. Membrane properties for 0.1 nm slices normal to the
membrane surface averaged over 25 ns for non-pore-forming
simulations are given for DPPC, POPC, and DOPC bilayers and octane
slabs without an external field and at the minimum field for poration
for lipids and at 0.50 V/nm for octane. Density profiles: lipid tails and
octane – magenta, water – blue, carbonyl oxygen – black, nitrogen –
green, phosphorus – red. Carbonyl groups are treated as a single united
atom. Lipid and water density values are divided by 10 to facilitate
comparison. All other profiles: solid blue – no field, dashed green –
minimum field for poration. Dipole moment and force are calculated
for water molecules.
102
103
In phospholipid systems the force asymmetry is not as pronounced, and the magnitude of
the driving force on water molecules is smaller than the force in octane systems. In
octane the water hydrocarbon interface is abrupt and water must reorient in a limited
amount of space resulting in a larger field gradient. The phospholipid water interface
104
provides a solvation buffer zone which reduces the field gradient by allowing water to
optimize hydrogen bonding over a larger volume, as indicated by the slope of the water
density profiles in Figures 4.3a-d (density panels). A smaller slope is associated with a
wider interface. Further insight may be gained by decomposing the field profile into its
constituent parts (Figure 4.4). When no external field is present, water and lipids balance
their dipoles, and no field exists in the membrane interior. In the presence of a porating
field, however, the interior membrane field is positive, a result of water dipole alignment.
Approximately 1 nm from the bilayer center the contribution to the field from the lipids is
zero; the total field tracks directly with water (Figure 4.4). Previous studies of dipole
potentials in lipids provide evidence that water molecules orient in response to the field
generated by the dipoles of the head groups (21, 53). A large enough external field,
however, will have a stronger influence over interfacial water orientation than the head
group dipoles. All of these factors combine to generate the local field and water dipole
moment profiles.
To confirm that the water dipoles indeed cause the positive field in the membrane
interior, an additional POPC simulation was performed in which the phosphorus and
nitrogen atoms of the lipid headgroups were frozen in their positions resulting from the
320 mV nm
-1
field (before pore formation). The field was then reduced to zero, and all
other molecules including water were allowed to move freely. The field in the membrane
interior dropped to zero within 100 ps after removing the external field. To determine
105
whether the lipid headgroups are responding to water or whether water is responding to
the lipid headgroups, the lipid phosphorus and nitrogen atoms were frozen in their
equilibrium, zero-field positions, and an external field of 320 mV nm
-1
was then applied
to the system. After 100 ps the field in the membrane interior increased to 0.51 V nm
-1
, a
sizeable fraction of the 0.84 V nm
-1
in the unrestricted simulations (Figure 4.3, POPC
field panel). Interface water dipoles oriented by the external electric field are the source
of the field in the membrane interior, and they also facilitate the small field-associated
rearrangement of the lipid headgroup dipoles.
Figure 4. The field profile for POPC without a field and at the
minimum field for poration (320 mV/nm) is decomposed into its lipid
and water components.
Figure 4.4. The field profile for POPC without a field and at the
minimum field for poration (320 mV/nm) is decomposed into its lipid
and water components.
106
Table 4.2. Phosphorus-nitrogen (phosphatidylcholine
head group dipole) vector without an external electric
field and at the minimum porating field
DOPC POPC DPPC
No Field
Anode
Cathode
22
23
12
11
12
12
Threshold Field
Anode
Cathode
18
24
8
15
8
13
Change
Anode
Cathode
4
1
4
4
4
1
4.4.4 Head group tilt not a factor in electroporation
It has been suggested that the tilting of the head group dipoles out of the plane of the
bilayer in an external electric field contributes to membrane defect formation and
poration (20, 77, 96, 187). Our simulations do not support this. To investigate the effect
of an external field on head group conformation in phospholipid bilayers we calculated
the head group dipole tilt angle distribution for DPPC, POPC, and DOPC in zero- and
minimum porating field conditions (Table 4.2). [Note: the POPC data in Table 4.2 and in
Figure 4.5 was presented in (209) and is included here for comparison.] At the threshold
field for poration the external field changes the mean head group dipole angle less than
about 4 degrees in the field direction. Despite the large applied field the head group
dipole angle distribution for each of the three phospholipids deviates minimally from the
zero-field case (Figure 4.5). The lack of any significant change in conformation and
107
mobility of the head group dipole suggests that head group dipole tilt is not a causative
component of the electropore formation mechanism.
Figure 5. The distribution of the phosphorus-nitrogen v ector (P →N)
angle is giv en for DPPC, POPC, and DOPC in the absence of an
external electric field and at the minimum porating field. The mean
P →N v ector angle is barely affected by the external field for all three
lipids. The most probable position for DPPC and POPC head group
dipoles is 90 degrees (co-planar with the membrane surface) while
DOPC head group P →N v ectors are tilted slightly into the water phase.
4.4.5 Nanoscale water dipole moment distribution drives electroporation
Figure 4.5. The distribution of the phosphorus-nitrogen vector (P →N)
angle is given for DPPC, POPC, and DOPC in the absence of an
external electric field and at the minimum porating field. The mean
P →N vector angle is barely affected by the external field for all three
lipids. The most probable position for DPPC and POPC head group
dipoles is 90 degrees (co-planar with the membrane surface) while
DOPC head group P →N vectors are tilted slightly into the water phase.
108
Water dipole moment profiles for these systems reveal several lipid-specific and shared
features. Common to DPPC, POPC, and DOPC is a small increase in the mean dipole
moment in the bulk in response to the electric field. Adjacent to the interface, where the
water density and therefore availability of hydrogen bonding partners is lower, there is a
larger response (Figures 4.3a-d dipole moment panels). (A positive dipole moment at the
anode leaflet represents the water dipole vector (O →H) pointing from the bulk toward
the interface. At the cathode-facing leaflet a positive water dipole moment points away
from the interface into the bulk.) For DOPC and POPC a secondary dipole peak closer to
the interior of the membrane is observed on the cathode side approximately 1 nm away
from the center. These water molecules are closer to the bilayer center, almost completely
in the hydrocarbon phase, and they do not hydrate the carbonyl oxygens, whose density
peak is 2 nm from the bilayer center (Figure 4.3a-c density panels). Typical hydrogen
bond lengths are 0.2 nm. The extra rigidity provided by the double bonds in DOPC and
POPC result in an increased amount of water in the hydrocarbon phase, a narrower
distribution of the carbonyl oxygen, and a larger cathode leaflet dipole moment (Figure
4.3a-c density and dipole panels). These structural changes when combined with an
external electric field result in significant changes in both the local electrostatic profile
and localized forces on water dipoles discussed further below.
The maximum force on a water dipole in each lipid system investigated here is between
15 and 20 pN, whereas the force on a water molecule in the octane slab is 8 pN even
though the field applied to octane was at least 120 mV nm
-1
greater than the fields applied
109
to the lipids. The force on a water molecule for both DOPC and octane is larger on the
anode leaflet (Figure 4.3b,d force panels). DOPC, however, initiates pores from the
cathode side (Figure 4.2) while octane initiates pores from the anode side
9
. One
explanation for the cathode-side poration preference observed with DOPC bilayers can be
derived from local electrostatics by carefully examining the field profiles of each lipid.
All lipids investigated here except DOPC exhibit a change in the polarity of the field at
both leaflets when the minimum porating field is applied. Because these field polarity
reversals are associated with a torque on interface water dipoles that tends to rotate them
relative to their zero-field orientation, we can reasonably expect that the free energy of
penetration of water into the membrane is dependent on the orientation of water entering
the membrane.
A second explanation for the exclusively cathodic poration of DOPC bilayers is the larger
area per lipid and number of water molecules per head group for DOPC membranes. The
two unsaturated DOPC tails result in a larger area per lipid, which allows more water into
the head group region than is the case for DPPC and POPC. As water penetrates toward
the bilayer center it becomes more susceptible to the external electric field, since the
availability of hydrogen bonding targets (and the local dielectric permittivity) is reduced
in the membrane interior. At the cathode leaflet the electric field works against water’s
natural orientation (by exerting a rotating torque on the water dipoles), destabilizing the
interface and aligning water with its oxygen closer to the membrane interior. The
electronegative oxygen atom of water is exposed directly to the positive field of the
110
membrane interior (water molecules from 0.6 nm–1.5 nm on the cathode leaflet (Figure
4.3b dipole and field panels)) and is drawn into the membrane interior as seen in Figure
4.2.
4.5 Conclusions
Small differences in the local electrostatic profiles and hydrophobic thickness have a
significant effect on the location of pore initiation and the minimum field for electric-
field-driven pore formation in lipid bilayers. DOPC, with two unsaturated fatty acid
residues, shows a strong preference for cathode-side pore initiation. Pores initiate from
both sides in DPPC (two saturated fatty acid tails) and POPC (one unsaturated residue).
Octane and an artificial test system of chargeless DOPC porate preferentially from the
anode side. These differences arise primarily from differences in the configuration of
water at the interface.
In a purely hydrocarbon system (octane and chargeless DOPC), where the hydrophobic
effect is the major determinant of water organization, a porating electric field reinforces
the water dipole orientation at the anode leaflet. At the cathode leaflet the electric field
works to disorganize the interfacial water configuration, eliminating the field gradient
and the resulting inwardly directed force on water molecules at the cathode leaflet. In the
absence of a driving force on the cathode leaflet, pores initiate on the anode side.
111
The interfacial electric field landscape is more complex for phospholipids, with their
polar head groups and the additional effects of variations in the hydrocarbon tails, which
changes the hydration level, the area per lipid (122), and the head group dipole angle
distribution. This leads to complexity in the poration behavior of PC lipids. MD
simulations show rapid changes in head group dipole tilt across a wide distribution of
angles, larger than might be suspected from a summary analysis of x-ray scattering or
NMR observations of PC lipids for which static values for the head group dipole angle
are reported. Despite these large and dynamic variations in head group dipole orientation,
the average value of the head group tilt angle is affected only slightly by a porating
electric field, and this does not seem to be a significant factor in electropore initiation.
We have not been able to identify a definitive pre-pore state from our simulations. Bumps
in the membrane surface consisting of aggregations of head groups and water molecules
may in fact constitute an essential pre-pore structure, but characterization of an
unambiguous set of structural prerequisites for electropore formation remains elusive.
The minimum porating field for electroporation provides several advantages as a
reference point for comparison of simulations of electroporation in lipid membranes and
should be used for this purpose, in our opinion, rather than the transmembrane potential.
The transmembrane potential in the simulations presented here is actually a dipole
112
potential generated from water and head group orientations. The absence of an ion
gradient between membrane compartments complicates comparisons of these simulated
transmembrane potentials to those obtained from experiments on living cells.
Furthermore, the simulated transmembrane potential is a calculated value derived from an
error-prone double integration of charge density. By using the minimum porating field,
an empirical property directly associated with a clearly defined simulation end point
which correlates well with actual values used in the laboratory, such problems are
avoided.
4.6 Acknowlegement
We thank Rainer Böckmann and Peter Tieleman for stimulating discussions. Computing
resources for this work were provided by the USC Center for High Performance
Computing and Communications. The authors are supported by MOSIS, Information
Sciences Institute, Viterbi School of Engineering, University of Southern California.
113
CHAPTER 5: STRUCTURE OF AN ELECTROPORE
Portions of this chapter are adapted from:
Vernier, P. T. and Ziegler, M. J. 2007. Nanosecond field alignment of head group and
water dipoles in electroporating phospholipid bilayers. Journal of Physical Chemistry B
111:12993-12996.
5.1 Abstract
To investigate the mechanism of biological cell membrane electroporation at the
nanosecond, nanometer scale, we track pore-forming lipids and water in molecular
dynamics simulations of a palmitoyloleoylphosphatidylcholine (POPC) bilayer in a
minimum porating electric field. Although the field-generated torque tilts the mean head
group dipole a few degrees away from its equilibrium, zero-field position relative to the
bilayer plane, this change in conformation does not appear to contribute directly to the
development of the pore-initiating aggregation of lipid head groups and water that leads
to the formation of a membrane-spanning hydrophilic pore. Field-directed rotation of the
head group dipoles in the plane of the incipient pore wall, in combination with water
114
dipole and solvation interactions at the aqueous-lipid interface, is one component in the
coordinated ensemble of electroporation events.
5.2 Introduction
Supraphysiological transmembrane electrical potentials cause rapid increases in
biological membrane conductance and permeability (92, 185). Although electrical
measurements (10) and fluorescence microscopy (206) indicate that permeabilization can
occur in less than 10 ns, the molecular details are not well established. Continuum
electrophysical models are consistent with nanosecond electropermeabilization (57, 87),
but visualization and real-time analysis are challenging at this time scale (48).
Molecular dynamics (MD) simulations of phospholipid bilayers provide a framework for
experimental design and interpretation in the absence of direct observation (61, 182, 194,
208). Results presented here show the cooperative, electric field-driven interactions of
water and head group dipoles in pre-pore lipid-water defects and in the subsequent
development of hydrophobic and hydrophilic pores (207).
A stabilizing role for the field-aligned head group P →N dipole, which lies nearly flat in
the plane of a hydrated, fluid-phase phospholipid bilayer, was proposed in early studies
of electroporation (125). More recently it has been suggested on the basis of coarse-
grained models(77, 96) that poration defects are initiated by substantial changes in the tilt
of the head group dipole relative to the bilayer plane, but in all-atom MD simulations
115
electroporation is observed even in a water-octane system (192, 193), indicating that lipid
head group electrostatic interactions are not an essential component of the poration
mechanism. Water dipole rearrangements at the membrane interface in a porating electric
field are striking, however, and may be a driving factor in defect initiation in at least
some systems (90, 194, 207).
We focus here on lipid and water molecules directly involved in pore formation.
Following methods previously described (208), we find that a minimum porating electric
field (320 mV/nm for this system) normal to a bilayer containing 128
palmitoyloleoylphosphatidylcholine (POPC) molecules and 4480 waters produces a mean
head group dipole tilt of less than 5° from the zero-field value everywhere in the bilayer
before pore initiation. During pore formation, as the head groups migrate into the bilayer
interior, the P →N dipoles rotate in the plane of the pore wall to align with the electric
field, which is normal to the plane of the unporated bilayer but tangential to the plane of
the pore wall.
5.3 Simulation Methods
5.3.1 Molecular dynamics simulations
The bilayer test system consists of 128 lipids and 4480 Simple Point Charge (SPC) water
molecules (11). Each of the two leaflets of the lipid bilayer includes 64
116
palmitoyloleoylphosphatidylcholine molecules (POPC). The phospholipid parameters
were derived from the OPLS-united atom based parameters of Berger et al. (14). CH
2
and
CH
3
groups are treated as a united atom. Test systems were energy minimized using a
steepest descent algorithm and then equilibrated for 25 ns. Production simulations with
applied electric fields were run for 25 ns with a time step of 2 fs using periodic boundary
conditions. Bond lengths were constrained using LINCS (68) for lipids and SETTLE
(117) for water. Electrostatic interactions were calculated with a Particle Mesh Ewald
algorithm (42) using fast Fourier transforms with a cutoff radius of 1 nm for coulombic
and Lennard-Jones interactions. Simulations were performed under an NPT ensemble.
Solvent and lipid were coupled to a 310 K temperature bath with a coupling constant of
100 fs (12). Pressure was maintained at 1 bar with a weak-coupling algorithm (12) with a
relaxation time of 1 ps. The system size was allowed to fluctuate semi-anisotropically
with a compressibility of 4.5 x 10
-5
/bar in the z direction and in the xy-plane to maintain
the bilayer shape throughout the simulation. Equilibrated area per lipid (POPC) averaged
over 25 ns is 0.651 nm
2
, comparable to measured values (0.683 nm
2
± 0.015 nm
2
) (97).
Applied external electric fields ranged from 150 mV/nm to 340 mV/nm. All simulations
were run using GROMACS version 3.2.1 (13, 105) or 3.3.1 on the University of Southern
California High Performance Computing and Communications Linux cluster
(http://www.usc.edu/hpcc/). Molecular graphics images were generated with VMD (80).
117
5.4 Results and Discussion
In order to avoid effects that might appear with excessively high electric fields we
defined a minimum porating electric field as the smallest external applied field necessary
to produce at least one porating event in three independent 25 ns simulations. By this
definition we found that 320 mV/nm is the minimum porating field for the POPC bilayer
system described here. No pores are formed within 25 ns with applied fields of 150, 220,
240, 260, 280, and 300 mV/nm.
Figure 5.1 shows the mean field-dipole angle with time for pore lipid and water
molecules in a representative simulation with a minimum porating electric field, where a
pore lipid is defined as a phospholipid molecule for which the phosphorus atom is at any
time during a simulation closer to the bilayer center than the interior minimum of the
phosphorus mass density distribution of the equilibrated, unporated bilayer. Before the
pore-initiating defect appears (22 ns after application of the field in this example), the
mean head group dipole angle remains roughly constant (tilted about 12° away from the
plane of the bilayer), and water intrusions into the membrane interior are sporadic. The
number of waters in the bilayer interior then begins to increase (Figure 5.1, open circles),
and, as a new, inwardly-directed hydrophilic surface of lipid head groups begins to form,
the P →N dipole vector for the pore-forming lipids shifts away from the plane of the
unporated region of the bilayer and toward the normal, aligned with the field in the plane
of the incipient pore wall. Figure 5.1 also shows the water dipoles of the hydrophobic
118
pore (207), highly oriented in the field, in the short-lived, membrane-spanning water
channel just before 24 ns (Figure 5.1, squares), after the first detectable rotation of the
lipid dipoles in the porating region (Figure 5.1, triangles).
0
20
40
60
80
100
21.0 22.0 23.0 24.0 25.0
Time (ns)
Dipole-Field Angle
0
50
100
150
200
250
300 Water Molecule Count
Figure 1. Dipole-electric field alignment during pore formation in a
POPC bilayer. Head group dipoles ( ) of pore-forming lipids (those
that become part of the pore during the simulation) remain in the
unporated bilayer plane (nearly perpendicular to the field) until a pre-
pore defect of lipids and penetrating (within 1 nm of the bilayer
center) water ( ) appears after 22 ns. P →N dipoles on the developing
pore surface rotate in the new pore wall plane to align with the
electric field, shown by the decrease in the head group dipole-field
angle. Deep interior (within 0.25 nm of bilayer center) water ( )
then appears for the first time, forms a transient, hydrophobic pore,
and solvates head groups in the pore wall. Open and closed circles are
the number of shallow ( ) and deep ( ) penetrating water molecules,
respectively.
Snapshots of this simulated membrane in cross-section (Figure 5.2) show how the
electroporating bilayer structure tracks the head group and water dipole alignment. No
significant disturbance of the membrane interface is apparent for the first 21 ns in this
Figure 5.1. Dipole-electric field alignment during pore formation in a
POPC bilayer. Head group dipoles ( ▲) of pore-forming lipids (those
that become part of the pore during the simulation) remain in the
unporated bilayer plane (nearly perpendicular to the field) until a pre-
pore defect of lipids and penetrating (within 1 nm of the bilayer center)
water ( □) appears after 22 ns. P →N dipoles on the developing pore
surface rotate in the new pore wall plane to align with the electric field,
shown by the decrease in the head group dipole-field angle. Deep
interior (within 0.25 nm of bilayer center) water ( ■) then appears for
the first time, forms a transient, hydrophobic pore, and solvates head
groups in the pore wall. Open and closed circles are the number of
shallow ( ○) and deep ( ●) penetrating water molecules, respectively.
119
simulation. (Pore initiation is a stochastic event; at higher fields pore initiation is more
frequent and tends to occur earlier in the simulations). At 22.6 ns an intrusive aggregation
of head groups and solvating and non-solvating waters appears, and by 23.8 ns a column
of water with no associated head groups (hydrophobic pore) spans the hydrocarbon
interior of the bilayer. A hydrophilic pore (lined with hydrated head groups) forms within
2 ns after the initial dimpling defect. Can we identify a causative role for head group
dipole tilt in this poration sequence?
120
Figure 2. Snapshots from the POPC bilayer electroporation simulation
used for Figure 1. (a) Membrane before poration. (b) Pre-pore intrusion
of water and lipid head groups into bilayer interior (lipid tails not
shown). (c) Hydrophobic pore. (d) Hydrophilic pore. Electric field is
directed normal to the bilayer, along the z-axis from left to right. P –
blue; N – red; H
2
O – pale red (O) and white (H); other phospholipid
atoms not shown.
21.1 ns 22.6 ns
23.8 ns 24.6 ns
a
b
c d
An electric field normal to a POPC bilayer exerts a torque on the P →N dipole, tilting the
N closer to the bilayer interior in one leaflet and farther away in the other, countered by
intra- and intermolecular bonding forces and steric factors. We calculated the field-driven
change in the head group dipole angle relative to the bilayer plane based on coordinates
from 128 POPC molecules at over 750 time points (at least 10 ps apart) in three separate
simulations. A 320 mV/nm field shifts the mean P →N angle from a zero-field value of
11.7° to 7.7° on the positive-facing leaflet and 14.3° on the negative side. No significant
change in the area per lipid — 0.65 nm
2
— is observed after the field is applied.
Although it is possible that this slightly modified conformation facilitates water intrusion
Figure 5.2. Snapshots from the POPC bilayer electroporation
simulation used for Figure 1. (a) Membrane before poration. (b) Pre-
pore intrusion of water and lipid head groups into bilayer interior (lipid
tails not shown). (c) Hydrophobic pore. (d) Hydrophilic pore. Electric
field is directed normal to the bilayer, along the z-axis from left to
right. P – blue; N – red; H
2
O – pale red (O) and white (H); other
phospholipid atoms not shown.
121
into the membrane interior, note from Figures 5.1 and 5.2 that even when a pre-pore
defect is clearly visible (22.6 ns) the angle of the head group dipoles for the lipids that
will form the pore 2 ns later remains relatively constant.
Figure 3. Distribution of head group P →N dipole angles relativ e to the
bilayer normal based on atomic coordinates from 128 POPC molecules
at ov er 750 time points (at least 10 ps apart) in non-porated sequences
from three separate simulations with and without a porating electric
field.
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0 20 40 60 80 100 120 140 160
Head Group Dipole Angle with Bilayer Normal
Relative Frequency
0 mV/nm
320 mV/nm
It is important to note that although these mean values are instructive, the head group
dipole angle varies rapidly and over a wide range, with or without an external electric
field. At any given time most P →N dipoles are within 10–15 degrees of the bilayer plane,
but substantial excursions from this mean, much larger than the shift induced by a
porating electric field, are common within a simulated population of lipids (Figure 5.3).
This population distribution reflects the dynamic Brownian jostling of the head group
Figure 5.3. Distribution of head group P →N dipole angles relative to
the bilayer normal based on atomic coordinates from 128 POPC
molecules at over 750 time points (at least 10 ps apart) in non-porated
sequences from three separate simulations with and without a porating
electric field.
122
dipole angle for each individual lipid in the bilayer, to which a porating electric field
adds only a small bias. We see no trends or characteristics in the dipole angle histories of
individual POPC molecules that would distinguish a pore-forming lipid from any other
(Figure 5.4).
Figure 4. P →N dipole angles with bilayer normal versus time for
representative lipid molecules from a 128 POPC, 4480 water simulation
with a minimum porating electric field (320 mV/nm) before pore
formation. The molecules in the upper two plots did not become part of
the pore; the molecules in the lower two plots are pore-forming lipids as
defined in the text. Angles are from atomic coordinates sampled every
100 ps. Solid line represents a 10-point moving average (1 ns).
0
20
40
60
80
100
120
140
160
180
05 10 15
Time (ns)
Angle with Bilayer Normal
Non-pore lipid 15
0
20
40
60
80
100
120
140
160
180
05 10 15
Time (ns)
Angle with Bilayer Normal
Non-pore lipid 21
0
20
40
60
80
100
120
140
160
180
05 10 15
Time (ns)
Angle with Bilayer Normal
0
20
40
60
80
100
120
140
160
180
05 10 15
Time (ns)
Angle with Bilayer Normal
Pore lipid 12 Pore lipid 13
Figure 5.4. P →N dipole angles with bilayer normal versus time for
representative lipid molecules from a 128 POPC, 4480 water
simulation with a minimum porating electric field (320 mV/nm) before
pore formation. The molecules in the upper two plots did not become
part of the pore; the molecules in the lower two plots are pore-forming
lipids as defined in the text. Angles are from atomic coordinates
sampled every 100 ps. Solid line represents a 10-point moving average
(1 ns).
123
To further test the idea that electric torque-driven twisting of head group dipoles
contributes only in a minor way, if at all, to electroporation of phospholipid bilayers, we
set all head group charges in an equilibrated POPC bilayer system to zero. Although a
larger electric field is required (1 V/nm), a "chargeless" POPC membrane forms water
defects and porates in the same way as a normal system, as previously reported for octane
(193).
Figure 5. Field-driven head group dipole alignment. (a) Projection of
P →N dipoles in a porated membrane on a plane parallel with a 320
mV/nm electric field directed along the pore axis (left to right). (b)
Projection of P →N dipoles in an unporated membrane on a plane
parallel with the bilayer surface with a 140 mV/nm electric field
applied in the same plane (left to right).
b a
P
-
Æ N
+
In contrast to the small average tilt of the P →N vector relative to the unporated bilayer
surface, we observe that head group dipoles rotate readily in the plane of the bilayer,
following a tangential field both in the pore wall and on an unporated membrane surface
Figure 5.5. Field-driven head group dipole alignment. (a) Projection of
P →N dipoles in a porated membrane on a plane parallel with a 320
mV/nm electric field directed along the pore axis (left to right). (b)
Projection of P →N dipoles in an unporated membrane on a plane
parallel with the bilayer surface with a 140 mV/nm electric field
applied in the same plane (left to right).
124
(Figure 5.5). In Figure 5.5a, which can be viewed as a cross section of a pore, the
alignment of the P →N vectors in the pore walls with the field from left to right is clear,
while at the same time the small field-induced tilt of the dipoles out of the plane of the
bilayer cannot be detected by a visual inspection of the image. As one would expect,
applying a field in the plane of the unporated bilayer, rather than normal to it, results in
an alignment of the head group dipoles, even when the field is 140 mV/nm, much smaller
than the magnitude of a threshold porating field (Figure 5.5b). In both cases — head
group dipole alignment in the pore walls and in the unporated bilayer plane — reversing
the field direction causes a rapid (within 1 ns) tracking rotation of the dipole vectors, and
removing the field results in immediate thermalization (complete within 400 ps) of the
mean P →N vector angle (not shown).
5.5 Conclusions
We see no evidence in these simulations for a significant contribution to pore initiation
arising from an electric field-driven tilt of head group dipoles out of the plane of the
bilayer. The results are consistent with the hypothesis that electroporation is the result of
a statistically infrequent, nanoscale, progressive alignment of a small population of water
and lipid dipoles, driven by electric force gradients at the membrane interface, for which
the energy barrier is lowered by an external electric field. Although the directed electrical
torque on head group dipoles in the pre-porated membrane surface does not seem to play
125
an essential causative role in poration (electroporation of lipid membranes occurs even in
the absence of charged head groups), the concerted association of polarized, pore-
forming water and lipid atomic groups during and after the porating electric field is
applied suggests that the magnitude and intramolecular equilibrium orientation of the
phospholipid head group dipole (that is, the P →N vector) could be one of several factors
that determine the frequency of pore initiation, the kinetics of pore growth, and the
stability (lifetime) of fully formed electropores. Longer and more area-extensive
simulations will be required to evaluate these possibilities.
In addition, our simulations predict an alignment of head group dipoles at any point in a
phospholipid bilayer system where there is a tangential external electric field of sufficient
magnitude. This will result in a localized segregation of aligned and non-aligned
molecular dipoles at the site of an electric field-induced pore, as described here, and also
around the unporated equator of a phospholipid vesicle or a biological cell between
parallel plate electrodes, where an applied external field will also be tangential
(electroporation occurs preferentially at the anode and cathode poles, where the field is
normal to the membrane, with very little permeabilization at points on the cell 90° from
the poles (51, 70)). These conformational alignments during exposure to pulsed electric
fields, which still require experimental confirmation, could have consequences for
membrane lipid-lipid and lipid-protein associations, and may contribute to the
perturbative effects of nanosecond electric pulses on biological systems (160, 206).
126
5.6 Acknowledgement
We thank D. Peter Tieleman for critical comments and valuable discussions. Computing
resources for this work were provided by the USC Center for High Performance
Computing and Communications (http://www.usc.edu/hpcc/). The authors are supported
by MOSIS, Information Sciences Institute, USC.
127
CHAPTER 6: EQUILIBRIUM PROPERTIES OF CALCIUM IN LIPID
BILAYERS
Portions of this chapter are adapted from:
Vernier, P. T., Ziegler, M. J., and Dimova, R. 2008. Calcium Binding and Head Group
Dipole Angle in Phosphatidylserine-Phosphatidylcholine Bilayers. Langmuir. Submitted.
6.1 Abstract
Manipulating the plasma membrane, gateway to the cell interior, with chemical and
physical agents for genetic and pharmacological therapy, and understanding the
interactions of lipid membrane components with proteins and other structural and
functional elements of the cell, requires a detailed biomolecular membrane model. We
report here progress along one path toward such a model — molecular dynamics
simulations of mixed, zwitterionic-anionic, asymmetric phospholipid bilayers with
monovalent and divalent cations. With phosphatidylcholine:phosphatidylserine systems
we identify temporal and concentration boundaries for equilibration of calcium with the
bilayer and saturation of the calcium capacity of the membrane, we demonstrate the
electrostatic- and entropic-driven associations of calcium and sodium ions with polar
groups in the bilayer interface region, expressed in spatial distribution profiles and in
changes in the orientation of the phospholipid head groups, and we describe for the first
time simulations of dynamic, calcium-mediated adjustments in the conformation of
128
mixed phospholipid species co-resident in the same leaflet of the bilayer. The results are
consistent with experimental observations and point the way to further refinement and
increased realism of these molecular models.
6.2 Introduction
The complex biophysical landscape of cell membranes reflects the functions of partition,
transport, and communication that define the instant physiological state of the cell.
Considerable effort has been devoted to understanding the detailed structure and dynamic
organization of this nanoscale barrier (84), using living cells, model membranes, and
computer simulations. Atomic-resolution molecular dynamics (MD) simulations of lipid
bilayers have contributed to our knowledge of the kinetics and energetics of molecular
interactions in the membrane and at the aqueous interface, and in recent years the
expansion of available computational power has enabled the practical analysis of larger
systems (hundreds of lipid and thousands of water molecules) for longer times (hundreds
of nanoseconds), so that we can now see beyond snapshots of conformation to events of
cell biological significance like the operation of ion channels (85) or the electroporation
of the membrane (194).
We present here the results of MD simulations of asymmetric phospholipid bilayers
composed of dioleoylphosphatidylcholine (DOPC) in one leaflet and a mixture of DOPC
129
and dioleoylphosphatidylserine (DOPS) in the other leaflet, with calcium, sodium, and
chloride ions in the aqueous medium, a simple representation of the configuration in
living cell membranes. The interaction of calcium with phospholipid bilayers in general,
and with anionic phospholipids like DOPS in particular, is of interest because of the role
of calcium ion gradients in a wide variety of cell physiological processes (16) and
because of the importance of anionic phospholipid localization in cell signaling (41, 43,
223) and structuring of the membrane (139, 181).
Experimental investigations of phospholipid bilayer stacks and vesicles in aqueous
electrolytes have shown that the properties of these systems are highly dependent on the
specific ions present and their concentration (113), and that measured values for the
calcium binding constant (2, 3, 17, 45, 46, 112, 113, 131, 167, 183),
calcium:phospholipid stoichiometry (3, 45, 67, 79, 112, 113), and head group dipole
angle (2, 18, 29, 166) can vary, in some cases over a wide range, for different
phospholipids and analytical techniques. Molecular dynamics simulations of
phospholipid bilayers have yielded results generally consonant with these observations
for binding constant (18, 19, 134, 139, 171, 224), stoichiometry(139), and head group
dipole angle (19, 74, 135, 153).
A recent experimental study (172) of calcium binding to PC:PS vesicles addresses the
spread of published values for the binding constant of calcium to phospholipid bilayers
130
with a systematic physical analysis under controlled conditions of ionic strength and
osmolarity. These results make clear the importance of the ion spatial distribution profile
and associated electrostatics for modeling calcium binding to phospholipid bilayers. The
ions that bind to the membrane interface partition from the local concentration adjacent to
the membrane, not from the bulk, and the differential can be quite large — a factor of 100
or more (167). The data of Sinn et al. suggest that it is this high surface calcium
concentration, driven by electrostatics, not a preferential attraction to PS, that accounts
for the greater observed binding of calcium to PC:PS vesicles relative to those composed
of pure PC. They also provide evidence that calcium binding to vesicles is entropy-driven
(and endothermic), and suggest that the entropy increase results from the loss of some
water from the calcium hydration shell for membrane-bound calcium ions and at the
same time a calcium-mediated dehydration of the lipid membrane (calcium replaces some
water in the membrane interface).
We place our simulations of mixed-composition phospholipid bilayers with calcium in
the context of these experimental results and of previously reported MD simulations of
lipid membranes with aqueous electrolytes. We report lipid area and head group dipole
angle effects, detailed atomic density profiles, radial distribution functions, and
equilibrium binding distributions of calcium and water in the membrane interfacial region
for bilayer systems containing: 128 phospholipids (DOPC:DOPS — 128:0, 118:10,
108:20), with all of the DOPS molecules in one leaflet of the bilayer in each case, to
131
represent the normal asymmetric distribution in cell membranes); 0, 10, and 100 ionized
CaCl
2
molecules; 0, 10, and 20 Na
+
counter-ions for the anionic DOPS residues; and
4180–4480 water molecules. For pure DOPC and for mixed, asymmetric DOPC:DOPS
systems we observe calcium-induced structural changes, binding kinetics and
distribution, and stoichiometry that are generally consistent with previous work, and we
describe for the first time the interacting and opposing effects of calcium and the anionic
phospholipid PS on the properties of DOPC:DOPS bilayers, including the PC head group
dipole angle.
6.3 Materials and Methods
6.3.1 Molecular dynamics simulations
All simulations were performed using the GroMACS set of programs version 3.3.1 (13,
105, 198) on the University of Southern California High Performance Computing and
Communications Linux cluster (http://www.usc.edu/hpcc/). Lipid parameters were
derived from the OPLS, united-atom parameters of Berger et al. (14). Additional head
group parameters and charges for phosphatidylserine were taken from Mukhopadhyay et
al. (118) or were provided by D. P. Tieleman (personal communication). The Simple
Point Charge (SPC) model (11) for water was chosen for compatibility with the lipid
models, computational advantages, and the ability to produce the correct area per lipid.
Systems were coupled to a temperature bath at 310 K with a relaxation time of 0.1 ps and
132
a pressure bath at 1 bar with a relaxation time of 1 ps, each using a weak coupling
algorithm (12). Pressure was scaled semi-isotropically with a compressibility of 4.5 x 10
-5
bar
-1
in the plane of the membrane and 4.5 x 10
-5
bar
-1
perpendicular to the membrane.
Bond lengths were constrained using LINCS (68) for lipids and SETTLE (117) for water.
Short-range electrostatics and Lennard-Jones interactions were cut off at 1.0 nm. Long-
range electrostatics were calculated by the PME algorithm (42) using fast Fourier
transforms and conductive boundary conditions. Real-space interactions were cut off at
1.0 nm, and reciprocal-space interactions were evaluated on a 0.12 nm grid with fourth
order B-spline interpolation. The parameter ewald_rtol, which controls the relative error
for the Ewald sum in the direct and reciprocal space, was set to 10
-5
. Periodic boundary
conditions were employed to mitigate system size effects.
6.3.2 Structures
A sample DOPC lipid structure was obtained from the Protein Data Bank
(http://www.rcsb.org/), and custom code was used to create a bilayer system with 128
DOPC and 4480 water molecules (35 water per lipid). The system was equilibrated until
the area per lipid was constant. Using custom code and the GroMACS function genion,
12 systems were created from the pure DOPC system, summarized in Table 6.1.
133
Table 6.1. System configurations
System Name
Outer (left)
Leaflet
Inner (right)
Leaflet
Ions
Water
DOPC DOPC DOPS Calcium Sodium Chloride
0 PS 0 Ca 64 64 0 0 0 0 4480
0 PS 1 Ca 64 64 0 1 0 2 4477
0 PS 10 Ca 64 64 0 10 0 20 4450
0 PS 100 Ca 64 64 0 100 0 200 4180
10 PS 0 Ca 64 54 10 0 10 0 4480
10 PS 1 Ca 64 54 10 1 10 2 4477
10 PS 10 Ca 64 54 10 10 10 20 4450
10 PS 100 Ca 64 54 10 100 10 200 4180
20 PS 0 Ca 64 44 20 0 20 0 4480
20 PS 1 Ca 64 44 20 1 20 2 4477
20 PS 10 Ca 64 44 20 10 20 20 4450
20 PS 100 Ca 64 44 20 100 20 200 4180
For systems with DOPS, custom code was used to substitute a DOPS for a DOPC by
modifying the head group and adding a corresponding sodium counter-ion. A short 10 ps
equilibration was performed to test each lipid substitution for bad contacts. For systems
with calcium the GroMACS function genion was used to place ions. genion intelligently
replaces a random solvent water molecule with an ion in a low energy environment and
does not require further equilibration. For example, in order to add CaCl
2
, three random
water molecules are removed and replaced by calcium and chloride ions.
6.3.3 Calcium binding criteria
Criteria for binding of calcium to other atom groups are similar to those used in
simulations of calcium binding to palmitoyloleoylphosphatidylcholine (POPC) (19). A
radial distribution function was created by calculating the distance between each calcium
134
ion and the phosphatidyl, carbonyl, and carboxyl oxygens in the system. Atoms closer
than 0.3 nm are considered bound. No binding was observed between calcium and
nitrogen or calcium and phosphorus.
Images. Molecular graphics images were generated with Visual Molecular Dynamics
(VMD) (80).
6.4 Results and Discussion
6.4.1 System equilibration time, area per lipid, bilayer thickness
Establishment of equilibrium configurations for the nine DOPC and DOPC:DOPS
systems investigated in this study, using area per lipid as an indicator of stability, occurs
within about 50 ns (Figure 6.1A). Calcium-lipid binding patterns continue to evolve,
however, until approximately 100 ns (Figure 6.1B). No further changes in these
properties were observed after an additional 100 ns (200 ns total). Area per lipid and the
corresponding bilayer thickness data for each system, time-averaged from 140 ns to 150
ns, are compiled in Table 6.2. Adding calcium ions to the solvent or PS to the bilayer
causes a decrease in area per lipid and a corresponding increase in bilayer thickness,
consistent with previous experimental observations (112) and simulations (19).
135
Figure 1. Time Evolution of Area per Lipid and Calcium Pairs (A) Area
per lipid as a function of time for a DOPC:DOPS (108:20) system.
Adding calcium to the system results in tighter lipid packing and a
decrease in the bilayer area. Data points every 10 ps. Heavy lines are 2.5
ns moving average. (B) Number of calcium-oxygen coordinated pairs
(calcium-phosphatidyl, calcium-carbonyl, calcium-carboxyl) as a
function of time for a DOPC:DOPS (108:20) system containing 100
calcium ions. Binding configuration stabilizes after about 100 ns.
0.56
0.58
0.60
0.62
0.64
0.66
0.68
0
20
40
60
0 50 100 150
Time (ns)
Number of Ca
2+
-
Coordinated Pairs
A
B
Area per Lipid
(nm
2
)
PC PO
4
108 DOPC, 20 DOPS
0Ca
10 Ca
100 Ca
PC C=O
PS COO
108 DOPC, 20 DOPS, 100 Ca
Figure 6.1. Time Evolution of Area per Lipid and Calcium Pairs (A)
Area per lipid as a function of time for a DOPC:DOPS (108:20)
system. Adding calcium to the system results in tighter lipid packing
and a decrease in the bilayer area. Data points every 10 ps. Heavy lines
are 2.5 ns moving average. (B) Number of calcium-oxygen coordinated
pairs (calcium-phosphatidyl, calcium-carbonyl, calcium-carboxyl) as a
function of time for a DOPC:DOPS (108:20) system containing 100
calcium ions. Binding configuration stabilizes after about 100 ns.
136
Table 6.2. Variation in bilayer area per lipid and thickness with number of calcium ions
and PS molecules.
Area per lipid (nm
2
)
Thickness, z
P-P
(nm)
PS 0 Ca 10 Ca 100 Ca 0 Ca 10 Ca 100 Ca
0 0.66 0.62 0.60 3.9 4.3 4.3
10 0.64 0.60 0.58 4.0 4.3 4.3
20 0.61 0.59 0.58 4.2 4.4 4.3
Mean values over the interval 140 ns to 150 ns of the simulations. Thickness is the
distance along the bilayer normal between the planes defined by the phosphorus atoms in
each leaflet of the bilayer. Standard deviations for the area per lipid values are < 0.01.
6.4.2 Calcium distribution in DOPC bilayers
Calcium distribution profiles in pure DOPC bilayer systems containing 10 (Figure 6.2A)
and 100 (Figure 6.2B) calcium ions show a single Ca
2+
peak located between the head
group phosphate and the backbone carbonyl peaks, consistent with neutron diffraction
measurements (67). In the 10-calcium system all of the calcium is found in this peak, but
in the 100-calcium system a significant concentration of calcium in the bulk is observed,
indicating that the calcium capacity for this 128-DOPC system saturates somewhere
between 10 and 100 calcium ions.
137
Figure 6.2 Number density (nm-3) profiles for selected components of
DOPC and DOPC:DOPS (108:20) bilayer systems containing 10 (A, C)
and 100 (B, D) calcium ions and corresponding counter-ions (Table 1).
For convenient scaling, calcium density is multiplied by 10, DOPS
carboxyl density is multiplied by 2, and water density is divided by 2.
The phosphate and carbonyl profiles are for the DOPC-associated
residues. To simplify the diagram the DOPS phosphates and carbonyls,
which coincide almost exactly with their DOPC counterparts, are not
shown. The plots are centered on the density minimum at the midpoint
of the bilayer.
138
0
5
10
15
-4 -3 -2 -1 0 1 2 3 4
Distance From Bilayer Center (nm)
0
5
10
15
0
5
10
15
128 DOPC, 0 DOPS, 10 Ca
128 DOPC, 0 DOPS, 100
Ca
0
5
10
15
108 DOPC, 20 DOPS, 10 Ca
108 DOPC, 20 DOPS, 100 Ca
Number Density (nm
-3
)
A
B
C
D
Water/2
DOPC C=O
DOPC PO4
Ca·10
DOPS COO·2
COO
-
PO
4
-
Ca
2+
Ca
2+
Ca
2+
C=O
H
2
O
139
6.4.3 Calcium distribution in DOPC:DOPS bilayers
Calcium distribution profiles in DOPC:DOPS (108:20) systems containing 10 (Figure
6.2C) and 100 (Figure 6.2D) calcium ions show a single Ca
2+
peak in the pure DOPC
leaflet of the bilayer (left side of diagrams in Figure 6.2), but the distribution is shifted to
the DOPS carboxyl group in the mixed DOPC:DOPS leaflet (right side of Figure 6.2C
and 6.2D). In the 10-calcium, mixed-lipid system, as with pure DOPC, all of the calcium
is localized in the bilayer. In the 100-calcium system again there is a significant calcium
concentration in the bulk, indicating that the calcium capacity for this mixed, 128-lipid
DOPC:DOPS bilayer, like that of the pure DOPC system, saturates between 10 and 100
calcium ions.
6.4.4 Calcium and chloride interfacial distribution in DOPC and DOPC:DOPS
In the calcium-saturated DOPC and DOPC:DOPS systems, the calcium concentration
reaches a maximum in the phosphate or carboxyl region of the bilayer interface, drops to
a minimum at the solvent boundary, then rises again in the bulk medium (Figure 6.2B,
6.2D). This MD-simulated manifestation of the electric double layer (18, 135, 171) can
be seen in more detail in Figure 6.3, which also elucidates the differences between the
pure DOPC leaflet and the mixed DOPC:DOPS leaflet. At the DOPC-only interface
(Figure 6.3A) the phosphate-localized calcium peak, the counteracting calcium and
140
chloride profiles, and the exclusion of sodium from the head group region are
straightforward. At the aqueous boundary of the mixed DOPC:DOPS leaflet on the other
side of the bilayer (Figure 6.3B) the calcium and chloride peaks are distinct, but the
calcium profile in the head group region traces a landscape complicated by the presence
of the serine carboxyl groups, and the calcium is accompanied by a significantly greater
amount of sodium than in the DOPC-only leaflet.
141
FIGURE 3 Number density (nm
-3
) profiles for the interface regions of
the DOPC-only (A) and DOPC:DOPS (B) leaflets of a DOPC:DOPS
(128:20) asymmetric bilayer. The DOPC-only leaflet contains 64
DOPC. The DOPC:DOPS leaflet contains 44 DOPC and 20 DOPS. (A)
and (B) are expanded views of the left side and right side, respectively,
of Fig. 2D. In addition to the scaling applied in Fig. 2, the chloride and
sodium densities are each multiplied by 10.
0
5
10
15
12 34
0
5
10
15
-4 -3 -2 -1
Ca
2+
Distance From Bilayer Center (nm)
A
B
Water/2
DOPC PO4
Ca·10
DOPS COO·2
Cl·10
Na·10
Cl
-
Ca
2+
Cl
-
COO
-
PO
4
-
Ca
2+
Number Density (nm
-3
)
Figure 6.3. Number density (nm-3) profiles for the interface regions of
the DOPC-only (A) and DOPC:DOPS (B) leaflets of a DOPC:DOPS
(128:20) asymmetric bilayer. The DOPC-only leaflet contains 64
DOPC. The DOPC:DOPS leaflet contains 44 DOPC and 20 DOPS. (A)
and (B) are expanded views of the left side and right side, respectively,
of Fig. 2D. In addition to the scaling applied in Fig. 2, the chloride and
sodium densities are each multiplied by 10.
142
6.4.5 Calcium density in DOPC and DOPC:DOPS systems
To understand the calcium distribution in the DOPC:DOPS system in a more quantitative
way, we consider next the concentration of calcium in the phosphate and carbonyl
regions of the bilayer, defined in Table 6.3 and shown graphically in Figure 6.2. As
shown in Table 6.3, the concentration of calcium in the phosphate and carbonyl regions
of the DOPC-only leaflet (left side of Figure 6.2) is not significantly affected by the
presence of PS in the opposite leaflet. In the DOPC:DOPS leaflet, however, increasing
the fraction of PS increases the concentration of calcium in both the phosphate and
carbonyl regions. (As shown in Figures 6.2 and 6.3, the PS carboxyl group lies almost
entirely inside the phosphate distribution, so it is not tabulated separately here.)
143
Table 6.3. Calcium density in phosphate and carbonyl regions of the phospholipid
interface in the DOPC and DOPC:DOPS leaflets.
Ca
2+
density (nm
-3
)
PO
4
PC
PO
4
PC:PS
PS 10 Ca 100 Ca 10 Ca 100 Ca
0 0.07 0.41 0.08 0.36
10 0.08 0.37 0.08 0.44
20 0.05 0.42 0.13 0.56
C=O
PC
C=O
PC:PS
PS 10 Ca 100 Ca 10 Ca 100 Ca
0 0.07 0.34 0.08 0.32
10 0.08 0.32 0.08 0.35
20 0.04 0.36 0.09 0.43
PO
4
+C=O
PC
PO
4
+C=O
PC:PS
PS 10 Ca 100 Ca 10 Ca 100 Ca
0 0.05 0.31 0.06 0.28
10 0.06 0.29 0.06 0.34
20 0.04 0.32 0.10 0.46
Densities are the time-averaged values (140 ns to 150 ns) of the number of calcium ions
in the rectangular volumes defined by the maximum and minimum position along the
bilayer normal of any of the oxygen atoms connected to phosphorus (PO
4
region) and any
of the carbonyl oxygens (C=O region) — the boundaries on either side of the peaks in
Figure 6.2. Phosphate and carbonyl regions overlap, and the PC phosphate, PS phosphate,
and PS carboxyl regions are essentially coincident (see Figure 6.2). PO
4
and C=O with
superscripts PC and PC:PS refer to the phosphate and carbonyl regions in the PC-only
leaflet (Figure 6.2, left side) and the mixed PC:PS leaflet (Figure 6.2, right side),
respectively. PO
4
+C=O is the combined phosphate and carbonyl region in each case,
defined by the extreme positions of phosphate and carbonyl oxygens.
Converting the atomic number densities in Table 6.3 to molar concentrations permits
comparison to physiological conditions and to previous analyses of calcium binding to
phospholipid membranes. For the 10-calcium systems presented here the concentration of
calcium in the bulk solvent is essentially zero (Figure 6.2A, 6.2C), so we may consider
144
the membrane calcium concentrations in these systems to lie near the lower limit of
physiological calcium concentrations in the aqueous-membrane interface. Simplifying by
taking values from the combined phosphate-carbonyl region in Table 6.3, we find
interface calcium concentrations in the range 60–100 mM in the DOPC-only leaflet and
110–170 mM in the DOPC:DOPS leaflet. These values, and the distribution of calcium in
the bilayer, are consistent with the model and measurements reported by Seelig (167) —
for [Ca
2+
]
bulk
= 1 mM, [Ca
2+
]
bound
(confined within a stratum of 1 nm thickness) = 167
mM, and with the simulation results reported by Pedersen et al. for charged phospholipid
bilayers (139).
Experimental studies suggest a surface calcium ion concentration that is 3 (167) to
30(172) times higher than the bulk concentration, and that this high surface
concentration, not preferential binding to PS, accounts for greater binding of Ca to PC:PS
relative to pure PC. The small size of our systems, in which a single calcium ion in the
bulk solvent region corresponds to a concentration of about 20 mM, does not permit
detection of this surface charge layer.
The 100-calcium systems are clearly well beyond membrane binding site saturation and
physiological conditions. The bulk calcium concentration in these systems is about 1.25
M, and the concentration in the phosphate-carbonyl region, again from Table 6.3, ranges
from 470 mM to 760 mM.
145
6.4.6 Water density in DOPC and DOPC:DOPS systems
Experimental studies indicate that calcium binding to phospholipid bilayers is
endothermic and entropy-driven, and that the entropy increase results at least in part from
the displacement of interfacial water by calcium (17, 172). In Table 6.4 we see this effect
most clearly by comparing the water densities in the DOPC-only systems with 0 and 10
calcium ions. In the 100-calcium systems the interfacial water concentration is higher
because the additional calcium in the membrane brings with it additional waters of
hydration. The effects of calcium and of the PS:PC ratio on water density distributions in
the mixed bilayer DOPC:DOPS systems are more complex because of the interactions of
packing configurations and electrostatics, and simulations at intermediate calcium levels
and PS:PC fractions will be required for the development of a more systematic analysis.
146
Table 6.4. Water density in phosphate and carbonyl regions of the phospholipid interface
H
2
O density (nm
-3
)
PO
4
PC
PO
4
PC:PS
PS 0 Ca 10 Ca
100
Ca 0 Ca 10 Ca
100
Ca
0 17.1 14.3 15.0 16.1 14.0 15.3
10 16.3 16.2 14.5 15.2 15.3 14.9
20 15.1 15.5 14.9 15.7 16.3 14.7
C=O
PC
C=O
PC:PS
PS 0 Ca 10 Ca
100
Ca 0 Ca 10 Ca
100
Ca
0 9.4 7.4 9.0 9.2 8.2 8.9
10 9.4 9.1 8.8 8.1 7.9 8.7
20 8.3 9.4 8.8 8.1 8.1 8.1
PO
4
+C=O
PC
PO
4
+C=O
PC:PS
PS 0 Ca 10 Ca
100
Ca 0 Ca 10 Ca
100
Ca
0 14.0 11.3 11.5 13.0 11.1 12.1
10 13.2 13.0 11.3 12.1 12.1 11.4
20 12.2 12.4 11.2 12.4 12.5 12.0
Densities and regions are defined as in Table 6.3, for water molecules rather than for
calcium ions.
6.4.7 Bound, interface, and bulk calcium
In addition to the spatial distributions for calcium described above, we also analyzed the
proximity of calcium to phospholipid binding sites (the electronegative phosphatidyl and
carbonyl oxygen atoms). Table 6.5 shows that for the 10-calcium systems almost all of
the calcium ions are bound to a lipid (located within 0.3 nm of a phosphatidyl or carbonyl
oxygen), as previously reported for PC bilayers (19) and that PS in the bilayer increases
147
the fraction of lipid-bound calcium. The attractive effect of PS is clearly evident in the
saturated, 100-calcium systems as well. These calcium-saturated bilayers also contain
more calcium ions in the interfacial region which are not bound to lipid, and which are
therefore more fully hydrated, consistent with the higher interfacial water densities in
100-calcium systems noted above.
Table 6.5. Distribution of calcium in DOPC and DOPC:DOPS bilayers
Number of calcium ions
System Bulk Interface Bound
0 PS, 10 Ca 1.0 0.7 8.3
10 PS, 10 Ca 0.0 0.3 9.7
20 PS, 10 Ca 0.0 0.1 9.9
0 PS, 100 Ca 60.5 2.6 36.9
10 PS, 100 Ca 58.3 4.1 37.6
20 PS, 100 Ca 52.0 2.6 45.4
Bound calcium is located within 0.3 nm of a lipid binding site (phosphatidyl, carbonyl, or
carboxyl oxygen). Interface calcium is located within the bilayer (defined by the
maximum excursions of phosphorus atoms from the bilayer center) but not bound to a
lipid. Bulk calcium is located outside the bilayer. Data are mean values averaged over the
interval 140 ns to 150 ns.
A detailed breakdown of the distribution of bound calcium in DOPC and DOPC:DOPS
bilayers is presented in Table 6.6. Not only does the PS carboxyl bind a significant
percentage of the total bound calcium, in agreement with previous simulations of anionic
phospholipid bilayers (139), it also dominates the DOPS intramolecular calcium binding
space. In the 10-calcium, 20 PS system, the PS carboxyls account for half of the total
calcium binding sites, and no calciums are found in binding proximity to the PS
148
phosphate or PS carbonyl sites. Low levels of PS phosphate and PS carbonyl binding are
observed in the 100-calcium system. The distribution of phosphate and carbonyl binding
in the DOPC-only system (0 PS) is consistent with infrared spectroscopic evidence for
PC vesicles at high calcium:lipid molar ratios (17).
Table 6.6. Equilibrium calcium binding in DOPC and DOPC:DOPS systems
Phosphate-, Carbonyl-, and Carboxyl-Bound Ca
2+
(< 0.3 nm), 140–150 ns
PC
Phosphate
PC
Carbonyl
PS
Phosphate
PS
Carbonyl
PS
Carboxyl
Total Ca
2+
Pairs
PS 10 Ca
100
Ca
10
Ca
100
Ca
10
Ca
100
Ca
10
Ca
100
Ca
10
Ca
100
Ca
10
Ca
100
Ca
0 12.8 96.0 11.9 46.7 – – – – – – 24.7 142.7
20 10.9 65.2 8.2 46.7 0.0 6.9 0.0 5.2 18.4 34.5 37.5 158.5
A bound calcium-lipid pair is recorded for each occurrence of a calcium ion located
within 0.3 nm of a lipid binding site (phosphatidyl, carbonyl, or carboxyl oxygen). One
calcium ion may be bound to multiple lipid binding sites and may thus be associated with
multiple bound pairs. Data are mean values of the number of bound pairs for each group
averaged over the interval 140 ns to 150 ns.
Calorimetry data for DOPC:DOPS vesicles with PS/PC = 0.11 suggests that calcium, at
calcium:lipid molar ratios less than 0.1, binds preferentially to PS, then, at higher calcium
concentrations, to PC(172). Since our 10-calcium, 20-PS DOPC:DOPS systems have
PS/(PC+PS) = 0.16 (but with an asymmetric distribution between the two leaflets of the
bilayer) and a calcium:lipid molar ratio of about 0.08, the calorimetric evidence would
lead us to expect that we are in the range where both PS and PC binding is occurring, and
that is what we see (Table 6.3, Table 6.6).
149
Graphical representations of the radial distribution functions for calcium and the
phospholipid binding sites (carbonyl, phosphatidyl, and carboxyl oxygens) in pure DOPC
(Figure 6.4A) and mixed DOPC:DOPS bilayers (Figure 6.4B) confirm the specificity of
the binding relationships presented in tabular form above. Note in particular the very
strong coordination with the PS carboxyl relative to PC carbonyl and phosphatidyl
oxygens and the absence of coordination with other PS oxygens (Figure 6.4B). When the
number of calcium ions is increased from 10 (Figure 6.4) to 100 (not shown), PS
carbonyl and phosphate coordination peaks appear, but at lower magnitudes than the
corresponding PC carbonyl and phosphate peaks, consistent with the data in Table 6.6.
150
Figure 6.4 Radial distribution functions from equilibrated pure DOPC
and mixed DOPC:DOPS (108:20) systems with 10 Ca
2+
, time-averaged
from 140 ns to 150 ns, for calcium and phospholipid binding sites
(C=O – glycerol backbone carbonyl oxygens, C-O – backbone ether
oxygens, PO4 – oxygens bonded to phosphorus, COO – serine carboxyl
oxygens). g(r) values are plotted on a log scale to facilitate
visualization of the secondary associations at radial separations greater
than 0.3 nm.
151
1
10
100
1,000
0.2 0.3 0.4 0.5 0.6
r (nm)
PC C=O
PC C-O
PC PO4
PS C=O
PS C-O
PS PO4
PS COO
g(r)
C=O
CO
O
-
1
10
100
1,000
PC C=O
PC C-O
PC PO4
152
6.4.8 Calcium-PS complex
Experimental evidence for a calcium:phosphatidylserine stoichiometric ratio of 1:2 (3,
45, 112) is supported by structures like the representative snapshot of a Ca(PS)
2
(H
2
O)
4
complex shown in Figure 6.5. In this ion-molecular assembly we see the replacement of
solvating water in the lipid head group by calcium, and the concomitant partial
dehydration of the calcium water shell — calcium coordinates 8–9 water molecules (115,
132) — which follow from the enthalpic and entropic considerations mentioned above in
the discussion of water density. A movie showing the formation of calcium-carboxyl
oxygen complexes like the one in Figure 6.5 is provided as Supplementary Material. The
calcium-carboxyl oxygen coordination numbers in Table 6.7 show the prevalence of 3-
and 4-oxygen structures — Ca(PS)
2
— in the 10-calcium systems. With 100 calcium ions
the binding affinity between the multiplicity of calcium ions and the limited number of
serine carboxyl oxygens dominates over the tendency to form the coordinated dimer.
Although no calcium-coordinated complexes involving both carboxyl oxygens and
phosphatidyl oxygens (from PC or PS) were observed during the 140 ns – 150 ns
sampling period, and no instances of calcium coordinating phosphatidyl oxygens from
two phospholipids, PC oxygens form a variety of associations with calcium — for
example, complexes with three carbonyl oxygens (from three different DOPC
molecules), with two carbonyl oxygens and a phosphatidyl oxygen (again from three
different DOPC molecules), and with carboxyl (DOPS) and carbonyl (DOPC) oxygens
153
detected (not shown). The lifetime of these complexes and their effect on the diffusion of
calcium and phospholipids in the membrane will be addressed in future work.
Table 6.7. Calcium coordination with DOPS carboxyl oxygens in DOPC:DOPS.
Calcium-Carboxyl Coordination Number Distribution, 140–150 ns
10 PS 10 Ca 10 PS 100 Ca 20 PS 10 Ca
[1] [2] [3] [4] [1] [2] [3] [4] [1] [2] [3] [4]
0.8 0.2 0.0 2.0 5.5 7.4 0.0 0.0 0.7 0.4 2.3 2.5
For each calcium the carboxyl coordination number is the number of carboxyl oxygens
within a radius of 0.3 nm. The values in the table are the number of calcium ions in each
of the three systems — DOPC:DOPS (118:10) with 10 Ca
2+
, DOPC:DOPS (118:10) with
100 Ca
2+
, and DOPC:DOPS (108:20) with 10 Ca
2+
— with coordination numbers 1
through 4, averaged over the period 140 ns to 150 ns. A calcium-carboxyl coordination
number of 5 was not detected over this interval. Calcium bound to 2 PS molecules will
have a coordination number of 3 or 4 (2 carboxyl oxygens from one PS and 1 or 2
carboxyl oxygens from the second PS).
154
Figure 5. Calcium-DOPS-water complex from an equilibrated
DOPC:DOPS (108:20) system. Calcium binds to DOPS and water with
a stoichiometric ratio of 1:2:4. Blue lines connect calcium and atoms
within 0.3 nm (oxygens from water and the DOPS carboxyl groups).
Violet – Ca, red – O, white – H, teal – C, blue – N, gold – P.
Figure 6.5. Calcium-DOPS-water complex from an equilibrated
DOPC:DOPS (108:20) system. Calcium binds to DOPS and water with
a stoichiometric ratio of 1:2:4. Blue lines connect calcium and atoms
within 0.3 nm (oxygens from water and the DOPS carboxyl groups).
Violet – Ca, red – O, white – H, teal – C, blue – N, gold – P.
155
6.4.9 Calcium and sodium in the bilayer interface
The sodium ions included as counter-ions for the PS carboxyl groups contribute, along
with calcium ions and the electronegative lipid binding sites, to the electrostatics and
energetics of the bilayer. Although the interactions of these components with sodium ions
are complex, some general patterns can be discerned from shifts in mean sodium
concentration (Table 6.8) and distribution peaks (Figure 6.6) in the membrane interface
As shown in Figure 6.6.A, and consistent with earlier studies (18, 60, 102), sodium
penetrates to the carbonyl level of a DOPC-only leaflet of a DOPC:DOPS (108:20)
bilayer. In the 10-calcium system, the presence of calcium ions in the phosphate region is
associated with an outward shift in the peak of the sodium ion distribution, away from the
carbonyl region (Figure 6.6B), while at the same time facilitating an increase in the
number of sodium ions in the interface. In the 100-calcium system this increased sodium
is expelled into the bulk solvent (Figure 6.6C), consistent with experimental
observations(112) and simulations of other phospholipid systems (139).
156
Figure 6.6. Number density (nm-3) profiles, including sodium, for the
DOPC-only leaflet of a DOPC:DOPS (128:20) asymmetric bilayer,
with some components not shown for clarity. (B) and (C) are
expanded views of the left sides of Figure. 6.2C and 6.2D, respectively.
As in Figure 6.3, calcium and sodium densities are each multiplied by
10, and water density is divided by 2.
157
0
5
10
15
-4 -3 -2 -1
0
5
10
15
0
5
10
15
Water/2
DOPC C=O
DOPC PO4
Ca·10
Na·10
Ca
2
+
Distance From Bilayer Center (nm)
A
B
Na
+
PO
4
-
Ca
2
+
Number Density (nm
-3
)
C
10 Ca
Na
+
100 Ca
0 Ca
Na
+
C=
O
158
On the DOPC:DOPS side of the bilayer, sodium is found in the phosphate region,
chaperoned by the carboxyl groups in the DOPC:DOPS (108:20) system, even with 100
calcium (Figure 6.3B). With only 10 PS in the leaflet instead of 20, however, sodium is
excluded from the interface of a 100-calcium system (not shown), just as it is from the
DOPC-only leaflet (Figure 6.6C).
Table 6.8. Sodium density in phosphate and carbonyl regions of the phospholipid
interface in the DOPC and DOPC:DOPS leaflets.
Na
+
density (nm
-3
)
PO
4
PC
PO
4
PC:PS
PS 0 Ca 10 Ca 100 Ca 0 Ca 10 Ca 100 Ca
10 0.029 0.015 0.010 0.122 0.148 0.019
20 0.029 0.062 0.030 0.292 0.286 0.154
C=O
PC
C=O
PC:PS
PS 0 Ca 10 Ca 100 Ca 0 Ca 10 Ca 100 Ca
10 0.029 0.015 0.002 0.121 0.145 0.008
20 0.029 0.058 0.020 0.278 0.253 0.139
PO
4
+C=O
PC
PO
4
+C=O
PC:PS
PS 0 Ca 10 Ca 100 Ca 0 Ca 10 Ca 100 Ca
10 0.023 0.012 0.008 0.098 0.116 0.014
20 0.024 0.049 0.023 0.231 0.219 0.125
Densities and regions are defined as in Table 6.3, for sodium ions rather than for calcium
ions.
159
6.4.10 Calcium and PC and PS head group dipole angle
The effects of monovalent and divalent ion binding on the phospholipid head group
dipole angle have been studied experimentally (166) and more recently with the tools of
molecular dynamics (19, 74, 135). Here we examine calcium-induced changes in the PC
and PS head group dipole angle distribution in DOPC and mixed DOPC:DOPS bilayers.
160
In the pure DOPC system the broad PC P Æ N dipole angle distribution peaks at
approximately 78 degrees from the bilayer normal (the choline nitrogen is on the average
lifted about 12 degrees out of the plane of the bilayer in the direction of the bulk solvent
(Figure 6.7A). In the same system after equilibration with 100 calcium ions, the head
group dipole angle peak shifts to about 43 degrees (Figure 6.7A), a very large change. In
a 10-calcium system (not shown), we observe an intermediate value, about 64 degrees,
for the peak of the dipole angle distribution.
161
Figure 6.7. Calcium-induced changes in PC and PS head group
orientation in DOPC and DOPC:DOPS bilayers. Distribution of head
group P →N dipole angles relative to the bilayer normal based on
atomic coordinates from (A) 128 DOPC molecules in a pure DOPC
bilayer, (B) 44 DOPC molecules in the DOPC-only leaflet of a
DOPC:DOPS (128:20) bilayer, (C) 20 DOPS molecules in the mixed
DOPC:DOPS leaflet of a DOPC:DOPS (108:20) bilayer. Each panel
shows distributions from simulations containing zero (diamond
symbol) and 100 (triangle symbol) calcium ions. Data was sampled
every 100 ps between 140 ns and 150 ns.
162
0.0
1.0
2.0
3.0
0 20 40 60 80 100 120 140 160 180
0.0
0.5
1.0
1.5
2.0
0.0
0.5
1.0
1.5
2.0
B
A
θ
P Æ N,PC
108 DOPC, 20
DOPS
44 PC: 20 PS
leaflet
θ
P Æ N,PC
128 DOPC, 0
DOPS
C
θ
P Æ N,PS
108 DOPC, 20
DOPS
Relative Frequency
Angle From Bilayer Normal (degrees)
0 Ca
100 Ca
163
In the PC:PS leaflet of a DOPC:DOPS (108:20) system, the PC dipole angle distribution
is bimodal, with a larger peak at about 76 degrees and a smaller peak at about 25 degrees
(Figure 6.7B). This “splitting” of the PC head group dipole angle arises from interactions
with neighboring PS head groups. The PC dipole angle distribution in the PC-only leaflet
of the mixed, asymmetric DOPC:DOPS (108:20) bilayer is essentially identical to the
distribution in the pure DOPC system. In the presence of 100 calcium ions a single PC
dipole angle peak is observed at about 44 degrees (Figure 6.7B). As might be expected,
we observe intermediate values for the PC dipole angle distribution in a 10-calcium
DOPC:DOPS (108:20) system, with peaks at 72 degrees and 44 degrees (not shown).
For comparison we looked also at the effect of calcium on the P Æ N vector component of
the head group dipole angle for PS in the DOPC:DOPS (108:20) system. We refer to this
as the PS head group dipole angle, although we are ignoring the contribution of the
negatively charged carboxyl group. A dominant peak in the distribution is observed at
about 85 degrees in the absence of calcium (Figure 6.7C). Interestingly, in a 100-calcium
system the PS head group dipole angle distribution becomes bimodal, with peaks at about
97 degrees and about 70 degrees (Figure 6.7C), suggesting, along with the results for the
PC dipole angle, complex and subtle adjustments in conformation and aggregation among
the phospholipids in a mixed zwitterionic:anionic bilayer in response to local changes in
monovalent and divalent cation concentrations. Again, the 10-calcium system dipole
164
angle distribution is intermediate, with a single peak shifted to about 93 degrees, and no
obvious second peak (not shown).
6.5 Conclusions
Molecular dynamics simulations of mixed, zwitterionic-anionic, asymmetrically arranged
DOPC:DOPS bilayers with monovalent and divalent cations are consistent with
experimental observations and with previously reported results for simulations of simpler
phospholipid systems. In addition to confirming for the DOPC:DOPS system the area-
decreasing, thickness-increasing, interface-dehydrating effects of calcium ions, the
greater affinity of phospholipid binding sites for calcium relative to sodium, the greater
affinity of calcium for PS carboxyl sites relative to phosphate and carbonyl sites, and the
calcium-induced rotation of the PC head group dipole toward the solvent, we have also
identified some new features of these model systems.
Area per lipid, often used as an indicator of attainment of a thermal and electrostatic
equilibrium in simulations of lipid bilayer systems, is misleading in the case of
DOPC:DOPS with calcium, for which the area per lipid stabilizes after about 50 ns, but
the establishment of a stable number of coordinated calcium-lipid pairs, and the
associated ion-lipid structural rearrangements, takes 100 ns.
165
For a 128-lipid DOPC:DOPS system, saturation of the calcium capacity of the bilayer
occurs between 10 and 100 calcium ions.
The PC head group dipole angle distribution in PC:PS bilayers is bimodal, an effect of
the interaction between the anionic PS and the zwitterionic PC, and this distribution can
be further modulated by varying the number of calcium ions in the system. Atomically
local binding profiles for anionic versus zwitterionic phospholipids and their
combinations, with monovalent and divalent counter-ions, may have significant
consequences for membrane lipid and protein aggregations.
The ion models from the GroMACS ffgmx topology database used here have limitations
arising from their non-polarizability (102), which we accept because of computational
resource limits. In particular, their performance in aqueous and non-aqueous systems may
differ markedly; adjusting the force field parameters for one environment can cause
unsatisfactory behavior in another (214). Because calcium in the interface region remains
at least partially hydrated (with some water oxygens replaced by carboxyl, carbonyl, and
phosphatidyl oxygens), and because the coordination number and geometry of this lipid-
bound calcium is similar to that of fully hydrated calcium (132), we think it is reasonable
to expect that the calcium model does not introduce unacceptable inaccuracies in these
phospholipid bilayer systems. This expectation is supported by the correspondence of our
results with experimental studies.
166
From another perspective, this work underlines the necessity for improvement of the
accuracy of the atomic-level interactions. Continuing and extending the development of
molecular models and mechanisms for membrane fusion, budding, invagination, ruffling,
pseudopod extension, and other remodeling processes, and for electroporation and
electric field-driven and mechanical restructuring of the membrane, will require taking
into account, as we have done here, the effects of monovalent and divalent cations,
calcium in particular, but not only calcium, in heterogeneous, asymmetric bilayers. But
beyond this, to begin to answer questions that approach the level of cell physiology,
functionally if not dimensionally, a simulation of ions at the phospholipid interface must
properly represent not only systems at equilibrium but also those with large spatial
(nanometer scale) and temporal (nanosecond scale) ion concentration gradients, such as
those encountered, for example, when plasma membrane ion channels are activated or
when the contents of intracellular calcium compartments are released. Continuum
properties like diffusion coefficients, binding constants, and stoichiometries, which guide
us at larger scales, are less helpful in this dynamic regime than the precision in local
electrostatics and van der Waals forces and in bond angles and dimensions which ensures
that molecular velocities and residence times reflect and respect the complexities of
hydration, hydrogen bonding, and hydrophobic interactions in the stereospecifically
complex, low-density (water), highly non-uniform electrical permittivity environment of
the membrane interface.
167
6.6 Acknowledgement
We thank Rainer Böckmann and Peter Tieleman for stimulating discussions. Computing
resources were provided by the USC Center for High Performance Computing and
Communications. This work was made possible in part by support from the Air Force
Office of Scientific Research. PTV and MJZ are supported by MOSIS, Information
Sciences Institute, Viterbi School of Engineering, University of Southern California.
168
CHAPTER 7: CALCIUM BINDING KINETICS
Portions of this chapter are adapted from:
Ziegler, M. J., Vernier, P. T. 2008. Nano-Kinetics of Calcium Binding to Mixed
Phosphatidylserine-Phosphatidylcholine Bilayers. Journal of Physical Chemistry B.
Submitted.
7.1 Abstract
Motivated by recent studies of calcium binding to zwitterionic and anionic lipid bilayers,
we performed molecular dynamics (MD) simulations of dioleoylphosphatidylcholine
(DOPC) bilayers with varying concentrations of dioleoylphosphatidylserine (DOPS) on
the inner leaflet in the presence of calcium ions to examine the effects of calcium binding
on membrane structural properties and the nanosecond time scale kinetics of the process.
The simulations presented here are consistent with experimentally determined binding
constants and diffusion rates, and provide molecular details of the calcium binding
process to lipid bilayers including selectivity to various electronegative centers.
Calcium binds to the carbonyl oxygen of the phospholipid tails, phosphatidyl oxygen, and
the carboxyl group of DOPS. When the amount of available calcium is limited, calcium
binds preferentially to the carboxyl oxygen. We observe the formation of many calcium
lipid complexes including calcium coordinated with two carboxyl oxygens, calcium
coordinated with up to three carbonyl oxygens, calcium coordinated with phosphatidyl
oxygens, and combinations of these. The diffusion rates of bound calcium and lipids are
169
similar, and complexes once formed do not dissociate for the remainder of the 150 ns
simulation. A given Ca
2+
, however, may coordinate with multiple lipids. The kinetics of
the coordination are best described by a set of first-order serial reactions.
7.2 Introduction
Calcium ions are of fundamental importance to the biochemistry of cells. Many cellular
processes including signal transduction (101), energy production, membrane fusion(138,
147, 180) and muscle contraction (83) depend on calcium. Calcium also affects the
electrical and structural properties of the membrane by modulating the transmembrane
and surface potentials (21, 29, 114, 140, 183), reducing the area per lipid (131, 139),
adjusting the head group dipole angle (2, 29, 166), dehydrating the membrane (142), and
binding to negatively charged (45, 113, 172) and zwitterionic molecules (147, 167, 172).
Calcium binding to lipids has been studied extensively experimentally utilizing a wide
variety of techniques including X-ray diffraction (33, 46), NMR (2, 3, 58, 81, 150, 151)
and IR spectroscopy (17), calcium-sensitive electrodes (148, 172), precipitation and
titration (40, 45, 172), and, to a limited extent, molecular simulations (19, 79, 139).
Measurements from these techniques and the properties derived from them vary
considerably. For example, reported values of the binding constant between calcium and
phosphatidylcholine bilayers differ by as many as four orders of magnitude (172).
170
The atomic resolution of molecular dynamics simulations may help to clarify the
experimental picture, but several factors limit the applicability of molecular dynamics
methods to this problem. One is the long time, 100–200 ns, necessary to equilibrate ions
with the membrane interface (19). Until recently the computational resources required for
simulations of atomistic models on this time scale were out of reach for most
investigators.
A second limitation is the calcium model itself. Non-polarizable models of divalent
cations, a necessity given current readily available computing power, can do an excellent
job of reproducing solvation free energies with appropriate parameter adjustments, but
these same adjusted parameters do not necessarily produce the correct free energies in
non-aqueous environments (214), a consequence of the symmetry of the solvation
analysis. Because an ion in bulk water is uniformly surrounded by water molecules, any
polarization-related effects can be lumped into the parameters describing the size of the
molecule. The amount of error introduced by a non-polarizable model in calcium binding
simulations depends on the geometry of the binding interactions.
Molecular simulations and experimental observations show that calcium binds to lipid
bilayers by forming coordinated complexes (19). These complexes, which involve
significant changes in head group and tail orientation, consist of a calcium molecule
surrounded by various oxygen atoms — carboxyl, carbonyl, phosphatidyl, glycerol ester,
and water oxygens — which can bridge multiple lipids (3, 19, 104). The oxygen shells
171
surrounding calcium in the bilayer interface — lipid oxygens in combination with a
partial water hydration shell — are similar to those in bulk water. On this basis it is
reasonable to expect that a non-polarizable model can provide a useful representation of
calcium binding in phospholipid bilayer simulations, and this expectation is supported by
the correlation of our results with experimental data.
7.3 Methods
7.3.1 Simulation Parameters
All simulations were performed using the GROMACS set of programs version 3.3.1 (13,
105, 198) on the University of Southern California High Performance Computing and
Communications Linux cluster (http://www.usc.edu/hpcc/). Lipid parameters were
derived from the OPLS, united-atom parameters of Berger et al. (14) Additional head
group parameters and charges for phosphatidylserine were borrowed from
Mukhopadhyay et al. (118) and Zu et al. (226). The Simple Point Charge (SPC) model
(11) for water was chosen for compatibility with the lipid models, computational
advantages, and the ability to produce the correct area per lipid. Systems were coupled to
a temperature bath at 310 K with a relaxation time of 0.1 ps and a pressure bath at 1 bar
with a relaxation time of 1 ps, each using a weak coupling algorithm (12). Pressure was
scaled semi-isotropically with a compressibility of 4.5 x 10
-5
bar
-1
in the plane of the
membrane and 4.5 x 10
-5
bar
-1
perpendicular to the membrane. Bond lengths were
constrained using LINCS (68) for lipids and SETTLE (117) for water. Short-range
172
electrostatics and Lennard-Jones interactions were cut off at 1 nm. Long-range
electrostatics were calculated by the PME algorithm (42) using fast Fourier transforms
and tin-foil boundary conditions. Real-space interactions were cut off at 1.0 nm, and
reciprocal-space interactions were evaluated on a 0.12 nm grid with fourth order B-spline
interpolation. The parameter ewald_rtol which controls the relative error for the Ewald
sum in the direct and reciprocal space was set to 10
-5
. Periodic boundary conditions were
employed to mitigate system size effects.
7.3.2 Structures
A sample DOPC lipid structure was obtained from the Protein Databank
(http://www.rcsb.org/), and custom code was used to create a bilayer system with 128
DOPC lipids and 4480 water molecules. The system was equilibrated until the area per
lipid was constant. 12 systems were created from the pure DOPC system, summarized in
Table 7.1.
173
Table 7.1. System configurations
System Name
Outer
Leaflet
Inner
Leaflet
Ions
Water
DOPC DOPC DOPS Calcium Sodium Chloride
0 PS 0 Ca 64 64 0 0 0 0 4480
0 PS 1 Ca 64 64 0 1 0 2 4477
0 PS 10 Ca 64 64 0 10 0 20 4450
0 PS 100 Ca 64 64 0 100 0 200 4180
10 PS 0 Ca 64 54 10 0 10 0 4480
10 PS 1 Ca 64 54 10 1 10 2 4477
10 PS 10 Ca 64 54 10 10 10 20 4450
10 PS 100 Ca 64 54 10 100 10 200 4180
20 PS 0 Ca 64 44 20 0 20 0 4480
20 PS 1 Ca 64 44 20 1 20 2 4477
20 PS 10 Ca 64 44 20 10 20 20 4450
20 PS 100 Ca 64 44 20 100 20 200 4180
For systems with DOPS, custom code was used to substitute a DOPS lipid for a DOPC
lipid by modifying the head group and adding a corresponding sodium counter-ion. A
short 10 ps equilibration was performed to test each lipid substitution for bad contacts.
For systems with calcium the GROMACS function genion was used to place ions. genion
randomly replaces a solvent water molecule with an ion. Molecular settling occurs on the
picosecond time scale and further equilibration is not necessary. For example, in order to
add CaCl
2
, 3 random water molecules are removed and replaced by either a calcium or
chloride ion.
174
0.21 0.25 0.30
0.00
0.02
0.04
0.06
0.08
0.10
Distance (nm)
Relative Frequency
Pair Correlation Function - 0 PS 100 Ca
DOPC O1
DOPC O4
DOPC O2 O3
DOPC O6 O8
DOPC N1
DOPC O5 O7
0.21 0.25 0.30
0.00
0.05
0.10
0.15
0.20
0.25
0.30
Distance (nm)
Relative Frequency
Pair Correlation Function - 20 PS 100 Ca
DOPS O3
DOPS O6
DOPS O4 O5
DOPS O8 O10
DOPS N1
DOPS O7 O9
DOPS O1 O2
Figure 1. Pair correlation functions of the minimum distance between
a calcium ion and any electronegative group in DOPC (top) and DOPS
(bottom) at equilibrium. Electronegative groups are described in Figure
2. Data are averaged from 140 – 150 ns in the simulation trajectory for
a system with 128 DOPC lipids and 100 calcium ions.
Figure 7.1. Pair correlation functions of the minimum distance between
a calcium ion and any electronegative group in DOPC (top) and DOPS
(bottom) at equilibrium. Electronegative groups are described in Figure
2. Data are averaged from 140 – 150 ns in the simulation trajectory for
a system with 128 DOPC lipids and 100 calcium ions.
175
7.3.3 Calcium Binding Criteria
Binding criteria for calcium to carbonyl atoms were similar to the criteria used in
simulations of calcium binding to palmitoyloleoylphosphatidylcholine (POPC) (19). A
pair correlation function (Figure 7.1) was created by calculating the minimum distance
between each calcium ion and all electronegative atoms in the system to screen for
potential binding sites. A summary of all electronegative atoms and their associated
partial charges is given in Figure 7.2. The pair correlation functions revealed peaks
between 0.23 nm and 0.30 nm. Binding interactions were observed for carboxyl,
carbonyl, and phosphatidyl oxygen. Binding to the ester oxygens of the glycerol
backbone or the amine group of the phosphatidylcholine head group was too infrequent
and sporadic for analysis on the time scale of these simulations.
176
N1 : -0.3
O1 : -0.82
O2 : -1.02
O3 : -1.02
O4 : -0.80
O6: -0.69
O5 : -0.66
O8 : -0.69
N1 : -0.60
O1 : -0.73
O2 : -0.73
O3 : -0.82
O5 : -1.02
O9 : -0.66
O10 : -0.69
O7 : -0.66
O8 : -0.69
O4 : -1.02
O6 : -0.80
O7 : -0.66
DOPC DOPS
Figure 2. Potential binding atoms for DOPC and DOPS. All negative
partial charges in the DOPC and DOPS models are shown Blue –
Nitrogen, Red–Oxygen, White–Hydrogen, Gold–Phosphorus, Teal
–Carbon
Figure 7.2. Potential binding atoms for DOPC and DOPS. All negative
partial charges in the DOPC and DOPS models are shown Blue –
Nitrogen, Red – Oxygen, White – Hydrogen, Gold – Phosphorus, Teal
– Carbon
177
7.3.4 Kinetic Models
Kinetic models were formed based on mass action rules with forward only reactions. The
rate constant matrix K was derived by discretizing the resulting differential equations
summarized below (Equations 7.1-7.3) and performing a least squares fit as described by
Equation 7.4 on the concentration C of each species with time. The superscript T
represents the transpose of the concentration matrix.
(7.1)
·
·
(7.2)
·
·1
·
(7.3)
d
dt
(7.4)
The notation Ca · nS represents a calcium molecule bound simultaneously to n binding
sites (S). The mathematics described above hold for any definition of binding site.
178
Discussion of specific site definitions and their effect on the rate constant are described
later.
7.3.5 Equilibrium Properties
Equilibrium properties of the system were determined from the time averaged simulation
data between 140 ns and 150 ns.
7.3.6 Images
Molecular graphics images were generated with Visual Molecular Dynamics (VMD)
(80).
179
50 ns
Figure 3. Snapshots from the 0
PS 100 Ca simulation. At the
beginning of the simulation the
calcium distribution is random.
As calcium enters the membrane
two populations of calcium form
inside and outside the bilayer with
a gap between regions.
0 ns
Figure 7.3. Snapshots from the 0 PS 100 Ca simulation. At the
beginning of the simulation the calcium distribution is random. As
calcium enters the membrane two populations of calcium form inside
and outside the bilayer with a gap between regions.
180
7.4 Results
7.4.1 Qualitative observations from the trajectory
The trajectories of MD simulations of fully hydrated DOPC bilayers in the presence of
excess calcium (100-calcium systems) reveal on visual inspection a story about how
calcium binds and interacts with lipid bilayers. Initially calcium is randomly distributed
in the system (Figure 7.3, 0 ns; Figure 7.4) and there is a large flux of calcium into the
membrane. An equilibrium between the amount of calcium in the membrane and in the
bulk solution is achieved after approximately 50 ns and a distinctive electrical double
layer is observed (Figure 7.3, 50 ns). As calcium enters the membrane, the area per lipid
slowly decreases while the thickness of the membrane increases. This process also
reaches equilibrium around 50 ns and is visualized in Figure 7.5. The calcium in the
bilayer interface settles between the phosphate and carbonyl groups and binds to
available oxygen from each species. Calcium in the interface begins to coordinate with
lipid binding sites, and this process takes longer to reach equilibrium than the
bulk:interface partition and the establishment of a new area per lipid. It is also dependent
on the amount of calcium in the system. The 100 Ca systems, representative of excess
calcium, reach equilibrium more quickly than the 10 Ca systems. Calcium molecules
remaining in the bulk have a larger diffusion rate than the bound population, and the flux
of calcium in and out of the membrane is small relative to the total number of calcium
ions in the bulk.
181
0 50 100 150
0
5
10
Number Calcium
0 50 100 150
20
40
60
80
Time (ns)
Number Calcium
Free
Bound
0 PS 10 Ca
0 PS 100 Ca
Figure 4. Progression of free to bound calcium.
When the amount of calcium is limited (10 Ca systems) it also takes approximately 50 ns
for the calcium to reach equilibrium concentrations in the bulk and in the interface
(Figure 7.6). In contrast to the 100-calcium systems, structural optimizations continue to
take place for the remainder of the simulation, reflected in the calcium-lipid coordination
Figure 7.4. Progression of free to bound calcium. In both the 10-
calcium (upper panel) and 100-calcium systems, the ratio of free to
bound calcium reaches an equilibrium value after 50 ns.
182
number (described below in the discussion of binding kinetics). Sample movies of the
trajectories are available as supplemental information online.
183
Figure 5. Thedependenceof areaper lipid on the number of DOPS
lipids on the inner leaflet calculated by dividing the number of lipids
on each leaflet by the product of the box dimensions that comprise the
bilayer plane. The number of calcium ions in the system is varied
between 0 calcium (blue) 10 calcium (red) and 100 calcium (green).
Figure 7.5. Time course of the equilibration of area per lipid in
asymmetric DOPC:DOPS bilayers with calcium ions. Upper panel:
pure DOPC; middle panel: one leaflet contains 64 DOPC, the other
contains 54 DOPC, 10 DOPS; lower panel: one leaflet contains 64
DOPC, the other contains 44 DOPC, 20 DOPS. The area per lipid is
calculated by dividing the number of lipids in each leaflet by the
product of the box dimensions that comprise the bilayer plane. Each
panel includes plots for systems containing 0 calcium ions (blue), 10
calcium ions (red), and 100 calcium ions (green).
184
Figure 6. Snapshots of the
progression of calcium from bulk
water to the membrane for the 0
PS 10 Ca system. At the
simulation onset calcium (purple
spheres) resides in the bulk water
phase. Within 50 ns all the
calcium is adsorbed by the
membrane. The complete
adsorption of calcium by the
membrane corresponds with the
equilibrium point in the time
evolution of the area per lipid,
however, further membrane
structural rearrangements
continue as calcium achieves an
optimal coordination number with
available lipid binding sites and
water in the interface.
a)
50 ns
0 ns
Figure 7.6. Snapshots showing the migration of calcium from bulk
water to the membrane for the 0 PS 10 Ca system. At the beginning of
the simulation calcium (purple spheres) is located primarily in the
solvent phase. Within 50 ns all the calcium is adsorbed by the
membrane.
185
7.4.2 Time evolution of the area per lipid
In simulations of lipid bilayers area per lipid is used as a guide in determining whether
the system is well equilibrated or suitable for analysis of equilibrium properties. Figure
7.5 shows the time evolution of the area per lipid for the systems studied here. At 0 ns
calcium and phosphatidylserine are added to the system as required to meet the desired
system composition. The area per lipid reaches equilibrium in 2 ns for the 10 PS 0 Ca
system and in 7 ns for the 20 PS 0 Ca system. For all systems with calcium the area per
lipid equilibration takes approximately 50 ns. This is consistent with the trajectory
observations (Figures 7.3, 7.6) and indicates that the bilayer has either adsorbed the entire
population of calcium (10 Ca systems) or reached a saturation point (100 Ca systems).
It should be noted that the Berger-based force field (14) used in these simulations
underestimates the experimental value of 0.724 nm
2
(99) for the area per lipid of pure
DOPC by 5%. It does capture, however, the general physical trends that are expected.
Both calcium and phosphatidylserine decrease the area per lipid experimentally, and that
is what we see in the simulations. Suitable alternative force fields for DOPC that
successfully reproduce the area per lipid without resorting to a fixed surface tension
(NP γT) or fixed volume (NVT) ensemble are not currently available.
186
7.4.3 Calcium binding to carbonyl oxygen and associated kinetics
DOPC and DOPS lipids have two and three centers of negative charge, respectively,
which provide binding sites for calcium. The phosphatidyl oxygens of the head group and
the carbonyl oxygens of the fatty acid tails are common to both lipids, while DOPS also
has a carboxyl component. Other electronegative atoms, including the saturated amine
nitrogen and the ester oxygens of the glycerol backbone, show much lower affinities for
calcium. In order to screen for binding between calcium and the electronegative centers,
we calculated a pair correlation function of the minimum distance between calcium and
any atom in the system with a negative charge (Figure 7.1). The lipid models do not
include any molecular polarizability and all charges remain constant throughout the
simulation. The partial charges on each lipid atom are summarized in Figure 7.2.
Figure 7.1 illustrates the binding interactions between each of the potential binding atoms
for 100-Ca DOPC systems with and without DOPS. Calcium is physisorbed (strongly
bound within 0.30 nm), to phosphatidyl oxygen, carbonyl oxygen and carboxyl oxygen.
The tightest associations, centered at 0.22 nm, are with the carboxyl, carbonyl, and
phosphoester oxygens (Figure 7.1: O1, O2 DOPS; O6, O8 DOPC; O1, O4 DOPC). Even
though the DOPC phosphoester oxygens are less electronegative in this model, they bind
more tightly than the non-ester phosphorus oxygens (Figure 7.1: O2, O3 DOPC). Despite
the similar partial charges of the four oxygens, calcium binds to carbonyl oxygens but not
to the glycerol backbone ester oxygens of DOPC (Figure 7.1: O5, O7 DOPC). Calcium is
187
most strongly attracted to the carboxyl group on DOPS, and this interaction causes a
reduction in the amount of calcium bound to the other binding sites. Unlike the glycerol
backbone oxygens of DOPC, the DOPS glycerol backbone ester oxygens (Figure 7.1: O7,
O9 DOPS) do bind to calcium.
Simulations of calcium and other cations have shown evidence of the formation of intra-
leaflet calcium complexes of multiple lipids bridged together by combinations of the
interactions described above (18, 19, 102, 118, 136, 224). . These calcium-facilitated
lipid aggregations, consisting of lipids on opposing leaflets, are precursors of membrane
fusion (126, 138, 147, 180). Böckmann et al. modeled the kinetics describing the
formation of these aggregations in POPC bilayers with a single-step, forward-only, mass-
action-based model (19). Here we apply the same model to simulations with DOPC and
asymmetric DOPC:DOPS bilayers. The system of reactions describing the production
rates of each species R is described by Equation 7.5 where the component matrix A and
stoichiometric ν matrix are defined in Equations 7.6 and 7.7 and the notation Ca · S refers
to a calcium molecule bound to a single carbonyl oxygen binding site S.
(7.5)
188
Ca
Ca · S
Ca · 2S
Ca · 3S
Ca · 4S
S
(7.6)
1 1 0 0 01
0 1 1 0 01
0 0 1 1 01
0 0 0 1 11
(7.7)
R
d Ca dt
,R
d Ca · S dt
,…R
d Ca · nS dt
(7.8)
Table 7.2. Calcium·Carbonyl Binding Rate Constants
System k
1
(ns
-1
) k
2
(ns
-1
) k
3
(ns
-1
) k
4
(ns
-1
)
0 PS 10 Ca 9.5e-5 3.1e-7 2.1e-5 1.2e-4
0 PS 100 Ca 7.2e-6 3.8e-5 2.7e-5 1.6e-4
20 PS 10 Ca 2.1e-5 7.1e-5 4.9e-5 7.2e-5
20 PS 100 Ca 8.2e-6 3.9e-5 1.7e-5 1.1e-4
Kinetic constants for the forward-only, mass-action-based, carbonyl-only binding model,
summarized in Table 7.2, were derived by performing a least squares fit of the system of
equations described by Equations 7.1-7.8 on the simulation data. Data for the 0 PS 10 Ca
189
system are plotted in Figure 7.7 along with curves generated from the fitted parameters in
Table 2. The model somewhat successfully reproduces [Ca , [Ca · S , and S values, as
illustrated by Figure 7.7. It fails, however, to model calcium coordinated with 2 or more
carbonyl oxygens. While this may be expected for [Ca · nS species for n ≥ 2 due to the
limited amounts of higher-order reactant species in the system, a direct consequence of
the small system size, it does not explain the poor fit for n < 2. The volatility in the data
suggests that back reactions should be included in the model. Figure 7.7 also illustrates a
considerable amount of direct Ca → Ca · 2S complexation that increases with DOPS
concentration. The separate addition of back reactions and direct Ca · 2S complexation
did not result in any noticeable improvement of the fits.
190
Figure 7. Fit of carbonyl only binding model - 0 PS 10 Ca. Ca · nS
refers to calcium bound to n carbonyl oxygen molecules. 128 DOPC
lipids provides a total of 256 sites.
Ca
free
Ca · S
Ca · 2S
Ca · 3S
Ca · 4S
Number Calcium
7.4.4 Generalized calcium binding
The unsatisfactory fit of the simulation data to the carbonyl-only binding model
suggested by Böckmann et al. (19) led us to further investigate the calcium complexes in
the system in order to refine the binding model. In contrast to sodium, calcium has a
stronger affinity for groups other than the carbonyl oxygen. Figures 7.8a and 7.8b show
some of the complexes and binding interactions observed in these simulations.
Figure 7.7. Fit of carbonyl-only binding model for the 0 PS 10 Ca,
system. Ca·nS refers to calcium bound to n carbonyl oxygen molecules.
128 DOPC lipids provide a total of 256 sites.
191
Figure 8. Calcium complexes A) Ca·2PS Complex. Calcium
coordinates two DOPS lipids by binding to both carboxyl oxygen
molecules of the headgroup along with a 4 water molecules for total
coordination number of 8. B) 2CO·Ca·PO
4
Complex. Calcium binds to
the phosphatidyl and sn-2 tail carbonyl moieties of one lipid and the
carbonyl group of another lipid along with 3 water molecules for a total
coordination number of 8.
Ca·CO
Ca·PO
4
2-
AB
Ca·COO
-
7.4.5 Calcium-carboxyl complex Ca·COO
-
The calcium-carboxyl complex is specific to DOPS and consists of a single calcium ion
bound to both ionized oxygens of the DOPS head group. Calcium can pair with up to two
DOPS molecules in this manner as illustrated by Figure 7.8a.
Figure 7.8. Calcium complexes A) Ca·2PS complex. Calcium
coordinates two DOPS lipids by binding to both carboxyl oxygen
molecules of the head group along with 4 water molecules for a total
coordination number of 8. B) 2CO·Ca·PO
4
complex. Calcium binds to
the phosphatidyl and sn-2 tail carbonyl moieties of one lipid and the
carbonyl group of another lipid along with 3 water molecules for a total
coordination number of 8.
192
7.4.6 Calcium-carbonyl complex Ca·CO
Calcium-carbonyl complexes (Figure 7.8b) are simple pair associations between calcium
and any carbonyl oxygen in the system.
7.4.7 Calcium-phosphatidyl oxygen complex Ca·PO
4
2-
The calcium-phosphatidyl oxygen complex (Figure 7.8b) consists of a calcium ion within
0.3 nm of all four phosphatidyl oxygens of a lipid molecule. These oxygen molecules
form two shells — one at 0.23 nm and one at 0.26 nm — which can be seen in the pair
correlation function (Figure 7.1). The square tetrahedral geometry of the phosphatidyl
group along with a carbonyl oxygen of the same lipid creates a pocket for calcium at the
base of the pyramid (Figure 7.8b).
7.4.8 Combination complexes
Several combination complexes existed in the simulation as well. A dual carboxyl
complex (Figure 7.8a) was observed, consistent with stoichiometric ratios obtained by
titration experiments (45). Multiple carbonyl complexes were also observed, consistent
with simulations of POPC systems (19). Mixed Ca·CO–Ca·PO
4
2-
and Ca·CO–Ca·COO
-
were observed, but no Ca·COO
-
–Ca·PO
4
2-
or multiple Ca·PO
4
2-
complexes appeared
193
during the 10 ns window of observation between 140 and 150 ns. We observed a
maximum calcium lipid coordination number of 4 in all simulations except for the 20 PS
100-Ca system, in which there were fleeting instances of calcium coordinated with 5
lipids.
194
7.4.9 Calcium and the head group dipole angle
The formation of specific complexes could yield clues explaining the changes in general
physical properties observed when calcium is added to model bilayer systems. For
example, the center of the distribution of the head group phosphorus nitrogen dipole
(P →N) angle with the bilayer normal shifts from near planar (~75º) to more
perpendicular with the bilayer (~40º) in the presence of calcium (210). Figure 7.8b
illustrates a possible mechanism for this shift. Calcium simultaneously binding to the
phosphatidyl group and a carbonyl oxygen of the same lipid tilts the P →N angle outward,
whereas a lone carbonyl binding interaction allows a wider range of P →N angles. As the
simulation progresses and the lipid coordination number of calcium increases, the P →N
angle of DOPC shifts accordingly. Figure 7.9 shows the P →N angle profile for DOPC
averaged and normalized over 3 intervals during the simulation. The continued evolution
of the P →N angle with time, long after the equilibration in area, illustrates that multiple
time scales are involved in the calcium binding process.
195
Figure 9. Time evolution of the P →N dipole averaged and normalized
over three simulation intervals with a +/- 5 degree moving average to
improve readability. The P →N dipole angle shifts from a more planar
orientation to a more bilayer normal orientation as the simulation
progresses and the calcium coordination increases. The shift in P →N
angle is caused by calcium binding to both the phosphatidyl and
carbonyl oxygen of a lipid molecule simultaneously (Figure 8).
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
0 50 100 150
Relative Frequency
PN Dipole Angle
Time Evolution of P →N Dipole Angle 0 PS 100 Ca
10‐30 ns
50‐70 ns
140‐150 ns
7.4.10 Updated kinetic model
The kinetic model was adjusted to be more consistent with the complexes found in the
simulation system, and the rate constants were recalculated using the methods described
above. Specifically, DOPC was treated as a 3-site adsorbent. Calcium was considered
bound to a binding site S when it was within 0.3 nm of all 4 phosphatidyl oxygen
Figure 7.9. Time evolution of the P →N dipole averaged and
normalized over three simulation intervals with a +/- 5 degree moving
average to improve readability. The P →N dipole angle shifts from a
more planar orientation toward the bilayer normal as the simulation
progresses and the calcium coordination number increases. The shift in
P →N angle is caused by calcium binding to both the phosphatidyl and
carbonyl oxygen of a lipid molecule simultaneously (Figure 7.8).
196
molecules or within 0.3 nm of either carbonyl oxygen. DOPS was treated as a 4-site
adsorbent. The additional site was considered occupied when calcium was within 0.3 nm
of both carboxyl oxygen molecules. The new rate constants for the expanded definition of
binding are summarized in Table 7.3.
Figure 10. Number of calcium complexes based on coordination
number for the 3-site DOPC kinetic model and the 0 PS 100 Ca system.
Calcium is rapidly adsorbed by the bilayer and the shift from lower
coordination numbers to higher coordination numbers is linear.
Ca
free
Ca · S
Ca · 2S
Ca · 3S
Ca · 4S
Figure 7.10. Number of calcium complexes based on coordination
number for the 3-site DOPC kinetic model and the 0 PS 100 Ca system.
Calcium is rapidly adsorbed by the bilayer and the shift from lower
coordination numbers to higher coordination numbers is linear.
197
Using this model, the 10 Ca and 100 Ca system behaviors were fundamentally different.
At the 100 Ca level (Figure 7.10), which corresponds to a bulk calcium concentration
greater than 1 M, calcium is rapidly taken into the membrane. A single calcium occupies
between 1 and 4 sites, and as the simulation progresses there is a shift in the calcium lipid
coordination number from lower to higher coordination numbers. The rate constant k
1
,
which describes the binding of free (non-bound) calcium with the membrane, increases
with phosphatidylserine concentration, indicating that DOPS enhances the initial uptake
of calcium.
198
Ca
free
Ca · S
Ca · 2S
Ca · 3S
Ca · 4S
Figure 11. Number of calcium complexes based on coordination
number for the 3 and 4-site DOPC and DOPS kinetic models for the 0
PS 10 Ca system (top) and the 20 PS 10 Ca system (bottom).
0 PS 10 Ca
20 PS 10 Ca
At
Figure 7.11. Number of calcium complexes based on coordination
number for the 3 and 4-site DOPC and DOPS kinetic models for the 0
PS 10 Ca system (top) and the 20 PS 10 Ca system (bottom).
199
the 10 Ca level ([Ca
2+
]
bulk
≈ 0) we observe a multi-step evolution in the coordination
number of the calcium complexes (Figure 7.11), as previously reported for POPC (19).
The initial step involves calcium entry into the membrane and coordination with a single
binding site. This process is enhanced by the presence of phosphatidylserine, and the
value of k
1
is increased by a factor of 3 (Table 7.3). The inflection point for the
concentration of free calcium , which coincides with the maximum concentration of Ca·S
(calcium coordinated with a single binding site) for the 0 PS 10 Ca system occurs near 50
ns, while for both the 10 PS 10 Ca and 20 PS 10 Ca systems it occurs much sooner —
around 15 ns (Figure 7.11). It is worth noting that the amount of time required for the
equilibration in area per lipid and the amount of time required to maximize the
concentration of Ca · S are similar. The rate constants associated with the migration from
a singly coordinated complex to a dual coordinated complex are smaller than the initial
binding step. This contrasts with the results of Böckmann et al., which suggested that k
2
> k
1
, (binding to a second group was faster than the initial binding) (19). For systems
with phosphatidylserine the maximum concentration of Ca·2S occurs between 80 and 130
ns. Additional longer simulations will be needed to examine the formation of calcium
complexes with coordination number of 3 and higher.
200
Table 7.3. Generalized Binding Rate Constants
System k
1
(ns
-1
) k
2
(ns
-1
) k
3
(ns
-1
) k
4
(ns
-1
)
0 PS 10 Ca 6.8e-5 3.0e-5 - -
10 PS 10 Ca 2.6e-4 8.1e-5 2.2e-5 4.2e-5
20 PS 10 Ca 1.9e-4 3.8e-5 2.1e-5 1.9e-5
0 PS 100 Ca 9.9e-6 2.8e-5 3.0e-5 3.7e-5
10 PS 100 Ca 8.2e-6 3.2e-5 1.9e-5 3.1e-5
20 PS 100 Ca 1.3e-5 3.6e-5 3.4e-5 3.6e-5
7.4.11 Equilibrium binding constants
·
1 ·
(7.9)
Equilibrium binding constants (K) were calculated using a Langmuir isotherm (Equation
7.9) where θ is the fraction of occupied sites. Literature values for the stoichiometry
between Ca
2+
and lipid bilayers vary between 1:1 and 1:2, and the association of calcium
in the membrane with different electronegative groups is difficult to extract from
experimental evidence (15). Simulations of lipids with Na
+
disagree on the location of
sodium binding as well, with some indicating only ester binding (89, 118), others
carbonyl and carboxyl binding (135), phosphate and carbonyl binding (136), carbonyl-
only binding (18), and still others no deep binding at all (153). In a review by Berkowitz
et al. the variation among simulations is attributed to varying simulation conditions and
models (15). Binding and area per lipid appear correlated, however, area is most strongly
dependent on the model and simulation conditions so this observation may be
confounded. To facilitate comparisons with literature values for the equilibrium binding
201
constant, we considered all possible scenarios by varying the number of binding sites per
lipid in our calculations. Seelig (167) and Sinn et al. (172) use a 1:2 binding model, while
Pandit et al. use a 1:1 model (136). Böckmann et al. use a 1:2 model in which calcium
binds only to carbonyl oxygen (19). We also present results using a 1:6,8 binding model
supported by the pair correlation function in Figure 7.2, in which we show binding
interactions for all 4 phosphatidyl oxygens, 2 carbonyl oxygens, and 2 carboxyl oxygens,
resulting in a 1:6 ratio for DOPC and a 1:8 ratio for DOPS. The intrinsic binding constant
of Ca
2+
to PC has been reported to be as low as 0.1 M
-1
and as high as 100 M
-1
(172).
With the exception of K
1:1
all the data presented in Table 7.4 fall within this range.
Table 7.4. Equilibrium binding constants
System C
bulk
(M) Free Ca
2+
Θ
1:6,8
K
1:6,8
(M
-1
) K
1:2
(M
-1
) K
1:1
(M
-1
)
0 PS 10 Ca 0.04 1.7 0.03 0.82 2.62 5.84
0 PS 100 Ca 1.54 63.1 0.19 0.15 0.82 -
10 PS 10 Ca 0.01 0.2 0.05 9.70 33.00 80.30
10 PS 100 Ca 1.52 62.5 0.18 0.14 0.76 -
20 PS 10 Ca 0.00 0.1 0.05 20.74 71.10 171.00
20 PS 100 Ca 1.20 54.6 0.20 0.21 1.35 -
Subscripts indicate the ratio of the number of calcium to the number binding atoms on a lipid in each
model. Θ is the fraction of bound atoms and K is the equilibrium constant.
Seelig suggested that the surface charge introduced by the phosphatidylserine carboxyl
group increases the availability of calcium to the membrane as the result of a general
electrostatic attraction apart from any effects of direct calcium-phosphatidylserine
binding (167). We find both — evidence of direct binding to the carboxyl group as well
as enhanced binding to the membrane, indicated by the fraction of occupied sites Θ and
increasing values of K with increasing concentrations of DOPS. Experiments by
202
Feigenson (45) also support calcium carboxyl complexation. The Ca·2PS stoichiometry
of the precipitate from these titration experiments of calcium and phosphatidylserine
bilayers is consistent with the 1:2 ratio of calcium to phosphatidylserine in complexes
observed in the simulation. Figure 7.8a shows a representative snapshot of calcium
coordinated with two carboxyl oxygens of nearby DOPS lipids.
The ratios of bound to added calcium for the 10 Ca data is in almost exact agreement
with Figure 7.7 of Sinn et al. (172). The 100 Ca results however are six times higher,
which suggests that the system is supersaturated. We performed additional simulations
with 100 Ca with three times more water and saw a decrease in the amount of calcium
bound to the membrane and a reduction in bulk calcium concentration from 1.3 M to 0.3
M. Further simulations will be necessary to reproduce more precisely the physiological
balance between membrane-bound and free calcium.
203
7.4.12 Water hydration shell
In these simulations we use a simple distance cut-off for determining binding, while
Pandit et al. suggest an additional criteria for binding — the loss of 50% of the water
hydration shell (136) . The preferred natural coordination number of calcium with water
is 8 (108). Experimentally observed coordination numbers vary with concentration and
are reported to be somewhere between 5 and 10, consistent with the results reported in
Table 7.5.
Figure 12. Snapshot of bulk calcium coordinated with 8 water
molecules.
Figure 7.12. Snapshot of bulk calcium coordinated with 8 water
molecules.
204
An example of the Ca · 8 H
O structure taken from the 0 PS 100 Ca system is shown in
Figure 7.12. Calcium-oxygen coordinated pairs with a separation of ≤ 0.3 nm are drawn
as solid lines. The cut-off of 0.3 nm is representative of the first hydration shell of
calcium (69). A summary of the hydration shell of calcium in bulk water is presented in
Table 7.5. Of the systems with 10 calcium, only the pure DOPC (0 PS 10 Ca) system had
sufficient calcium in the bulk for the determination of coordination numbers. The average
coordination numbers, which decrease with increasing calcium concentration, agree well
with those reported by Hewish et al. (69) and Palinkas and Heinzinger (132). For all
systems the most probable coordination number is 8, and at least 84% of calcium ions
have a coordination number of 5 or more. A calcium ion with a water coordination
number of 4 or less can be considered to have lost its hydration shell (136).
Table 7.5. Distribution of the coordination number of bulk calcium with water
System C
bulk
(M) [1]
*
[2] [3] [4] [5] [6] [7] [8] AVG
**
0 PS 10 Ca 0.04 - - - 0.10 - 0.20 - 0.70 7.20
0 PS 100 Ca 1.54 0.01 0.01 0.05 0.09 0.12 0.15 0.10 0.47 6.50
10 PS 100 Ca 1.52 - 0.01 0.04 0.08 0.14 0.15 0.12 0.48 6.65
20 PS 100 Ca 1.20 - 0.01 0.03 0.04 0.17 0.10 0.12 0.53 6.80
*[n] is the fraction of calcium with water coordination number n. **Average coordination number of
calcium with water.
If we adjust the equilibrium binding constants in Table 7.4 by treating calcium with water
coordination numbers of 5 or more as unbound, we obtain the values in Table 7.6, which
are in better agreement with values reported by Sinn et al. (172). It is readily apparent
205
that the choice of binding definitions (based on the types of complexes, availability and
type of binding sites, and inclusion of the water hydration shell) changes the value of the
equilibrium constant considerably.
Table 7.6. Equilibrium binding constants (Loss of 50% Hydration Shell)
System Θ
1:6,8
K
1:6,8
(M
-1
) K
1:2
(M
-1
) K
1:1
(M
-1
)
0 PS 10 Ca 0.003 0.084 0.256 0.518
0 PS 100 Ca 0.025 0.017 0.053 0.117
10 PS 10 Ca 0.010 1.942 6.031 12.468
10 PS 100 Ca 0.028 0.019 0.061 0.134
20 PS 10 Ca 0.008 3.582 11.225 23.062
20 PS 100 Ca 0.039 0.034 0.114 0.264
7.4.13 Water organization at the lipid interface
In the absence of calcium the dipole potential of the lipid headgroups is balanced by
water molecules aligning at the interface. The water molecular dipole points toward the
membrane interface in the absence of calcium. When calcium is present, it displaces
water from the membrane (172), and water reorganizes in response to the new
electrostatic environment. Water dipole moment profiles (Figure 7.13) reveal that
calcium causes a more random water orientation at the interface, indicated by a decrease
in dipole moment on the bulk side of the interface and an increase on the side adjoining
the membrane interior. Phosphatidylserine inhibits water intrusion into the membrane on
the opposing leaflet (Figure 7.13, bottom panel), indicated by the extension of the zero-
value of the dipole moment beyond the cutoff line.
206
Figure 13. The bilayer normal component of the water dipole moment.
In order to reduce noise and improve readability of the plot, values of
the dipole moment where the density is less than 2.5 % bulk density
and between +/- 0.75 nm are set to zero. The cut-off is indicated by
additional tick marks.
-3 -2 -1 0 1 2 3
-1.0
-0.5
0.0
Average Water Dipole Moment (D)
64 DOPC : 64 DOPC
-3 -2 -1 0 1 2 3
-1.0
-0.5
0.0
Distance From Bilayer Center (nm)
64 DOPC : 44 DOPC / 20 DOPS
0 Ca
2+
10 Ca
2+
100 Ca
2+
7.4.14 Diffusion rates of bound and free calcium
In order to verify and quantify observations from the trajectory visualizations described
above indicating the existence of two populations of calcium with differing diffusion
constants (faster in the bulk, slower in the membrane interface), we calculated the
diffusion constant of both membrane bound and free calcium using the Einstein diffusion
equation given below.
Figure 7.13. The bilayer normal component of the water dipole
moment. In order to reduce noise and improve readability of the plot,
values of the dipole moment where the density is less than 2.5 % bulk
density and between +/- 0.75 nm are set to zero. The cut-off is indicated
by additional tick marks.
207
6 (7.10)
The diffusion constant D is related to ratio of the mean square deviation of the distance
r(t) of a particle from its initial location r(t
0
) to the total elapsed time t. When calculating
the diffusion constant with systems that use periodic boundary conditions, it is critical
that jumps over the box edges are taken into account and that the time step selected is
such that the distance traveled by any molecules between sampling is less than half the
box dimension in the direction of travel.
Table 7.7. Diffusion Coefficients for Bulk and Bound Calcium
Bulk Calcium Bound Calcium
Average
Standard
Deviation
Median Average
Standard
Deviation
Median
0 PS 0 Ca NA NA NA NA Na Na
0 PS 10 Ca 1.7E-05 no data 1.7E-05 1.3E-06 7.6E-07 7.7E-07
0 PS 100 Ca 1.8E-06 1.8E-06 1.1E-06 3.4E-08 2.8E-08 2.8E-08
10 PS 0 Ca NA NA NA NA NA NA
10 PS 10 Ca * * * 6.6E-07 2.7E-07 6.2E-07
10 PS 100 Ca 2.2E-06 1.8E-06 1.5E-06 2.6E-08 9.1E-09 2.4E-08
20 PS 0 Ca NA NA NA NA NA NA
20 PS 10 Ca * * * 5.7E-07 2.3E-07 5.2E-07
20 PS 100 Ca 2.2E-06 1.7E-06 1.3E-06 6.4E-08 1.6E-08 5.9E-08
Units reported in cm2/s *All calcium was bound to the bilayer
Calcium diffusion rates in bulk water (Table 7.7) are within the range of values reported
experimentally (64, 211). Increasing concentrations of calcium appear to decrease the
diffusion coefficient of calcium in water which agrees with the theoretical predictions of
Onsager and Fuoss (130). The presence of phosphatidylserine at a molar fraction of 0.3 in
208
one leaflet (20 PS systems) does not affect the bulk calcium diffusion coefficient within
statistical bounds.
The diffusion coefficient of bound calcium closely matches the diffusion coefficient of
DOPC (Table 7.8), an indication that the binding interactions are strong and long-lasting,
consistent with observations from the trajectory described above.
7.4.15 Diffusion rates of lipids
DOPS diffusion coefficients in these simulations are greater than those for DOPC (Table
7.8), probably an artifact resulting from the strong affinity of calcium for DOPS carboxyl
oxygens. The largest DOPS diffusion coefficient is found for the system in which the
average starting distance between a calcium and a DOPS molecule is greatest, the 10 PS
10 Ca system. This indicates that the diffusion constant for DOPS is reflective of
electrophoretic motion rather than solely brownian motion and, along with the evolution
of the head group dipole angle shown in Figure 7.9, suggests that longer simulation times
are needed to fully equilibrate DOPC:DOPS·Ca
2+
systems.
209
Table 7.8. Diffusion Coefficients for DOPS and DOPS
Average
DOPC
Standard
Deviation Average
DOPS
Standard
Deviation
0 PS 0 Ca 1.2E-06 6.9E-07 NA NA
0 PS 10 Ca 9.8E-07 6.0E-07 NA NA
0 PS 100 Ca 3.6E-07 2.8E-07 NA NA
10 PS 0 Ca 8.6E-07 5.3E-07 5.8E-06 2.1E-06
10 PS 10 Ca 7.2E-07 6.2E-07 1.0E-05 4.9E-06
10 PS 100 Ca 2.8E-07 1.9E-07 2.8E-06 3.1E-07
20 PS 0 Ca 7.0E-07 5.5E-07 2.6E-06 8.2E-07
20 PS 10 Ca 7.7E-07 4.5E-07 3.3E-06 2.6E-06
20 PS 100 Ca 6.6E-07 3.5E-07 4.2E-06 2.9E-06
Calcium reduces the diffusion constant of DOPC, consistent with the formation of slower
moving complexes of multiple lipids. Although more simulations are needed to reduce
the statistical uncertainty, the DOPC diffusion coefficient appears to be reduced in
bilayers containing DOPS. Decomposition of the diffusion data into separate leaflets did
not reveal differences between the pure DOPC side and the mixed DOPC:DOPS side
(data not shown). Both Ca
2+
and DOPS reduce the area per lipid (210) leading to a more
rigid membrane (131). The diffusion coefficients for DOPC in homogeneous bilayers
reported here are roughly an order of magnitude faster than experimentally determined
values (47) and those reported for other simulated systems(173). Molecular dynamics
diffusion rates for pure lipids are dependent on the water model (173); faster diffusing
water models result in higher diffusion coefficients for lipids. The SPC water model used
here overestimates the diffusion rate of pure water, which accounts for the fast diffusion
constants reported here. Although the absolute values of the diffusion coefficients are not
accurate, the validity of the results presented here rests on the relative effects of the
210
various components of these bilayer systems, which we have shown to be consistent with
experimental observations.
7.5 Conclusions
The molecular dynamics simulations presented here illustrate several key features of the
calcium-phospholipid binding process and are consistent with both experimental findings
and other simulations of calcium and lipid bilayers. The simulations show that calcium
binds to carbonyl, phosphatidyl, and carboxyl oxygen; phosphatidylserine enhances
binding kinetics by drawing more calcium into the membrane; and the binding
interactions are long-lasting on the 150-ns time scale of the simulations. Calcium bound
to the membrane has a diffusion constant similar to the lipids. Calcium binds to multiple
lipids simultaneously, forming coordinated complexes. The formation of the complexes is
consistent with first-order, mass-action-rate kinetics. The equilibrium binding constants
in the simulation are consistent with experimental findings, although more precision is
needed in these measurements. Calculations of the binding constant are strongly
dependent on the assumptions which are used define binding. By adjusting the definition
of a binding site or by incorporating the loss of the water hydration shell as a criterion for
binding, the binding constant can change by an order of magnitude. This finding is
consistent with and helps to explain the range of values found in the literature. Although
divalent cations are difficult to model with non-polarizable models, and lipid properties
211
can be dependent on the choice of simulation ensemble (173), the simulations presented
here reproduce most physical properties and are consistent with experimental trends.
7.6 Acknowledgement
We thank Rainer Böckmann, Peter Tieleman, Arieh Warshel, and Rumiana Dimova for
stimulating discussions. Computing resources were provided by the USC Center for High
Performance Computing and Communications. MJZ and PTV are supported by the Air
Force Office of Scientific Research and MOSIS, Information Sciences Institute, Viterbi
School of Engineering, University of Southern California.
212
REFERENCES
1. Abidor, I. G., V. B. Arakelyan, L. V. Chernomordik, Y. A. Chizmadzhev, V. F.
Pastushenko, and M. R. Tarasevich. 1979. Electric Breakdown of Bilayer Lipid-
Membranes .1. Main Experimental Facts and Their Qualitative Discussion.
Bioelectrochemistry and Bioenergetics 6:37-52.
2. Akutsu, H., and J. Seelig. 1981. Interaction of metal ions with
phosphatidylcholine bilayer membranes. Biochemistry 20:7366-7373.
3. Altenbach, C., and J. Seelig. 1984. Ca2+ binding to phosphatidylcholine bilayers
as studied by deuterium magnetic resonance. Evidence for the formation of a Ca2+
complex with two phospholipid molecules. Biochemistry 23:3913-3920.
4. Anezo, C., A. H. de Vries, H. D. Holtje, D. P. Tieleman, and S. J. Marrink. 2003.
Methodological issues in lipid bilayer simulations. Journal of Physical Chemistry B
107:9424-9433.
5. Balasubramanian, K., and A. J. Schroit. 2003. Aminophospholipid asymmetry: A
matter of life and death. Annual Review of Physiology 65:701-734.
6. Barnett, A., and J. C. Weaver. 1991. Electroporation - a Unified, Quantitative
Theory of Reversible Electrical Breakdown and Mechanical Rupture in Artificial Planar
Bilayer-Membranes. Bioelectrochemistry and Bioenergetics 25:163-182.
7. Basse, F., J. G. Stout, P. J. Sims, and T. Wiedmer. 1996. Isolation of an
erythrocyte membrane protein that mediates Ca2+-dependent transbilayer movement of
phospholipid. Journal of Biological Chemistry 271:17205-17210.
8. Beebe, S. J., P. M. Fox, L. J. Rec, K. Somers, R. H. Stark, and K. H. Schoenbach.
2002. Nanosecond pulsed electric field (nsPEF) effects on cells and tissues: Apoptosis
induction and tumor growth inhibition. Ieee Transactions on Plasma Science 30:286-292.
9. Beebe, S. J., P. M. Fox, L. J. Rec, E. L. Willis, and K. H. Schoenbach. 2003.
Nanosecond, high-intensity pulsed electric fields induce apoptosis in human cells.
FASEB J 17:1493-1495.
10. Benz, R., and U. Zimmermann. 1980. Pulse-Length Dependence of the Electrical
Breakdown in Lipid Bilayer-Membranes. Biochimica et biophysica acta 597:637-642.
11. Berendsen, H. J. C., J. P. M. Postma, W. F. van Gunsteren, and J. Hermans. 1981.
Intermolecular Forces, Reidel, Dordrecht.
213
12. Berendsen, H. J. C., J. P. M. Postma, W. F. Vangunsteren, A. Dinola, and J. R.
Haak. 1984. Molecular-Dynamics with Coupling to an External Bath. J Chem Phys
81:3684-3690.
13. Berendsen, H. J. C., D. Vanderspoel, and R. Vandrunen. 1995. Gromacs - a
Message-Passing Parallel Molecular-Dynamics Implementation. Comput Phys Commun
91:43-56.
14. Berger, O., O. Edholm, and F. Jahnig. 1997. Molecular dynamics simulations of a
fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and
constant temperature. Biophysical journal 72:2002-2013.
15. Berkowitz, M. L., D. L. Bostick, and S. Pandit. 2006. Aqueous solutions next to
phospholipid membrane surfaces: insights from simulations. Chem Rev 106:1527-1539.
16. Berridge, M. J., M. D. Bootman, and H. L. Roderick. 2003. Calcium signalling:
dynamics, homeostasis and remodelling. Nat Rev Mol Cell Biol 4:517-529.
17. Binder, H., and O. Zschornig. 2002. The effect of metal cations on the phase
behavior and hydration characteristics of phospholipid membranes. Chemistry and
physics of lipids 115:39-61.
18. Bockmann, R. A., A. Hac, T. Heimburg, and H. Grubmuller. 2003. Effect of
sodium chloride on a lipid bilayer. Biophysical journal 85:1647-1655.
19. Bockmann, R. A., and H. Grubmuller. 2004. Multistep binding of divalent cations
to phospholipid bilayers: a molecular dynamics study. Angewandte Chemie (International
ed 43:1021-1024.
20. Bockmann, R. A., B. L. de Groot, S. Kakorin, E. Neumann, and H. Grubmuller.
2008. Kinetics, Statistics, and Energetics of Lipid Membrane Electroporation Studied by
Molecular Dynamics Simulations. Biophysical journal.
21. Brockman, H. 1994. Dipole potential of lipid membranes. Chemistry and physics
of lipids 73:57-79.
22. Brooks, B. R., R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and
M. Karplus. 1983. Charmm - a Program for Macromolecular Energy, Minimization, and
Dynamics Calculations. J Comput Chem 4:187-217.
23. Buchner, R., J. Barthel, and J. Stauber. 1999. The dielectric relaxation of water
between 0 degrees C and 35 degrees C. Chemical Physics Letters 306:57-63.
214
24. Case, D. A., T. E. Cheatham, 3rd, T. Darden, H. Gohlke, R. Luo, K. M. Merz, Jr.,
A. Onufriev, C. Simmerling, B. Wang, and R. J. Woods. 2005. The Amber biomolecular
simulation programs. J Comput Chem 26:1668-1688.
25. Chandler, D. 2005. Interfaces and the driving force of hydrophobic assembly.
Nature 437:640-647.
26. Chaplin, M. 2008. Water Structure and Science.
27. Chernomordik, L. V., S. I. Sukharev, I. G. Abidor, and Y. A. Chizmadzhev. 1983.
Breakdown of Lipid Bilayer-Membranes in an Electric-Field. Biochimica et biophysica
acta 736:203-213.
28. Chizmadzhev, Y. A., and I. G. Abidor. 1980. Bilayer Lipid-Membranes in Strong
Electric-Fields. Bioelectrochemistry and Bioenergetics 7:83-100.
29. Clarke, R. J., and C. Lupfert. 1999. Influence of anions and cations on the dipole
potential of phosphatidylcholine vesicles: a basis for the Hofmeister effect. Biophysical
journal 76:2614-2624.
30. Clarke, R. J. 2001. The dipole potential of phospholipid membranes and methods
for its detection. Adv Colloid Interface Sci 89-90:263-281.
31. Cole, K. S. 1937. Electric impedance of marine egg membranes. Transactions of
the Faraday Society 33:0966-0971.
32. Comfurius, P., P. Williamson, E. F. Smeets, R. A. Schlegel, E. M. Bevers, and R.
F. A. Zwaal. 1996. Reconstitution of phospholipid scramblase activity from human blood
platelets. Biochemistry 35:7631-7634.
33. Coorssen, J. R., and R. P. Rand. 1995. Structural effects of neutral lipids on
divalent cation-induced interactions of phosphatidylserine-containing bilayers.
Biophysical journal 68:1009-1018.
34. DeBruin, K. A., and W. Krassowska. 1998. Electroporation and shock-induced
transmembrane potential in a cardiac fiber during defibrillation strength shocks. Ann
Biomed Eng 26:584-596.
35. DeBruin, K. A., and W. Krassowska. 1999. Modeling electroporation in a single
cell. II. Effects Of ionic concentrations. Biophysical journal 77:1225-1233.
36. DeBruin, K. A., and W. Krassowska. 1999. Modeling electroporation in a single
cell. I. Effects Of field strength and rest potential. Biophysical journal 77:1213-1224.
215
37. Diederich, A., G. Bahr, and M. Winterhalter. 1998. Influence of surface charges
on the rupture of black lipid membranes. Physical Review E 58:4883-4889.
38. Drago, G. P., M. Marchesi, and S. Ridella. 1984. The Frequency-Dependence of
an Analytical Model of an Electrically Stimulated Biological Structure.
Bioelectromagnetics 5:47-62.
39. Dzubiella, J., R. J. Allen, and J. P. Hansen. 2004. Electric field-controlled water
permeation coupled to ion transport through a nanopore. J Chem Phys 120:5001-5004.
40. Ekerdt, R., and D. Papahadjopoulos. 1982. Intermembrane contact affects calcium
binding to phospholipid vesicles. Proceedings of the National Academy of Sciences of
the United States of America 79:2273-2277.
41. Elliott, J. I., A. Surprenant, F. M. Marelli-Berg, J. C. Cooper, R. L. Cassady-Cain,
C. Wooding, K. Linton, D. R. Alexander, and C. F. Higgins. 2005. Membrane
phosphatidylserine distribution as a non-apoptotic signalling mechanism in lymphocytes.
Nature Cell Biology 7:808-U876.
42. Essmann, U., L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen.
1995. A Smooth Particle Mesh Ewald Method. J Chem Phys 103:8577-8593.
43. Fadok, V. A., A. de Cathelineau, D. L. Daleke, P. M. Henson, and D. L. Bratton.
2001. Loss of phospholipid asymmetry and surface exposure of phosphatidylserine is
required for phagocytosis of apoptotic cells by macrophages and fibroblasts. Journal of
Biological Chemistry 276:1071-1077.
44. Falck, E., M. Patra, M. Karttunen, M. T. Hyvonen, and I. Vattulainen. 2004.
Impact of cholesterol on voids in phospholipid membranes. J Chem Phys 121:12676-
12689.
45. Feigenson, G. W. 1986. On the nature of calcium ion binding between
phosphatidylserine lamellae. Biochemistry 25:5819-5825.
46. Feigenson, G. W. 1989. Calcium ion binding between lipid bilayers: the four-
component system of phosphatidylserine, phosphatidylcholine, calcium chloride, and
water. Biochemistry 28:1270-1278.
47. Filippov, A., G. Oradd, and G. Lindblom. 2003. Influence of cholesterol and
water content on phospholipid lateral diffusion in bilayers. Langmuir 19:6397-6400.
48. Frey, W., J. A. White, R. O. Price, P. F. Blackmore, R. P. Joshi, R. Nuccitelli, S.
J. Beebe, K. H. Schoenbach, and J. F. Kolb. 2006. Plasma membrane voltage changes
during nanosecond pulsed electric field exposure. Biophysical journal 90:3608-3615.
216
49. Fricke, H. 1925. A mathematical treatment of the electric conductivity and
capacity of disperse systems II. The capacity of a suspension of conducting spheroids
surrounded by a non-conducting membrane for a current of low frequency. Physical
review 26:0678-0681.
50. Fricke, H. 1925. The electric capacity of suspensions with special reference to
blood. Journal of General Physiology 9:137-152.
51. Gabriel, B., and J. Teissie. 1997. Direct observation in the millisecond time range
of fluorescent molecule asymmetrical interaction with the electropermeabilized cell
membrane. Biophysical journal 73:2630-2637.
52. Garon, E. B., D. Sawcer, P. T. Vernier, T. Tang, Y. H. Sun, L. Marcu, M. A.
Gundersen, and H. P. Koeffler. 2007. In vitro and in vivo evaluation and a case report of
intense nanosecond pulsed electric field as a local therapy for human malignancies.
International Journal of Cancer 121:675-682.
53. Gawrisch, K., D. Ruston, J. Zimmerberg, V. A. Parsegian, R. P. Rand, and N.
Fuller. 1992. Membrane dipole potentials, hydration forces, and the ordering of water at
membrane surfaces. Biophysical journal 61:1213-1223.
54. Glaser, R. W., S. L. Leikin, L. V. Chernomordik, V. F. Pastushenko, and A. I.
Sokirko. 1988. Reversible Electrical Breakdown of Lipid Bilayers - Formation and
Evolution of Pores. Biochimica et biophysica acta 940:275-287.
55. Goldman, D. E. 1944. Potential, impedance, and rectification in membranes.
Journal of General Physiology 27:37-60.
56. Gowrishankar, T. R., and J. C. Weaver. 2003. An approach to electrical modeling
of single and multiple cells. Proceedings of the National Academy of Sciences of the
United States of America 100:3203-3208.
57. Gowrishankar, T. R., A. T. Esser, Z. Vasilkoski, K. C. Smith, and J. C. Weaver.
2006. Microdosimetry for conventional and supra-electroporation in cells with
organelles. Biochemical and Biophysical Research Communications 341:1266-1276.
58. Grasdalen, H., L. E. Goran Eriksson, J. Westman, and A. Ehrenberg. 1977.
Surface potential effects on metal ion binding to phosphatidylcholine membranes 31P
NMR study of lanthanide and calcium ion binding to egg-yolk lecithin vesicles.
Biochimica et biophysica acta 469:151-162.
59. Grosse, C., and H. P. Schwan. 1992. Cellular Membrane-Potentials Induced by
Alternating-Fields. Biophysical journal 63:1632-1642.
217
60. Gurtovenko, A. A. 2005. Asymmetry of lipid bilayers induced by monovalent
salt: Atomistic molecular-dynamics study. J Chem Phys 122:-.
61. Gurtovenko, A. A., and I. Vattulainen. 2005. Pore formation coupled to ion
transport through lipid membranes as induced by transmembrane ionic charge imbalance:
Atomistic molecular dynamics study. Journal of the American Chemical Society
127:17570-17571.
62. Gurtovenko, A. A., and I. Vattulainen. 2007. Ion leakage through transient water
pores in protein-free lipid membranes driven by transmembrane ionic charge imbalance.
Biophysical journal 92:1878-1890.
63. Hamilton, W. A., and A. J. H. Sale. 1967. Effects of High Electric Fields on
Microorganisms .2. Mechanism of Action of Lethal Effect. Biochimica et biophysica acta
148:789-&.
64. Harned, H. S., and A. L. Levy. 1949. The Differential Diffusion Coefficient of
Calcium Chloride in Dilute Aqueous Solutions at 25-Degrees. Journal of the American
Chemical Society 71:2781-2783.
65. Harrison, R. L., B. J. Byrne, and L. Tung. 1998. Electroporation-mediated gene
transfer in cardiac tissue. Febs Letters 435:1-5.
66. Heinz, T. N., W. F. van Gunsteren, and P. H. Hunenberger. 2001. Comparison of
four methods to compute the dielectric permittivity of liquids from molecular dynamics
simulations. J Chem Phys 115:1125-1136.
67. Herbette, L., C. A. Napolitano, and R. V. McDaniel. 1984. Direct determination
of the calcium profile structure for dipalmitoyllecithin multilayers using neutron
diffraction. Biophysical journal 46:677-685.
68. Hess, B., H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije. 1997. LINCS: A
linear constraint solver for molecular simulations. J Comput Chem 18:1463-1472.
69. Hewish, N. A., G. W. Neilson, and J. E. Enderby. 1982. Environment of Ca-2+
Ions in Aqueous Solvent. Nature 297:138-139.
70. Hibino, M., M. Shigemori, H. Itoh, K. Nagayama, and K. Kinosita. 1991.
Membrane Conductance of an Electroporated Cell Analyzed by Submicrosecond Imaging
of Transmembrane Potential. Biophysical journal 59:209-220.
71. Hibino, M., H. Itoh, and K. Kinosita, Jr. 1993. Time courses of cell
electroporation as revealed by submicrosecond imaging of transmembrane potential.
Biophysical journal 64:1789-1800.
218
72. Hickman, S. E., J. Elkhoury, S. Greenberg, I. Schieren, and S. C. Silverstein.
1994. P2z Adenosine-Triphosphate Receptor Activity in Cultured Human Monocyte-
Derived Macrophages. Blood 84:2452-2456.
73. Hofsass, C., E. Lindahl, and O. Edholm. 2003. Molecular dynamics simulations
of phospholipid bilayers with cholesterol. Biophysical journal 84:2192-2206.
74. Hogberg, C. J., and A. P. Lyubartsev. 2008. Effect of local anesthetic lidocaine on
electrostatic properties of a lipid bilayer. Biophysical journal 94:525-531.
75. Homan, R., and H. J. Pownall. 1988. Transbilayer Diffusion of Phospholipids -
Dependence on Headgroup Structure and Acyl Chain-Length. Biochimica et biophysica
acta 938:155-166.
76. Hu, Q., R. P. Joshi, and K. H. Schoenbach. 2005. Simulations of nanopore
formation and phosphatidylserine externalization in lipid membranes subjected to a high-
intensity, ultrashort electric pulse. Physical review 72:031902.
77. Hu, Q., S. Viswanadham, R. P. Joshi, K. H. Schoenbach, S. J. Beebe, and P. F.
Blackmore. 2005. Simulations of transient membrane behavior in cells subjected to a
high-intensity ultrashort electric pulse. Physical Review E 71:-.
78. Hu, Q., V. Sridhara, R. P. Joshi, J. F. Kolb, and K. H. Schoenbach. 2006.
Molecular dynamics analysis of high electric pulse effects on bilayer membranes
containing DPPC and DPPS. Ieee Transactions on Plasma Science 34:1405-1411.
79. Huang, J., J. E. Swanson, A. R. Dibble, A. K. Hinderliter, and G. W. Feigenson.
1993. Nonideal mixing of phosphatidylserine and phosphatidylcholine in the fluid
lamellar phase. Biophysical journal 64:413-425.
80. Humphrey, W., A. Dalke, and K. Schulten. 1996. VMD: visual molecular
dynamics. Journal of molecular graphics 14:33-38, 27-38.
81. Huster, D., K. Arnold, and K. Gawrisch. 2000. Strength of Ca(2+) binding to
retinal lipid membranes: consequences for lipid organization. Biophysical journal
78:3011-3018.
82. Idziorek, T., J. Estaquier, F. Debels, and J. C. Ameisen. 1995. Yopro-1 Permits
Cytofluorometric Analysis of Programmed Cell-Death (Apoptosis) without Interfering
with Cell Viability. Journal of Immunological Methods 185:249-258.
83. Imagawa, T., J. S. Smith, R. Coronado, and K. P. Campbell. 1987. Purified
Ryanodine Receptor from Skeletal-Muscle Sarcoplasmic-Reticulum Is the Ca-2+-
219
Permeable Pore of the Calcium Release Channel. Journal of Biological Chemistry
262:16636-16643.
84. Jacobson, K., O. G. Mouritsen, and R. G. Anderson. 2007. Lipid rafts: at a
crossroad between cell biology and physics. Nat Cell Biol 9:7-14.
85. Jogini, V., and B. Roux. 2007. Dynamics of the Kv1.2 voltage-gated K+ channel
in a membrane environment. Biophysical journal 93:3070-3082.
86. Jorgensen, W. L., and J. Tiradorives. 1988. The Opls Potential Functions for
Proteins - Energy Minimizations for Crystals of Cyclic-Peptides and Crambin. Journal of
the American Chemical Society 110:1657-1666.
87. Joshi, R. P., Q. Hu, R. Aly, K. H. Schoenbach, and H. P. Hjalmarson. 2001. Self-
consistent simulations of electroporation dynamics in biological cells subjected to
ultrashort electrical pulses. Physical Review E 6401:-.
88. Kandasamy, S. K., and R. G. Larson. 2006. Molecular dynamics simulations of
model trans-membrane peptides in lipid bilayers: a systematic investigation of
hydrophobic mismatch. Biophysical journal 90:2326-2343.
89. Kandasamy, S. K., and R. G. Larson. 2006. Effect of salt on the interactions of
antimicrobial peptides with zwitterionic lipid bilayers. Biochimica et biophysica acta
1758:1274-1284.
90. Kandasamy, S. K., and R. G. Larson. 2006. Cation and anion transport through
hydrophilic pores in lipid bilayers. J Chem Phys 125:-.
91. Kierzenka, J., and L. F. Shampine. 2001. A BVP solver based on residual control
and the MATLAB PSE. ACM Transactions on Mathematical Software 27:299-316.
92. Kinosita, K., and T. Y. Tsong. 1977. Voltage-Induced Pore Formation and
Hemolysis of Human Erythrocytes. Biochimica et biophysica acta 471:227-242.
93. Knecht, V., H. J. Risselada, A. E. Mark, and S. J. Marrink. 2008. Electrophoretic
mobility does not always reflect the charge on an oil droplet. Journal of Colloid and
Interface Science 318:477-486.
94. Koopman, G., C. P. M. Reutelingsperger, G. A. M. Kuijten, R. M. J. Keehnen, S.
T. Pals, and M. H. J. Vanoers. 1994. Annexin-V for Flow Cytometric Detection of
Phosphatidylserine Expression on B-Cells Undergoing Apoptosis. Blood 84:1415-1420.
220
95. Kotnik, T., and D. Miklavcic. 2000. Second-order model of membrane electric
field induced by alternating external electric fields. Ieee Transactions on Biomedical
Engineering 47:1074-1081.
96. Kotulska, M., K. Kubica, S. Koronkiewicz, and S. Kalinowski. 2007. Modeling
the induction of lipid membrane electropermeabilization. Bioelectrochemistry 70:64-70.
97. Kucerka, N., S. Tristram-Nagle, and J. F. Nagle. 2005. Structure of fully hydrated
fluid phase lipid bilayers with monounsaturated chains. The Journal of membrane biology
208:193-202.
98. Kucerka, N., S. Tristram-Nagle, and J. F. Nagle. 2006. Closer look at structure of
fully hydrated fluid phase DPPC bilayers. Biophysical journal 90:L83-85.
99. Kucerka, N., S. Tristram-Nagle, and J. F. Nagle. 2006. Structure of fully hydrated
fluid phase lipid bilayers with monounsaturated chains. Journal of Membrane Biology
208:193-202.
100. Kuthi, A., P. Gabrielsson, M. R. Behrend, P. T. Vernier, and M. A. Gundersen.
2005. Nanosecond pulse generator, using fast recovery diodes for cell
electromanipulation. Ieee Transactions on Plasma Science 33:1192-1197.
101. Lee, H. C., R. Aarhus, and T. F. Walseth. 1993. Calcium mobilization by dual
receptors during fertilization of sea urchin eggs. Science 261:352-355.
102. Lee, S. J., Y. Song, and N. A. Baker. 2008. Molecular dynamics simulations of
asymmetric NaCl and KCl solutions separated by phosphatidylcholine bilayers: potential
drops and structural changes induced by strong Na+-lipid interactions and finite size
effects. Biophysical journal 94:3565-3576.
103. Leekumjorn, S., and A. K. Sum. 2007. Molecular studies of the gel to liquid-
crystalline phase transition for fully hydrated DPPC and DPPE bilayers. Biochimica Et
Biophysica Acta-Biomembranes 1768:354-365.
104. Liao, M. J., and J. H. Prestegard. 1981. Structural properties of a Ca2+-
phosphatidic acid complex: small angle X-ray scattering and calorimetric results.
Biochimica et biophysica acta 645:149-156.
105. Lindahl, E., B. Hess, and D. van der Spoel. 2001. GROMACS 3.0: a package for
molecular simulation and trajectory analysis. J Mol Model 7:306-317.
106. MacKenzie, A., H. L. Wilson, E. Kiss-Toth, S. K. Dower, R. A. North, and A.
Surprenant. 2001. Rapid secretion of interleukin-1 beta by microvesicle shedding.
Immunity 15:825-835.
221
107. Mackenzie, A. B., M. T. Young, E. Adinolfi, and A. Surprenant. 2005.
Pseudoapoptosis induced by brief activation of ATP-gated P2X(7) receptors. Journal of
Biological Chemistry 280:33968-33976.
108. Malenkov, G. G., and M. M. FrankKamenetskii. 1996. Computer simulation of
the structure of the hydration shells of calcium and calcium-water-zeolite A. Journal of
Structural Chemistry 37:76-83.
109. Marrink, S. J., M. Berkowitz, and H. J. C. Berendsen. 1993. Molecular-Dynamics
Simulation of a Membrane Water Interface - the Ordering of Water and Its Relation to
the Hydration Force. Langmuir 9:3122-3131.
110. Marrink, S. J., J. Risselada, and A. E. Mark. 2005. Simulation of gel phase
formation and melting in lipid bilayers using a coarse grained model. Chemistry and
physics of lipids 135:223-244.
111. Marszalek, P., D. S. Liu, and T. Y. Tsong. 1990. Schwan Equation and
Transmembrane Potential Induced by Alternating Electric-Field. Biophysical journal
58:1053-1058.
112. Mattai, J., H. Hauser, R. A. Demel, and G. G. Shipley. 1989. Interactions of metal
ions with phosphatidylserine bilayer membranes: effect of hydrocarbon chain
unsaturation. Biochemistry 28:2322-2330.
113. McLaughlin, S., N. Mulrine, T. Gresalfi, G. Vaio, and A. McLaughlin. 1981.
Adsorption of divalent cations to bilayer membranes containing phosphatidylserine. The
Journal of general physiology 77:445-473.
114. McLaughlin, S. G., G. Szabo, and G. Eisenman. 1971. Divalent ions and the
surface potential of charged phospholipid membranes. The Journal of general physiology
58:667-687.
115. Megyes, T., I. Bako, S. Balint, T. Grosz, and T. Radnai. 2006. Ion pairing in
aqueous calcium chloride solution: Molecular dynamics simulation and diffraction
studies. Journal of Molecular Liquids 129:63-74.
116. Mir, L. M., M. F. Bureau, J. Gehl, R. Rangara, D. Rouy, J. M. Caillaud, P.
Delaere, D. Branellec, B. Schwartz, and D. Scherman. 1999. High-efficiency gene
transfer into skeletal muscle mediated by electric pulses. Proceedings of the National
Academy of Sciences of the United States of America 96:4262-4267.
117. Miyamoto, S., and P. A. Kollman. 1992. Settle - an Analytical Version of the
Shake and Rattle Algorithm for Rigid Water Models. J Comput Chem 13:952-962.
222
118. Mukhopadhyay, P., L. Monticelli, and D. P. Tieleman. 2004. Molecular dynamics
simulation of a palmitoyl-oleoyl phosphatidylserine bilayer with Na+ counterions and
NaCl. Biophysical journal 86:1601-1609.
119. Muller, K. J., V. L. Sukhorukov, and U. Zimmermann. 2001. Reversible
electropermeabilization of mammalian cells by high-intensity, ultra-short pulses of
submicrosecond duration. Journal of Membrane Biology 184:161-170.
120. Murtola, T., E. Falck, M. Patra, M. Karttunen, and I. Vattulainen. 2004. Coarse-
grained model for phospholipid/cholesterol bilayer. J Chem Phys 121:9156-9165.
121. Nagle, J. F., R. Zhang, S. Tristram-Nagle, W. Sun, H. I. Petrache, and R. M.
Suter. 1996. X-ray structure determination of fully hydrated L alpha phase
dipalmitoylphosphatidylcholine bilayers. Biophysical journal 70:1419-1431.
122. Nagle, J. F., and S. Tristram-Nagle. 2000. Structure of lipid bilayers. Biochimica
Et Biophysica Acta-Reviews on Biomembranes 1469:159-195.
123. Neu, J. C., K. C. Smith, and W. Krassowska. 2003. Electrical energy required to
form large conducting pores. Bioelectrochemistry 60:107-114.
124. Neumann, E., and K. Rosenheck. 1972. Permeability changes induced by electric
impulses in vesicular membranes. The Journal of membrane biology 10:279-290.
125. Neumann, E., M. Schaeferridder, Y. Wang, and P. H. Hofschneider. 1982. Gene-
Transfer into Mouse Lyoma Cells by Electroporation in High Electric-Fields. Embo
Journal 1:841-845.
126. Nir, S., J. Bentz, and J. Wilschut. 1980. Mass-Action Kinetics of
Phosphatidylserine Vesicle Fusion as Monitored by Coalescence of Internal Vesicle
Volumes. Biochemistry 19:6030-6036.
127. Nuccitelli, R., U. Pliquett, X. H. Chen, W. Ford, R. J. Swanson, S. J. Beebe, J. F.
Kolb, and K. H. Schoenbach. 2006. Nanosecond pulsed electric fields cause melanomas
to self-destruct. Biochemical and Biophysical Research Communications 343:351-360.
128. Nymeyer, H., and H. X. Zhou. 2008. A method to determine dielectric constants
in nonhomogeneous systems: Application to biological membranes. Biophysical journal
94:1185-1193.
129. O'Shea, P. 2003. Intermolecular interactions with/within cell membranes and the
trinity of membrane potentials: kinetics and imaging. Biochem Soc Trans 31:990-996.
223
130. Onsager, L., and R. M. Fuoss. 1932. Irreversible processes in electrolytes
diffusion, conductance, and viscous flow in arbitrary mixtures of strong electrolytes.
Journal of Physical Chemistry 36:2689-2778.
131. Pabst, G., A. Hodzic, J. Strancar, S. Danner, M. Rappolt, and P. Laggner. 2007.
Rigidification of neutral lipid bilayers in the presence of salts. Biophysical journal
93:2688-2696.
132. Palinkas, G., and K. Heinzinger. 1986. Hydration Shell Structure of the Calcium-
Ion. Chemical Physics Letters 126:251-254.
133. Pan, J., S. Tristram-Nagle, N. Kucerka, and J. F. Nagle. 2008. Temperature
dependence of structure, bending rigidity, and bilayer interactions of
dioleoylphosphatidylcholine bilayers. Biophysical journal 94:117-124.
134. Pandit, S. A., and M. L. Berkowitz. 2002. Molecular dynamics simulation of
dipalmitoylphosphatidylserine bilayer with Na+ counterions. Biophysical journal
82:1818-1827.
135. Pandit, S. A., D. Bostick, and M. L. Berkowitz. 2003. Mixed bilayer containing
dipalmitoylphosphatidylcholine and dipalmitoylphosphatidylserine: lipid complexation,
ion binding, and electrostatics. Biophysical journal 85:3120-3131.
136. Pandit, S. A., D. Bostick, and M. L. Berkowitz. 2003. Molecular dynamics
simulation of a dipalmitoylphosphatidylcholine bilayer with NaCl. Biophysical journal
84:3743-3750.
137. Pandit, S. A., D. Bostick, and M. L. Berkowitz. 2004. Complexation of
phosphatidylcholine lipids with cholesterol. Biophysical journal 86:1345-1356.
138. Papahadjopoulos, D., S. Nir, and N. Duzgunes. 1990. Molecular Mechanisms of
Calcium-Induced Membrane-Fusion. Journal of Bioenergetics and Biomembranes
22:157-179.
139. Pedersen, U. R., C. Leidy, P. Westh, and G. H. Peters. 2006. The effect of
calcium on the properties of charged phospholipid bilayers. Biochimica et biophysica
acta 1758:573-582.
140. Peitzsch, R. M., M. Eisenberg, K. A. Sharp, and S. McLaughlin. 1995.
Calculations of the electrostatic potential adjacent to model phospholipid bilayers.
Biophysical journal 68:729-738.
224
141. Petrache, H. I., S. Tristram-Nagle, K. Gawrisch, D. Harries, V. A. Parsegian, and
J. F. Nagle. 2004. Structure and fluctuations of charged phosphatidylserine bilayers in the
absence of salt. Biophysical journal 86:1574-1586.
142. Petrache, H. I., I. Kimchi, D. Harries, and V. A. Parsegian. 2005. Measured
depletion of ions at the biomembrane interface. Journal of the American Chemical
Society 127:11546-11547.
143. Phillips, J. C., R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C.
Chipot, R. D. Skeel, L. Kale, and K. Schulten. 2005. Scalable molecular dynamics with
NAMD. J Comput Chem 26:1781-1802.
144. Plonsey, R., and K. W. Altman. 1988. Electrical-Stimulation of Excitable Cells - a
Model Approach. Proceedings of the Ieee 76:1122-1129.
145. Popescu, D., C. Rucareanu, and G. Victor. 1991. A Model for the Appearance of
Statistical Pores in Membranes Due to Selfoscillations. Bioelectrochemistry and
Bioenergetics 25:91-103.
146. Popescu, D., and G. Victor. 1991. The Transversal Diffusion-Coefficient of
Phospholipid Molecules through Black Lipid-Membranes. Bioelectrochemistry and
Bioenergetics 25:105-108.
147. Portis, A., C. Newton, W. Pangborn, and D. Papahadjopoulos. 1979. Studies on
the mechanism of membrane fusion: evidence for an intermembrane Ca2+-phospholipid
complex, synergism with Mg2+, and inhibition by spectrin. Biochemistry 18:780-790.
148. Rehfeld, S. J., N. Duzgunes, C. Newton, D. Papahadjopoulos, and D. J. Eatough.
1981. The exothermic reaction of calcium with unilamellar phosphatidylserine vesicles:
titration microcalorimetry. FEBS Lett 123:249-251.
149. Rols, M. P. 2006. Electropermeabilization, a physical method for the delivery of
therapeutic molecules into cells. Biochimica Et Biophysica Acta-Biomembranes
1758:423-428.
150. Roux, M., and M. Bloom. 1990. Ca2+, Mg2+, Li+, Na+, and K+ distributions in
the headgroup region of binary membranes of phosphatidylcholine and
phosphatidylserine as seen by deuterium NMR. Biochemistry 29:7077-7089.
151. Roux, M., and M. Bloom. 1991. Calcium binding by phosphatidylserine
headgroups. Deuterium NMR study. Biophysical journal 60:38-44.
152. Sachs, J. N., P. S. Crozier, and T. B. Woolf. 2004. Atomistic simulations of
biologically realistic transmembrane potential gradients. J Chem Phys 121:10847-10851.
225
153. Sachs, J. N., H. Nanda, H. I. Petrache, and T. B. Woolf. 2004. Changes in
phosphatidylcholine headgroup tilt and water order induced by monovalent salts:
molecular dynamics simulations. Biophysical journal 86:3772-3782.
154. Saiz, L., and M. L. Klein. 2002. Electrostatic interactions in a neutral model
phospholipid bilayer by molecular dynamics simulations. J Chem Phys 116:3052-3057.
155. Sale, A. J. H., and W. A. Hamilton. 1967. Effects of High Electric Fields on
Microorganisms .I. Killing of Bacteria and Yeasts. Biochimica et biophysica acta
148:781-&.
156. Sale, A. J. H., and W. A. Hamilton. 1968. Effects of High Electric Fields on
Micro-Organisms .3. Lysis of Erythrocytes and Protoplasts. Biochimica et biophysica
acta 163:37-&.
157. Sansom, M. S. P., J. Breed, H. J. C. Berendsen, and P. Tieleman. 1998.
Alamethicin channels - Simulation studies. Biophysical journal 74:A232-A232.
158. Sansom, M. S. P., D. P. Tieleman, L. R. Forrest, and H. J. C. Berendsen. 1998.
Molecular dynamics simulations of membranes with embedded proteins and peptides:
porin, alamethicin and influenza virus M2. Biochemical Society Transactions 26:438-
443.
159. Schoenbach, K. H., S. J. Beebe, and E. S. Buescher. 2001. Intracellular effect of
ultrashort electrical pulses. Bioelectromagnetics 22:440-448.
160. Schoenbach, K. H., R. P. Joshi, J. F. Kolb, N. Y. Chen, M. Stacey, P. F.
Blackmore, E. S. Buescher, and S. J. Beebe. 2004. Ultrashort electrical pulses open a new
gateway into biological cells. Proceedings of the Ieee 92:1122-1137.
161. Schote, U., and J. Seelig. 1998. Interaction of the neuronal marker dye FM1-43
with lipid membranes - Thermodynamics and lipid ordering. Biochimica Et Biophysica
Acta-Biomembranes 1415:135-146.
162. Schwan, H. P. 1957. Electrical properties of tissue and cell suspensions. In
Advances in Biological and Medical Physics. J. H. Lawrence, and C. A. Tobias, editors.
Academic Press, New York.
163. Schwan, H. P., and K. R. Foster. 1977. Microwave Dielectric Properties of Tissue
- Some Comments on Rotational Mobility of Tissue Water. Biophysical journal 17:193-
197.
226
164. Seelig, A., and J. Seelig. 1974. The dynamic structure of fatty acyl chains in a
phospholipid bilayer measured by deuterium magnetic resonance. Biochemistry 13:4839-
4845.
165. Seelig, J. 1970. Spin Label Studies of Oriented Smectic Liquid Crystals - a Model
System for Bilayer Membranes. Journal of the American Chemical Society 92:3881-&.
166. Seelig, J., P. M. Macdonald, and P. G. Scherer. 1987. Phospholipid Head Groups
as Sensors of Electric Charge in Membranes. Biochemistry 26:7535-7541.
167. Seelig, J. 1990. Interaction of phospholipids with Ca2+ ions. On the role of the
phospholipid head groups. Cell biology international reports 14:353-360.
168. Seigneuret, M., and P. F. Devaux. 1984. Atp-Dependent Asymmetric Distribution
of Spin-Labeled Phospholipids in the Erythrocyte-Membrane - Relation to Shape
Changes. Proceedings of the National Academy of Sciences of the United States of
America-Biological Sciences 81:3751-3755.
169. Sham, Y. Y., Z. T. Chu, H. Tao, and A. Warshel. 2000. Examining methods for
calculations of binding free energies: LRA, LIE, PDLD-LRA, and PDLD/S-LRA
calculations of ligands binding to an HIV protease. Proteins 39:393-407.
170. Sher, L. D., E. Kresch, and H. P. Schwan. 1970. On Possibility of Nonthermal
Biological Effects of Pulsed Electromagnetic Radiation. Biophysical journal 10:970-&.
171. Shinoda, K., W. Shinoda, and M. Mikami. 2007. Molecular dynamics simulation
of an archaeal lipid bilayer with sodium chloride. Phys Chem Chem Phys 9:643-650.
172. Sinn, C. G., M. Antonietti, and R. Dimova. 2006. Binding of calcium to
phosphatidylcholine-phosphatidylserine membranes. Colloid Surface A 282:410-419.
173. Siu, S. W., R. Vacha, P. Jungwirth, and R. A. Bockmann. 2008. Biomolecular
simulations of membranes: physical properties from different force fields. J Chem Phys
128:125103.
174. Smith, K. C., T. R. Gowrishankar, A. T. Esser, D. A. Stewart, and J. C. Weaver.
2006. The spatially distributed dynamic transmembrane voltage of cells and organelles
due to 10-ns pulses: Meshed transport networks. Ieee Transactions on Plasma Science
34:1394-1404.
175. Spohr, E. 1997. Effect of electrostatic boundary conditions and system size on the
interfacial properties of water and aqueous solutions. J Chem Phys 107:6342-6348.
227
176. Stewart, D. A., I. R. Gowrishankar, and J. C. Weaver. 2004. Transport lattice
approach to describing cell electroporation: Use of a local asymptotic model. Ieee
Transactions on Plasma Science 32:1696-1708.
177. Sugar, I. P. 1981. The Effects of External Fields on the Structure of Lipid
Bilayers. Journal De Physiologie 77:1035-1042.
178. Sugar, I. P., and E. Neumann. 1984. Stochastic model for electric field-induced
membrane pores. Electroporation. Biophysical chemistry 19:211-225.
179. Sun, Y. H., P. T. Vernier, M. Behrend, L. Marcu, and M. A. Gundersen. 2005.
Electrode microchamber for noninvasive perturbation of mammalian cells with
nanosecond pulsed electric fields. Ieee Transactions on Nanobioscience 4:277-283.
180. Szule, J. A., S. E. Jarvis, J. E. Hibbert, J. D. Spafford, J. E. A. Braun, G. W.
Zamponni, G. M. Wessel, and J. R. Coorssen. 2003. Calcium-triggered membrane fusion
proceeds independently of specific presynaptic proteins. Journal of Biological Chemistry
278:24251-24254.
181. Tanaka, Y., and A. J. Schroit. 1986. Calcium/phosphate-induced immobilization
of fluorescent phosphatidylserine in synthetic bilayer membranes: inhibition of lipid
transfer between vesicles. Biochemistry 25:2141-2148.
182. Tarek, M. 2005. Membrane electroporation: a molecular dynamics simulation.
Biophysical journal 88:4045-4053.
183. Tatulian, S. A. 1987. Binding of alkaline-earth metal cations and some anions to
phosphatidylcholine liposomes. European journal of biochemistry / FEBS 170:413-420.
184. Taupin, C., M. Dvolaitzky, and C. Sauterey. 1975. Osmotic-Pressure Induced
Pores in Phospholipid Vesicles. Biochemistry 14:4771-4775.
185. Teissie, J., and T. Y. Tsong. 1981. Electric field induced transient pores in
phospholipid bilayer vesicles. Biochemistry 20:1548-1554.
186. Teissie, J., and M. P. Rols. 1993. An Experimental Evaluation of the Critical
Potential Difference Inducing Cell-Membrane Electropermeabilization. Biophysical
journal 65:409-413.
187. Teissie, J. 2007. Biophysical effects of electric fields on membrane water
interfaces: a mini review. Eur Biophys J 36:967-972.
228
188. Tekle, E., H. Oubrahim, S. M. Dzekunov, J. F. Kolb, K. H. Schoenbach, and P. B.
Chock. 2005. Selective field effects on intracellular vacuoles and vesicle membranes with
nanosecond electric pulses. Biophysical journal 89:274-284.
189. Tieleman, D. P., and H. J. C. Berendsen. 1996. Molecular dynamics simulations
of a fully hydrated dipalmitoyl phosphatidylcholine bilayer with different macroscopic
boundary conditions and parameters. J Chem Phys 105:4871-4880.
190. Tieleman, D. P., and H. J. C. Berendsen. 1998. A molecular dynamics study of
the pores formed by Escherichia coli OmpF porin in a fully hydrated
palmitoyloleoylphosphatidylcholine bilayer. Biophysical journal 74:2786-2801.
191. Tieleman, D. P., J. Breed, H. J. C. Berendsen, and M. S. P. Sansom. 1998.
Alamethicin channels in a membrane: molecular dynamics simulations. Faraday
Discussions:209-223.
192. Tieleman, D. P., H. J. C. Berendsen, and M. S. P. Sansom. 2001. Voltage-
dependent insertion of alamethicin at phospholipid/water and octane/water interfaces.
Biophysical journal 80:331-346.
193. Tieleman, D. P., H. Leontiadou, A. E. Mark, and S. J. Marrink. 2003. Simulation
of pore formation in lipid bilayers by mechanical stress and electric fields. Journal of the
American Chemical Society 125:6382-6383.
194. Tieleman, D. P. 2004. The molecular basis of electroporation. BMC biochemistry
5:10.
195. Tieleman, D. P., and S. J. Marrink. 2006. Lipids out of equilibrium: Energetics of
desorption and pore mediated flip-flop. Journal of the American Chemical Society
128:12462-12467.
196. Timoshkin, I. V., S. J. MacGregor, R. A. Fouracre, B. H. Crichton, and J. G.
Anderson. 2006. Transient electrical field across cellular membranes: pulsed electric field
treatment of microbial cells. Journal of Physics D-Applied Physics 39:596-603.
197. van der Ploeg, P., and H. J. C. Berendsen. 1982. Molecular-Dynamics Simulation
of a Bilayer-Membrane. J Chem Phys 76:3271-3276.
198. van der Spoel, D., E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. C.
Berendsen. 2005. GROMACS: Fast, flexible, and free. J Comput Chem 26:1701-1718.
199. Vasilkoski, Z., A. T. Esser, T. R. Gowrishankar, and J. C. Weaver. 2006.
Membrane electroporation: The absolute rate equation and nanosecond time scale pore
creation. Physical Review E 74:-.
229
200. Verhoven, B., R. A. Schlegel, and P. Williamson. 1995. Mechanisms of
Phosphatidylserine Exposure, a Phagocyte Recognition Signal, on Apoptotic T-
Lymphocytes. Journal of Experimental Medicine 182:1597-1601.
201. Vernier, P. T., A. M. Li, L. Marcu, C. M. Craft, and M. A. Gundersen. 2003.
Ultrashort pulsed electric fields induce membrane phospholipid translocation and caspase
activation: Differential sensitivities of Jurkat T lymphoblasts and rat glioma C6 cells.
Ieee Transactions on Dielectrics and Electrical Insulation 10:795-809.
202. Vernier, P. T., Y. H. Sun, L. Marcu, S. Salemi, C. M. Craft, and M. A.
Gundersen. 2003. Calcium bursts induced by nanosecond electric pulses. Biochemical
and Biophysical Research Communications 310:286-295.
203. Vernier, P. T., Y. Sun, L. Marcu, C. M. Craft, and M. A. Gundersen. 2004.
Nanoelectropulse-induced phosphatidylserine translocation. Biophysical journal 86:4040-
4048.
204. Vernier, P. T., Y. H. Sun, L. Marcu, C. M. Craft, and M. A. Gundersen. 2004.
Nanosecond pulsed electric fields perturb membrane phospholipids in T lymphoblasts.
Febs Letters 572:103-108.
205. Vernier, P. T., Y. Sun, J. Wang, M. M. Thu, E. Garon, L. Valderrabano, L.
Marcu, H. P. Koeffler, and M. A. Gundersen. 2005. Nanoelectropulse intracellular
perturbation and electropermeabilization technology: phospholipid translocation, calcium
bursts, chromatin rearrangement, cardiomyocyte activation, and tumor cell sensitivity. In
IEEE Engineering in Medicine and Biology Society 27th International Conference,
Shanghai.
206. Vernier, P. T., Y. H. Sun, and M. A. Gundersen. 2006. Nanoelectropulse-driven
membrane perturbation and small molecule permeabilization. Bmc Cell Biology 7:-.
207. Vernier, P. T., M. J. Ziegler, Y. Sun, M. A. Gundersen, and D. P. Tieleman. 2006.
Nanopore-facilitated, voltage-driven phosphatidylserine translocation in lipid bilayers--in
cells and in silico. Phys Biol 3:233-247.
208. Vernier, P. T., M. J. Ziegler, Y. H. Sun, W. V. Chang, M. A. Gundersen, and D.
P. Tieleman. 2006. Nanopore formation and phosphatidylserine externalization in a
phospholipid bilayer at high transmembrane potential. Journal of the American Chemical
Society 128:6288-6289.
209. Vernier, P. T., and M. J. Ziegler. 2007. Nanosecond field alignment of head group
and water dipoles in electroporating phospholipid bilayers. The journal of physical
chemistry 111:12993-12996.
230
210. Vernier, P. T., M. J. Ziegler, and R. Dimova. Submitted. Calcium Binding and
Head Group Dipole Angle in Phosphatidylserine-Phosphatidylcholine Bilayers.
Biophysical journal.
211. Wang, J. H. 1953. Tracer-Diffusion in Liquids .4. Self-Diffusion of Calcium Ion
and Chloride Ion in Aqueous Calcium Chloride Solutions. Journal of the American
Chemical Society 75:1769-1770.
212. Wang, Q. F., L. Q. Wang, Y. H. Feng, X. Li, R. Zeng, and G. I. Gorodeski. 2004.
P2X(7) receptor-mediated apoptosis of human cervical epithelial cells. American Journal
of Physiology-Cell Physiology 287:C1349-C1358.
213. Warshel, A., P. K. Sharma, M. Kato, and W. W. Parson. 2006. Modeling
electrostatic effects in proteins. Biochimica et biophysica acta 1764:1647-1676.
214. Warshel, A., M. Kato, and A. V. Pisliakov. 2007. Polarizable force fields:
History, test cases, and prospects. J Chem Theory Comput 3:2034-2045.
215. Weaver, J. C., and R. A. Mintzer. 1981. Decreased Bilayer Stability Due to
Transmembrane Potentials. Physics Letters A 86:57-59.
216. Weaver, J. C., K. T. Powell, R. A. Mintzer, H. Ling, and S. R. Sloan. 1984. The
Electrical Capacitance of Bilayer-Membranes the Contribution of Transient Aqueous
Pores. Bioelectrochemistry and Bioenergetics 12:393-404.
217. Weaver, J. C., K. T. Powell, R. A. Mintzer, S. R. Sloan, and H. Ling. 1984. Bis -
the Diffusive Permeability of Bilayer-Membranes the Contribution of Transient Aqueous
Pores. Bioelectrochemistry and Bioenergetics 12:405-412.
218. Weaver, J. C., and Y. A. Chizmadzhev. 1996. Theory of electroporation: A
review. Bioelectrochemistry and Bioenergetics 41:135-160.
219. White, J. A., P. F. Blackmore, K. H. Schoenbach, and S. J. Beebe. 2004.
Stimulation of capacitative calcium entry in HL-60 cells by nanosecond pulsed electric
fields. J Biol Chem 279:22964-22972.
220. Williamson, P., A. Christie, T. Kohlin, R. A. Schlegel, P. Comfurius, M.
Harmsma, R. F. A. Zwaal, and E. M. Bevers. 2001. Phospholipid scramblase activation
pathways in lymphocytes. Biochemistry 40:8065-8072.
221. Wohlert, J., and O. Edholm. 2004. The range and shielding of dipole-dipole
interactions in phospholipid bilayers. Biophysical journal 87:2433-2445.
231
222. Wurth, G. A., and A. Zweifach. 2002. Evidence that cytosolic calcium increases
are not sufficient to stimulate phospholipid scrambling in human T-lymphocytes.
Biochemical Journal 362:701-708.
223. Yeung, T., G. E. Gilbert, J. Shi, J. Silvius, A. Kapus, and S. Grinstein. 2008.
Membrane phosphatidylserine regulates surface charge and protein localization. Science
319:210-213.
224. Zhao, W., T. Rog, A. A. Gurtovenko, I. Vattulainen, and M. Karttunen. 2007.
Atomic-scale structure and electrostatics of anionic palmitoyloleoylphosphatidylglycerol
lipid bilayers with Na+ counterions. Biophysical journal 92:1114-1124.
225. Zheng, C., and G. Vanderkooi. 1992. Molecular origin of the internal dipole
potential in lipid bilayers: calculation of the electrostatic potential. Biophysical journal
63:935-941.
226. Zu, Z. 2008. Personal Communication.
227. Zweifach, A. 2000. FM1-43 reports plasma membrane phospholipid scrambling
in T-lymphocytes. Biochemical Journal 349:255-260.
228. Zweifach, A. 2000. Target-cell contact activates a highly selective capacitative
calcium entry pathway in cytotoxic T lymphocytes. J Cell Biol 148:603-614.
Abstract (if available)
Abstract
Recent advances in computing technology have facilitated the application of simulations to studying biological systems at the atomic level. In particular atomistic molecular dynamics provide an opportunity to model systems that are unobservable through conventional experimental methods as well as supplement understanding of observations. In this thesis molecular dynamics were applied to study biological cell membranes, specifically lipid bilayers, the primary constituent of the cell membrane, in electric fields, and to understand the mechanism and events associated with electroporation. Electroporation is a widely used experimental and commercial technique for introducing normally excluded compounds such as DNA, RNA, ions, drugs, and other chemicals into cells. Traditional electroporation utilizes kilovolt-per-meter electric fields applied on the order of microseconds that disrupt and scramble the cell membrane and allow diffusive entry, however, ultra-short nanosecond pulses at megavolt-per-meter fields produce different effects in cells which remain largely uncharacterized. One such effect, the migration of the negatively charged lipid phosphatidylserine from the inner to outer leaflet of the cell, is of particular biological interest because of its association with cell apoptosis, or programmable cell death. Control of such an event could be useful in developing a targeted treatment for removing unwanted cells such as tumors or melanoma. In this thesis I begin by introducing electroporation and its history and explain how molecular dynamics and its techniques can help advance our understanding of the field. In the following chapters, consisting of peer reviewed and submitted journal articles loosely tied together, I explore the mechanism of phosphatidylserine translocation induced by nanosecond pulses in megavolt-per-meter electric fields, and correlate experimental data and anecdotal evidence of phosphatidylserine translocation in vitro with a detailed molecular
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations
PDF
Molecular simulations of water and monovalent ion dynamics in the electroporation of phospholipid bilayers
PDF
Biophysical studies of passive transport across synthetic lipid bilayers
PDF
Biophysical studies of nanosecond pulsed electric field induced cell membrane permeabilization
PDF
The effects of oxidation on bilayer membranes studied using giant unilamellar vesicles
PDF
Simulating the helicase motor of SV40 large tumor antigen
PDF
Molecular dynamics simulation study of initial protein unfolding induced by the photo-responsive surfactants, azoTAB
PDF
Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations
PDF
Probing the effects of transmembrane domains on the continuum mechanics of lipid bilayers
PDF
Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations
PDF
Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials
PDF
A molecular dynamics study of interactions between the enzyme lysozyme and the photo-responsive surfactant azobenzene trimethylammonium bromide (azoTAB)
PDF
High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks
PDF
Physical principles of membrane mechanics, membrane domain formation, and cellular signal transduction
PDF
Molecular dynamics studies of protein aggregation in unbounded and confined media
PDF
Multiscale and multiresolution approach to characterization and modeling of porous media: From pore to field scale
PDF
Neural network for molecular dynamics simulation and design of 2D materials
PDF
Dynamics of water in nanotubes: liquid below freezing point and ice-like near boiling point
PDF
Structure-based computational analysis and prediction of TCR CDR3 loops in the TCR-peptide-MHC complex using solvation parameters and peptide molecular dynamics.
PDF
Photoexcitation and nonradiative relaxation in molecular systems: methodology, optical properties, and dynamics
Asset Metadata
Creator
Ziegler, Matthew James
(author)
Core Title
Molecular dynamics simulations of lipid bilayers in megavolt per meter electric fields
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Chemical Engineering
Publication Date
12/12/2008
Defense Date
10/06/2008
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
calcium binding,electroporation,GROMACS,headgroup dipole angle,kinetics,Lipids,molecular dynamics,OAI-PMH Harvest,phosphatidylcholine,phosphatidylserine,simulation,translocation
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Gundersen, Martin A. (
committee chair
), Vernier, P. Thomas (
committee chair
), Shing, Katherine S. (
committee member
), Warshel, Arieh (
committee member
)
Creator Email
mjzglr@gmail.com,mjziegle@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m1920
Unique identifier
UC1132233
Identifier
etd-Ziegler-2481 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-150169 (legacy record id),usctheses-m1920 (legacy record id)
Legacy Identifier
etd-Ziegler-2481.pdf
Dmrecord
150169
Document Type
Dissertation
Rights
Ziegler, Matthew James
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
calcium binding
electroporation
GROMACS
headgroup dipole angle
kinetics
molecular dynamics
phosphatidylcholine
phosphatidylserine
simulation
translocation