University of Southern California Dissertations and Theses

Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations

(USC Thesis Other)

Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content
THEORETICAL STUDIES OF LIPID BILAYER ELECTROPORATION USING
MOLECULAR DYNAMICS SIMULATIONS

by

Zachary Alan Levine

A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(PHYSICS)

August 2013

Copyright 2013 Zachary Alan Levine
ii

Dedication
To my wife, Morgan.

iii

Acknowledgements
I would like to express my deepest thanks and gratitude to my adviser, mentor, and friend
Paul Thomas Vernier for his constant support and guidance over the past five and a half years.
Your diligent guidance through the worlds of chemistry, biology, and molecular dynamics
simulations has undoubtedly made me a more well-rounded scientist. I also thank you for being
uncompromising when it came to making sure each and every detail was accounted for. I’ve
learned that this separates the good things from the truly great things, and I’m proud of what we
were able to accomplish in just a few short years.
There are also many others I would like to thank - Martin Gundersen for his invaluable
advice and mentorship throughout my time at USC; Stephan Haas for his friendship, sense of
humor, and trust in my abilities from day one; Aiichiro Nakano and the Collaboratory for
Advanced Computing and Simulations for their ongoing encouragement and support, especially
while I was chasing a computer science diploma; Noah Malmstadt for participating as my
outside committee member; and Moh El-Naggar and Hubert Saleur for stimulating conversations
and ongoing friendships.
I was also privileged to have worked with giants in the field such as Peter Tieleman,
James Weaver, Lluis Mir, Justin Teissié, Andrei Pakhomov, and Damijan Miklavčič to name a
few. On the computing side I want to thank USC’s High Performance Computing and
Communications Linux Cluster for allowing me to run our simulations on a world class parallel
computing cluster. I’m also grateful to the people at MOSIS and ISI for providing me with
countless opportunities, funding, and friendships throughout the years, in addition to a beautiful
iv

office in Marina del Rey that no graduate student deserves. I also want to thank Brad Cleaver
who has continuously helped me navigate through all of the treacherous technical experiments
I’ve encountered, and for helping me take baby steps into the UNIX and Linux multiverse.
I would also like to give special thanks to my family who has had to endure vague
descriptions of what I research for the past few years. I hope you understand that if I decide to
remain in academia, that description will only become longer and more difficult to pronounce.
You can direct your frustrations at all of the people mentioned above. Jokes aside, your support
has meant more to me than I could possibly convey on pieces of paper, and I am forever grateful
to each and every one of you.
I am also indebted to a number of people who are typically not thanked enough,
undergraduate professors and high school teachers that have put in extraordinary efforts to see
me succeed. To Barbara Neuhauser at San Francisco State University, your confidence in me as
a young scientist has helped shape the rest of my life and I am genuinely grateful for your words
and your time. To Jeff Ramstedt and Gary Reinstein of Santa Susana High School, thank you for
the inspiring books and enthusiasm you gave to me and to all of your students about math and
science on a daily basis. Your commitment to my education, and the education of others, has
helped shape the world in immeasurable ways.
And to my wife Morgan who has never stopped believing in me and supporting me, thank
you for enduring the countless nights where I was away writing or doing research. I’m excited to
start the next chapter of our lives outside of the confines of graduate school where the world
awaits our impact. May our adventures continue to be memorable and full of joy.
v

Table of Contents

Dedication ii
Acknowledgements iii
List of Tables viii
List of Figures ix
Abstract xii
Chapter 1: Introduction 1
1.1 Phospholipids and Phospholipid Bilayers ...........................................................................1
1.2 Electroporation ....................................................................................................................5
1.2.1 Overview ..................................................................................................................5
1.2.2 Continuum Models...................................................................................................7
1.3 Molecular Dynamics Simulations .....................................................................................11
1.3.1 Overview ................................................................................................................11
1.3.2 Implementation ......................................................................................................15
1.3.3 Water Models .........................................................................................................16
1.3.4 Phospholipid Models .............................................................................................18
1.3.5 Energy Minimization .............................................................................................21
1.3.6 Electrostatics ..........................................................................................................22
1.3.7 Thermodynamic Ensembles ...................................................................................24
1.3.8 Restraints and Holonomic Constraints...................................................................28
1.3.9 Analysis..................................................................................................................30

Chapter 2: Effect of Oxidized Lipids on Electroporation 31
2.1 Abstract .............................................................................................................................31
2.2 Introduction .......................................................................................................................32
2.3 Methods.............................................................................................................................33
2.3.1 Molecular Dynamics Parameters ...........................................................................33
2.3.2 Electroporation Simulations...................................................................................34
2.3.3 Assembly of Larger Systems .................................................................................35
2.3.4 Cell lines and Preparation ......................................................................................36
2.3.5 Pulse Generator and Pulse Exposures ....................................................................37
2.3.6 Fluorescence Microscopy and Microphotometry ..................................................37
2.4 Results ...............................................................................................................................38
2.4.1 Simulations .............................................................................................................38
2.4.2 Observations with Living Cells ..............................................................................44
2.5 Discussion .........................................................................................................................47
2.5.1 Electrotransfection and Electropermeabilization Enhancement .............................48
vi

2.5.2 Oxidative Stress and Apoptosis ..............................................................................49
2.5.3 Permeabilizing Defects and Membrane Boundaries ...............................................50
2.6 Acknowledgements ...........................................................................................................50

Chapter 3: Life Cycle of an Electropore 51
3.1 Abstract .............................................................................................................................51
3.2 Introduction .......................................................................................................................52
3.3 Materials and Methods ......................................................................................................53
3.3.1 Simulation .............................................................................................................53
3.3.2 Systems and Structures ..........................................................................................54
3.3.3 Electropore Life Cycle ...........................................................................................55
3.3.4 Pore Creation .........................................................................................................56
3.3.5 Pore Annihilation ...................................................................................................57
3.4 Results ...............................................................................................................................60
3.4.1 Pore Creation and Annihilation – POPC ................................................................60
3.4.2 Pore Creation – POPC versus DOPC......................................................................63
3.4.3 Pore Creation – Old DOPC models versus New DOPC models ............................65
3.5 Discussion .........................................................................................................................66
3.5.1 Stages of Electroporation ........................................................................................67
3.5.2 DOPC versus POPC ................................................................................................69
3.5.3 Deficiencies in the Electropore Life Cycle Scheme ...............................................70
3.6 Acknowledgements ...........................................................................................................72

Chapter 4: Calcium and Phosphatidylserine affect Pore Life Cycle 73
4.1 Abstract .............................................................................................................................73
4.2 Introduction .......................................................................................................................74
4.3 Materials and Methods ......................................................................................................76
4.3.1 Structures and Parameters ......................................................................................76
4.3.2 Updates to the Electropore Life Cycle ...................................................................78
4.4 Results ...............................................................................................................................81
4.4.1 Pore Creation – POPC and POPC:POPS ................................................................87
4.4.2 Pore Creation – POPC and Calcium .......................................................................87
4.4.3 Pore Creation – Calcium and POPC:POPS ............................................................88
4.4.4 Pore Annihilation – POPC and POPC:POPS..........................................................88
4.4.5 Pore Annihilation –POPC and Calcium ..................................................................90
4.4.6 Pore Annihilation – Calcium and POPC:POPS ......................................................90
4.4.7 Calcium Binding Isotherms ....................................................................................91
4.5 Discussion .........................................................................................................................92
4.5.1 Electropore Life Cycle for POPC:POPS Bilayers ..................................................92
4.5.2 Effects of Calcium on Pore Life Cycle ...................................................................93
4.5.3 Molecular and Continuum Models of Lipid Electropores ......................................95
4.6 Acknowledgements ...........................................................................................................96

vii

Chapter 5: The Role of Water in Electroporation 98
5.1 Abstract .............................................................................................................................98
5.2 Introduction .......................................................................................................................99
5.3 Water-Vacuum Systems..................................................................................................104
5.3.1 Materials and Methods .........................................................................................105
5.3.2 Results and Discussion ........................................................................................108
5.4 Water Bridges and Osmotic Swelling .............................................................................117
5.4.1 Simulation Methods .............................................................................................118
5.4.2 Experimental Methods .........................................................................................121
5.4.3 Results and Discussion ........................................................................................122
5.5 Conclusions and Summary .............................................................................................141
5.6 Acknowledgements .........................................................................................................143

Chapter 6: The Role of Lipids in Electroporation 144
6.1 Abstract ...........................................................................................................................144
6.2 Introduction .....................................................................................................................145
6.3 Lipid Temperature affects Pore Annihilation .................................................................148
6.3.1 Methods................................................................................................................148
6.3.2 Results and Discussion ........................................................................................149
6.4 Modulating Pore Radius with Alternating Electric Fields ..............................................162
6.4.1 Methods................................................................................................................162
6.4.2 Results and Discussion ........................................................................................163
6.5 Conclusions .....................................................................................................................175
6.6 Acknowledgements .........................................................................................................176

Conclusions and Outlook ..........................................................................................................177
Bibliography ...............................................................................................................................181
Appendix A: Analytical Derivations
A.1 A derivation of the Langevin-Debye dipole moment in an electric field ........................199
A.2 A derivation of the diffusional flow of water across an intact membrane .......................202

viii

List of Tables

Table 2.1 Effect of oxidized lipid concentration on electropore formation time in oxPLPC:PLPC
bilayers ...........................................................................................................................................40
Table 2.2 Effect of electric field magnitude on electropore formation time in 50%
oxPLPC:PLPC bilayers ..................................................................................................................41
Table 3.1 POPC pore creation times .............................................................................................62
Table 3.2 POPC pore annihilation times .......................................................................................62
Table 3.3 DOPC (2005) pore creation times .................................................................................64
Table 3.4 DOPC (2009) pore creation times .................................................................................64
Table 3.5 Normalized POPC pore creation times .........................................................................65
Table 4.1 POPC pore creation times without PS and without Ca
2
................................................83
Table 4.2 POPC pore creation times without PS and with 100 Ca
2+
............................................83
Table 4.3 POPC pore creation times with PS and without Ca
2+
...................................................84
Table 4.4 POPC pore creation times with PS and with 100 ..........................................................84
Table 4.5 POPC pore annihilation times without PS and without Ca
2+
........................................85
Table 4.6 POPC pore annihilation times without PS and with 100 Ca
2+
......................................85
Table 4.7 POPC pore annihilation times with PS and without Ca
2+
.............................................86
Table 4.8 POPC pore annihilation times with PS and with 100 Ca
2+
...........................................86

ix

List of Figures

Figure 1.1 Diagram of a typical zwitterionic phospholipid ............................................................2
Figure 1.2 Electrical properties of a phospholipid bilayer ..............................................................4
Figure 1.3 Electrochemotherapy utilizes membrane electropermeabilization and increases the
efficacy of traditional chemotherapy by enhancing the uptake of chemotoxic drugs into normally
impermeant tumor cells....................................................................................................................7
Figure 2.1 Oxidized and unoxidized phospholipid conformations change over time ..................39
Figure 2.2 Poration of an oxidized phospholipid bilayer ..............................................................43
Figure 2.3 Quilted PLPC:oxPLPC bilayer ....................................................................................43
Figure 2.4 Electropermeabilization enhancement by treatment with peroxidation agents ...........45
Figure 2.5 Enhanced electropermeabilization of peroxidized cells with both ultra-short and
conventional electric pulse treatments ...........................................................................................47
Figure 3.1 Life cycle of an electropore .........................................................................................56
Figure 3.2 Phosphorus and water distribution during pore annihilation .......................................58
Figure 3.3 Electric field profile for a POPC bilayer .....................................................................60
Figure 3.4 Log (pore creation time) versus membrane internal electric field for POPC, DOPC
(2005), and DOPC (2009) ..............................................................................................................66
Figure 4.1 Pore creation times for three different porating electric fields (400, 500, 600 MV/m)
for bilayers consisting of 128 POPC (0PS:0Ca), 128 POPC saturated with calcium (0PS:100Ca),
108 POPC and 20 POPS on a single leaflet without calcium present (20PS:0Ca), and bilayers
containing both PS and Ca
2+
(20PS:100Ca) ..................................................................................81
Figure 4.2 Pore annihilation times after pore formation in the same systems shown in Fig. 4.1. 82
Figure 4.3 Evolution of pore radius after removal of the porating electric field for the same
systems shown in Figure 4.1 ..........................................................................................................89
Figure 4.4 Calcium binding curves show rough correspondence between experimental and
simulated results.............................................................................................................................92
Figure 4.5 Lipid electropore creation time as a function of bilayer internal electric field for the
systems shown in Figure 4.1 ..........................................................................................................95
x

Figure 5.1 Formation of a water bridge across a phospholipid bilayer .......................................109
Figure 5.2 Formation of a water bridge across a vacuum gap separating two regions of liquid
water .............................................................................................................................................110
Figure 5.3 Pore initiation times for water-vacuum and water-lipid systems ..............................111
Figure 5.4 Classification of bulk water and interface water in a typical simulation...................112
Figure 5.5 Anti-correlation of protrusion height and total interaction energy ............................115
Figure 5.6 Comparison of constituent terms of the total interaction energy per protrusion
molecule between WLW (red curve) and WVW (blue curve) simulations .................................117
Figure 5.7 A one-dimensional charge density profile after an external electric field is applied for
a water-lipid and water-vacuum system ......................................................................................124
Figure 5.8 Measurements of interfacial water orientation for water-lipid and water-vacuum
systems after an external electric field was applied .....................................................................125
Figure 5.9 Three-dimensional plots of the mean water dipole moment (z-component) per voxel
during water bridge formation under an externally applied electric field ...................................126
Figure 5.10 Electropore creation sequence .................................................................................128
Figure 5.11 Water bridge precursor (pre-pore pedestal) in a POPC bilayer ...............................129
Figure 5.12 Random permeation event (I) in a phospholipid bilayer .........................................131
Figure 5.13 Random permeation event (II) .................................................................................132
Figure 5.14 Detail of initiation of permeation event shown in Figure 5.13 ................................133
Figure 5.15 Osmotic swelling after nano-electropulse exposure ................................................136
Figure 5.16 A simple exponential function describes the time course of nano-electropulse-
induced cell swelling....................................................................................................................140
Figure 6.1 Area per lipid and lipid temperature in a POPC bilayer simulation as a function of
time while temperature is continuously varied ............................................................................153
Figure 6.2 Area per lipid versus temperature of a POPC bilayer when the system is allowed to
equilibrate for 40 ns at a target temperature ................................................................................154
Figure 6.3 Order parameter S
z
for the sn-1 palmitoyl lipid tail at various temperatures ............155
Figure 6.4 Order parameter S
z
for the sn-2 oleoyl lipid tail at various temperatures .................155
Figure 6.5 Heat capacity at constant pressure was measured by observing the change in enthalpy
per unit temperature .....................................................................................................................156
Figure 6.6 POPC lateral diffusion decreases as temperature is lowered ....................................156
xi

Figure 6.7 Total and constituent pore creation times for sub-zero lipid temperatures ...............158
Figure 6.8 Total and constituent pore creation times for lipid temperatures which vary from 0 to
50 degrees Celsius........................................................................................................................158
Figure 6.9 At high temperatures pores reseal quickly, but at low temperatures pores appear to be
stable out to 200 ns ......................................................................................................................160
Figure 6.10 ‘Cold’ electropores take longer to minimize than warm ones, and they also contain
more waters in their interior .........................................................................................................161
Figure 6.11 The creation of ‘stabilized electropores’ in MD ......................................................164
Figure 6.12 A typical time-dependent bipolar electric field implemented in GROMACS ........166
Figure 6.13 The net water-dipole response of a POPC bilayer simulation during application of a
time-dependent electric field drawn in Figure 6.12 .....................................................................166
Figure 6.14 A typical time-dependent unipolar electric pulse implemented in GROMACS .....167
Figure 6.15 The net water-dipole response of a POPC bilayer simulation during application of a
time-dependent electric field drawn in Figure 6.14 .....................................................................168
Figure 6.16 Measurements of the MD time-dependent electropore radius in a POPC membrane
as a function of pulse amplitude and pulse shape (bipolar (A), unipolar (B)) .............................169
Figure 6.17 A simple DC circuit with a diode attached which results in output signals similar to
the pore radii observed in Figure 6.16 .........................................................................................170
Figure 6.18 A side view of a ‘breathing’ electropore as a time-dependent electric field is applied
......................................................................................................................................................171
Figure 6.19 A top-down view of a ‘breathing’ electropore through the lipid bilayer as a time-
dependent electric field is applied ................................................................................................172
Figure 6.20 A measure of the simulated electropore radius in a POPC bilayer as a function of
time while an external, time-dependent electric field is applied .................................................173
Figure 6.21 A measure of the rate of simulated electropore expansion in a POPC bilayer as a
function of time while an external, time-dependent electric field is applied ...............................174

xii

Abstract
Computer simulations of physical, chemical, and biological systems have improved
tremendously over the past five decades. From simple studies of liquid argon in the 1960s to
fully atomistic simulations of entire viruses in the past few years, recent advances in high-
performance computing have continuously enabled simulations to bridge the gap between
scientific theory and experiment. Molecular dynamics simulations in particular have allowed for
the direct observation of spatial and temporal events which are at present inaccessible to
experiments. For this dissertation I employ all-atom molecular dynamics simulations to study the
transient, electric field-induced poration (or electroporation) of phospholipid bilayers at MV/m
electric fields.
Phospholipid bilayers are the dominant constituents of cell membranes and act as both a
barrier and gatekeeper to the cell interior. This makes their structural integrity and susceptibility
to external perturbations an important topic for study, especially as the density of
electromagnetic radiation in our environment is increasing steadily. The primary goal of this
dissertation is to understand the specific physical and biological mechanisms which facilitate
electroporation, and to connect our simulated observations to experiments with live cells and to
continuum models which seek to describe the underlying biological processes of electroporation.
In Chapter 1 I begin with a brief introduction to phospholipids and phospholipid bilayers,
followed by an extensive overview of electroporation and atomistic molecular dynamics
simulations. The following chapters will then focus on peer-reviewed and published work we
performed, or on existing projects which are currently being prepared for submission. Chapter 2
xiii

looks at how external electric fields affect both oxidized and unoxidized lipid bilayers as a
function of oxidation concentration and oxidized lipid type. Oxidative damage to cell membranes
represents a physiologically relevant system where lipids can become damaged or severely
impacted from interacting with reactive oxygen species, and these events become more frequent
with age. The results are then compared to experiments where we show agreement between our
simulations, theoretical models, and experiments with peroxidized cells in our lab. In Chapter 3 I
outline a set of unique metrics which can be used to quantitatively measure the life cycle of a
discrete electropore for the first time, across multiple lipid species, and I compare these results to
analytical models where we find good agreement with theory. In Chapter 4 I use the life cycle of
an electropore as a tool to measure the effects of electrolyte and lipid headgroup charge on
electroporation compared to electrolyte-free and zwitterionic systems, in addition to presenting
ion binding isotherms to determine the validity of our simulated electrolyte models. Chapters 5
and 6 focus on the roles of water and lipid respectively on electroporation using simplified
water:vacuum systems, osmotic swelling simulations, systems at varying temperature, and
systems where we successfully modulated the electropore radius using customized time-
dependent electric fields. I conclude this dissertation with a brief summary of these studies
followed by a short outlook on the future of electroporation simulations as a whole.

1

Chapter 1: Introduction

1.1 Phospholipids and Phospholipid Bilayers

In cell biology, the macromolecules of life can be neatly subdivided into four broad
molecular classes - proteins, carbohydrates, nucleic acids, and lipids [131]. For the most part,
proteins act as the catalyzers for the chemical reactions required for life, carbohydrates store
energy for later utilization, nucleic acids store genetic information for overall cell function and
division, and lipids form protective structural barriers inside and around the perimeter of a cell to
maintain organelle compartmentalization and chemical regulation. Although there exists a large
number of studies which analyze the function or conformation of proteins and nucleic acids, this
dissertation will focus instead on the roles of lipids and their interactions with external electric
fields.
There are hundreds of structurally distinct lipid species that are relevant in biology such
as triglycerides, sphingolipids, and saccharolipids to name a few, but cell membranes consist
primarily of cholesterol and a group of lipids called ‘phospholipids’, molecules which are
characterized by the presence of a phosphate headgroup and a glycerol backbone that joins two
fatty acid tails together (Figure 1.1). Furthermore the phosphate in the lipid headgroup region is
joined together with an additional organic molecule which can create positive, negative, or
zwitterionic (dipolar) headgroup charge densities [110]. In contrast to the polar lipid headgroup
region, lipid tails are composed of cis- or trans- oriented hydrocarbon chains which are
electrostatically neutral and characterized primarily by their length or degree of saturation (where
single bonds join saturated atoms and double bonds join unsaturated atoms).
2

Lipids are also ampiphilic molecules which mean that their uncharged tails are
hydrophobic (or ‘water-hating’) while their polar headgroups are hydrophilic (or ‘water-loving’).
This dichotomous property results in unique aggregate structures when multiple lipids are mixed
together in an aqueous or cytosolic solution. These structures can take on a variety of forms
based on the solvent temperature or the average lipid size. Examples include micelles which are
spherical lipid aggregate structures with tails facing radially inward and the headgroups facing
radially outward; bilayers which are two-dimensional lateral sheets composed of opposing lipids
with tails enclosed; and vesicles (or liposomes) which are aggregate structures formed from a
bilayer which curves around on itself and joins together to create a hollow spherical shell.
Vesicles or micelles can be used to transport foreign molecules within their hydrophobic or

Figure 1.1 Diagram of a typical zwitterionic phospholipid.
3

hydrophilic interiors, however lipid bilayers are the geometric foundations of cellular
membranes, so I will focus primarily on these structures for the remainder of this dissertation.
As a continuous biological unit, the cell membrane is an integral and dynamic part of the
cell that is quickly able to adapt to minute changes in extracellular and intracellular activity.
Through a complex network of membrane-embedded proteins, water channels, electrolyte
interactions, mechanically- and electrically- sensitive ion channels, cytoskeletal scaffolds, and
other secondary structures, membranes are able to facilitate a wide variety of complex tasks.
These tasks range from the regulation of osmotic pressures to the mediation and externalization
of intracellular molecules for apoptosis signaling [195]. Although lipid bilayers can be thought of
as a ‘passive glue’ which link secondary membrane structures together, using them as a starting
point to model real biological membranes offers a simple and intuitive basis for researchers to
use which they can iteratively build new models from.
The mechanical and electrical properties of simple lipid bilayers closely resemble the
properties of real cell membranes. Lipid bilayers typically span 4-5 nm in thickness depending
on the lipid type, and they hold a net resting potential of roughly 70 mV [153]. This is due to a
number of factors such as surface potentials which stem from differences between bulk
electrolyte charge and the membrane surface charge, dipolar potentials which inhabit the lipid
headgroup region, and other various trans-membrane potentials based on differences in
electrolyte concentration from one side of the bilayer to the other (Figure 1.2). This implies that a
single charge would need only 0.07 eV of energy to cross the bilayer, however this is not the
case. The lipid headgroup region contains a very large 280 mV potential difference across it,
although that potential is hidden because of the equal and opposite headgroup dipoles located on
opposing leaflets. A single charge would actually require about 0.3 eV of energy just to cross the
4

headgroup region (consider a charge 1 e traveling 1 nm against a ~300 MV/m electric field).
However, a single charge needs only to traverse halfway across the bilayer because the
hydrophobic lipid interior will push the particle outward with little to no additional energy. For
symmetric bilayers one can imagine the headgroup polarization spanning at most 2.5 nm (half-
way across the bilayer) in which case the translocation energy approaches 1 eV, though this can
vary quite a lot if there are asymmetries present.

Approximate Translocation Energy = 𝑞𝐸 𝑑 = (1 e) (300 MV/m) (2.5 nm) = 0.75 eV (1.1)

Figure 1.2 Electrical properties of a phospholipid bilayer [31, 153]. ∆Ψ
m
, ∆Ψ
s
,

and

∆Ψ
d
stand
for the transmembrane, surface, and membrane dipole potentials respectively whereas Ε
m
, Ε
s
,

and

Ε
d

stand for the corresponding electric field values.
5

1.2 Electroporation
1.2.1 Overview
One way we can study cells and cell membranes is to observe their behavior under a
variety of physiological and non-physiological conditions to understand how they function or
interact with the outside world. In medicine, one of the central problems researchers have is
finding an efficient way to target and deliver therapeutic molecules to specific cells, and
furthermore how to reliably get molecules across the cellular membrane and into the cell interior.
A large number of potentially viable drugs become underutilized because it becomes difficult to
target precise biological locations without costly treatment plans that are often difficult or
impossible to implement. Biology has, after all, evolved for millions of years to provide natural
defenses from external tampering, however life itself has created its own mechanisms for
transporting material into cells through viruses, peptides, ion channels, and proteins, but these
delicate biological processes are not easily utilized and are still not well understood. As a result
we turn to the interactions of biological systems with physical systems where we are free to
apply non-physiological forces such as external electric fields or external surface tensions with
the intent of reliably manipulating the cell membrane.
Although the structural integrity of the cell membrane is robust and difficult to
disassemble, it can be disrupted under the influence of large MV/m external electric fields
through a process called electroporation [175]. Electroporation (sometimes referred to as
‘electropermeabilization’) was developed in the late 1970’s and early 1980’s as a technique to
introduce external DNA into mouse lyoma cells [119, 204] by applying quick, high-voltage
pulses (8 kV/cm, 5 microsecond) once every three seconds to cell suspensions where DNA had
6

been added and allowed to adsorb onto the cell surface. This facilitated the entry of linear and
circular plasmid DNA from the cell surface into the cell interior which subsequently encouraged
expression of normally absent DNA in electropermeabilized cells. Although pulsed electric fields
had been used previously to stimulate the release of ATP, catecholamine, and proteins from
chromaffin granules of bovine adrenal medullae [118], the use of electric fields to transfer
biologically-active agents into the cell interior is a relatively new technique that could
significantly affect the development of novel therapeutic drugs.
When a cell becomes permeabilized there is a measurable increase in its electrical
conductance along with an influx and efflux of normally impermeant small and large molecules
[108, 119, 142] across its membrane. Small changes in conductance can be measured within a
few nanoseconds after the application of a pulsed electric field [11], and fluorescent dyes can be
used as molecular markers to indicate the degree of membrane permeabilization, but the detailed
mechanisms of electropermeabilization are difficult to access experimentally [173], even with
nanoscale instrumentation such as atomic force microscopy. Other analytical methods, including
continuum and atomic-level molecular simulations, provide windows for investigating
electropore formation at the nanoscale, where there are at present no direct observational tools.
Electroporation is also used extensively in the lab for introducing polar molecules such as
dyes into cells for fluorescence tagging, or in clinical settings where it can aid in the treatment of
tumors by facilitating the transfer of chemotherapy drugs into normally resistant tumors using
electrochemotherapy [2, 5, 55, 140] (Figure 1.3). It is hypothesized that aqueous pores lined with
phospholipid headgroups (so-called ‘hydrophilic pores’) are created in the membrane interior
[50, 106, 200, 203]. Although several analytical models characterize aspects of electroporation
7

such as the size and distributions of pores, the discrete atomic processes of electric field-driven
pore creation and annihilation remain poorly understood.

1.2.2 Continuum Models

Early models of membrane permeabilization in the 1950’s were based on experiments
with cells where the conductance and dielectric properties of the membrane varied as a function
of the external electric field frequency [32, 155]. Around the same time researchers also
discovered a change in membrane impedance could be triggered in artificial membranes using
strong currents generated from ion gradients across the membrane [53]. It was not until
Neumann and Rosenheck demonstrated in the early 1970’s [118] that permeability, and thus
conductance could be a transient event if the electric field was removed quickly (or pulsed with a
short enough period). This type of electroporation was termed ‘reversible electroporation’
Figure 1.3 Electrochemotherapy utilizes membrane electropermeabilization and increases
the efficacy of traditional chemotherapies by enhancing the uptake of chemotoxic drugs
into normally impermeant tumor cells. Original image courtesy of Wikipedia.
8

because its effects were temporary if the external field was removed in a sufficient amount of
time. This is in contrast to ’irreversible electroporation’ which describes instead systems that are
permanently modified from prolonged exposure to electric fields, sometimes leading to eventual
cell death if the field is maintained.
Irreversible electrical breakdown of membranes was first described by Hamilton and Sale
[60, 147, 148] in a set of three publications which described the non-thermal killing of yeast,
bacteria, and erythrocytes under kV/cm electric fields. More specifically membrane lysis was
observed at transmembrane voltages of 0.3 volts to 1.15 volts. Reversible electrical breakdown
and permeabilization of cell membranes however appeared to be more closely correlated with the
initial studies on changes in membrane conductance [1, 11, 27, 81] in which the membrane was
modeled as a viscoelastic medium. Additionally, because reversible electroporation studies were
transient and generally non-lethal to cells, it became an attractive system for further study and
subsequently laid the groundwork for modern continuum models to build upon.
Observing electroporation directly was not (and to the author’s knowledge, still not)
possible experimentally, thus a plethora of analytical models [36, 37, 57, 165, 200] and indirect
experimental observations [47, 65, 112, 119, 138] were combined to try and explain the
molecular mechanisms of electroporation using continuum electrodynamic equations. The
simplest analytical models treat cell membranes or lipid bilayers as a smooth hydrophobic
surface and apply Poisson's equation to determine the potential at the surface. Double layer
models (which include the effects of adsorbed surface charges built up at the water/lipid
interface) such as the Gouy-Chapman-Stern model [164] use more sophisticated methods, but
lipid headgroups in these models are coarsely-discretized and have their own modified charges
which greatly complicates the overall electrostatic environment.
9

Models of electroporation fall into two broad categories which are deemed ‘surface
potential models’, or models which calculate a static potential profile for a bilayer and combine it
with stability conditions (e.g. surface tension bounds), and ’stochastic pore models’ which
assume the existence of a distribution of transient pores with some associated probability.
Chizmadzhev and colleagues pioneered models of the first type [1, 29] which use circuit analogs
to analyze the electrical permeation of lipid bilayers (e.g. Pore Energy = 𝑒 𝑑𝑔𝑒 𝑒𝑛 𝑒𝑟𝑔𝑦 −
𝑠 𝑢 𝑟 𝑓𝑎 𝑐𝑒 𝑒𝑛 𝑒𝑟𝑔 𝑦 − 𝑓 ( 𝑚𝑒 𝑚𝑏 𝑟 𝑎 𝑛 𝑒 𝑐𝑎𝑝 𝑎𝑐𝑖𝑡 𝑎𝑛 𝑐𝑒 , 𝑟 𝑒 𝑠 𝑖 𝑠 𝑡 𝑎 𝑛𝑐𝑒 , 𝑒 𝑡 𝑐 )) while Sugar and Neumann
developed models of the second type [166, 167] using bilayer transition probabilities. Modern
models use hybrids of these such as the asymptotic model of electroporation [117] which include
both the electrical properties of the membrane and probabilistic estimates for the number of open
pores N per unit time (Equation 1.2) where α and β are normalization constants and N
eq
is the
equilibrium number of pores. Monte Carlo models can also be used to study lipid bilayers, but
because electroporation is a transient non-equilibrium event, molecular dynamics (MD)
simulations are much more suitable for this endeavor.

𝑑𝑁 𝑑𝑡 = 𝛼 𝑒 𝛽 𝑉 2
(1 −
𝑁 𝑁 𝑒𝑞
) (1.2)

Much of the interest in (historically microsecond) electroporation stems from its ability to
reliably transport molecules across a cellular membrane without permanently disrupting the
function or structure of the cell, however recent observations may indicate that short, nanosecond
MV/m electric pulses induce different types of cellular responses. These studies have shown
[192, 195, 202] that ‘ultrashort’ nanosecond MV/m electric fields can trigger intracellular
calcium bursts and the externalization of phosphatidylserine, a lipid headgroup which can help
10

signal apoptosis, from the plasma membrane’s inner leaflet to the outer leaflet. This suggests that
the frequency of the external electric field is an important factor to study, in addition to how
frequency affects other physical properties such as the radius of a discrete electropore which can
affect how readily ions, dyes, and other cytosolic molecules travel through the membrane [191].
Later in this dissertation I will show using numerical simulations what the effects of modulating
the pulse frequency are on lipid bilayer electroporation and relate extracted parameters from our
observations in molecular dynamics to those in continuum models for direct comparison.
The results obtained through ultra-short nanosecond electroporation differ considerably
from traditional microsecond electroporation, and much work remains to be done in presenting a
coherent, atomic-level description of membrane electropermeabilization in cells which
encompass major experimental observations. In order to efficiently study these systems in
atomic-level detail over nanosecond timescales we can use computer simulations of lipid bilayer
electroporation to explain experimental results using ab-initio physical and chemical models.
Our goal then is to take continuum electrodynamic models and experimental observations and
combine them through all-atom molecular dynamics simulations in order to better explain the
molecular mechanisms of electroporation while also providing simple, testable predictions for
future experiments.

11

1.3 Molecular Dynamics Simulations

1.3.1 Overview

The first molecular dynamics (MD) simulations of a lipid bilayer system were performed
in 1982 by Herman J.C. Berendsen using two layers of sixteen decane molecules held together
by a simple dihedral potential [181] at a rate of 80 picoseconds of simulation time per 10 hours
of computation time. From these simulations the authors were able to correctly extract the
carbon-deuterium order parameter S
cd
relative to the bilayer normal using a second-order
Legendre polynomial as a geometric measurement of hydrocarbon orientation and rigidity.
Furthermore Seelig confirmed through similar simulations [157, 158] that lipid tail rigidity
decreases as one traverses down the hydrocarbon chain, a result which was consistent with
experimental observations. These studies had in fact broken new ground proving that it was
possible to use atomic simulations as a means of explaining complex biological observations
using simple, classical physics models as a basis, however because of the small system
geometries and short timescales available for simulation back then, only a limited number of
measurements were possible.
As computers became faster and more robust, the timescales and system sizes which
could be simulated grew as well which enabled larger and more complex studies to take place.
Eventually user-defined external potentials were replaced with all-atom solvent models of water,
and advances in long-range electrostatic algorithms allowed for new force fields and integration
techniques to be implemented. Eventually the accuracy and precision of simulations became so
high that direct experimental observations could not validate them, and so a paradigm shift
occurred where there was a demand for higher-resolution experiments to confirm the predictions
12

of simulations, rather than the other way around. Simulations now, for instance, can show in
atomic detail the ion transport and current flow of electrolyte passing through a permeabilized
membrane, in addition to extracting the energy barriers required for lipid flip-flop [151, 177],
however these quantities remain difficult to study experimentally and often must be observed
indirectly. One must also check that models and results obtained in MD are at least somewhat
consistent with experimental observations, a feat which is difficult to accomplish when working
in regimes that are not well documented or studied.
There are also interactions that cannot be directly resolved with classical mechanical
models. An example is the quantum behavior of hydrogen since electrons can sometimes
quantum tunnel through large potential barriers. Another example is when liquid helium atoms
turn superfluidic due to BCS theory. The quantum threshold is generally located near frequencies
ν where the energy hν is greater than kT (where h is Planck’s constant, k is Boltzmann’s
constant, and T is temperature) and this is the case for most hydrogen and helium bonds. In this
case we can either add static quantum corrections to each classical harmonic bond, or we can use
united-atom models where atom-to-hydrogen links are treated as single unified atoms. In the
united case, existing bonds are rigidly constrained at a fixed distance which describes quite well
the quantum oscillator with electrons in the ground state. Similarly MD models lack explicit
electron polarizability which makes simulations of dielectric mediums difficult, although some
have introduced explicit Drude-model corrections to compensate for this [89]. Computing
limitations also limit the total number of atoms that can be modeled in a single system and
require interaction cutoff distances to be defined.
Despite these limitations, advances in high performance computing over the past few
decades have contributed much to the study of electroporation, and molecular dynamics
13

simulations of phospholipid bilayers in porating electric fields have revealed striking details
about pore creation such as the importance of interface water orientation [171, 175, 209]. MD
has also confirmed that electroporation can facilitate various physiological functions such as the
externalization of phosphatidylserine (PS) [195], a process which is important for cell signaling
and which was first observed in experiments [189], further validating simulation models. Some
of our preliminary MD simulations have also suggested that pore creation may be related to the
electroevaporation of water molecules in a vacuum [122], a process where highly ordered
structures or ‘pedestals’ of water form at a low-dielectric interface in order to minimize their
total potential energies. These pedestals serve as launchpads for surrounding waters to free
themselves from the hydrogen bonds present in bulk solution, but few systematic studies exist
which thoroughly compare electroevaporation to lipid bilayer electroporation. Many simulations
of electroporation are cursory and are often incomplete with only a single trial of data presented,
thus there is an explicit need for additional studies to thoroughly and systematically characterize
the atomic mechanisms of electroporation. This would allow existing analytical models to be
validated by large, statistical data sets with fine granularity, further bridging the gap between
theory and experiment.
The focus of this dissertation is to systematically characterize and analyze the structures
which give rise to lipid bilayer electroporation using measurements of pore formation times, lipid
tail order parameters, binding affinities, radial distribution functions, and other kinetic properties
associated with the restructuring of lipid bilayers in electric fields. We accomplish this by using
the open-source molecular dynamics package GROMACS (GROningen MAchine for Chemical
Simulations) [64] which is a fast, scalable, and powerful leap-frog integrator developed at the
University of Groningen for iteratively solving Newton’s equations of motion. Other molecular
14

dynamics packages exist such as NAMD [115] and AMBER [129], however GROMACS is best-
optimized for parallel computation using the Message Passing Interface (MPI) scheme with
efficient electrostatics algorithms to attain a computational complexity of O(NLog(N)) through
the use of binary tree summation (where N is the total number of simulated atoms). The
GROMACS integrator produces trajectories which are identical to the Velocity-Verlet algorithm,
first taking all positions and potentials into account, then calculating the resulting force on each
atom

𝐹 𝑖 = −
𝜕𝑉
𝜕 𝑟 𝑖 (1.3)

which determines all subsequent atom accelerations and velocities.

𝜕 2
𝑟 𝑖 𝜕 𝑡 𝑡 =
𝐹 𝑖 𝑚 𝑖
𝜕 𝑟 𝑖 𝜕𝑡
= 𝑣 𝑖 (1.4)

Once new positions are calculated, updated potentials are generated and the algorithm repeats
itself. The only required input quantities are initial atom positions, atom topologies or force
fields which label static properties of the system (e.g. masses, Lennard-Jones coefficients,
charges, etc), and a configuration file which specifies the ensemble to be used, the thermostat,
the barostat, the field to apply (if any), the integration time step, and the duration of the
simulation. GROMACS is distributed under the GNU General Public License (GPL) which
makes it free to operate and modify, along with a set of built-in analysis tools included with the
source code.
15

1.3.2 Implementation

The potentials represented in GROMACS are Coulomb, Lennard-Jones, bond length,
bond angle, and bond dihedral/torsion, having the functional form:

𝑈 = ∑
𝑞 𝑖 𝑞 𝑗 4 𝜋 𝜀 0
𝑟 𝑖𝑗
+
𝑖 < 𝑗 ∑ 4 𝜀 𝑖 𝑗 �
𝜎 𝑖𝑗
1 2
𝑟 𝑖𝑗
1 2
+
𝜎 𝑖𝑗
6
𝑟 𝑖𝑗
6
� +
𝑖 < 𝑗 ∑
1
2
𝑘 𝑖𝑗 𝑏 � 𝑟 𝑖𝑗 − 𝑏 𝑖𝑗 0
�
𝑏 𝑜 𝑛𝑑𝑠
(1.5)
+ �
1
2
𝑘 𝑖𝑗𝑘
𝜃 ( 𝜃 𝑖𝑗𝑘
− 𝜃 𝑖𝑗𝑘
0
)
2
+
𝑎𝑛𝑔𝑙𝑒 𝑠 � 𝑘 𝜑 [1 + cos� 𝑛 ( 𝜑 − 𝜑 0
) �
𝑑𝑖 ℎ𝑒 𝑑𝑟 𝑎𝑙𝑠

Cutoff (or truncation) distances are implemented around each atom (approximately 1 nm) for
Coulomb and Lennard-Jones interactions, however electrostatics calculations past the cutoff
radius are still considered using the Particle Mesh Ewald (PME) method [41]. While short-range
electrostatics quickly converge in real space, PME performs Fast Fourier Transforms (FFTs) on
long-range atoms which allow those contributions to converge quickly in reciprocal Fourier
space. Although spring constant analogs can be assigned to bond lengths and angles, one finds
that the classical ground state frequencies of many simple bonds is in the Terahertz range which
implies that the timescales of bond oscillations are considerably longer than our integration
timestep of 2 fs. In response, we (and many researchers) have chosen to use rigid constraints
where bonds are held at fixed distances, similar to maintaining bonds in their quantum ground
state. This also allows for the use of larger integration time steps which enables us to simulate
over longer, physiological time scales.
Most of our lipid bilayer systems contain about 128 lipids consisting of (for example) 1-
palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine (POPC) or 1-palmitoyl-2-oleoyl-sn-
16

glycero-3-phosphatidylserine (POPS) — and about 4480 to 9000 water molecules (35 to 70
waters per lipid). Typical box sizes are approximately 7 nm x 7 nm x 7 nm or 7 nm x 7 nm x 10
nm in x, y, and z respectively. Simulations with multiple trials are run in parallel, starting from
the same initial positions. To ensure that each trial is independent, every atom is assigned a
randomized velocity based on a Maxwell-Boltzmann distribution at the beginning of the
simulation with a mean velocity such that
1
2
𝑚 𝑣 2
= 𝑘 𝑏 𝑇 . All systems are equilibrated to constant
area per lipid. Periodic boundary conditions are employed to mitigate system size effects. All
simulations are run on the USC High Performance Computing and Communications Linux
cluster (http://www.usc.edu/hpcc/) with an average of 72 cores per simulation at a rate of about
50 ns of simulation time per day.

1.3.3 Water Models

The interactions between biological molecules and aqueous solvents are absolutely
crucial for life to exist, so one must pay careful attention to ensure that the water models being
implemented are realistic and experimentally consistent. Failure to do so can dramatically
change the intended function or conformation of each molecule we introduce into an aqueous
system, and invalid water models could quickly negate decades of research in force field
optimization by providing misfolded lipid and protein structures during hydration. Many
different explicit and implicit water models exist with some models utilizing three atoms while
others utilize up to five atoms, however each model is optimized to exhibit specific
characteristics that may be desirable during different types of simulations. This implies that there
is no single water model which is universally better than others, only water models which best
deliver the characteristics one is trying to model.
17

The most common water model is an explicit model called the Simple Point Charge
(SPC) model [12] which treats each water atom as a classical point charge. The partial charge
assignment for SPC water is 0.41 e for hydrogen and -0.82 e for oxygen (where e is the electron
charge) however this results in an incorrect net dipole moment, so instead of using the
experimentally verified H-O-H angle of 104.45º, SPC water utilizes an ideal tetrahedral shape
with an H-O-H angle of 109.47º, resulting in a dipole moment of 2.23 Debye. Because SPC
water is so simplistic it also tends to diffuse somewhat quicker at low temperatures compared to
real water because SPC’s solid-phase structure is affected by the lack of electron polarization.
This effect however is small and the computational speedup gained from utilizing SPC water
tends to outweigh the small perturbations at low temperatures, making SPC water a top candidate
for most general simulations compared to other more complicated models. More recent studies
have updated the SPC model with partial charge assignments to four significant digits rather than
two, and they also include additional polarization correction terms to the potential energy [13].
This extended SPC model (termed SPC/E) results in more realistic diffusion behavior at low
temperatures which creates ice densities that are more accurate than the original model, however
the additional correction term to the potential energy function makes it more costly to implement
than SPC water while the improvements at human body temperature (310 K) are minimal.
Other three-site water models exist such as the TIP3P model [79] which use the
experimental H-O-H angle of 104.45 º and the more realistic OH bond length of 0.95 Å, however
due to lack of polarization this model also requires costly corrections that add up when hundreds
of thousands of water molecules are simulated at once. Other models try and work around this
using explicit correction terms by inserting virtual dummy atoms in tetrahedral geometries to
mimic non-local charge densities such as the TIP4P model which is a four-site model and TIP5P
18

which is a five-site model [98]. These models are capable of obtaining water densities that differ
from experiment by only 6 mg cm
-3
at various physiological temperatures and pressures,
however the modifications can sometimes lead to unwanted interactions or conformations
between dummy atoms on the four- and five-site water models with other biological molecules
present in the system.
For the purpose of studying water:phospholipid interfaces (e.g. hydrated phospholipid
bilayers) and their relationship to electroporation, SPC water provides the most experimentally
consistent results while also being the most computationally efficient model. The efficiency of
this model is realized through the use of SETTLE constraints (an analytical version of the
SHAKE algorithm) which quickly constrains all rigid OH bonds in two simple steps, thus
increasing performance by factors of 3 or 4. More information about these techniques can be
found in the ‘Holonomic Constraints’ section.

1.3.4 Phospholipid Models

Once we are satisfied with the accuracy of our solvent force fields we need next to utilize
valid molecular models of the smallest biological structures such as lipids or proteins in order to
begin simulating the macro biological processes which occur in aqueous environments.
Simulations on the order and complexity of a phospholipid bilayer must be built upon
experimental observations whenever they become available to ensure that the simulation results
are meaningful and grounded in reality. These observations often include x-ray diffraction
measurements to obtain the bilayer area per lipid; nuclear magnetic resonance (NMR) on
deuterated lipids to observe order parameters or relaxation times for tail dihedrals; black lipid
19

membrane observations to determine permeability; IR spectroscopy for trans-gauche defects; and
fluorescence microscopy, ESR spectroscopy, measurements of membrane action potentials,
electrophoretic mobility, atomic force microscopy, and differential scanning calorimetry for
determining other membrane properties such as the phase transition temperature. In the first lipid
simulations there were no obvious a priori length or time scales to use and there are countless
initial lipid geometries, electrostatic cutoffs, and boundary conditions to choose from. How does
one begin tackling a problem of such high complexity?
Early lipid models began with a crystalline lipid geometry, however because biological
membranes are not usually in their crystalline phase, artificial defects were added for
stochasticity [186] so long as the models still agreed with experimental order parameters
obtained through NMR. The most common models differed predominantly by their
parameterization methods and their short-range electrostatic cut-offs, however every model
contains a similar potential function which is comprised of a Lennard-Jones contribution, a
Coulomb contribution, a quadratic contribution for bonds and angles, and a cosinusoidal dihedral
series. Some models implement ‘all-atom’ descriptions which explicitly include every atom
including hydrogen, whereas other models use ‘united-atom’ descriptions which map the
properties of one or two small atoms (usually hydrogen) onto a larger, adjacent atom (e.g.
carbon) for simplicity. There also exist ‘course-grained’ models which combine a much larger
set of atoms into a single course atom for faster computation times [101] however these models
often ignore the quick picosecond or nanosecond motion of small atoms and should only be used
when small interactions can be adequately ignored. Although some models have recently begun
to incorporate atomic polarization [89] a majority of the existing force fields completely lack
explicit electron polarization, thus most of the lipid (and water) force fields are incomplete,
20

however they are sufficient for studying first order effects in biological systems, especially when
multiple trials are run in order to collect statistically significant observations.
The major atomic force fields used for lipid bilayer simulations are AMBER, CHARMM,
GROMOS, OPLS, and the Berger model. AMBER, a molecular dynamics package and a
molecular force field (sometimes referred to as the Generalized Amber Force Field (GAFF)
[78]), minimizes the chemical bond stretching energy and uses a united-atom approach for CH
2

atoms in the lipid tail region. Although AMBER has been used in many lipid [73, 130] and lipid-
protein simulations [72] it is slowly becoming deprecated due to a lack of long-term support
from its creators at UC San Francisco. The CHARMM force field on the other hand, originally
implemented at Harvard University in the early 1970s by Martin Karplus and updated every few
years [23], has seen considerable use in lipid bilayer simulations and was originally intended to
be used with TIP3P water and explicit electrostatic cut-off schemes which vary the short-range
cut-off radius for various atomic interactions. Subsequent revisions to CHARMM have added
explicit all-atom hydrogen atoms which make simulations more-closely resemble experiments,
although this is often a computationally expensive option.
The lipid force fields used primarily throughout this dissertation are comprised of the
remaining three model types not yet discussed - GROMOS, OPLS, and Berger. The GROMOS
force field [39] (originally the default force field of GROMACS) is intended for use with SPC
water and utilizes united-atom CH
2
atoms as AMBER does, however the initial GROMOS
parameters resulted in a lipid density which was too high (a property of most united-atom
models). Berger et al. [15] thus re-parameterized DPPC (dipalmitooleoylphosphatidylcholine)
lipids using both GROMOS force fields for bonded parameters and an existing force field called
OPLS (optimized potentials for liquid simulations) for non-bonded interactions, and these new
21

simulations exhibited lipid bilayer densities that were within 1% of experimental observations.
The authors also parameterized Lennard-Jones interactions on the lipid tail by using parameters
obtained from previous simulations of pentadecane, thus making this model one of the most
widely-used and robust force fields for lipid bilayer simulations.

1.3.5 Energy Minimization

GROMACS provides a number of tools and methods to energy minimize a given set of
atomic coordinates. The potential energy landscape of a macromolecular system is generally
complex and difficult to minimize, however there always exists at least one global minimum
which we are after where the derivatives of the potential energy are zero and all second
derivatives (i.e. the eigenvalues of the Hessian Matrix) are positive, although many local
minimas have similar characteristics which make identifying a global minimum challenging. To
accommodate this GROMACS must easily transition from local minimum to local minimum by
traversing through ‘mountain paths’ or saddle points and paths characterized by Hessian matrices
containing only one negative eigenvalue. Because the dimensionality of the entire landscape is
very high, complete knowledge of saddle points and minima is difficult in practice to obtain, and
no one minimization method guarantees convergence, however a small number of techniques
have been implemented historically which deliver reasonable accuracy.
Three general energy minimization techniques are: those that require only functional
evaluations (i.e. simplex methods) where a step is made based on the success of previous
evaluations, those that utilize first derivatives of the potential, and those that utilize both first and
second derivatives of the potential. GROMACS is able to implement all three of these methods
22

using different algorithms. It can evaluate the nearest local minimum using the steepest descent
technique which is of the second type mentioned above, simulated annealing which is of the first
type mentioned above, and conjugate gradient techniques which is of the third type which use
energy gradient information from previous steps to find configuration minimas. In addition to
these, GROMACS also supports Replica-Exchange Molecular Dynamics (REMD) where
multiple replicas of a system are run at different temperatures. In REMD two systems will
periodically swap their temperatures if the switch further minimizes their free energy, or if a
random number chooses to do so, in order to adequately sample all possible confirmations. This
technique is somewhat analogous to Monte Carlo sampling and other Metropolis algorithms
which perturb some property of the system based on both a random number and the evaluation of
some system characteristic which must be optimized.

1.3.6 Electrostatics

Electrostatic pair calculations account for almost all of the computation time in molecular
dynamics simulations, so the algorithms used to calculate these interactions needs be quick and
efficient since its implementation can dramatically slow down or speed up simulations by
multiple orders of magnitude. Because the functional form of the Coulomb potential scales as a
function of inverse distance, close forces result in large electrostatic contributions and forces
which are far result in small electrostatic contributions. As a result GROMACS treats short-range
electrostatics differently from long-range electrostatics by how often it calculates slowly-varying
long-range interactions versus quickly-varying short-range interactions. This technique is
implemented through a twin-range cutoff scheme where a cutoff radius determines when short-
range interactions transition to long-range interactions. Short range interactions are classified as
23

all electrostatic interactions that lie within 1 nanometer of the area we are interested in and these
forces are summed every 2 fs. Long-range interactions (i.e. all electrostatic interactions beyond 1
nm) are summed instead every 20 fs (or ten time-steps) in reciprocal Fourier-space using a
Particle Mesh Ewald (PME) method since the inverse sum converges quickly over large
distances.
While PME is in fact a robust and efficient electrostatic summing algorithm for long-
range interactions, it has its roots in crystallography where the structure in question, generally a
crystal, was assumed to be periodic, charge-neutral, and isotropic. When we model membranes
we use periodic boundary conditions to mirror our system in the x, y, and z directions, however
long-range anisotropies can affect the accuracy of PME if a net system charge is present. To
account for any net charges PME will create an equal and opposite charge density which is
evenly distributed around the system in the form of a surface dipole moment so that the resulting
system is net neutral, thus allowing for the creation of mirror images through periodic boundary
conditions without an infinite number of charges being introduced [198]. Although our
simulations often contain counter-ions which balance local charge asymmetries, the application
of a global external electric field results in a net system dipole moment which can introduce error
when interacting with the equal and opposite dipole moment created with PME. To account for
this error and the magnitude of the additional dipole term in PME we can adjust a referential
dielectric constant at the unit cell surface to one (to achieve vacuum boundary conditions), to
infinity (to achieve tin-foil boundary conditions), or to zero (to regain the traditional PME
algorithm).

24

1.3.7 Thermodynamic Ensembles

If an ensemble is not specified in GROMACS then by default it will implement a simple
microcanonical ensemble where the number of particles, the system volume, and the total system
energy is conserved in time (i.e. an NVE ensemble). Real systems however are not well
represented by microcanonical ensembles, so various canonical ensembles can be maintained by
holding volume and temperature constant (an NVT ensemble), pressure and temperature constant
(an NPT ensemble), or pressure and surface-tension constant (an NPγT ensemble) which require
the continual addition of thermodynamic energy to the system. Although NVT ensembles are
used to first create lipid bilayer systems with experimental area per lipid values, most of the
simulations contained in this dissertation will utilize an NPT canonical ensemble after
equilibration with an NVT ensemble to mimic the 1 bar of pressure real membranes experience.
There are many different thermostats and barostats GROMACS can take advantage of to
maintain specific thermodynamic ensembles. One of the most popular thermostats for
maintaining temperature is the Berendsen thermostat [14] named after Herman J.C. Berendsen at
the University of Groningen which rescales atomic velocity vectors over a time constant τ, based
on the deviation of the system temperature from some reference temperature (Equations 1.6 to
1.8)

𝑑𝑇
𝑑𝑡
=
𝑇 0
− 𝑇 ( 𝑡 )
𝜏 (1.6)
𝑣 ′
= 𝜆𝑣 (1.7)
25

𝜆 = [1 +
𝛥 𝑡 𝜏 �
𝑇 0
𝑇 ( 𝑡 −
𝛥 𝑡 2
)
− 1 � ]
1
2
(1.8)
Here T is the current temperature, T
o
is the reference temperature, t is time, τ is a time constant,
𝑣 ’ is the rescaled velocity, 𝑣 the initial velocity, λ the rescaling factor, and ∆t the integration
time step. Although the Berendsen thermostat is widely-used in lipid bilayer simulations, other
researchers such as Shuichi Nosé and William Hoover took alternate approaches by introducing
a unitless interaction variable s which takes into account the interaction of the external bath with
the system rather than continually correcting with a weak coupling algorithm. This alternate
approach became known as the Nosé-Hoover thermostat and it produces a fundamentally
different Hamiltonian from weak-coupling methods, producing instead equations of motions that
contain the interaction variable s and its associated conjugate momentum p
s
, thus allowing the
user to fine-tune the degree of canonical interaction from a simple microcanonical system.
More recently in 2007 Michele Parrinello’s group developed what is known as a velocity-
rescaling algorithm (24) which uses stochastically-determined velocity rescaling values rather
than fixed ones. Unlike Berendsen’s weak coupling algorithm to a fixed external bath or Nosé-
Hoover’s method of using fixed rescaling variables which generate hybrid interaction
coordinates and conjugate momenta, Parrinello’s group coupled atom velocities to a dynamic
factor which is calculated from the instantaneous kinetic energy deviation from the average
kinetic energy as calculated from thermodynamic equipartition (Equation 1.9) at a given
temperature.
𝐾 �
=
𝑁 𝐹 2
𝑘 𝑏 𝑇 (1.9)
26

𝐾 �
is the average kinetic energy at temperature T, 𝑘 𝑏 is Boltzmann’s constant, and 𝑁 𝐹 is the number
of degrees of freedom per atom. Initial systems presented in this dissertation were coupled to a
temperature bath at 310 K (body temperature) with a relaxation time of 0.1 ps using the
Berendsen weak coupling algorithm, although later studies utilized the velocity-rescaling
thermostat which produces more realistic canonical ensembles.
As for pressure there are two dominant barostats in use with GROMACS which are the
Berendsen barostat [14] and the Parrinello-Rahman barostat, both developed by the same
individuals who contributed to the thermostat algorithms. The Berendsen barostat once again
uses a weak-coupling algorithm to slowly introduce pressure into a system from an outside bath
by scaling atom positions by a constant which depends on the pressure differential (Equations
1.10 to 1.12).
𝑑𝑃
𝑑𝑡
=
𝑃 𝑜 − 𝑃 ( 𝑡 )
𝜏 (1.10)
𝑞 ′
= 𝜇𝑞 (1.11)
𝜇 = 𝛿 𝑖 𝑗 −
𝛥𝑡
3 𝜏 𝛽 𝑖𝑗
( 𝑃 𝑜 , 𝑖𝑗
− 𝑃 𝑖𝑗
( 𝑡 )) (1.12)
Here P is the current pressure, P
o
the reference pressure, t is time, τ is a time constant, 𝑞 ’ is the
rescaled spatial coordinate, 𝑞 the initial spatial coordinate, µ the rescaling factor, ∆t the
integration time step, 𝛽 𝑖 𝑗 the isothermal compressibility in the i
th
to j
th
direction, and 𝛿 𝑖 𝑗 the
Kroenecker delta. Alternatively Michele Parrinello and Annesur Rahman (one of the founders of
the molecular dynamics method) produced another barostat in a way similar to how the Nosé-
Hoover thermostat was developed. They introduced a unitless coupling parameter h which
defined the degree of volume distortion (e.g. from a square to a parallelepiped geometry) and a
mass-like inertial variable W which defined the degree of resistance to box deformations [128].
27

This resulted in a new Hamiltonian which produced equations of motion that included the
pressure p and the scaling factors h and W along with their conjugate momentum, and became
known as the Parrinello-Rahman barostat. The simulations in this dissertation however are
coupled semi-isotropically with the Berendsen barostat (using a compressibility of 4.5 x 10
-5

bar
-1
) normal to and in the plane of the membrane with a pressure bath of 1 bar and a relaxation
time of 1 ps.
GROMACS also calculates a Virial which helps measure the instantaneous pressure and
temperature at any given time. If we define the Virial and kinetic energy as:
Ξ = −
1
2
∑ 𝑟 𝑖𝑗 ⨂𝐹 𝑖𝑗 𝑖 < 𝑗 (1.13)
K
�
=
1
2
∑ 𝑚 𝑖 𝑣 𝑖 ⨂𝑣 𝑖 𝑁 𝑖 (1.14)
then we can determine the temperature T using Equation 1.9 or the pressure tensor P using the following
relationship:

1
2
𝑷 𝑉 = K
�
− Ξ (1.15)
where V is volume, r
ij
is the atomic position in the i
th
to j
th
axis, F
ij
is the force in the i
th
to j
th
axis,
and m
i
is the mass of the i
th
atom. If we need to obtain the average scalar pressure for use with
the isotropic pressure coupling algorithms mentioned previously we can simply take the average
of the diagonal eigenvalues of the pressure tensor P to obtain:
𝑃 𝑠𝑐𝑎𝑙𝑎𝑟
=
𝑇𝑟𝑎 𝑐 𝑒 ( 𝑷 )
3
(1.16)

28

1.3.8 Restraints and Holonomic Constraints

In GROMACS, "constraints" and "restraints" refer to fundamentally different simulation
processes. Constraints are algorithms specifically designed to hold bond lengths or angles
constant after the integration of Newton’s equations whereas restraint algorithms modify the
Hamiltonian in such a way so that simulations line-up with some preconceived theory.
Restraints are often created by delivering some sort of energy penalty based on the degree of
divergence from a predefined value (e.g. coordinates, atom spacing, bond dihedrals). An
example of a restraint is the ability to artificially freeze a set of atoms such that their atomic
positions do not change over time, however this can sometime interfere with system drift
velocities or available free energies because the thermostats and barostats use metrics like atom
positions to sample the box pressures and temperatures which can subsequently effect the
rescaling of the box size and thus future atom velocities. Most of the time there will not be any
external restraints applied to our lipid bilayer systems, however there are a few circumstances
which require restraints in order to proceed.
Constraints however will be present in every simulation because we need to ensure that
bond lengths, angles, and dihedrals (both proper and improper) are maintained after each time
step. The time steps in classical molecular dynamics are limited by the period of bond
oscillations, however most of these high frequency oscillations are quantum in nature and
diverge from classical oscillator frequencies. As mentioned previously we can deal with such
inconsistences between classical and quantum mechanics by constraining each bond as a rigid
rod to mimic the ground-state quantum oscillator. In GROMACS we can do this using the
LINCS (LINear Constraint Solver) algorithm [63] or the SHAKE algorithm [145].
29

SHAKE was the first holonomic constraint algorithm used with GROMACS and all it
does is change an unconstrained set of coordinates into a constrained set of coordinates based on
a predefined list of holonomic constraints. SHAKE does this by method of Lagrange multipliers
which it obtains through differentiation of each holonomic constraint equation. Then each
multiplier is added in to the equations of motion until a specific tolerance is achieved by
minimizing the associated constraint force, however because each bond is intertwined with other
bonds this often becomes a nonlinear problem to solve. When a constraint force, and thus a
constraint acceleration is calculated, SHAKE adjusts each atom’s position over a distance in
which the constraint acceleration would push it over a single integration timestep. Interestingly,
constraining water OH bonds often takes up more than 80% of the simulation time, so a special
technique is used for water called SETTLE [109] which uses a solved analytical version of
SHAKE for fast convergence.
The LINCS algorithm however is a more recent constraint technique and it is often
preferred in GROMACS because it is 3-4 times faster than SHAKE while converging under a
wider variety of conditions than its predecessors. Essentially LINCS resets bond lengths after
unconstrained atomic motion occurs in just two fixed steps whereas SHAKE iteratively solves
systems of Lagrange multipliers for an undetermined number of steps without the guarantee of
success. LINCS is able to work so efficiently because it creates a constraint coupling matrix
which is symmetric and sparse, thus dramatically reducing the time required to perform matrix
multiplication.

30

1.3.9 Analysis

Many of the post-simulation analyses found in this dissertation are a result of custom-built Perl
scripts which vary in form and function, however a majority of these scripts typically read output
data from the molecular dynamics integrator (e.g. a GROMACS trajectory) and produce specific
results or measurements based on the contents of that trajectory. Trajectories typically include
atom positions, atom charges, atom masses, water dipole moments, system energetics, or other
physical characteristics like temperature or pressure. GROMACS also contains a set of its own
analysis tools which can directly read trajectory information and produce charge density profiles,
local electric field profiles, local electric potential profiles, radial distribution functions, and
other useful quantities which may appear throughout this dissertation, however each tool will be
discussed in greater detail as they are used in later chapters. Visual Molecular Dynamics (VMD)
[75] and GNUPlot [137] were also used in later chapters to generate graphics and movies, and
the COMSOL finite-element modeling package was used in a handful of circumstances to model
the electrostatic properties around fixed geometries.

31

Chapter 2: Effect of Oxidized Lipids on
Electroporation

Portions of this chapter are adapted from:
P. T. Vernier, Z. A. Levine, Y. H. Wu, V. Joubert, M. J. Ziegler, L. M. Mir, and D. P. Tieleman,
"Electroporating Fields Target Oxidatively Damaged Areas in the Cell Membrane," PLoS One 4
(11) (2009).

2.1 Abstract

Reversible electropermeabilization (electroporation) is widely used to facilitate the
introduction of genetic material and pharmaceutical agents into living cells. Although
considerable knowledge has been gained from the study of real and simulated model membranes
in electric fields, efforts to optimize electroporation protocols are limited by a lack of detailed
understanding of the molecular basis for the electropermeabilization of the complex
biomolecular assembly that forms the plasma membrane. We show here, with results from both
molecular dynamics simulations and experiments with living cells, that the oxidation of
membrane components enhances the susceptibility of the membrane to electropermeabilization.
Manipulation of the level of oxidative stress in cell suspensions and in tissues may lead to more
efficient permeabilization procedures in the laboratory and in clinical applications such as
electrochemotherapy and electrotransfection-mediated gene therapy.

32

2.2 Introduction
External electric fields of sufficient strength and duration cause a rapid increase in the
electrical conductance of biological membranes, with an associated increased permeability to
ions and small and large molecules [60, 81, 119]. If the dose is limited, cells can survive the
treatment. Electropermeabilization (also called electroporation) technology is widely used in
laboratories to facilitate transfection in cells and tissues [61, 108, 138] and recently has appeared
in the clinic as a component of systems for electrochemotherapy [102] and tumor killing and
ablation [49, 120, 121, 144]. Although the physical and electrochemical fundamentals of electric
field-induced permeabilization of lipid membranes are well known [1, 29, 210], the mechanistic
details of the membrane restructuring that follows electric field exposure in living cells have not
been definitively established. Electrical measurements

[11], flow cytometry [112], and
fluorescence microscopy [191] indicate that permeabilization can occur in less than 10 ns,
implying a direct rearrangement of membrane components, but real-time analysis of any kind at
this time scale is difficult [45, 136]. Observations with living cells [48], artificial membranes and
lipid vesicles [38, 40, 85, 174], continuum electrophysical models [80, 161], and molecular
dynamics (MD) simulations of phospholipid bilayers [16, 58, 171, 176, 195] provide a valuable
perspective, but significant gaps remain between these model systems and the complexity of the
tissue in an electropermeabilized tumor.
Optimizing the efficiency of electroporation methods requires, in addition to knowledge
of the effects of varying the relevant physical parameters (electric field strength and duration,
temperature, cell concentration, composition of the medium, electrode arrangement) [107, 168,
173], an understanding of the relation between the physiological state of cells and their
susceptibility to electropermeabilization. It may be, for example, that starved cells can be most
33

effectively treated under conditions that are very different from those that are optimal for
actively respiring cells.
Because oxidative stress is readily imposed and commonly encountered in cultured cells
and in whole organisms under a variety of conditions, because it impacts cellular activities across
the metabolic spectrum, and because it directly affects the physical properties of the cell
membrane [51, 163, 172], we speculated that this might be an important factor in the
effectiveness of electroporation methods. Studies of the peroxidation of membrane lipids after
electroporation have been reported [20, 46, 97, 208], but the effects of pre-treatment oxidative
stress have received little experimental attention, and so we investigated the effect of membrane
oxidative damage on the sensitivity of cells to permeabilizing electrical pulses. Molecular
dynamics simulations have recently shown that incorporating oxidized lipids into phospholipid
bilayers increases the water permeability of these membranes [205], suggesting that bilayers
containing oxidized lipids will also electroporate more readily (the formation of membrane-
spanning water defects is one of the initial steps in molecular dynamics representations of
electroporation). That is the hypothesis that we tested in the work reported here, in simulations
and in living cells.

2.3 Methods

2.3.1 Molecular Dynamics Parameters

All simulations were performed using the GROMACS set of programs version 3.3.1
[182] on the University of Southern California High Performance Computing and
Communications Linux cluster (http://www.usc.edu/hpcc/). Lipid parameters were derived from
34

OPLS, united-atom parameters [15] modified for PLPC and oxPLPC [205]. We used the Simple
Point Charge (SPC) model [12] for water. 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 [14]. 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 [63] for lipids and SETTLE [109]
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 [41] using fast Fourier
transforms and conductive boundary conditions. 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.

2.3.2 Electroporation Simulations

To determine a baseline porating electric field [209], simulations of equilibrated (constant
area per lipid), fully hydrated (40 water/lipid) PLPC bilayer systems containing 72 lipid
molecules and 2880 water molecules were run with applied electric fields ranging from 300 to
500 mV/nm. The lowest field which results in pore formation within 25 ns in at least one of three
parallel PLPC simulations is 360 mV/nm, and that value was used for comparisons of pore
formation times in oxidized (PLPC:oxPLPC) and non-oxidized (PLPC) systems. Oxidized lipid
systems included 8 (11.1%), 18 (25%), or 36 (50%) randomly distributed molecules of oxPLPC,
either 12-al or 13-tc [205].
35

Electroporation times were calculated by measuring at each time step the number of
phosphorus atom groups in a system, where a group is a cluster of phosphorus atoms each no
more than 1.2 nm from another phosphorus atom [159]. The two phosphorus atom groups in an
intact bilayer (one in each of the two leaflets) merge when the bilayer interior is bridged by a
hydrophilic pore. We call the time at which this occurs the electroporation time. In some
simulations the merged groups split again, in each case in less than 400 ps after merger. This was
not considered pore formation.

2.3.3 Assembly of larger systems.

Systems with four times the area were created by doubling pure PLPC systems in x and y
using the GROMACS function genconf. Using custom code, oxidized groups were inserted on
the sn-2 linoleate tails in two opposite quadrants to create quilted systems where two quadrants
contained pure PLPC and two quadrants contained 50% oxPLPC, either 12-al or 13-tc. For 12-al
an aldehyde group was added at C12. For 13-tc, a hydroperoxy group was added at C13, and the
double bond at C12 was shifted to C11. The system was energy-minimized for 1 ps using the
GROMACS steepest descent method to test for bad contacts each time a new oxidized group was
inserted. As with the smaller systems, the oxPLPC lipids were distributed equally between the
two leaflets of the bilayer, and the water-to-lipid ratio was 40. Each system was then equilibrated
long enough (11.5 ns for 12-al, 32.5 ns for 13-tc) to bring the area per lipid within 2% of its final
equilibrated mean value while maintaining the four-quadrant distribution of the lipids. Molecular
graphics images were generated with Visual Molecular Dynamics (VMD) [75].

36

2.3.4 Cell Lines and Preparation

Jurkat T lymphoblasts (ATCC TIB-152) were grown in RPMI 1640 medium (Mediatech)
containing 10% heat-inactivated fetal bovine serum (FBS; Mediatech), 2 mM L-glutamine
(Invitrogen), 50 units/mL penicillin (Gibco), and 50 μg/mL streptomycin (Gibco) at 37 °C in a
humidified, 5% carbon dioxide atmosphere. DC-3F adherent cells (Chinese hamster lung
fibroblast cells) were grown in Minimum Essential Medium (Invitrogen, France) containing 10%
heat-inactivated FBS (Invitrogen), 500 U/ml penicillin, 500 µg/ml streptomycin (Invitrogen) at
37 °C in a humidified, 5% carbon dioxide atmosphere.
Jurkat cells were exposed to peroxidizing conditions — H
2
O
2
(200 and 500 µM) and
FeSO
4
(400 and 1000 µM) in RPMI 1640 — for 15 minutes before pulse exposure. After
peroxidation, 5 μM YO-PRO-1 (Molecular Probes, Invitrogen) was added as a permeabilization
indicator, and the cells were immediately pulsed. DC-3F cells were incubated with 1 µM calcein-
AM (Sigma-Aldrich, France) in culture medium for 1 hour at 37 °C, then rinsed with PBS,
trypsinized, and centrifuged at 1000 rpm for 10 minutes. For peroxidation, DC-3F cells were
incubated with H
2
O
2
and FeSO
4
(each 1000 µM for long pulse experiments, 1500 µM for short
pulse experiments) for 1 hour at 37 °C, rinsed with PBS and centrifuged at 1000 rpm for 10
minutes. The pellets were suspended in sucrose buffer (250 mM sucrose, 10 mM Tris, 1 mM
MgCl
2
— Sigma), pH 7, containing 2% low-melting agarose (Tebu-Bio, France) at about 50 000
cells/mL. Cell suspensions were kept at 37 °C until electric pulse exposure.

37

2.3.5 Pulse Generator and Pulse Exposures

For the Jurkat cell experiments 30 ns, 3 MV/m pulses were delivered at a 50 Hz
repetition rate to cell suspensions in commercial electroporation cuvettes (VWR) with a 1 mm
electrode spacing in ambient atmosphere at room temperature from a USC pulse generator based
on a magnetic compression, diode opening-switch architecture [170]. For the DC-3F experiments
the cell suspension was placed within the 2 mm gap between two copper electrodes fixed to a
microscope slide. Cells were exposed to 1 long pulse (100 µs, 50 or 60 kV/m) delivered by a
micropulse generator (Cliniporator, IGEA, Italy) or to 1000 short pulses (10 ns, 2.5 MV/m, 10
Hz) delivered by a FID Technology (Russia) pulse generator.

2.3.6 Fluorescence Miscroscopy and Microphotometry

Jurkat cell images were captured and analyzed with a Zeiss AxioCam MRm and
AxioVision 3.1 software (Carl Zeiss Goettingen, Germany) on a Zeiss Axiovert 200M
epifluorescence microscope. Intracellular YO-PRO-1 (λ
ex
= 480 nm, λ
em
= 535 nm) fluorescence
was used as an indicator of membrane permeabilization [66, 76]. Membrane lipid peroxidation
was monitored qualitatively with the fluorescent dye C11-BODIPY
581/591
. Upon oxidation, the
dye shifts from a red-emitting form (595 nm) to a green-emitting form (520 nm) [127]. Cells
were incubated in RPMI 1640 medium containing 4 µM C11-BODIPY
581/591
for 30 minutes at 37
°C before exposure to peroxidizing reagents. Images of calcein-loaded DC-3F cells were taken
immediately before and 10 min after pulse delivery at λ
ex
= 496 and λ
em
= 516 nm, using a Zeiss
AxioCam HRC camera coupled to a Zeiss Axiovert S100 microscope. Intracellular fluorescence
was averaged over more than 300 cells. Each experiment was repeated three times. For each
38

repetition of the calcein efflux experiments the intracellular fluorescence (F) was measured on
the same cells (each time, more than 300 cells) before (F
b
) and 10 minutes (F
10
) after the pulses
delivery. The normalized decrease in fluorescence (F
b
-F
10
)/F
b
was averaged (+/- SD) over the 3
repeats. * = p < 0.05; ** = p < 0.01; NS = no statistically significant difference.

2.4 Results

2.4.1 Simulations

Pore formation time for PLPC and oxidized PLPC bilayers. Two oxidized variants of 1-
palmitoyl-2-linoleoyl-sn-glycero-3-phosphatidylcholine (PLPC) were chosen for this study based
on the increased water permeability of PLPC bilayers containing these species in molecular
dynamics simulations [205]. Construction of the oxidized PLPC (oxPLPC) molecules containing
the modified linoleoyl residues 12-oxo-cis-9-dodecenoate (12-al) and 13-hydroperoxy-trans-11,
cis-9-octadecadienoate (13-tc) has been previously described [205]. Composite snapshots of
PLPC and the 12-al and 13-tc variants taken from the simulations with external electric fields
reported here (Figure 2.1) show the location of the oxidized molecular modifications and how
the conformation of the oxidized lipid tails contributes to an increase in area per lipid when these
species are incorporated into a PLPC bilayer (Table 2.1). The stabilization associated with
hydration of the tail oxygens in the oxidized lipids results in a tendency for the oxidized tails to
spend more time near the aqueous interface than the hydrocarbon tails of PLPC. In systems
containing 11% oxPLPC the mean distance from the phosphorus plane (membrane interface
39

region) to the C-12 oxygen atom in 12-al (0.6 nm) and to the C-13 oxygens in 13-tc (0.3 nm) is
much less than that for the sn-2 C-13 in PLPC (1.6 nm).

Figure 2.1 Oxidized and unoxidized phospholipid conformations change over time. Composite
snapshots (21 images captured at 0.5 ns intervals over a 10 ns period) of PLPC and two
oxidized variants, oxPLPC (12-al) and oxPLPC (13-tc), showing their conformations in
molecular dynamics simulations of PLPC with 11% oxPLPC bilayers in a 360 mV/nm field.
The spheres near the end of the lipid tails mark the location of the introduced oxygens or C-13
of PLPC. Structures of the individual lipid molecules are shown below the corresponding
composite. Teal – C, red – O, gold – P, blue – N, gray – C-13.

40

Table 2.1 Effect of oxidized lipid concentration on electropore formation time in oxPLPC:PLPC
bilayers.

Bilayer thickness (distance between the two mean phosphorus planes) and area per lipid are
mean values taken from three independent simulations for each condition over 100 ps time steps
from the beginning of the simulation until 1 ns before pore formation. For systems where
poration occurred in less than 1.5 ns, the cutoff is 0.5 ns before pore formation. Bilayer thickness
values from independent simulations vary by less than 10%. Area per lipid is the area of the
bilayer divided by the number of lipids in one leaflet of the bilayer - 36. Individual and mean
times to poration are shown for three independent simulations of each system. Systems contain
2880 water molecules and 72 total lipids. The applied electric field is 360 mV/nm.

To observe the effect of oxidized lipids on the sensitivity of a bilayer to electroporating
fields, equilibrated systems composed of 72 lipids and 2880 waters with varying percentages of
PLPC and one of the two oxPLPC species were assembled and subjected to an external electric
field of 360 mV/nm. The results are shown in Table 2.1. (Trial simulations showed that for pure
oxPLPC Thickness Area/Lipid Trial Poration Time (ns)
(nm) (nm
2
) Individual Mean (S.D.)
1 > 25
None 0% 3.8 0.68 2 > 25 > 22.6
3 22.6
1 16.5
11% 3.4 0.68 2 6.3 8.4 ± 7.3
3 2.3
1 12.5
12-al 25% 3.5 0.69 2 6.2 7.6 ± 4.4
3 4.0
1 2.9
50% 3.3 0.72 2 1.1 1.8 ± 1.0
3 1.4
1 > 55
11% 3.8 0.69 2 43.0 > 22.9
3 2.7
1 22.8
13-tc 25% 3.7 0.70 2 13.4 14.4 ± 7.9
3 7.1
1 5.4
50% 3.6 0.72 2 4.9 4.7± 0.9
3 3.7
41

PLPC bilayers this field strength leads to formation of a pore — a membrane-spanning column
of water approximately 1 nm diameter surrounded by phospholipid head groups [176] — in less
than 25 ns in at least one of three independent trajectories.) Increasing the oxidized lipid content
decreases the bilayer thickness, increases the average area per lipid, and reduces the time to
poration in an external electric field.
Additional 25 ns simulations with systems containing 50% oxPLPC show that even when
the field is reduced far below the minimal porating value for pure PLPC, the oxidized bilayers
still form pores, and that membranes containing 12-al porate more readily than 13-tc bilayers
(Table 2.2). Although 12-al bilayers have a slightly smaller area per lipid (which might be
expected to restrict the entry of water into the membrane interior), they are thinner than 13-tc
bilayers, which lowers the energy barrier for water penetration. In addition, the 12-al aldehyde
group spends more time in intermediate locations between the aqueous interface and the low
dielectric membrane mid-plane than is the case for the 13-tc hydroperoxy oxygens (Figure 2.1).
As we demonstrate below, the partially hydrated aldehyde group facilitates formation of the
membrane-spanning water column that marks the initiation of poration.

Table 2.2 Effect of electric field magnitude on electropore formation time in 50%
oxPLPC:PLPC bilayers.
Mean Poration Time (ns)
150 mV/nm 200 mV/nm 250 mV/nm 300 mV/nm
12-al 10.4 5.7 3.7 3.8
13-tc > 25 14.1 10.2 7.3

Field-dependent poration occurs in oxPLPC:PLPC bilayers even at fields where no poration is
observed in pure PLPC membranes during 25 ns simulations. 12-al:PLPC bilayers porate more
readily than 13-tc:PLPC bilayers. Systems contain 2880 water molecules, 36 PLPC, and 36
oxPLPC (12-al or 13-tc).
42

Pore formation occurs at local concentrations of oxidized lipids. Because the initial steps
leading to electroporation involve the intrusion of water into the bilayer interior, and because
bilayers containing oxidized species are more permeable to water, we expected that the site of
pore initiation in our simulations might be specifically associated with single or aggregated
oxPLPC molecules. This is what we observe in every case. The snapshots in Figure 2.2 are
frames taken from a representative simulation (11% 12-al:PLPC in a 360 mV/nm electric field)
just before the appearance of a membrane spanning water column (Figure 2.2a) and just after a
pore has formed (Figure 2.2b), showing clearly the association of the intruding water with the
oxidized residues. The initial water defect in bilayers containing 12-al or 13-tc is always
associated with the aldehyde or peroxy oxygens in the oxidized tail. To visualize this more
clearly we created “quilted” PLPC bilayers, with some sections containing 50% oxPLPC and
others containing only PLPC, and subjected them to a porating electric field (Figure 2.3). In all
of our simulations electropores invariably form in immediate association with one or more
oxidized lipids.

43

Figure 2.2 Poration of an oxidized phospholipid bilayer. Snapshots of a PLPC system with an
11% concentration of 12-al oxPLPC before (a) and after (b) an electropore is formed (separated
by about 2 ns.) Only water (small red and white "v's") and the head groups and sn-2 tails of the
oxPLPC molecules are shown. The large red spheres near the ends of the tails are aldehyde
oxygens, which appear to facilitate the entry of water into the bilayer interior. Dimensions of
the simulation box are approximately 5 nm x 5 nm x 7 nm.
Figure 2.3 Quilted PLPC:oxPLPC bilayer. The simulated system, bounded by the black
square in panel A, is divided (dashed lines) into two regions of approximately 100% PLPC
(light gray) and two regions of approximately 50% oxPLPC (12-al) (dark blue) and 50%
PLPC, as described in the text. To show more clearly where poration occurs after application
of an external electric field, copies of the simulated system are tiled to make a periodic 2 x 2
system (approximately 10 nm x 10 nm x 7 nm). Preferential electroporation of the PLPC
bilayer in regions of high oxidized lipid content is demonstrated in panel B, which shows the
system 1 ns after applying a 360 mV/nm field normal to the bilayer. The bilayer is in the plane
of the page.
44

2.4.2 Observations with Living Cells

To determine whether these simulations are indicative of the responses of living cells, we
treated Jurkat T lymphoblasts with sublethal doses of peroxidizing agents (hydrogen peroxide
and ferrous sulfate) and then monitored the influx of the fluorescent dye YO-PRO-1, a sensitive
indicator of membrane permeabilization [191], after exposure to ultra-short (30 ns), high-
intensity (3 MV/m) electric pulses. Conditions for the peroxidation treatment were developed by
maximizing the green:red fluorescence emission ratio from the membrane-staining fluorescent
dye C11-BODIPY
581/591
[114] while at the same time minimizing changes in cell morphology
and membrane integrity. Results are shown in Figure 2.4. Pulse-induced YO-PRO-1 uptake in
peroxidized Jurkat cells is significantly higher than in non-oxidized cells, as predicted by the
simulations.
45

In a parallel set of experiments with calcein-loaded DC-3F (Chinese hamster lung
fibroblast) cells, we monitored the decrease in intracellular calcein fluorescence caused by two
kinds of permeabilizing electric pulses: ultra-short (10 ns), high-intensity (2.5 MV/m)
nanoelectropulses; and long (100 µs), low-field (50 kV/m) pulses typical of those used in
Figure 2.4 Electropermeabilization enhancement by treatment with peroxidation agents. (A).
Fluorescence images of Jurkat cells showing pulse-induced (30 ns, 3 MV/m, 50 Hz) YO-
PRO-1 influx into control and peroxidized cells after 10 min exposure to 500 μM H
2
O
2
+
1000 μM FeSO
4
. (B). Integrated YO-PRO-1 fluorescence intensity from more than 300
individual cells from three independent experiments for each condition. Error bars are
standard error of the mean.
46

electrotransfection and electroporation protocols. This method for detecting permeabilization,
which takes advantage of the outward flux of normally impermeant calcein molecules across a
steep concentration gradient and into the large reservoir represented by the external medium may
be more sensitive than standard procedures for measuring electroporation, which rely on the
influx of dye molecules across a decreasing concentration gradient into the relatively small
intracellular volume.
The loss of fluorescence from calcein-loaded DC-3F cells measured 10 minutes after
exposure to one thousand 10 ns, 2.5 MV/m pulses delivered at 10 Hz is much greater for
peroxidized cells than for the controls, consistent with the Jurkat cell results. The same is true for
the response to a single 100 µs, 50 kV/m pulse (Figure 2.5). Figure 2.5 also shows that
increasing the pulse amplitude to 60 kV/m produces significant permeabilization of untreated
cells, indicating that the 100 µs, 50 kV/m dose is near the threshold for detectable
permeabilization for these cells, and that membrane peroxidization reduces the value of this
threshold.

47

2.5 Discussion

Molecular dynamics simulations of lipid bilayers and laboratory studies of model membranes
and cells in electric fields are consistent with the stochastic pore hypothesis for
electropermeabilization [117, 133, 167, 199]. Here we have shown that it is possible to bias pore
formation at the molecular level by oxidatively modifying the properties of the membrane in situ,
so that we have now a new method for lowering the energy barrier to poration and a mechanism
that might be used to localize where in the membrane poration occurs. This demonstration of the
Figure 2.5 Enhanced electropermeabilization of peroxidized cells with both ultra-short and
conventional electric pulse treatments. (A). The fluorescence of calcein-loaded DC-3F
(Chinese hamster lung fibroblast) cells exposed to 1000 10 ns, 2.5 MV/m pulses with a 10 Hz
repetition rate decreases after 10 min, indicating membrane permeabilization. The effect is
much greater in peroxidized cells than in untreated control cells. (B). Peroxidized cells treated
with a single 100 µs, 50 kV/m pulse show a similar increased susceptibility to
electropermeabilization. Increasing the 100 µs pulse amplitude to 60 kV/m results in
significant permeabilization of control cells, indicating that these doses are near the threshold
for a detectable response under these conditions. Peroxidized cells treated with a single 100
µs, 60 kV/m pulse show a strong decrease of fluorescence, indicating substantial membrane
permeabilization. Bars are the mean, and error bars are the standard deviation. * p < 0.05; **
p < 0.01; *** p < 0.001.
48

increased sensitivity of oxidatively damaged cells to electropermeabilization has practical
implications for the laboratory and the clinic.

2.5.1 Electrotransfection and Electropermeabilization Enhancement

The results suggest that an appropriately controlled peroxidizing regimen, that is, one
which enhances permeabilization without significantly affecting viability, might increase the
efficiency of electrotransfection protocols either indirectly by enabling the use of lower porating
voltages (higher voltages are associated with lower cell survival rates) or directly by increasing
the amount of genetic material that enters the cells for a given series of electrical pulses.
Alternatively, a reducing environment would be expected to protect cells against
electropermeabilization. Extending this line of thinking, electrochemotherapy and direct ablation
and killing of tumor cells using electrical pulse therapy may likewise be enhanced by procedures
which promote oxidative stress in the tumor tissue before pulse delivery — and inhibited when
the environment in or around the tumor is reducing [33]. In any mixed population of cells, the
different native sensitivities of various cell types to permeabilizing electric fields [25, 154, 189]
and oxidative stress might be exploited by adjusting electrical, physical, and chemical parameters
to selectively transfect subsets of cells, in vitro and in vivo.
Because the simulations show the molecular-scale localization of electroporation
enhancement by oxidized membrane lipids, it is not unreasonable to infer that improper handling
of cell suspensions, resulting in uncontrolled and unmonitored degrees of oxidative stress, may
account for inconsistent and variable results of electrotransfection and other methods dependent
on reversible electropermeabilization. Cells maintained at high concentrations in rich medium, or
49

for long periods in non-nutritive buffers may be unexpectedly and unpredictably sensitive to
porating electric pulses. Likewise, adding reducing agents to cell suspensions to protect labile
compounds or the cells themselves [22] may unintentionally result in decreased
electropermeabilization and electrotransfection yields. Investigations of these possibilities,
although they will of necessity be protocol- and cell type-specific, have the potential to lead to
more consistent and more efficient membrane permeabilization procedures.

2.5.2 Oxidative Stress and Apoptosis

The relationship between oxidative stress and electropermeabilization sensitivity may
have another dimension. It is known that exposure to nanosecond, megavolt-per-meter electric
pulses can not only permeabilize membranes [126, 191] but also induce apoptosis [9, 189]. The
role of reactive oxygen species (ROS) in apoptosis has been investigated [28, 185, 207], and we
report here that oxidative stress also modulates electropermeabilization. A more extensive
analysis of oxidative stress and short and long pulse electroporation may lead to a more
sophisticated approach to controlling electrotransfection yield, including both the efficiency of
the incorporation of the genetic material and the subsequent viability of the cells, which may be
decreased by pulse-induced apoptosis. Similar considerations apply to potential improvements in
electroporation protocols. Thus the knowledge gained from an exploration of the interplay
between nanosecond membrane permeabilization [124, 191] and conventional electroporation
[70] in the presence of excess reactive oxygen species and the initiation of apoptosis may result
in improved electroporation and electrotransfection procedures, with benefits from cell science to
cancer therapeutics.
50

2.5.3 Permeabilizing Defects and Membrane Boundaries

Facilitation of electric field-driven water defect formation by the aldehyde and
hydroperoxy oxygens of the oxPLPC species reported here may in fact be a relatively simple
example of a more general tendency for water intrusion, permeabilization, and other membrane
restructuring events to occur at membrane phase or domain boundaries, especially where these
discontinuities have charged or hydrophilic components that extend into the membrane interior
[82, 96, 105, 167]. In this context it will no doubt prove instructive to examine the combined
effects of lipid peroxidation and cholesterol on membrane electropermeabilization. Adding
cholesterol to a bilayer is likely to render the membrane more difficult to electropermeabilize
[10], but the consequences of incorporating peroxidized lipids, which promote the formation of
cholesterol domains and lipid rafts [7, 77], into the system are difficult to predict. Increasing
computational power makes it possible now to address these more complex (and more realistic)
systems with an approach that ties together atomically detailed simulations and experimental cell
biology.

2.6 Acknowledgements

We thank Martin Gundersen for stimulating discussions and support. Computing resources
were provided by the USC Center for High Performance Computing and Communications. This
work was made possible in part by MOSIS, Information Sciences Institute, Viterbi School of
Engineering, University of Southern California.

51

Chapter 3: Life Cycle of an Electropore

Portions of this chapter are adapted from:
Z. A. Levine and P. T. Vernier, "Life Cycle of an Electropore: Field-Dependent and Field-
Independent Steps in Pore Creation and Annihilation," Journal of Membrane Biology 236 (1),
27-36 (2010).

3.1 Abstract

Electropermeabilization, an electric field-induced modification of the barrier functions of
the cell membrane, is widely used in laboratories and increasingly in the clinic, but the
mechanisms and physical structures associated with the electro-manipulation of membrane
permeability have not been definitively characterized. Indirect experimental observations of
electrical conductance and small molecule transport and molecular dynamics simulations have
led to models in which hydrophilic pores form in phospholipid bilayers with increased
probability in the presence of an electric field. Presently available methods do not permit the
direct, nanoscale examination of electroporated membranes that would confirm the existence of
these structures. To facilitate the reconciliation of poration models with the observed properties
of electropermeabilized lipid bilayers and cell membranes, we propose a scheme for
characterizing the stages of electropore formation and re-sealing. This electropore life cycle,
based on molecular dynamics simulations of phospholipid bilayers, defines a sequence of
discrete steps in the electric field-driven restructuring of the membrane that leads to the
52

formation of a head group-lined, aqueous pore and then, after the field is removed, to the
dismantling of the pore and reassembly of the intact bilayer. Utilizing this scheme we can
systematically analyze the interactions between the electric field and the bilayer components
involved in pore initiation, construction, and re-sealing. We find that the pore creation time
depends strongly on the electric field gradient across the membrane interface and that the pore
annihilation time is at least weakly dependent on the magnitude of the pore-initiating electric
field and is in general much longer than the pore creation time.

3.2 Introduction

The structural integrity of the cell membrane, essential for compartmentalization of
extracellular and intracellular materials, can be disrupted with electrically-induced
permeabilization (electroporation) [60, 141], which results in increased electrical conductance
and the influx and efflux of normally impermeant small and large molecules [108, 119, 139].
Small changes in conductance can be measured within a few nanoseconds after the application of
a pulsed electric field [11], and fluorescent dyes can be used as molecular markers to indicate the
degree of membrane permeabilization, but the detailed mechanisms of electropermeabilization
are difficult to access experimentally [173], even with nanoscale instrumentation such as atomic
force microscopy (AFM). Other analytical methods, including continuum and atomic-level
molecular simulations, provide windows for investigating electropore formation at the nanoscale,
where there are at present no direct observational tools.
In molecular dynamics (MD) simulations of electroporation [171, 176] the electric field-
driven formation of a water column across the bilayer interior is followed by the construction of
53

a bridge of hydrophilic lipid head groups and additional hydrating water molecules. Although it
has been shown that the pore creation time decreases as the pore-initiating electric field increases
[209], the detailed time course of electropore construction and annihilation in MD simulations
has not been quantitatively characterized, though several pioneering studies have described the
key steps in the electroporation process [171, 175, 176]. Here we propose a scheme for
objectively classifying and delineating the stages in the life cycle of an electropore.
The establishment of objective criteria for analyzing the appearance and disappearance of
membrane pores in electric fields permits quantifying and tying more closely to molecular
mechanisms the dependence of poration-associated processes on the electric field strength and
other variables, and facilitates comparisons of electropermeabilization in different model
systems. Because electroporation is strongly dependent on the interactions between interface
water molecules and lipid headgroups [171, 175, 209], even small differences in force fields and
other simulation parameters may be expected to produce significant differences in electropore
creation and annihilation dynamics. Examples are provided.

3.3 Materials and Methods

3.3.1 Simulation Conditions

All simulations were performed using the GROMACS set of programs version 3.3.3
[182] on the University of Southern California High Performance Computing and
Communications Linux cluster (http://www.usc.edu/hpcc/). Lipid topologies were derived from
OPLS united-atom parameters [15] and were obtained from Peter Tieleman of the University of
54

Calgary (http://moose.bio.ucalgary.ca). The Simple Point Charge (SPC) water model was used
[12]. 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 [14].
Pressure was coupled semi-isotropically (using a compressibility of 4.5 x 10
-5
bar
-1
) normal to
and in the plane of the membrane. Bond lengths were constrained using the LINCS algorithm
[63] for lipids and SETTLE [109] for water. Short-range electrostatic and Lennard-Jones
interactions were cut off at 1.0 nm. Long-range electrostatics were calculated by the PME
algorithm [41] using fast Fourier transforms and conductive boundary conditions. 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.

3.3.2 Systems and Structures

All systems contained 128 lipid molecules — 1-palmitoyl-2-oleoyl-sn-glycero-3-
phosphatidylcholine (POPC) or 1,2-dioleoyl-sn-glycero-3-phosphatidylcholine (DOPC) — and
at least 4480 water molecules (~35 waters/lipid), which resulted in a system box size of
approximately 7 nm x 7 nm x 7 nm. POPC systems used for tracking the complete life cycle,
both pore creation and annihilation, contained more water (~70 waters per lipid) and thus were
about 10 nm in the z dimension. Simulations with multiple trials were run in parallel, starting
from the same initial positions. To ensure that each trial was independent, every atom was
55

assigned a randomized velocity from a Maxwell distribution at the beginning of the simulation (a
built-in function of GROMACS). All systems were equilibrated until they attained a constant
area per lipid: POPC — 0.66 nm
2
; DOPC — 0.65 nm
2
; second DOPC model (see below) — 0.67
nm
2
. Pore creation times were measured at four external electric fields — 300 MV/m, 400
MV/m, 500 MV/m, and 600 MV/m — with three independent trials at each electric field
strength. (300 MV/m in vacuum corresponds to an effective electric field in aqueous media of
about 4 MV/m.) POPC pore annihilation times were taken from a separate set of simulations in
which the pore-initiating external electric fields were 320 MV/m, 400 MV/m, and 600 MV/m. At
least three trials per field were used for annihilation times. We tested two DOPC models, an
older version which was obtained from Peter Tieleman in 2005 (personal communication) and a
newer version, currently available at http://moose.bio.ucalgary.ca/. We refer to these as DOPC
(2005) and DOPC (2009), respectively. These models differ slightly in the charge distribution for
the atoms in the head group and glycerol backbone.

3.3.3 Electropore Life Cycle

We divide the electropore life cycle broadly into pore creation and pore annihilation. Pore
creation is comprised of three stages: initiation, construction, and maturation (defined below).
Pore annihilation begins when the external electric field is removed from a mature pore structure
and proceeds through destabilization, degradation, deconstruction, and dissolution (Figure 3.1).
56

3.3.4 Pore Creation

Pore initiation begins with the application of an external electric field and ends when the
two groups of water molecules, initially separated by the bilayer, merge to become a single
group (non-zero water density in every 100 pm slice between the two mean planes of phosphorus
atoms). In some simulations the combined water groups split again, in less than 400 ps after the
Figure 3.1 Life cycle of an electropore. Only water and phosphorus atoms are shown for
simplicity. Pore creation in an electric field begins with the introduction of a water defect into
the bilayer interior (pore initiation), followed by the reorganization of phospholipid head
groups in each leaflet around the defect (pore construction). Migration of additional water and
head groups into the pore continues until an arbitrarily defined mature pore structure is
formed (pore maturation). Pore annihilation begins with the removal of the porating electric
field. The pore structure is quasi-stable at this time (pore destabilization) but soon there is a
decrease in pore size as head groups and water begin to migrate out of the membrane interior
(pore degradation). The head groups separate again into two groups (pore deconstruction).
Water quickly follows (pore dissolution), and the intact structure of the bilayer is restored.
57

joining in each case. These transient events were not counted as mergers. Pore construction
begins with the formation of the membrane-spanning water column (hydrophobic pore [1]) that
marks the merger of the water groups and ends when the phosphorus groups that are initially
found on the two leaflets of the bilayer follow the water into the membrane interior and merge
into a single phosphorus group. A phosphorus group is defined as a set of atoms, each separated
by a maximum distance of 1.2 nm [159]. Because water and the charged phospholipid head
groups now bridge the membrane interior, this structure is comparable to what has been called in
several contexts a hydrophilic pore [52, 92, 201]. Continued application of the porating electric
field results in an evolution, or maturation, of the hydrophilic pore. We define a mature pore
arbitrarily as a hydrophilic pore in which at least 10 phosphorus atoms from the initial anodic
leaflet are found within 1.2 nm of phosphorus atoms from the cathodic leaflet.

3.3.5 Pore Annihilation

Pore destabilization, the first stage of pore annihilation, is the quasi-stable period after the
field is removed during which the number of anode to cathode phosphorus connections fluctuates
around the mature pore criterion (10 connections). Pore degradation begins when the number of
anode to cathode phosphorus connections drops below 10 and ends when there is only one anode
to cathode phosphorus connection (Figure 3.2). The pore diameter decreases during this time to a
minimum, about 0.4–0.6 nm. Pore deconstruction, another quasi-stable period, terminates when
the single phosphorus group of the porated bilayer splits into two groups which remain separate
for the remainder of the simulation. At the end of pore deconstruction only the water column
remains. We call the disassembly of the water column pore dissolution.

58

Individual pore life cycle times were calculated using a custom Perl program which
codifies the stage boundaries described above. The local electric field profile was extracted using
the GROMACS function ‘g_potential’, which integrates the charge density found in bins (0.1 nm
in width) parallel to the plane of the bilayer. Charges are rearranged due to the application of a
user-defined external electric field which is constant throughout the system. Figure 3.3 shows a
typical electric field profile for a POPC bilayer with an external field of 300 MV/m. The average
value of the local electric field within 0.5 nm of the bilayer mid-plane is referred to here as the
‘membrane internal electric field’, and was calculated from the electric field profile averaged
Figure 3.2 Phosphorus and water distribution during pore annihilation. The number of
inter-leaflet phosphorus connections (anodic phosphorus atoms located within 1.2 nm of
cathodic phosphorus atoms) remains relatively constant for some time (destabilization
stage) after the external electric field is removed from a system containing a mature pore,
even though a significant amount of water leaves the pore during this period. During pore
degradation the number of phosphorus connections falls from 10 to 1. Pore deconstruction
ends when all phosphorus connections have been broken. Pore dissolution, the final step in
pore annihilation, terminates when water is no longer present in the bilayer interior.
59

from 100 ps after the external field was applied to 200 ps before pore construction began.
Because the net electric field in bulk water is much smaller than the field in the membrane
interior (as a consequence of the high permittivity of water and the low permittivity of the
phospholipid hydrocarbon tails), the magnitude of the membrane internal electric field may be
regarded as essentially equivalent to the electric field gradient across the membrane interface.
Additionally, we calculated the ‘box potential’ in each system by multiplying the applied electric
field times the box length in the direction of the electric field (E·L
z
). This value is comparable to
the transmembrane potential reported in other papers [171, 175] where the membrane thickness
is a large fraction of the box length. Molecular graphics images were generated with Visual
Molecular Dynamics (VMD) [75].

60

3.4 Results

3.4.1 Pore Creation and Annihilation - POPC

As models [200] and intuition predict, the pore creation time for POPC bilayers decreases
with increasing electric fields (Table 3.1). This is almost entirely a result of the time required for
pore initiation, the first stage of pore creation. Pore construction time and pore maturation time
(an end point defined somewhat arbitrarily for this initial study) are only slightly dependent on
the electric field strength.
Figure 3.3 Electric field profile for a POPC bilayer. When an external electric field is applied
to a POPC bilayer, the local electric field can be calculated by integrating the charge density in
discrete planes (parallel to the plane of the bilayer and 0.1 nm in width). This value is
determined by the rearrangement of charge in the system. The ‘membrane internal electric
field’ is calculated by averaging the local electric field at the bilayer mid-plane (± 0.5 nm) from
100 ps after the external electric field is applied to 200 ps before pore construction begins.
61

Pore annihilation times are significantly longer than pore creation times (Table 3.2). In
pore creation, once the initial membrane-spanning water column is formed, it takes less than 1 ns
for the head groups to follow the water into the bilayer interior. When the electric field is
removed however, tens of nanoseconds elapse before the head groups migrate back to the bilayer
leaflets. During pore destabilization, the period immediately following removal of the electric
field, the number of water molecules in the pore (within 0.5 nm of the bilayer mid-plane)
decreases substantially in about a nanosecond, and the pore diameter decreases to about 0.6 nm.
In a small percentage of cases pore annihilation continues to the next stage with no delay, but
usually the pore remains intact, with the number of inter-leaflet phosphorus connections
fluctuating around 10, for a time that ranges from less than one to over 20 ns (Figure 3.2). At
some point (the end of the destabilization stage) the number of anode to cathode phosphorus
atom connections permanently declines from 10 or more to one — the reverse evolution of the
developed pore to a minimum, single phosphorus group structure. The time required for this
decline, the pore degradation time, is typically several nanoseconds, significantly longer than
pore maturation, the corresponding pore creation stage. Destabilization and degradation times
have large standard deviations, indicating the stochastic nature of the processes underlying these
stages of pore annihilation. At the highest porating field, however, there is a significant reduction
in the destabilization time, suggesting that pore destabilization has a field-dependent component.
Pore deconstruction (separation into two phosphorus groups) and dissolution (expulsion of water
from the membrane interior) are apparently independent of field strength. Deconstruction times
are tens of nanoseconds, with large standard deviations; the minimal hydrophilic pore is quasi-
stable. Pore dissolution on the other hand, is rapid; the hydrophobic pore is highly unstable.

62

Table 3.1 POPC pore creation times.
Initiation Construction Maturation Creation
E
ext

V
ext

E
int

Time

Mean

Time

Mean

Time

Mean

Time

Mean

(MV/m) (Volts) (MV/m) (ns) (ns) (ns) (ns)
11.2 0.2 0.8 12.2
300 3.2 1253 10.0 9.5 ± 2.0 0.2 0.3 ± 0.1 0.4 0.6 ± 0.2 10.6 10.4 ± 1.9
7.3 0.4 0.7 8.4

2.8 0.0 0.5 3.3
400 4.2 1630 2.3 2.0 ± 1.0 0.3 0.2 ± 0.2 1.1 0.8 ± 0.3 3.7 2.9 ± 1.1
0.8 0.2 0.7 1.7

0.9 0.1 1.0 2.0
500 5.3 2180 0.6 0.6 ± 0.3 0.2 0.2 ± 0.1 0.5 0.8 ± 0.3 1.3 1.6 ± 0.4
0.4 0.3 0.9 1.6

0.5 0.1 0.1 0.7
600 6.3 2347 0.5 0.4 ± 0.2 0.0 0.1 ± 0.1 0.3 0.2 ± 0.1 0.8 0.7 ± 0.1
0.2 0.2 0.2 0.6

Table 3.2 POPC pore annihilation times.
Destabilization Degradation Deconstruction Dissolution Annihilation
E
pore

Time

Mean

Time

Mean

Time

Mean

Time

Mean

Time

Mean

(MV/m) (ns) (ns) (ns) (ns) (ns)
0.1 0.6 31.1 0.8 32.6
320 4.7 12.9 ± 12.4 6.8 5.3 ± 5.0 37.7 27.6 ± 11.2 2.9 1.3 ± 1.1 52.1 47.1 ± 15.9
21.6 2.1 11.5 1.1 36.3
25.2 11.6 30.0 0.4 67.2

7.3 3.0 24.6 2.9 37.8
400 15.5 13.0 ± 12.5 1.0 2.2 ± 1.3 7.7 18.5 ± 8.5 2.5 1.9 ± 1.0 26.7 35.6 ± 18.1
29.2 3.7 25.9 1.1 59.9
0.1 1.2 15.6 1.0 17.9

11.4 1.0 8.1 1.5 22.0
600 3.8 5.2 ± 5.7 6.4 3.4 ± 2.7 28.5 26.6 ± 17.6 0.7 0.9 ± 0.5 39.4 36.1 ± 12.8
0.3 2.9 43.1 0.6 46.9

63

3.4.2 Pore Creation – POPC versus DOPC

POPC and DOPC bilayers have similar minimum porating electric fields [209], so we
might expect similarities in their electropore life cycles. Here we report individual pore creation
stage times for DOPC (2005) and DOPC (2009) bilayers (Tables 3.3 and 3.4) and compare these
patterns with those found for POPC (Table 3.5). To simplify the comparisons, the POPC
simulations were carried out with the same system size as the DOPC simulations. DOPC pore
initiation is strongly field-dependent, similar to POPC. DOPC pore construction and maturation
times, like those for POPC, exhibit only a small dependence on electric field strength. Each of
the DOPC stage times is comparable to the corresponding POPC time. All three systems appear
to have a smaller variance in pore creation time as the strength of the external electric field is
increased, implying that stochastic effects are less dominant at higher fields.
Pore creation times in all three systems are inversely dependent on the magnitude of the
bilayer internal electric field (Figure 3.4). Membrane internal electric fields (and thus the field
gradient across the interface) in DOPC and POPC systems of the same size are similar (Tables
3.3, 3.4, 3.5). In the POPC simulations with a longer box z-dimension (Table 3.1) the membrane
internal electric field for a given external field is higher, a consequence of how the external
electric field is implemented in GROMACS [59].

64

Table 3.3 DOPC (2005) pore creation times.
Initiation Construction Maturation Creation
E
ext

V
ext

E
int

Time

Mean

Time

Mean

Time

Mean

Time

Mean

(MV/m) (Volts) (MV/m) (ns) (ns) (ns) (ns)
>25 N/A N/A >25
300 2.2 820 >25 N/A N/A N/A N/A N/A >25 >25
>25 N/A N/A >25

>25 N/A N/A >25
400 3.0 1100 21.3 >15.4 ± 8.3 0.2 0.2 ± 0.0 0.8 0.7 ± 0.1 22.3 >16.3 ± 8.5
9.5 0.2 0.6 10.3

6.1 0.3 0.4 6.8
500 3.7 1333 4.6 3.9 ± 2.6 0.1 0.2 ± 0.1 0.3 0.4 ± 0.1 5.0 4.5 ± 2.6
1.0 0.2 0.5 1.7

1.4 0.1 0.2 1.7
600 4.4 1747 1.3 1.1 ± 0.5 0.1 0.2 ± 0.1 0.3 0.2 ± 0.1 1.7 1.5 ± 0.4
0.5 0.3 0.2 1.0

Table 3.4 DOPC (2009) pore creation times.
Initiation Construction Maturation Creation
E
ext

V
ext

E
int

Time

Mean

Time

Mean

Time

Mean

Time

Mean

(MV/m) (Volts) (MV/m) (ns) (ns) (ns) (ns)
>25 N/A N/A >25
300 2.1 793 >25 N/A N/A N/A N/A N/A >25 N/A
>25 N/A N/A >25

>25 N/A N/A >25
400 2.9 1147 15.6 >14.8 ± 1.2 0.2 0.3 ± 0.1 0.7 0.5 ± 0.2 16.5 >15.6 ± 1.3
13.9 0.4 0.4 14.7

7.0 0.2 0.3 7.5
500 3.6 1437 1.6 3.1 ± 3.4 0.1 0.2 ± 0.1 0.5 0.4 ± 0.1 2.2 3.7 ± 3.3
0.8 0.3 0.4 1.5

2.4 0.1 0.3 2.8
600 4.3 1837 1.6 1.8 ± 0.6 0.1 0.1 ± 0.1 0.2 0.3 ± 0.1 1.9 2.2 ± 0.6
1.3 0.2 0.3 1.8
65

Table 3.5 Normalized POPC pore creation times
*
.
Initiation Construction Maturation Creation
E
ext

V
ext

E
int

Time

Mean

Time

Mean

Time

Mean

Time

Mean

(MV/m
)
(Volts) (MV/m
)
(ns) (ns) (ns) (ns)
>25 >25 >25 >25
300 2.1 861 >25 >25 >25 >25 >25 >25 >25 >25
>25 >25 >25 >25

15.2 0.3 0.5 16.0
400 2.9 1090 12.8 11.2 ± 5.1 0.1 0.1 ± 0.2 0.6 0.5 ± 0.1 13.5 11.8 ± 5.3
5.5 0.0 0.4 5.9

4.8 0.0 0.3 5.1
500 3.6 1420 3.3 2.9 ± 2.1 0.2 0.1 ± 0.1 0.3 0.4 ± 0.2 3.8 3.4 ± 2.0
0.6 0.0 0.6 1.2

1.8 0.3 0.4 2.5
600 4.3 1783 1.2 1.3 ± 0.5 0.1 0.2 ± 0.1 0.2 0.3 ± 0.1 1.5 1.8 ± 0.6
0.9 0.1 0.3 1.3
*
System geometry similar to the DOPC (2005) and DOPC (2009) systems, as described in ‘Materials and
Methods’.

3.4.3 Pore Creation – Old DOPC Models versus New DOPC Models

Even though DOPC (2005) and DOPC (2009) have somewhat different charge
distributions, the differences in the electric field profile across the membrane interface for these
two systems do not appreciably affect electropore creation dynamics. Pore creation times (and
the times for the initiation, construction, and maturation stages) are similar, as shown in Tables
3.3 and 3.4 and in Figure 3.4. The common feature in the two DOPC systems (and the
corresponding POPC system) is the magnitude of the electric field in the membrane interior and
the resulting net gradient across the interface.
66

3.5 Discussion

We have presented a scheme for dividing the electroporation of lipid bilayers, as
represented in molecular dynamics simulations, into discrete stages, using the evolving
distributions of water and phospholipid head groups (phosphorus atoms) as indicators of the
stage boundaries. This objective model for the life cycle of an electropore facilitates comparisons
of bilayer systems of different compositions, subjected to different conditions, and it is our
expectation that it will contribute, along with other efforts along these lines [16], to the
Figure 3.4 Log (pore creation time) versus membrane internal electric field for POPC (▲),
DOPC (2005) (■), and DOPC (2009) (♦). Error bars (one standard deviation) are shown for
both the internal electric field values and pore creation times. Values for POPC were taken
from Table 3.5. These systems have dimensions comparable to those in the DOPC systems.
67

development of a coherent mechanism for membrane electroporation that is consistent with
molecular and continuum models, with experimental systems of artificial bilayers, and with
living cell membranes.
Our results show that pore creation and annihilation times are non-linearly dependent on
the porating electric field, and that some stages of electroporation are strongly field-dependent
while others are insensitive to the magnitude of the field. This molecular dynamics model of
electroporation is at least broadly consistent with the stochastic pore hypothesis for
electropermeabilization [117, 133, 167, 199], and that model may now be enhanced by the
incorporation of multiple, field-dependent and field-independent, steps in pore creation and
annihilation.

3.5.1 Stages of Electroporation

Pore initiation, the formation of a water column across the membrane interior, is the
major field-dependent step in pore creation. Pores initiate more quickly in higher electric fields.
This is not surprising in light of previous work showing that the electric field gradient in the
bilayers interface region drives water dipoles into the membrane [175]. Pore construction and
maturation, the stages following initiation, are not strongly dependent on the magnitude of the
porating electric field over the range investigated here, indicating that the reorganization and
inward migration of lipid headgroups during pore formation is driven at least as much by
intermolecular forces promoting association with the water in the bilayer interior as by the
applied electric field. If the external field is removed immediately after pore initiation, however,
pore construction does not proceed, and the bilayer-bridging water column collapses.
68

Pore annihilation is broadly analogous to pore creation, although it is by no means simply
pore creation run in reverse. During destabilization, the first step in annihilation, we observe
shrinkage but no other obvious changes in the gross structure of the pore. Water and head group
dipoles are relaxing during this phase, and without the electric field there is no energy to
maintain the intermolecular arrangements associated with the hydrophilic pore [209]. Although
more data is needed, our results show that the pore destabilization time (and the overall pore
annihilation time) decreases slightly at the highest field tested. We hypothesize that these “high-
field pores” are less ordered than pores formed at lower fields, and thus easier to disassemble,
since the water and phospholipid head group components will not have settled so far into local
energy minima. Still, this correlation is small, and stochastic effects may be dominating our
results. Additionally, the structure of the pore, independent of how it was created, may be the
dependent variable which affects pore annihilation, and in this case we are observing the effects
of modulating the pore-initiating electric field on the pore structure itself.
A measure of the stability of a pore, once formed, may be taken from the observation that
in these simulations pore deconstruction on average takes about 75 times longer than pore
construction. Pore deconstruction and dissolution (like construction and maturation) show little
field dependence. Overall, pore annihilation appears to consist of both deterministic and
stochastic elements, with the primarily stochastic stages accounting for most of the annihilation
time.

69

3.5.2 DOPC versus POPC

Pore creation times for the POPC and DOPC systems reported here are similar when
simulations with similar internal electric fields are compared. Because electropore formation
involves primarily a reorganization of water and head groups at the interface, it is not surprising
that the saturation or unsaturation of the second fatty acid residue does not strongly affect the
pore creation time of a homogeneous bilayer, despite other differences in physical properties
[178]. It may be, however, that in more complex systems containing multiple phospholipids and
other components like cholesterol [88], sphingolipids, and cytoskeletal attachments [142], the
extent of saturation in the phospholipids may be significant.
Pore creation times for POPC and DOPC are inversely dependent on the membrane
internal electric field, and, regardless of the lipid type, systems with similar internal electric
fields have similar pore creation times, suggesting that it may be useful in molecular dynamics
studies to regard the electric field differential across the membrane interface as a fundamental
determinant of the poration process. We note that the y-intercept in Figure 3.4 can be interpreted
as the projected pore creation time when the electric field is zero, a value which is also
represented in continuum models for pore creation and which is related to the basal permeability
of artificial lipid bilayers and cell membranes. This kind of intersection may contribute to the
bridging of gaps between atomically detailed simulations, mathematical and physical models,
and experiment.

70

3.5.3 Deficiencies in the Electropore Life Cycle Scheme

The electropore life cycle model presented here has at least three significant deficiencies.
First, although our definition of a “mature” pore is convenient and not altogether arbitrary, it is
not based on the physical properties of pores in real systems, because these are not known at this
time at the molecular level. We would like to define a “mature” pore that corresponds to the
quasi-stable structures that form in artificial lipid bilayers and living cell membranes exposed to
electric fields, but at present no such definition exists [173]. We hope to investigate this
unexplored territory through simulations in which a balance is struck between pore diameter,
electrical conductance, and applied electric field.
Second, the simulation times are very short, much less than a microsecond. Although this
molecular dynamics regime is well suited for comparisons to data from experiments using
nanosecond pulsed electric fields, it is possible that the much longer rise times, lower fields, and
longer pulses (microseconds to milliseconds) associated with conventional electroporation
technology result in a different energy landscape for pore formation, and porated membrane
structures that are significantly different from those we observe in the systems reported here.
More powerful computing hardware and smarter simulation software combined with well-
designed experiments will tell us whether this is so.
Third, an important stage may be missing. It is possible that what we have called pore
initiation actually includes two steps, somewhat analogous to the first two steps in pore
annihilation. In the first of these two steps, the true initiating step, thermal jostling of water
dipoles and head groups in the interface, driven more along one path than another by the applied
electric field and the membrane internal electric field, results in the pre-pore bumps described
71

previously [209]. A hypothetical structure which we have not yet identified develops in some
small percentage of these bumps, either as the bump is formed or as a consequence of the
continued interaction of the bump with the electric field. This structure represents the low energy
path to formation of the membrane-spanning water column, a clear signature of commitment to
pore formation. If there is a missing stage, it would represent the conversion of a pre-pore bump
to a structure from which water is launched into the membrane interior — a pore pedestal. This
pedestal structure, at least the water components of it, may have features in common with the
pyramidal or conical features associated with the electric field-enhanced evaporation of water
[122].
It is also possible that the missing stage is not really a separate step at all, but rather a
statistical selection of a pore-promoting structure from all of the randomly formed pre-pore
bump configurations. Our efforts to identify a signature for this putative pore pedestal have been
unsuccessful to date [209]. It has been proposed that lipid protrusions, dragging water into the
bilayer interior and promoting a dielectric avalanche, are the primary actors in electropore
formation [16]. But we observe many lipid protrusions (bumps of water and lipid) in our
simulations. Only very few of them result in pore formation, and we have been unable to detect
any head group dipole configuration or density fluctuation associated with the lipids in these
bumps that correlates with pore formation [193]. Although cooperative associations of water and
lipid head groups undoubtedly play an important role in pore formation, the fact that we observe
electroporation in simulations of phospholipid bilayers with chargeless head groups [193] and
even in octane systems [176] leads us to believe that rather than look at lipids we should
cherchez l’eau.
72

3.6 Acknowledgements

We thank D. Peter Tieleman for stimulating discussions and insightful input. Computing
resources were provided by the USC Center for High Performance Computing and
Communications. This work was made possible in part by the Air Force Office of Scientific
Research and by MOSIS, Information Sciences Institute, Viterbi School of Engineering,
University of Southern California.

73

Chapter 4: Calcium and Phosphatidylserine
affect Pore Life Cycle

Portions of this chapter are adapted from:
Z. A. Levine and P. T. Vernier, "Calcium and phosphatidylserine inhibit lipid electropore
formation and reduce pore lifetime," J Membr Biol 245 (10), 599-610 (2012).

4.1 Abstract

Molecular dynamics simulations of electroporation of homogeneous phospholipid bilayers show
that the pore creation time is strongly dependent on the magnitude of the applied electric field.
Here we investigate whether heterogeneous bilayers containing phospholipids with zwitterionic
and anionic headgroups exhibit a similar dependence. To facilitate this analysis we divide the life
cycle of an electropore into several stages marking the sequence of steps for pore creation and
pore annihilation (restoration of the bilayer after removal of the electric field). We also report
simulations of calcium binding isotherms and the effects of calcium ions on the electroporation
of heterogeneous lipid bilayers. Calcium binding simulations are consistent with experimental
data using a 1:2 Langmuir binding isotherm. We find that calcium ions and phosphatidylserine
(PS) increase pore creation time and decrease pore annihilation time. For all systems tested, pore
creation time was inversely proportional to the bilayer internal electric field.

74

4.2 Introduction

Lipid membranes play many fundamental roles in cell biology, one of which is the
partitioning of interior cellular components from the outside world. This barrier function can be
disrupted with the application of a sufficiently high external electric field, which permeabilizes
the membrane, a process referred to as electropermeabilization, or electroporation [60, 141].
Within a few nanoseconds of applying such a field, increases in membrane electrical
conductance can be detected [11], and normally impermeant molecules which were previously
excluded from the cell are able to penetrate the membrane [108, 119, 139]. Fluorescent dyes
which interact only with intracellular material can be used as indicators of the extent of
permeabilization, but the mechanisms and physical structures associated with
electropermeabilization are far from being completely understood and are not easily accessible
by experiment [173].
Furthermore many biological processes are significantly affected by small changes in the
local energy landscape, as is the case when an asymmetric distribution of charge is present, thus
a mechanistic understanding of electropermeabilization cannot be restricted to simple,
homogeneous lipid bilayer systems. It becomes necessary then to take into account also the
effects of inorganic ions which may, for instance, cause membrane lipids to aggregate [18, 209],
possibly modifying local surface tensions. Similarly the behavior of heterogeneous phospholipid
bilayers, which may have not only zwitterionic, but also anionic and cationic headgroups, must
be considered since these additions further perturb the electrochemical landscape and greatly
increase the required complexity of analytical models.
75

Recent advances in high-performance computing have enhanced studies of
electropermeabilization using molecular dynamics (MD) simulations of phospholipid bilayers.
These simulations suggest that electropermeabilization results at least in part from the formation
of discrete, nanoscale electropores which develop in a characteristic sequence: (1) appearance of
an electric field-driven water column across the bilayer interior, a process dominated by the
energy minimization of interfacial water molecules under an applied electric field [175, 209]; (2)
construction of a bridge of hydrophilic lipid head groups and additional hydrating water
molecules; (3) expansion of the pore while the external field continues to be applied [16, 171,
175, 176]. Previously we proposed definitions for these individual stages in pore formation,
which we called pore initiation, pore construction, and pore maturation [93]. In this paper we
refine the nomenclature for these steps.
We have also observed that homogeneous, zwitterionic (containing phosphatidylcholine
(PC)) bilayers with varying degrees of hydrocarbon saturation require different ‘minimum’
external electric fields to form pores within a finite period of time [209] which indicates, not
surprisingly, that pore formation is affected by lipid properties. Because anionic lipids such as
phosphatidylserine (PS) can form complexes in the presence of calcium ions [194], it is
reasonable to expect that the incorporation of PS into phospholipid bilayers and the inclusion of
Ca
2+
in the system will also affect pore formation and annihilation. We present here a description
of the electropore life cycles of PC:PS bilayers with and without calcium, extending the analysis
we developed for homogeneous PC bilayers [93].
MD simulations of electropermeabilization must be verified by alignment with existing
continuum theories [200] and experiments [160]. For validating MD representations of calcium
binding to phospholipid bilayers, stoichiometric or coordination complex measurements have
76

been used as a metric [17, 134, 194]. Here we propose using binding isotherms as an additional
metric which describes the concentration of adsorbed ions at an interface relative to the total ion
concentration. At low calcium concentrations one would expect to see linear 1:1 binding,
implying that every added calcium ion binds to the interface. At higher calcium concentrations
however one expects to observe a Langmuir binding isotherm which is characterized by a
transition from linear ion binding to a fixed amount of ion binding when the interface becomes
saturated, as others have reported [17, 160]. We compare our simulated calcium binding
isotherms with existing experimental and theoretical binding constants for PC:PS bilayers and
Ca
2+
.

4.3 Materials and Methods

4.3.1 Structures and Parameters

All simulations were performed using the GROMACS set of programs version 4.0.5 [64]
on the University of Southern California High Performance Computing and Communications
Linux cluster (http://www.usc.edu/hpcc/). Lipid topologies were derived from OPLS united-atom
parameters [15], and the Simple Point Charge (SPC) water model was implemented [12]. Each
system was 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 [14].
Pressure was coupled semi-isotropically (using a compressibility of 4.5 x 10
-5
bar
-1
) normal to
and in the plane of the membrane. Bond lengths were constrained using the LINCS algorithm
[63] for lipids and SETTLE [109] for water. Short-range electrostatic and Lennard-Jones
interactions were cut off at 1.0 nm. Long-range electrostatics were calculated by the PME
77

algorithm [41] using fast Fourier transforms and conductive boundary conditions. Reciprocal-
space interactions were evaluated on a 0.12 nm grid with fourth order B-spline interpolation.
Periodic boundary conditions were employed to mitigate system size effects.
All systems contain a total of 128 lipids — either 1-palmitoyl-2-oleoyl-sn-glycero-3-
phosphatidylcholine (POPC) or 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylserine (POPS)
— and about 9000 water molecules (~70 waters/lipid), which resulted in an initial system box
size of approximately 7 nm x 7 nm x 10 nm. Homogeneous bilayers consist only of POPC.
Heterogeneous bilayers were obtained by replacing 20 POPC molecules on a single leaflet with
20 (anionic) POPS molecules and 20 sodium counter ions in the bulk water, followed by
equilibration until the total area per lipid became constant. Simulations with multiple trials were
run in parallel, starting from the same initial positions. To ensure that each trial was independent,
every atom was assigned a randomized velocity with a Maxwell distribution at the beginning of
the simulation (a built-in function of GROMACS). All systems were equilibrated for a constant
area per lipid. Pore creation times were measured for three values of applied external electric
field — 400 MV/m, 500 MV/m, and 600 MV/m — with three independent trials at each electric
field strength. Note that 600 MV/m in vacuum corresponds to an effective electric field in
aqueous media of about 8 MV/m due to their relative dielectric permittivities. Following pore
formation the external field was removed from a randomly selected simulation, and three
independent trials were run to track the annihilation of that single pore. Additionally, we
extracted calcium binding coefficients and pore creation times in the presence of calcium from
systems where the GROMACS function ‘genion’ was used to replace bulk water molecules with
one calcium ion and two chloride counter ions. Calcium binding coefficients were extracted after
150 ns, to allow the system to reach equilibrium after calcium addition. Binding coefficients
78

were determined using a 1:2 Langmuir isotherm. Pore radius measurements were obtained by
first identifying a pore axis which passes directly through the pore in the z-dimension (based on
perpendicular x and y water density profiles). Then, bins were assigned between mean lipid
phosphorus planes where, for each bin, maximally-distant water molecules were identified
relative to the pore axis in x and y to obtain local semi-major and semi-minor pore diameters.
Following this the semi-major and semi-minor pore diameters were averaged together in each
bin, and finally all local pore diameters were averaged together to obtain an average pore
diameter or pore radius. Pore creation times for systems with calcium present were observed
after an electric field was applied to the equilibrated system.

4.3.2 Updates to the Electropore Life Cycle

All simulations The life cycle of an electropore can be divided broadly into a pore
creation step and a pore annihilation step, as described previously [93]. Pore creation consists of
three stages: initiation, construction, and expansion (defined below). Pore annihilation begins
when the external electric field is removed from an expanded pore and proceeds through settling,
stabilization, deconstruction, and dissolution.
Pore creation. Pore initiation begins with the application of an external electric field and
ends when the two groups of water molecules initially separated by the bilayer merge to become
a single group (non-zero water density in every 0.1 nm slice between the two mean planes of
phosphorus atoms). In some simulations the combined water groups split again in less than 400
ps after the joining. These transient events are not counted as mergers. Pore construction begins
with the formation of the membrane-spanning water column (hydrophobic pore [1]) that marks
the merger of the water groups and ends when the phosphorus groups that are initially found on
79

the two leaflets of the bilayer follow the water into the membrane interior and merge into a
single phosphorus group. (A phosphorus group is defined as a set of atoms, each separated by a
maximum distance of 1.2 nm [159].) Because water and the charged phospholipid head groups
now bridge the membrane interior, this structure is comparable to what is sometimes called a
hydrophilic pore [52, 92, 201]. Continued application of the porating electric field results in an
evolution, or expansion, of the hydrophilic pore. We define an expanded pore arbitrarily as a
hydrophilic pore in which at least 10 phosphorus atoms from the initial anodic leaflet are found
within 1.2 nm of phosphorus atoms from the cathodic leaflet.
Pore annihilation. Pore settling, the first stage of pore annihilation, is the quasi-stable
period after the field is removed during which the number of anode-to-cathode phosphorus
connections fluctuates around the expanded pore criterion (10 connections). Pore stabilization
begins when the number of anode-to-cathode phosphorus connections drops below 10 and ends
when there is only one anode-to-cathode phosphorus connection. The pore radius decreases
around this time to a minimum, about 0.4–0.6 nm. We note that these so-called ‘minimal pores’
remain hydrophilic and facilitate the conduction of sodium ions when very small external electric
fields remain, consistent with theoretical models such as the asymptotic model of electroporation
[117] which also identify a minimum hydrophilic pore radius of about 0.5 nm. Pore
deconstruction, another quasi-stable period, terminates when the single phosphorus group of the
porated bilayer splits into two groups which remain separate for the remainder of the simulation.
At the end of pore deconstruction only the water column remains. We call the disassembly of the
water column pore dissolution.
Individual pore life cycle times were calculated using a custom Perl program which
codifies the stage boundaries described above. The average value of the electric field at the
80

bilayer mid-plane (membrane internal electric field) was extracted using the GROMACS
function ‘g_potential’ over the period from 100 ps after the external electric field was applied to
100 ps before pore construction began. As in our previous work, internal electric field was used
as a normalizing term [93]. Molecular graphics images were generated with Visual Molecular
Dynamics (VMD) [75].

81

4.4 Results

A complete summary of the pore life cycle results reported in this study can be found in
Figures 4.1 and 4.2 and in Tables 4.1-4.8. Details are described below.

Figure 4.1 Pore creation times for three different porating electric fields (400, 500, 600
MV/m) for bilayers consisting of 128 POPC (0PS:0Ca), 128 POPC saturated with calcium
(0PS:100Ca), 108 POPC and 20 POPS on a single leaflet without calcium present
(20PS:0Ca), and bilayers containing both PS and Ca
2+
(20PS:100Ca). Systems containing 20
PS also contain 20 Na
+
as counter ions, and systems containing Ca
2+
contain two chloride
counter ions for every calcium ion. Systems which contain both PS and Ca
2+
have both
sodium and chloride counter ions present. Ca
2+
and POPS in the bilayer increase the pore
initiation time.
82

Figure 4.2 Pore annihilation times after pore formation in the same systems shown in Figure
4.1. Ca
2+
and POPS in the bilayer decrease the pore annihilation time.
83

Table 4.1 POPC pore creation times without PS and without Ca
2+
[93]
Initiation Construction Expansion Creation

Average

2.8 0.0 0.5 3.3
400 1630 2.3 2.0 ± 1.0 0.3 0.2 ± 0.2 1.1 0.8 ± 0.3 3.7 2.9 ± 1.1
0.8 0.2 0.7 1.7

0.9 0.1 1.0 2.0
500 2180 0.6 0.6 ± 0.3 0.2 0.2 ± 0.1 0.5 0.8 ± 0.3 1.3 1.6 ± 0.4
0.4 0.3 0.9 1.6

0.5 0.1 0.1 0.7
600 2347 0.5 0.4 ± 0.2 0.0 0.1 ± 0.1 0.3 0.2 ± 0.1 0.8 0.7 ± 0.1
0.2 0.2 0.2 0.6

Table 4.2 POPC pore creation times without PS and with 100 Ca
2+
Initiation Construction Expansion Creation

Average

6.2 0.1 0.5 6.8
400 1373 4.6 5.0 ± 1.0 0.1 0.2 ± 0.1 0.4 0.4 ± 0.2 5.1 5.6 ± 1.1
4.3 0.3 0.2 4.8

0.9 0.1 0.2 1.2
500 1837 0.8 0.8 ± 0.2 0.1 0.1 ± 0.0 0.2 0.3 ± 0.1 1.1 1.1 ± 0.1
0.6 0.1 0.4 1.1

0.4 0.2 0.2 0.8
600 2292 0.4 0.4 ± 0.0 0.1 0.1 ± 0.1 0.2 0.2 ± 0.0 0.7 0.7 ± 0.1
0.4 0.1 0.2 0.7

84

Table 4.3 POPC pore creation times with PS and without Ca
2+

Initiation Construction Expansion Creation

Average

2.2 0.3 0.4 2.9
400 1707 3.2 3.1 ± 0.9 0.1 0.2 ± 0.1 0.4 0.4 ± 0.0 3.7 3.7 ± 0.8
3.9 0.2 0.4 4.5

0.5 0.2 0.2 0.9
500 1930 1.4 0.9 ± 0.5 0.1 0.1 ± 0.1 0.5 0.4 ± 0.2 2.0 1.4 ± 0.6
0.9 0.1 0.4 1.4

0.2 0.2 0.2 0.6
600 2197 0.3 0.4 ± 0.3 0.1 0.1 ± 0.1 0.1 0.1 ± 0.1 0.5 0.7 ± 0.2
0.7 0.1 0.1 0.9

Table 4.4 POPC pore creation times with PS and with 100 Ca
2+

Initiation Construction Expansion Creation

Average

8.9 0.0 0.6 9.5
400 1417 5.1 17.7 ± 18.6 0.1 0.1 ± 0.1 0.9 0.7 ± 0.2 6.0 18.4 ± 18.5
39.0 0.0 0.6 39.6

1.5 0.1 0.3 1.9
500 1833 1.8 1.7 ± 0.1 0.2 0.1 ± 0.1 0.3 0.4 ± 0.1 2.3 2.2 ± 0.2
1.7 0.1 0.5 2.3

1.0 0.4 0.4 1.8
600 2230 0.8 0.9 ± 0.1 0.0 0.2 ± 0.2 0.1 0.3 ± 0.2 0.9 1.3 ± 0.5
0.8 0.1 0.3 1.2

85

Table 4.5 POPC pore annihilation times without PS and without Ca
2+
[93]
Settling Stabilization Deconstruction Dissolution Annihilation

7.3 3.0 24.6 2.9 37.8
400 29.2 19.4 ± 11.1 3.7 2.5 ± 1.5 25.9 25.0 ± 0.8 1.1 1.7 ± 1.0 59.9 48.6 ± 11.1
21.6 0.8 24.6 1.1 48.1

84.4 7.5 17.3 1.4 110.6
500 17.3 53.7 ± 33.9 4.0 4.7 ± 2.5 49.8 29.3 ± 17.8 4.2 2.9 ± 1.4 75.3 90.7 ± 18.1
59.5 2.7 20.9 3.0 86.1

11.4 1.0 8.1 1.5 22.0
600 3.8 7.1 ± 3.9 6.4 5.4 ± 4.0 28.5 25.1 ± 15.5 0.7 0.9 ± 0.6 39.4 38.4 ± 16.0
6.2 8.7 38.6 0.4 53.9

Table 4.6 POPC pore annihilation times without PS and with 100 Ca
2+

Settling Stabilization Deconstruction Dissolution Annihilation

0.0 0.4 4.6 1.0 6.0
400 0.0 0.6 ± 1.0 2.3 1.4 ± 1.0 3.3 3.3 ± 1.3 0.7 1.0 ± 0.3 6.3 6.3 ± 0.3
1.7 1.6 2.1 1.2 6.6

0.1 0.1 1.2 0.5 1.9
500 0.0 0.1 ± 0.1 0.1 0.1 ± 0.0 0.1 0.5 ± 0.6 0.2 0.9 ± 0.9 0.4 1.5 ± 0.9
0.0 0.1 0.1 1.9 2.1

0.0 0.5 2.3 2.7 5.5
600 9.0 3.0 ± 5.2 5.4 2.4 ± 2.6 28.8 10.6 ± 15.8 4.4 2.4 ± 2.1 47.
6
18.5 ± 25.3
0.0 1.3 0.8 0.2 2.3

86

Table 4.7 POPC pore annihilation times with PS and without Ca
2+

Settling Stabilization Deconstruction Dissolution Annihilation

0.8 2.8 4.9 0.1 8.6
400 0.3 0.6 ± 0.3 3.7 2.7 ± 1.1 16.7 11.6 ± 6.0 0.7 0.3 ± 0.3 21.4 15.2 ± 6.4
0.8 1.6 13.1 0.1 15.6

3.2 0.6 2.4 3.1 9.3
500 0.9 1.5 ± 1.5 1.8 6.1 ± 8.5 8.7 14.2 ± 15.3 0.1 1.6 ± 1.5 11.5 23.4 ± 22.5
0.3 15.9 31.4 1.7 49.3

1.1 1.5 13.8 3.6 20.0
600 1.0 1.7 ± 1.1 1.2 2.3 ± 1.7 7.2 8.3 ± 5.1 9.5 4.7 ± 4.4 18.9 16.9 ± 4.5
2.9 4.2 3.8 0.9 11.8

Table 4.8 POPC pore annihilation times with PS and with 100 Ca
2+

Settling Stabilization Deconstructio

Dissolution Annihilation

0.0 7.1 2.5 0.2 9.8
400 0.1 0.1 ± 0.1 2.3 4.7 ± 2.4 2.0 1.7 ± 1.0 0.6 0.8 ± 0.7 5.0 7.2 ± 2.4
0.0 4.8 0.5 1.5 6.8

0.0 1.7 1.5 1.6 4.8
500 0.1 0.1 ± 0.1 1.9 2.2 ± 0.8 0.4 2.3 ± 2.4 1.0 1.5 ± 0.4 3.4 6.1 ± 3.5
0.1 3.1 5.0 1.8 10.
0

0.1 4.5 1.2 3.7 9.5
600 0.2 0.1 ± 0.1 3.0 2.8 ± 1.8 0.2 0.8 ± 0.5 1.1 1.9 ± 1.6 4.5 5.6 ± 3.5
0.0 0.9 0.9 0.9 2.7

87

4.4.1 Pore Creation – POPC and POPC:POPS

At smaller external electric fields, POPC:POPS bilayers have pore creation times that are,
on average, slightly longer than homogeneous POPC bilayers (Figure 4.1), primarily because of
an increase in the pore initiation time, the time it takes water to bridge the membrane interior.
The area per lipid — which is strongly correlated with membrane permeability [103] — in
POPC:POPS systems was about 0.60 nm
2
compared to 0.66 nm
2
for POPC bilayers. At higher
applied electric fields the differences in pore creation time become minimal, and mixed bilayers
have pore creation times which are not significantly different from those for homogenous
bilayers. Pore initiation times for POPC:POPS are inversely related to the externally applied
electric field, and pore construction times remain constant over all electric fields sampled,
similar to what we previously reported for homogeneous POPC bilayers [93]. Pore expansion
times decrease slightly as higher external electric fields are applied.

4.4.2 Pore Creation – POPC and Calcium
POPC bilayer systems containing calcium ions have, at small external electric fields, pore
creation times which are about twice as long as systems without calcium, and about one and a
half times as long as those in mixed (POPC:POPS) bilayers (Figure 4.1). As with systems
containing PS, when the external electric field is increased, pore creation times of pure POPC
systems with calcium and those without calcium are not significantly different. The convergent
area per lipid of POPC bilayers with calcium is about 0.56 nm
2
, a value significantly smaller than
the area per lipid without calcium reported above. Again as with systems containing PS, we
observe an inverse relationship between externally applied electric field and pore initiation time
88

and pore expansion time for POPC bilayer systems containing calcium. Pore construction times
remain constant at all electric fields applied.

4.4.3 Pore Creation –Calcium and POPC:POPS
POPC:POPS bilayers containing calcium have pore creation times which are significantly
longer than those for pure POPC systems, with or without calcium, and also longer than
POPC:POPS systems without calcium (Figure 4.1). This is true for all electric fields applied,
even very high fields where for all of the systems described above the pore creation times were
not significantly different. The average area per lipid of POPC:POPS bilayers with calcium was
about 0.55 nm
2
before a field was applied, slightly smaller than the area per lipid of POPC
bilayers with calcium. As with all other systems reported here, both pore initiation times and
pore expansion times decrease as the external electric field increases, while pore construction
times are similar for all values of electric field applied in these simulations.

4.4.4 Pore Annihilation – POPC and POPC:POPS
POPC:POPS bilayers have pore annihilation times about three times smaller than pure
POPC bilayers (Figure 4.2). All pore annihilation stages except stabilization are shortened
relative to pure POPC, and this is true regardless of the value of the electric field used to create
the pore. The average pore radius at the start of the annihilation step was about 2.3 nm for
POPC:POPS systems, similar to the average pore radius for pure PC bilayers, 2.2 nm. Also, pure
POPC bilayer pore annihilation times are dominated by pore settling and pore deconstruction.
89

POPC:POPS bilayers exhibit very short pore settling times and significantly reduced pore
deconstruction times compared to POPC systems, reducing the overall time required to
annihilate POPC:POPS electropores. An initial decrease in pore radius occurs immediately after
the field is removed in both POPC and POPC:POPS bilayers, from about 2.3 nm to about 0.5 nm
over the first few nanoseconds. This initial reduction of the pore radius does not appear to be
correlated with pore settling time (Figure 4.3).

Figure 4.3 Evolution of pore radius after removal of the porating electric field for the same
systems shown in Figure 4.1. This example, from a simulation in which E
porating
= 500 MV/m,
is representative of pore annihilation behavior for all values of E
porating
examined in this work.
Average initial pore radii are 2.2 nm for 0PS:0Ca systems, 2.4 nm for 20PS:0Ca systems, 1.7
nm for 0PS:100Ca systems, and 1.9 nm for 20PS:100Ca systems.
90

4.4.5 Pore Annihilation – POPC and Calcium
Pure POPC systems containing calcium have dramatically reduced pore annihilation
times compared to POPC systems without calcium (Figure 4.2). For pores created at 400 MV/m
and 500 MV/m, POPC–Ca
2+
systems exhibit virtually no pore settling, while all remaining stages
take no more than a few nanoseconds to complete. For pores created at 600 MV/m, calcium still
significantly reduces the pore annihilation time, but the variation from simulation to simulation is
large. As with POPC:POPS and pure POPC bilayers without calcium, the pore radius is reduced
in POPC systems containing calcium to 0.5 nm after only a few nanoseconds, but with calcium
present the initial pore radius immediately after the external field is removed is 1.7 nm, about 0.6
nm smaller than the pores in the annihilation simulations without calcium.

4.4.6 Pore Annihilation – Calcium and POPC:POPS
POPC:POPS systems with calcium have pore annihilation times which are also
significantly shorter than those for pure POPC bilayers and POPC:POPS systems without
calcium. POPC:POPS–Ca
2+
annihilation times are comparable to those for POPC–Ca
2+
and
POPC:POPS (no calcium) systems. For pores created at the highest electric field, POPC:POPS–
Ca
2+
systems have the shortest pore annihilation times of all systems sampled. Pore settling time
for these systems is too short to measure in about half of the trials; the longest settling time
measured for POPC:POPS–Ca
2+
is 200 ps. Stabilization is the dominant step in pore annihilation
for these systems, and the variance in pore annihilation time is much smaller than for pure POPC
systems without calcium. As with POPC–Ca
2+
systems, the initial pore radius for POPC:POPS–
Ca
2+
systems after the porating field is removed is about 1.8 nm. As soon as the pore radius
91

decreases to about 0.5 nm in the POPC:POPS–Ca
2+
systems, the pores dissipate quickly, in
contrast with pure POPC bilayers without calcium, where pores remain open with a radius
around 0.5 nm for many tens of nanoseconds.

4.4.7 Calcium Binding Isotherms
To assess the validity of our calcium ion models in POPC:POPS systems, we constructed
a binding isotherm (Figure 4.4), which plots the amount of bound calcium ions against the bulk
calcium ion concentration. Our data can be described by a 1:2 Langmuir binding isotherm
between calcium and phospholipid [3] and with binding isotherms taken from experiments with
pure and mixed vesicles [160]. Additionally, for small calcium concentrations we see linear 1:1
binding. Calcium appears to bind preferentially to PS carboxyl and PC:PS phosphoryl oxygens.
From systems equilibrated for 150 ns [194] we extracted a calcium binding coefficient K = 2.56
M
-1
.

92

4.5 Discussion

4.5.1 Electropore Life Cycle for POPC:POPS Bilayers

We have shown that the incorporation of the anionic phospholipid POPS into a POPC
(zwitterionic) bilayer slightly increases the time required for pore creation at lower external
electric fields and drastically decreases the time required for pore annihilation for electropores
created at all external electric fields. PS bilayers have a smaller area per lipid [111] compared to
PC bilayers, and changes in surface tension, which may be affected by locally varying area per
lipid, have been shown to affect pore formation [95, 175, 176].
Figure 4.4 Calcium binding curves show rough correspondence between experimental and
simulated results. The data can be described by a 1:2 Langmuir binding isotherm, consistent
with formation of Ca:PS
2
complexes.
93

At higher electric fields pore creation times for POPC:POPS bilayers are similar to those
for pure POPC bilayers. As previous studies have shown, interfacial water is a dominant
component of pore formation [171, 175, 209], and we can speculate that the formation of bilayer-
spanning water bridges in the strong interfacial electric field gradients that result from the
application of large external electric fields is only weakly influenced by PS or PC interactions.
This is consistent with the stochastic pore hypothesis for electropermeabilization [117, 133, 167,
199] and previous observations [209], in which pore creation time is inversely correlated with
the magnitude of the applied electric field. In this scheme our highest external electric field —
600 MV/m — is at or near a saturating value for porating fields for POPC and POPC:POPS
bilayers. In this saturation range electropores are created in similar, asymptotically convergent
times. We also observe comparable pore radii for POPC:POPS and pure POPC bilayers at the
end of the pore creation step, despite small differences in area per lipid, suggesting that the
radius of an electropore is only partially dependent on phospholipid areal density.

4.5.2 Effects of Calcium on Pore Life Cycle

Calcium in pure POPC bilayer systems delays pore creation at lower external electric fields
and greatly reduces the time required for pore annihilation. POPC:POPS systems containing
calcium have even longer pore creation times at lower electric fields. The pore radius after
creation is lower in both POPC and POPC:POPS systems containing calcium systems than for
systems without calcium. As indicated above, this may be partially associated with significantly
smaller area per lipid values found both experimentally [104] and in simulations [17] when
calcium is present.
94

Even though the area per lipid is similar for POPC–Ca
2+
and POPC:POPS–Ca
2+

systems, it has been shown experimentally and in MD simulations that calcium binds with more
affinity to PC:PS vesicles than to pure PC vesicles [160, 194], forming Ca
2+
–PS complexes
involving PS carboxyl and phosphoryl oxygens. By increasing the total number of calcium-lipid
complexes, we are effectively decreasing the area per lipid of our system and thus increasing
membrane surface tension. Until the molecular-level details of water intrusion and bridge
construction across the bilayer have been fully explained, we can only presume that increases in
surface tension associated with the presence of Ca
2+
and PS lead to increases in the energy
barrier for interfacial water entry into the bilayer interior, and an associated increase in pore
creation times. Alternatively when the external field is removed, a change in surface tension
would result in a change in the pore edge tension since the two quantities are related [146], thus
we would expect to see modified pore annihilation times for systems with PS and/or Ca
2+

present, compared to systems without PS and Ca
2+
.”
In addition, POPC:POPS–Ca
2+
systems exhibit very short pore annihilation times (a few
nanoseconds), similar to the time scale of pore creation. This is an order of magnitude faster than
pure POPC systems, which have pore annihilation times that are tens or even more than 100
nanoseconds in simulations [93]. These timescales are similar to the resealing timescales
reported in other molecular dynamics simulations [171], however a large discrepancy still exists
between simulated pore annihilation times and experimental resealing times which occur on
timescales of milliseconds [106] to hundreds of seconds [83]. These long-lasting pores occur in
living cell membranes, not simple phospholipid bilayers, suggesting that some membrane
restructuring beyond lipid nanopore formation occurs, and that the permeabilizing structures in
cell membranes have significantly stabilizing features.
95

4.5.3 Molecular and Continuum Models of Lipid Electropores

Our results agree (Figure 4.5) with numerical models which predict an exponential
relationship between transmembrane potential and pore creation [117] for all configurations
tested. (Note that there is a direct relationship between the internal electric field and the
transmembrane potential.) The electric field is the dominant term for pore creation, and although
we see variations in pore creation time at smaller fields for POPC versus POPC:POPS or with
the introduction of calcium ions, pore creation occurs at similar times for very high electric
fields. Pore annihilation is a stochastic process in some models [200], and indeed we see large
variances in the pore annihilation time, which suggests that local fluctuations and dynamics
strongly influence restructuring of the bilayer.

Figure 4.5 Lipid electropore creation time as a function of bilayer internal electric field for
the systems shown in Figure 4.1. The log scale reflects the exponential dependence of the
probability of pore formation on the applied electric field. Internal electric field is used here
as a normalizing quantity.
96

Our simulations of calcium binding to POPC and POPC:POPS bilayers provide support
for the validity of the calcium ion model we are using. Our binding isotherm for POPC:POPS
bilayers is similar to experimental data for dioleoyl-phosphatidycholine:phosphosphatidylserine
(DOPC:DOPS) vesicles [160] using a 1:2 Langmuir isotherm. The binding coefficient extracted
from our simulations is 2.56 M
-1
, in the range of the measured value of 13.8 M
-1
[3], although
experiments at different NaCl and Ca
2+
concentrations are difficult to compare directly. Better
models may lead to increased accuracy in determining proper binding coefficients.
Although our simple, largely intuitive characterization of the stages in the lipid
electropore life cycle provides a useful scheme for analysis of pore creation and pore
annihilation, a more sophisticated approach will build on this foundation to incorporate
systematic measurements of the pore radius and pore energies [117], and to include the atomic-
scale electric field landscape and the interactions of water oxygen and hydrogen with the
electron-dense acyl oxygens deep in the phospholipid bilayer interface. Once these key
components of the molecular structure of electropores are accurately represented, the additional
complexities of lipid heterogeneity, membrane proteins, and cytoskeletal and glycocalyx
attachments can be added to the model one by one, until we approach a useful representation of
the living cell membrane in a porating electric field.

4.6 Acknowledgements

We thank Rumiana Dimova for stimulating discussions and insightful input on calcium
binding. Computing resources were provided by the USC Center for High Performance
Computing and Communications (http://www.usc.edu/hpcc/). This work was made possible in
97

part by the Air Force Office of Scientific Research and by MOSIS, Information Sciences
Institute, Viterbi School of Engineering, University of Southern California.

98

Chapter 5: The Role of Water in
Electroporation

Portions of this chapter are adapted from:
M. Tokman, J. H. Lee, Z. A. Levine, M. C. Ho, M. E. Colvin, and P. T. Vernier, "Electric Field-
Driven Water Dipoles: Nanoscale Architecture of Electroporation," PLoS One 8 (4) (2013).
and
P. T. Vernier, Z. A. Levine, and M. A. Gundersen, "Water Bridges in Electropermeabilized
Phospholipid Bilayers," Proceedings of the IEEE 101 (2), 494-504 (2013).
and
S. Romeo, Y. Wu, Z.A. Levine, M. A. Gundersen, P. T. Vernier, "Water influx and cell swelling
after nanosecond electropermeabilization," Biochimica et Biophysica Acta (BBA) -
Biomembranes 1828 (8), 1715-1722 (2013)

5.1 Abstract

Pulsed electric-field permeabilization of living cell membranes forms a novel basis for
new biotechnological protocols and therapeutic modalities; however there exists considerable
disagreement as to which dominant physical structures are responsible for facilitating
electroporation. Experimental observations of artificial membranes and whole cells provide
99

evidence that a membrane-spanning, hydrophilic, conductive pore can form in nanoseconds. An
external electric field lowers the energy barrier for this stochastic process, reducing the mean
time for pore formation and increasing the pore areal density. Molecular dynamic simulations
have recently suggested that perhaps interfacial water, or water present at the aqueous:lipid
interface, plays a key role in facilitating electropermeabilization. Here we present comparisons of
simulations with and without a lipid bilayer present to show that the atomic process of
electroporation is driven by the field-induced reorganization of water dipoles at the water-lipid or
water-vacuum interface into energetically-favorable configurations with molecular dipoles
oriented towards the external electric field. Although the contributing role of water in
electroporation has been noted previously, we propose that interfacial water molecules are the
primary players in the process. The role of the lipid layer, to a first-order approximation, is then
reduced to a relatively passive barrier. This new view of electroporation simplifies the study of
the problem and opens up new opportunities in both theoretical modeling of the process and
experimental research to better control or to use it in new, innovative ways. These model systems
then can deepen our understanding of the factors governing pore initiation, construction, and
lifetime, knowledge that will translate to enhanced utilization of this method in biomedicine and
bioengineering.

5.2 Introduction

Electroporation is the breaching of the integrity of the cell membrane that follows the
application of an external electric field of sufficient magnitude and duration. In this process
permeabilizing structures (pores) appear in the membrane, allowing molecular transport across
100

the normally impermeable barrier [119, 138, 211]. Electroporation has a broad range of
applications in biology, biotechnology, and medicine, from drug and gene delivery into cells to
tumor therapy [62, 102, 120, 144]. Despite the wide laboratory use of electroporation, the details
of the effects of electric fields on biological membranes, and particularly the molecular
mechanisms of pore creation in living cells, are not well understood [173]. Our limited
knowledge of this phenomenon causes difficulties in controlling the process in clinical
applications and limits development of new technologies.
Since the first reports of reversible and irreversible modifications of membrane
conductance by electric fields [34, 60, 162], steady progress has been made toward a
phenomenological understanding of the nature of the permeabilized membrane and the processes
that restructure the phospholipid bilayer in an externally applied electric field [1, 29, 210].
Observational studies have been enhanced, and guided by, the development of continuum
electrophysical models of electroporation [35, 117, 133, 167, 184, 199, 200], which in turn have
contributed to experimental designs and perspectives. These models have been corroborated over
the past decade with molecular dynamics (MD) simulations [42, 58, 71, 171, 175, 196], which
have revealed plausible molecular structures and a very nearly complete, but still hypothetical,
physical mechanism for the poration of phospholipid bilayers if not yet for cell membranes.
These MD systems represent only the simplest cases, with one or two or at most three
phospholipid and one or two cation species in bilayers with areas less than 200 nm
2
, and which
contain no membrane proteins or other intra- and extracellular associations.
Permeabilization can be monitored by tracking the transport of normally impermeant
materials across the cell membrane or by measuring changes in the electrical properties of the
membrane, but direct experimental observation is difficult because of the small spatial and fast
101

temporal scales of this process. Theoretical models have been developed to facilitate the
interpretation of experimental data and the understanding of the mechanism of electroporation.
Although continuum models [84, 200] can predict large-scale features of electroporation, they
lack details regarding pore initiation, growth, and decay and contain empirically fitted
parameters. Molecular dynamics (MD) simulations provide access to the microscopic structure
of a membrane and its interaction with the surrounding solvent and ions in atomic detail. Since
we are interested in understanding the molecular mechanism of electropore creation and
evolution, we employ MD to study this problem.
Acknowledging that we are setting aside the complications of more realistic mixtures of
lipids and electrolytes, membrane proteins, cytoskeletal attachments, and mechanical fluctuations
over larger areas (in other words the compositional and structural details of a living cell
membrane), we have proposed from the results of MD simulations an electropore life cycle [93],
a scheme for the step-by-step development (and dissolution) of the electrically conductive
phospholipid bilayer defects that contribute at least in part to what we call a permeabilized
membrane.
Even within the boundaries of our simple systems many questions remain. Why does a
pore form here in the bilayer and not there? Why does pore initiation take this long and not that
long? Why does a pore have this diameter and not that diameter? Why does a pore have this
lifetime and not that lifetime? I present in this chapter some of what we have learned from our
efforts to answer these questions and to determine, in atomic detail, the specific electric field-
driven molecular interactions that lead to the formation of ion-conducting, bilayer-bridging water
defects. A better understanding of the physics of the single, isolated nanopores in our model
systems provide an underpinning for understanding and controlling populations of electric field-
102

generated pores in real micro- (cellular) and macro-scale (cell suspensions and tissues) systems
[173].
Nanoscale permeabilization in models and experiments. Single electric pulses as short as
3 ns can permeabilize cell membranes [11, 190, 191, 197]. Experimental data and MD
simulations are consistent with stochastic pore formation models, in which the electric field-
facilitated construction of water bridges across the low-dielectric membrane interior is followed
by head group migration along the water column and the establishment of a quasi-stable,
nanometer-diameter, hydrophilic pore. In this model the time to initiation of individual pore
formation is variable, but strongly dependent on the applied electric field. Once initial water
intrusion occurs, the hydrophilic pore is fully formed within about 2 ns, even at the lowest
porating field magnitudes studied [93]. Molecular dynamics studies indicate that interfacial water
dipoles are the primary transducers for converting electric field energy into membrane structural
rearrangements [193, 195, 209].
Pore initiation: does water lead or follow? In one proposed mechanism for electropore
initiation, water entry into the bilayer interior is promoted by the tilting of phospholipid head
group dipoles in a porating electric field [71] or by lipid protrusions into the bilayer interior [16],
or a combination of these. In our simulations, with electric field magnitudes close to the
minimum required to produce pores within a few tens of nanoseconds, very few lipid protrusion
events (temporary bumps of water and lipid) result in pore formation, and in those cases the lipid
head groups follow water into the membrane. In addition, we have been unable to detect any
head group dipole configuration or lipid density fluctuation in these bumps that correlates with
pore formation [193]. Although cooperative associations of water and lipid head groups
unquestionably play an important role in pore formation, the fact that electroporation is observed
103

in simulations of phospholipid bilayers with chargeless head groups [193], in water/octane/water
systems [176], and even in water-vacuum systems, supports the view that although there are
many variations in the molecular configurations that precede the appearance of membrane-
spanning defects, the primary actor in electropore formation is water.
At higher electric fields, like those used in [16], which are not likely to be attained under
the conditions of most electroporation protocols, we observe larger fluctuations in the structure
of the phospholipid-water interface. These deformations, which can involve deep intrusions of
phospholipid head groups into the bilayer interior, become increasingly dominant in the
permeabilizing mechanism when the magnitude of the electric field exceeds the minimum
porating value.
The molecular details of water defect formation in lipid bilayers have not been critically
characterized. One possibility is that a thermal jostling of water dipoles and head groups in the
interface, driven more along one path than another by the applied electric field and associated
interfacial field gradients, results in the pre-pore bumps described previously [193]. An
unidentified hypothetical water structure then develops in some of these bumps, either as the
defect is formed or as a consequence of the continued interaction of the bump components with
the electric field. This pedestal structure provides a low-energy departure point for the molecule-
by-molecule extension of a chain of water into a membrane-spanning water column. As the
water bridge assembles, the lipid head groups follow, and a hydrophilic pore structure is formed
— a clear signature of commitment to pore formation. The conversion of a pre-pore bump to a
structure from which water is launched into the membrane interior — a pore pedestal — may be
the actual pore-initiating event, if one can be defined. This pedestal structure, at least the water
components of it, may have features in common with the pyramidal or conical features
104

associated with the electric field-enhanced evaporation of water [122]. On the other hand, there
may be no separate pedestal-forming step, but rather a statistical selection of a pore-promoting
structure from all of the randomly formed pre-pore bump configurations. Efforts to identify a
signature for this putative pore pedestal have been unsuccessful in the past [193, 209], but I will
try and describe recent progress below.

5.3 Water-Vacuum Systems

Here I will present molecular dynamics simulations of simple water-vacuum (or more
precisely water-vacuum-water) systems which are analogous to lipid bilayer systems containing
layers of water-lipid-water. Both of these systems contain two high-dielectric slabs of water
separated by a low-dielectric slab of vacuum or lipid in the center. The purpose of comparing
these systems to one another is to see if a relationship exists between water-vacuum
electroevaporation and lipid-bilayer electroporation. If it can be shown that electroporation is a
special-case of water-vacuum electroevaporation which depends primarily on water transport
through a passive barrier then we could potentially reduce the complexity of electroporation
continuum models that utilize both water and lipids as first-order electropore mediators into
models which utilize only the properties of water, thus simplifying the overall picture of cell
membrane electropermeailization.

105

5.3.1 Materials and Methods

Model Systems. In order to clearly illustrate the primary role of water dipoles in
electropore formation we study two configurations using MD: water-vacuum-water (WVW) and
water-lipid-water (WLW).
WLW systems contain 128 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC)
lipids and 4480 water molecules (35 waters/lipid), which result in a system box size of
approximately 6.4 nm x 6.4 nm x 7.2 nm. The two directions tangential to the POPC bilayer are
defined as the 𝑋 and 𝑌 directions, with 𝑍 perpendicular to the plane of the membrane. To ensure
that replicated simulations are independent, each atom was assigned a randomized velocity from
a Maxwell distribution at the beginning of a simulation. POPC systems were equilibrated before
the application of an external electric field by allowing the simulations to proceed until a
constant area per lipid (approximately 0.66 nm
2
) was reached (typically in 10–30 ns).
WVW systems were comprised of 6877 water molecules arranged in two slabs of
thickness 4.2 nm separated by a 2.8 nm vacuum gap. This configuration was constructed by
generating a 7 nm x 7 nm x 7 nm periodic water box with the GROMACS utility ‘genbox’, then
using custom Perl scripts to remove a 2.8 nm slice of water molecules from the center of the
system, followed by 300 ps of equilibration at constant volume. The dimensions of this box and
the gap size have been chosen carefully in order to produce a realistic equilibrated initial system
in which the surface area of the water-vacuum interface is minimized and to ensure that the
magnitude of the electric field in the gap is comparable to that of the WLW simulations. Under
non-periodic conditions, the introduction of a vacuum gap into a constant volume system will
lead to formation of a spherical bubble that minimizes the area of the water-vacuum interface.
106

After equilibration of the initial WVW and WLW systems we impose an electric field in the 𝑍
direction (perpendicular to the membrane surface) and then observed the dynamics of the water
in the presence of a constant field. We performed simulations with external electric field values
between 450 MV/m and 1000 MV/m and a vacuum gap size of 2.8 nm.
Molecular Dynamics Simulations Protocols. All simulations were performed using the
GROMACS set of programs version 4.5.3, as previously described [93]. The Extended Simple
Point Charge (SPC/E) water model [13] was used for all simulations presented here, although we
obtained similar results using SPC [12] and SPC/E flexible [42] water models. Lipids are
parameterized with OPLS headgroups and Berger hydrocarbon tails [15]. Any atoms which are
not explicitly modeled with OPLS or Berger parameters use the native GROMOS force field
built into GROMACS. Lipid topologies were obtained from Peter Tieleman’s group
(http://moose.bio.ucalgary.ca).
All simulations 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 [14]. For lipid systems, pressure was coupled semi-isotropically (using a
compressibility of 4.5 x 10
-5
bar
-1
) normal to and in the plane of the membrane (NPT). No
pressure coupling was used for water-vacuum systems where volume was held constant (NVT).
WLW systems were simulated in the NPT ensemble at 1 bar to maintain a dynamically
controlled area per lipid before pore formation. WVW systems were run in the NVT ensemble to
constrain the box dimensions and water slab separation. Note that despite the different ensemble
used for the WLW system, over the time interval of interest the size of the computational box did
not fluctuate by more than 0.1 nm without external electric field and 0.5 nm in the presence of
electric field.
107

All simulations were performed with a time step of 2 fs and with Bussi’s stochastic
velocity rescaling algorithm [24] as a temperature coupling method. Bond lengths were
constrained using the LINCS algorithm [63] for lipids and for water. All bond lengths were fixed
using constraints after the integration of forces. Following Essmann, et al. [41], the Particle-
Mesh-Ewald (PME) method with tinfoil boundary conditions was used to handle long-range
electrostatic forces and cut-offs were employed for calculating van der Waal’s interactions. In
the simulations presented in the next section all of the relevant cut-off distances were set to 1.4
nm. However, we performed simulations with electrostatic cutoffs set to 1.0 nm, 1.2 nm, 1.4 nm,
1.6 nm, and 1.8 nm and found that the results were robust with respect to this parameter and that
the dynamics of the pore initiation were qualitatively unchanged.
For systems with an external electric field the orientation of polar molecules, such as
water, in the field induces a net dipole on the system which has been reported to cause spurious
dipole orientations [21]. For the WVW system we investigated several alternative treatments of
the long-range electrostatic forces, including PME with non-conducting infinite boundaries, very
long cut-offs (3.4 nm), and the reaction field approach, and found that the process of pore
initiation is similar, albeit at slightly different electric field strengths. After the pore forms and
begins to occupy a significant portion of the computational box, the average dipole saturation
differs for different choices of long-range electrostatics approximations and affects the growth
rate and stability of the pore. These post-initiation stages in pore development require different
models and will have to be addressed in later studies; here we focus on the pore initiation
process.

108

5.3.2 Results and Discussion

Dynamics of the Systems. Here I will concentrate on the analysis of 60 characteristic
replicate simulations of WLW systems and WVW systems at various external electric field
strengths, however similar results were obtained in additional simulations (hundreds in total)
with different thermodynamic ensembles, water models, electrostatic cutoff distances, and
applied electric fields, as described above. The dynamics described in this subsection were
observed in all simulations; I present in detail the results of one representative simulation to
illustrate characteristic behavior. Figure 5.1 shows the time evolution of a WLW configuration
over a few nanoseconds while an external electric field is applied. Similarly Figure 5.2 shows a
similar trajectory for a WVW system with an external electric field present. Both WVW and
WLW systems exhibit similar structural changes before electropore formation occurs. Similar to
previous studies where our group codified the discrete stages of lipid bilayer electroporation a
deformation directed towards the interior of the membrane forms at the water-lipid and water-
vacuum interfaces. This bump grows, and eventually water molecules from one side of the
system join water from the opposite side to form a bridge, closely followed by lipid head groups
in the WLW system.
One major difference however between WVW and WLW systems are the diameters of
the initial water channel. Figure 5.1 shows that in WLW systems there is usually a small, roughly
one-molecule-thick chain of water that precedes full hydrophilic pore formation. WVW systems
however exhibited wider water bridges on average spanning two-to-three water molecules in
diameter. Not every WLW system exhibited channels that were one-water-thick and not every
WVW system exhibited channels that were two-to-three waters-thick, but in general the channel
diameters created in WVW systems were larger than those created in lipid systems. In both
109

systems the newly formed column of water expanded to become an energetically minimized
electropore — a structure with water molecules in the middle and in WLW systems, lipid head
groups lining the periphery. A detailed description of this process can be found in earlier
publications [171, 175, 209].

Figure 5.1 Formation of a water bridge across a phospholipid (POPC) bilayer, showing water
only (A, B, C) and, for the same frames, water plus acyl oxygens (gray spheres), phosphorus
(gold), and nitrogen (blue) (D, E, F). Lipid tails are shown in cyan. Snapshots at 1585 ps (A,
D), 1761 ps (B, E), and 1830 ps (C, F) after application of a 320 MV/m electric field.
110

The dynamics of pore (water column) formation and the similarity between WLW and
WVW simulations are invariant across simulations over a wide range of parameters (e.g. external
electric field between 450 MV/m and 1000 MV/m, vacuum gap width ranging from 2.6 nm to
4.0 nm). WLW and WVW simulations differ mainly in the time scale over which the formation
of the water bridge occurs (i.e. the pore initiation time [93]). To compare initiation times
Figure 5.2 Formation of a water bridge across a vacuum gap separating two regions of liquid
water. Snapshots from 621 ps (A), 645 ps (B), 661 ps (D), 673 ps (E), and 683 ps (F) after
application of a 900 MV/m electric field. A magnified view of the water bridge precursor
structure (C), shows the alignment of the water dipoles in the electric field (directed
downward).
111

between WLW and WVW systems we selected a gap size for WVW systems that resulted in the
same local internal electric field magnitude in the vacuum gap as in the lipid bilayer interior.
This ensured that the interfacial water molecules are exposed to similar electric fields in the
WLW and WVW systems, permitting a fair comparison of pore initiation times. With equivalent
external and internal electric fields, WVW systems porated much quicker than WLW systems
(Figure 5.3). It appears as if the lipid bilayer acts as an inertial barrier to the rearrangement of
interfacial water dipoles, retarding the formation of the water bridge connecting the two water
layers. Similarly the rate of pore formation observed in simulations agrees well with the
stochastic pore hypothesis (Equation 1.2) which states that the time per unit pore (i.e. the inverse
of Equation 1.2) goes as exp(-voltage
2
).

Figure 5.3 Pore initiation times for water-vacuum and water-lipid systems. Average pore
initiation times for WLW and WVW systems calculated with three sets of simulations for
each configuration.
112

Energetics Analysis of the Systems. We now discuss the energetics of the simulated
systems. We can demonstrate that the pore formation process described above is driven by the
collective tendency of the interfacial water dipoles to minimize their electrostatic interactions
while adopting an orientation that minimizes the energy of the water dipole in the external
electric field, reflected in the steady drop in the per-molecule energy of waters in the nascent
pore as the protrusion develops. Below we demonstrate that this energetic behavior is present in
both WVW and WLW simulations.
In order to carry out this analysis we developed the following tools to examine the
energetics of water protrusion formation in ten WLW and ten WVW simulations. First we
defined the interface region of the water as follows. The computational box is split in the 𝑍
dimension into rectangular slices of thickness 0.1 nm. In each slice the average density of water
molecules is calculated, and the highest value is defined as bulk density. Slices with water
density not exceeding 80% of the bulk value are considered interfacial, and those with water
density equal to or greater than 80% of the bulk value are treated as bulk (Figure 5.4).

Figure 5.4 Classification of bulk water and interface water in a typical simulation.

113

Then in order to compute the time history of the protrusion energetics we identified the
water molecules directly involved in protrusion formation. Note that the trajectory of an
individual water molecule is very noisy so that selecting protrusion molecules for each time
frame is both impractical and arbitrary. Instead we selected a rectangular box that includes all of
the water molecules we consider to be in the protrusion at the end of our time interval of interest,
i.e. at the time when pore initiation is complete, then work backwards to track all water
molecules located in this box in each preceding time frame. To define the height of the box
containing water protrusion molecules we take the difference of the Z coordinate for the top
water molecule in the protrusion box to the Z coordinate at the boundary of the interfacial region
(i.e. the point above which the density of water does not exceed 80% of the bulk water density).
Once we isolated the water molecules comprising the protrusion we computed the time
history of the average potential energy per protrusion molecule and its constituent terms,
specifically: (i) the electrostatic interaction energy between water molecules in the protrusion
and all other water molecules, (ii) the Lennard-Jones approximation to the Van der Waal’s
interaction between protrusion and bulk water, and (iii) the interaction energy between the
protrusion water dipoles and the external electric field. For WLW simulations we also computed
the electrostatic and Lennard-Jones interaction energies between the protrusion water molecules
and the lipids. Note that all energetic terms are computed as averages per protrusion molecule.
That is, at each time frame we determine how many water molecules are in the protrusion box
and divide each of the interaction energy values by the number of protrusion molecules.
This analysis revealed a drop in the per-molecule energy of the waters in the protrusion
as the height of the protrusion grows for both WLW and WVW (Figure 5.5). To show a
correlation between this drop in energy and the protrusion growth for each simulation we
114

computed the Pearson’s correlation coefficient between the smoothed data for protrusion height
and each of the energies per protrusion water molecule over the time interval described above
(mid-point of protrusion growth). Specifically we examined the correlation coefficient between
the protrusion height and the sum of all the energetic interaction terms (i), (ii) and (iii) between
protrusion water and bulk water. Thirty independent simulations clearly demonstrated that the
growth of the protrusion is coupled with a drop in the average interaction energy for the water
molecules in the protrusion. On average the protrusion water molecule energy in WLW
simulations is reduced by 2.8 kJ/mol while for the WVW protrusion waters the total potential
energy decreases by 1.4 kJ/mol. Most of the change in potential energy arises from a change in
the dipole-dipole interaction energies between protrusion water and bulk water. The results for
the WVW simulations are noisier since water molecules move more freely into and out of the
interface and protrusion regions than they do in WLW systems, where the interface water
mobility is tempered by the presence of lipid headgroups.
These results present a scenario of protrusion evolution in which the interfacial waters,
constrained less than the molecules in the bulk, align with the external field and are restructured
into a protruding column where water dipoles are arranged vertically in a configuration that is
energetically more favorable than a horizontal arrangement in the plane of the membrane. The
external electric field drives the water molecules to overcome their interfacial and bulk
interactions and to form the protrusion. While it is not possible to consider the protrusion
molecules as an isolated subsystem with self-contained energetics and clearly decreasing total
energies, our analysis reveals similarities between the WVW and WLW systems at both
structural and energetic levels. Based on these results we argue that the same mechanism of
electrostatic energy minimization in the presence of external electric field drives formation of
115

pores for both WLW and WVW configurations, but the presence of lipid bilayer slows down this
process.

The following simple theoretical model illustrates the results of these simulations and the
energetic changes which can occur as a result of protrusion creation. Consider two configurations
of 7 dipoles: (I) a horizontal (i.e. perpendicular to the external electric field direction) sheet of
equidistant dipoles, where the mean dipole component < 𝜇 𝐸 > in the direction of the external
electric field 𝐸 is set according to the theoretically established dependency of < 𝜇 𝐸 > on the
magnitude of 𝐸 in bulk water (which can be computed analytically (Appendix A.1) for a system
of dipoles from a Langevin-Debye formula [183], < 𝜇 𝐸 >= 𝜇 �coth �
𝜇 𝐸 𝑘𝑇
� −
𝑘𝑇
𝜇 𝐸 �, or deduc ed
Figure 5.5 Anti-correlation of protrusion height and total interaction energy. Graphs
demonstrating anti-correlation between the increase of the protrusion height (black curve)
and the decrease of the total interaction energy per protrusion molecule of protrusion waters
with all other water molecules (both in the protrusion and in bulk) for (a) WLW (red curve)
and (b) WVW (blue curve) simulations.
116

from MD simulations [16]); (II) a vertical (i.e. parallel to the external electric field) single-file
chain of equidistant dipoles aligned with the external electric field. The distance between dipoles
is set to 0.31 nm in agreement with the average spacing of water molecules in the SPC/E model.
We now compute and compare the average energies of configurations (I) and (II) as follows. The
total energy for (I) is comprised of two terms: the dipole-dipole interaction, which represents the
dominant term of the electrostatic interaction between all of the dipoles, and the energy of the
dipole in the external electric field. In addition to these two terms the total energy for
configuration (II) includes an estimate of the desolvation energy it would require to remove the
dipoles from the bulk water [56].
For each of the two configurations we computed these total energies for different values
of an external electric field ranging from 0 to 1000 MV/m. Figure 5.6 shows the dependence of
the total energy of each of the configurations on the external electric field values. The results
presented in the figure demonstrate that we can expect the existence of a critical field value 𝐸 𝑐𝑟 𝑖𝑡
such that creation of configuration (II) becomes more energetically favorable then aligning the
dipoles in the bulk water. While obviously we cannot expect this basic model to provide a
precise estimate of the 𝐸 𝑐𝑟 𝑖 𝑡
value, it illustrates our theory of formation of protrusions and
bridges at the water-vacuum interface as a result of electrostatic energy minimization.

117

5.4 Water Bridges and Osmotic Swelling

In the previous section we explored how phospholipid bilayer electroporation could be
compared to simple water-vacuum systems which, under the influence of an external electric
field, can be transformed into a more energetically-minimized state. In this next section we will
explore the nature of these aqueous protrusions or ‘water bridges’ in water-lipid systems, water-
vacuum-systems, and in live cells to investigate at the atomic-scale explicitly how these systems
transform from one state to another. In the case of live cells we will also explore the role that
electropores may play in the osmotic swelling of cells immediately following the application of
electric pulses.

Figure 5.6 Comparison of constituent terms of the total interaction energy per protrusion
molecule between WLW (red curve) and WVW (blue curve) simulations. Interactions
consist of dipole–external electric field interactions, electrostatic interactions, and Lennard-
Jones interactions.
118

5.4.1 Simulation Methods

All simulations were performed using the GROMACS set of programs version 4.0.5 [64]
on the University of Southern California High Performance Computing and Communications
Linux cluster (http://www.usc.edu/hpcc/). Lipid topologies derived from OPLS united-atom
parameters [15] were obtained from Peter Tieleman (http://moose.bio.ucalgary.ca). The Simple
Point Charge (SPC) water model [12] was used for systems containing lipids; the Extended
Simple Point Charge (SPC/E) water model [13] was used for water-vacuum systems. All
simulations 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 [14].
For lipid systems, pressure was coupled semi-isotropically (using a compressibility of 4.5 x 10
-5

bar
-1
) normal to and in the plane of the membrane (NPT). No pressure coupling was used for
water-vacuum systems where volume was held constant (NVT). Bond lengths were constrained
using the LINCS algorithm [63] for lipids and SETTLE [109] for water. Short-range electrostatic
and Lennard-Jones interactions were cut off at 1.0 nm. Long-range electrostatics were calculated
with the PME algorithm [41] using fast Fourier transforms and conductive boundary conditions.
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. Diffusion of water was calculated from mean square displacements
over 50 ps using the ‘g_msd’ tool in GROMACS.
Systems and structures. All phospholipid systems contain 128 1-palmitoyl-2-oleoyl-sn-
glycero-3-phosphatidylcholine (POPC) lipids and 8960 water molecules (70 waters/lipid), which
results in a system box size of approximately 7 nm x 7 nm x 10 nm. To ensure that replicated
119

simulations are independent, each atom is assigned a randomized velocity from a Maxwell
distribution at the beginning of a simulation (a built-in function of GROMACS). POPC systems
are equilibrated before electroporation by waiting until they reach a constant area per lipid —
0.66 nm
2
. Water-vacuum systems are constructed by generating a 7 nm x 7 nm x 7 nm water box
with the GROMACS utility ‘genbox’, then using custom Perl scripts to remove a 4 nm slice of
water molecules from the center of the system, followed by 1 ns of equilibration at constant
volume. The final system contains 4650 SPC/E water molecules. Molecular graphics images
were generated with Visual Molecular Dynamics (VMD) [75].
For three-dimensional water dipole analyses, custom Perl scripts were used to generate
and normalize the data from single-precision molecular trajectories while MATLAB was used to
plot these trends on a 3D grid with dimensions equal to the simulation box and a resolution size
of 0.1 nm
3
per voxel. One-dimensional water dipole analyses was performed using the
GROMACS tool ‘g_h2order’. We also used MATLAB’s ‘equation and systems solver’ to create
custom trend lines for the cell swelling data in order to verify that our osmotic swelling
observations agreed with theoretical models.
A note on electric fields. The magnitude of the electric fields used in these simulations
may seem to some readers to be very high relative to the fields encountered in experimental
situations. However, consideration of the conditions being represented in the simulations, and a
closer look at the actual effective fields and potentials in experimental and simulated systems,
has led to a realization that these fields are not so high after all.
This can be understood first by referring to the well-established experimental fact that a
transmembrane potential as low as 400 mV, with 1 V considered to be a typical value, produces
membrane permeabilization. 1 V across the 4 nm thickness of the bilayer produces an electric
120

field of 250 MV/m, not far from the minimum porating electric field of 320 MV/m that we have
identified for POPC bilayers in previous publications [209]. Transmembrane potentials extracted
directly from molecular dynamics simulations are indeed in this range. And, as we and others
have shown, pore formation occurs in a similar fashion in our simulations when this porating
transmembrane potential is generated with an imbalance of ionic charge on the two sides of the
bilayer [196], instead of an externally applied field. In other words, the fields that we use in our
simulations are comparable to those that porating membranes actually experience in
experimental systems.
We can look at this from another perspective. The value of the electric field applied to a
molecular dynamics simulation volume (an electric field vector applied to every charge in the
system) is analogous to the field that would be measured across a vacuum of similar dimensions.
If we fill the volume with a dielectric (water in the case of our systems), the resulting effective
electric field is a smaller value, reduced roughly by a factor of the relative permittivity of the
dielectric (80 for water), and this is the field we see experimentally when we measure the
potential across a cell suspension between parallel-plate electrodes, for example, and divide it by
the electrode separation distance. Thus 320 MV/m in our molecular dynamics simulations of
phospholipid bilayers corresponds to a field in an electroporation cuvette of approximately 4
MV/m (320 / 80), well within the range of published experimental values for nanosecond pulse
exposures.

121

5.4.2 Experimental Methods

Cell Culture and Microscopy. Jurkat T lymphoblasts (ATCC TIB-152) were grown in
RPMI 1640 (Mediatech) containing 10% heat-inactivated fetal bovine serum (Gibco), 2 mM L-
glutamine (Gibco), 50 units/mL penicillin (Gibco), and 50 µg/mL streptomycin (Gibco). Cells
were cultured at 37 °C in a humidified, 5% CO
2
atmosphere and concentrated to 2 x 10
7
cells/mL
for pulse treatment.
Cell Preparation. For pulse treatment, cells were concentrated to 2 × 10
7
cells/mL and in
some cases incubated for 15 min with 0.5 μM calcein-AM (acetoxymethyl ester, Molecular
Probes, Eugene, OR), which enhances visualization of the cell outline using fluorescence
microscopy. After loading of the dye, cells were centrifuged and re-suspended in fresh RPMI
1640 (or in some cases 150 mM NaCl) at 2 × 10
7
cells/mL. For exposure in presence of
lanthanide or mercuric ions, GdCl
3
or LaCl
3
(Aldrich Chem. Co, Milwaukee, WI), 100 μM or 1
mM, or HgCl
2
, 50 μM, was added to the cell suspension 5 min before pulse treatment.
Pulsed Electric Field Exposures. For microscopic observation, cells were placed in a
microchamber 100 µm wide, 30 µm deep, and 15 mm long, with platinum electrode walls on a
glass microscope slide. A resonant-charged, solid-state Marx bank-driven, hybrid-core
compression, diode-opening switch pulse generator designed and assembled at the University of
Southern California [149] delivered 5 ns, 10 MV/m electrical pulses at 1 kHz repetition rate (1
Hz where noted) to the microchamber electrodes mounted on the microscope stage in ambient
atmosphere at room temperature.
DIC Microscopy. Observations of live cells during and after pulse exposure were made
with a Zeiss (Göttingen, Germany) Axiovert 200 epifluorescence microscope with 63× water
immersion objective and Hamamatsu (Higashi-ku, Hamamatsu City, Japan) ImageEM EM-CCD
122

camera. Captured images were analyzed with Hamamatsu SimplePCI and ImageJ
(http://imagej.nih.gov/ij/) software. The cell perimeter was tracked using a freehand selection
function, and then the area defined by the drawn perimeter was measured. To reduce variability,
cells in the center of the exposure chamber, not adjacent to the electrode surfaces, were selected
in the images captured immediately before pulsing, and then the same cells were analyzed in the
post-pulse images. For each pulsing condition, at least 3 experiments were carried out, and a total
of at least 30 cells were analyzed.

5.4.3 Results and Discussion

Although we observed in section 5.3 that water-vacuum systems experienced pore
initiation, or water bridge construction more quickly than in water-phosphilipid bilayer systems,
we did not present an atomic mechanism to explicitly account for this other than the presence of
a generalized ‘passive’ lipid bilayer slab which acted as an energetic barrier to pore formation.
Here however we will present a more detailed analysis on the differences between water-lipid
systems and water-vacuum systems to more completely define the relationship between
interfacial water and lipid bilayer electroporation.
Water orientation and pre-pore pedestals. When an external electric field of 700 MV/m
was applied across a 10 nm simulation box (V = E*L
z
= 0.7 V/nm*10 nm = 7 V) we observed
changes in the z-dipole moment (or orientation) of interface water such that they try to align
themselves in the direction of the field, similar to an ideal dipole. This occurred in both water-
vacuum and water-lipid systems. In bulk water the presence of numerous hydrogen bonds from
adjacent waters keep individual water molecules fixed from rotating unless they gain ~44 kJ/mol
123

of energy to escape [8]. In the interface region however low-density water is better able to orient
itself because it is further away from the isotropic hydrogen bonds present in bulk, and here we
see an interesting difference between lipid and vacuum systems. Although we, and others [209],
do not observe a cathode or anode preference for pore initiation to occur on, specific systems
appear to indicate that the external electric field perturbs the water-vacuum cathode interface
much more than in water-lipid systems [Figure 5.7]. Local charge densities reveal that although
interface water dipoles are initially directed inwards towards the low-dielectric interior (a direct
result of lipid headgroup hydration for water-lipid systems) the presence of a strong electric field
can negate, or sometimes reverse the polarization of interface water molecules on the cathode
end of water-vacuum systems. This is perhaps not a surprising result since one might imagine an
electric field strengthening water dipoles that are already oriented in the direction of the field,
and weakening those that are oppositely-aligned, however these simple dynamics may support
the observations that water-vacuum systems form water bridges more quickly because they are
systems which experience larger instabilities in the energetic requirements for bulk water to
transition into interface water.

124

To explore this idea further we measured the average water orientation on both the
cathode and anode facing interfaces to observe how interface water changes as a function of time
after a field has been applied. These measurements [Figure 5.8] confirm our previous
observations that water-vacuum systems experience a sudden change in the cathode-facing
interfacial water orientation as a field is applied, so much so that the average orientation on the
cathode end is net positive, or rotated from opposite-the-field to parallel-the-field. Lipid systems
on the other hand maintain their inward facing water dipole moments when an electric field is
present, likely because interface water is strongly associated with the lipid headgroup phosphate
and choline groups. As a result of preserving its initial polarization, lipid systems may maintain
much of its interfacial energy barriers in order to slow the large-scale delivery of bulk water
molecules into the interface as raw building materials for the construction of stable water
bridges. Although water-vacuum systems were run in an NVT ensemble and water-lipid systems
in an NPT ensemble, these results were still present when similar thermodynamics ensembles
Figure 5.7 A one-dimensional charge density profile after an external electric field is applied
for a water-lipid and water-vacuum system.
125

were used, however water-vacuum systems were much more stable at constant volume, which is
why we elected to run water-vacuum simulations at constant volume.

Figure 5.8 Measurements of interfacial water orientation for water-lipid and water-vacuum
systems after an external electric field was applied. The angle theta is measured with
respective to the direction of the applied electric field. Note that these two systems are un-
normalized and thus their pore creation times are not directly comparable.
126

To obtain additional information about the three-dimensional geometry of interface water
under an external electric field we also generated 3D voxel profiles of the water dipole moment
which illustrate the spatial orientation of all water molecules in the system (Figure 5.9). Just as
we saw in the one-dimensional case, water-vacuum interfaces are less ordered than water-lipid
interfaces, and the magnitude of interface polarization is much stronger in lipid systems which
suggest that charged interfaces may take longer to permeabilize compared to uncharged ones.
These results are also supported by the data in Figure 5.3 where lipid systems consistently
required more time to permeabilize when the local internal electric field was normalized to be
the same as the internal electric field found in lipid simulations.

Figure 5.9 Three-dimensional plots of the mean water dipole moment (z-component) per
voxel during water bridge formation under an externally applied electric field. Water-vacuum
interfaces are highly disordered, on average, while lipid interfaces are highly ordered.
127

Our simulations also reveal key details about how water bridges are explicitly formed at
the interface when the conditions are right. As we discussed in Chapter 3, the first identifiable
event in MD lipid nanopore formation is the appearance of a water intrusion into the lipid tail
region of the bilayer (Figure 5.10B). These intrusions are common (they occur even in the
absence of an external electric field, albeit less frequently), and they usually disappear, but
occasionally (more frequently at higher electric fields) they are extended, one, two, three
molecules at a time, until they form a bridge, spanning the low permittivity bilayer interior
(Figure 5.10C), an essentially irreversible stage in the process of formation of a hydrophilic pore,
a water column lined by phospholipid head groups. Water molecules lead the process of bridge
construction, followed closely by the phospholipid head groups as membrane structural
realignments accommodate the growing defect in the bilayer. The question we seek now to
answer is “What drives water, a highly polar compound, to form these chains across the apolar
environment of the membrane interior?”
In our explorations of the dynamics and energetics of water bridge formation we noticed
similarities between this process and electroevaporation, a phenomenon that is important in
several technological applications. Indeed, the nanoscale structures that nucleate the electric
field-enhanced evaporation of water [122] look like the cones of intruding water that we observe
in our simulations (Figure 5.10B). Figure 5.11 shows a magnified view of what may be called an
Okuno-Tanioka cone, after the similar structures reported by those authors in their physical
chemical analysis of water electroevaporation, where they demonstrate that this particular
conical structure lowers the energy needed to remove the water at the peak from the interface.
Note the alignment of the hydrogen-bonded water molecules in the electric field, which is
effectively higher in the vacuum than it is in the high-permittivity environment of the bulk water.
128

We suggest that this structure facilitates the formation of a membrane-spanning water column
just as it enhances the escape of a water molecule into vacuum or air.

Figure 5.10 Electropore creation sequence. (A) Molecular dynamics representation of a
POPC lipid bilayer system. Small red and white spheres at the top and bottom of the panel are
water oxygen and hydrogen atoms. Gold and blue spheres are head group phosphorus and
nitrogen, respectively, and large gray spheres are phospholipid acyl oxygens. For clarity,
atoms of the hydrocarbon chains in the interior of the bilayer are not shown. In the presence
of a porating electric field, a water intrusion appears (B), and extends across the bilayer (C).
Head groups follow the water to form a hydrophilic pore (D). The pore formation sequence,
from the initiation of the water bridge to the formation of the head-group-lined pore takes less
than 5 ns.
129

The water bridges formed in both water/vacuum/water and water/phospholipid/water
systems are highly dynamic and unique from one independent simulation to the next, and a small
set of snapshots cannot capture all of the variations and details of this phenomenon. The
formation of an Okuno-Tanioka cone-like structure, however, is an early, distinctly recognizable,
common feature of all of the thousands of systems we have inspected, and the similarities in the
electric field-driven rearrangements of water molecules, whether the gap is vacuum or
Figure 5.11 Water bridge precursor (pre-pore pedestal) in a POPC bilayer. Note similarity to
structure in Figure 5.2C. Acyl oxygens, relatively electron-dense sites deep in the membrane
interface, attract water and are frequently structurally associated with the water cone that
penetrates into the bilayer interior. Water molecules are red (oxygen) and white (hydrogen).
Isolated gray spheres are the acyl oxygens on the phospholipid backbone. Head group
nitrogen is blue; phosphorus is gold.

130

hydrocarbon, are striking.
Membrane permeation by ballistic water. Water can cross a lipid bilayer without the
assistance of an intruding cone, and without the formation of a water bridge. Occasionally, about
once per nanosecond in our POPC simulations at 310 K, an isolated water molecule crosses from
one leaflet of the bilayer to the other (Figures 5.12 and 5.13). Application of a porating electric
field has little effect on either the frequency or the direction of these random crossing events.
The primary determinant of the frequency (expressed experimentally as the permeability
coefficient) is the area per lipid [103]. Systems with a smaller area per lipid have smaller
permeability coefficients [180]. Consideration of this phenomenon can enhance our
understanding of the mechanism of water bridge construction, from pre-pore pedestal to
extension of the water chain across the bilayer.

131

Figure 5.12 Random permeation event (I) in a phospholipid bilayer. Single H
2
O molecule
(yellow) crosses the phospholipid bilayer interface without an external electric field present,
presumably due to thermal jostling. (A) 0 ps (B) 50 ps (C) 105 ps (D) 110 ps .
132

Figure 5.13 Random permeation event (II). Single H
2
O molecule crosses the phospholipid
bilayer interface. Applied electric field is 320 MV/m downward in the diagram. (A) 0 ps, (B)
40 ps, (C) 80 ps, (D) 120 ps. Only acyl oxygens and water molecules in the interface are
shown. The crossing water molecule is enlarged and colored yellow. Gray spheres are acyl
oxygens; red and white atoms are water. From a system containing 128 POPC, 4480 H2O, at
310 K, 1 bar.
133

Some of the detailed intermolecular interactions that “kick” a water molecule into and
across the membrane are shown in Figure 5.14, where the crossing water in Figure 5.13 is
enlarged. Thermal rotation of the arrowed green water molecule brings its hydrogens into close
and repulsive proximity to the hydrogen on the soon-to-be crossing water molecule, which is
already poised on the edge of the low-permittivity hydrocarbon region of the lipid tails that make
up the membrane interior. Note that there are acyl oxygens surrounding this event, an important
component of three-dimensional electrophysical contours of the site, noted also by Vacha et al.
[180]. Similar but less energetic interactions lead to the assembly of intruding Okuno-Tanioka
water cones, with their associated acyl oxygen and other lipid atom neighbors, and to the
molecule-by-molecule (sometimes plus one, minus two,…) growth of bilayer-spanning water
columns. On a macro scale, this is statistical mechanics. On the atomic scale in a phospholipid
bilayer the sum of these events is a new structure, the chain of waters through the membrane
Figure 5.14 Detail of initiation of permeation event shown in Figure 5.13. Repulsive
interactions between the hydrogen atoms of the crossing water (yellow, enlarged) and another
interface water molecule (green, arrowed) “kick” the crossing water out of the interface and
into the membrane interior. (A) 6 ps, (B) 10 ps, (C) 14 ps, (D) 18 ps. Times are relative to 0
ps in Figure 5.13.
134

interior (Figure 5.10C) that results in the formation of a hydrophilic pore (Figure 5.10D), a
conduction pathway through a normally impermeant barrier.
Figure 5.11 shows in detail a pre-pore pedestal in a POPC bilayer, what we have also
referred to here, focusing on the water component, as an Okuno-Tanioka cone. For simplicity we
display only water, acyl oxygens, and the phosphorus and nitrogen atoms of the head groups.
The surrounding environment is similar to that in Figure 5.1D. Structures like the one shown in
Figure 5.11 may be transient or they may lead to the formation of a lipid nanopore, as described
above, and shown in Figure 5.10. An external electric field, normal to the plane of the bilayer,
increases the probability of pore formation. The strong, aligning influence of the applied electric
field is evident in the orientation of the water dipoles nearest the tip of the cone. Toward the base
of the cone the water molecules are also oriented by the charged phosphate and amine groups
which they are solvating. Note also the presence of an acyl oxygen in the intruding structure.
Although we have shown here only cones intruding from the cathode side of the bilayer,
with the oxygen of the water at the tip pointing across the membrane interior toward the anode,
structures of the opposite polarity are also formed, at a lower frequency. In that case the cone
intrudes from the anode-facing leaflet, with the hydrogens of the water at the tip pointing across
the membrane interior at the cathode. It is not unusual for the appearance of a pedestal on one
leaflet to induce, through field enhancement and electrostatic attraction, the formation of a cone
on the opposite leaflet.

135

Osmotic Swelling and Electroporation. Water, primary player in electroporation, also
appears to affect live cells after pulsed electric fields are applied. More specifically, Jurkat T
lymphoblasts exposed to 5 ns, 10 MV/m electric pulses in physiological mediums (RPMI 1640)
exhibit an increase in cross sectional area (Figure 5.15) by as much as a factor of 1.5 (1.8-fold
volume increase, if the cells are perfectly spherical) within 3 min after pulses. The rate of
swelling is greatest immediately following pulse exposure, declining exponentially over tens of
seconds, but not falling to zero after 3 min. Cells were not observed to recover their initial
volume, under the conditions in these experiments, even after 7 min, the longest monitoring time
permitted by our apparatus. Cells do not respond uniformly, as shown in Figure 5.15. The
responses of individual cells vary not only in the magnitude of the volume increase, but also in
the visual appearance of the cells during the first minute after pulsing. Some cells expand more
or less uniformly, maintaining their roughly spherical shape. Others initially form one or more
blebs of various sizes which are subsumed as the cell continues to swell and recover its
spheroidal geometry. An increase in cell area can be detected even after a single 5 ns, 10 MV/m
pulse.

136

The phenomenon of pulse-induced osmotic imbalance and resulting cell swelling are
consistent with the mechanistic hypothesis that nanosecond pulses open a large number of small
pores on the plasma membrane, with a pore creation rate larger than the pore expansion one.
Inorganic ions and small molecules, but not larger solutes, pass through these defects [184].
From our presented results, and building on conclusions supported by previous studies [30, 54,
68, 123], we propose the following hypothetical sequence of events in Jurkat cells responding to
5 ns, 10 MV/m electric pulses delivered at 1 kHz. All of the observations we describe are
consistent with this scheme.
1. Pulse exposure results in the immediate formation of a dose-dependent number of lipid
nanopores of a dose-dependent total area, large enough to allow passage of small inorganic ions
[123]. These pores are formed during the 5 ns duration of a pulse. 2. Equilibration of intracellular
Figure 5.15 Osmotic swelling after nano-electropulse exposure. (A) Jurkat cells before pulse
exposure. (B) Same cells 60 s after exposure to 30, 5 ns, 10 MV/m pulses at 1 kHz. (C), (D).

137

and extracellular [Na+], [K+], and [Cl−] by diffusive ion transport through the permeabilized
membrane (presumably through the lipid nanopores but perhaps also through other structures)
proceeds for several minutes [152, 156]. 3. During [Na+], [K+], and [Cl−] equilibration, water
flux (bidirectional under balanced conditions) becomes net positive inward, driven by the
osmotic gradient caused by the reduction of the mobile ion concentration gradients and the
remaining intracellular colloidal anions and other large molecules. This osmotically driven water
influx is the primary cause of pulse-induced cell swelling [116]. 4. When the [Na+], [K+], and
[Cl−] concentration gradients across the membrane reach zero (or some minimum value), water
influx continues at a continually decreasing rate, as the internal osmolytes are diluted, until the
ion-permeant structures (pores) close, allowing re-establishment of small ion concentration
gradients, or until mechanical constraints become significant, or until the cell bursts.
There is also a relative cell area increase of 3% over the 10 s immediately following a
single 5 ns, 10 MV/m pulse. For a cell radius of 5 μm, that is a volume change of 2.3 × 10
−17
m
3
,
which corresponds to 7.7 × 10
11
water molecules over 10 s. If the expansion were linear this
would be a net 77 water molecules entering the cell per nanosecond. If the swelling is described
by a decaying exponential, then the initial swelling rate is higher. For an initial approximation
we assume 77 H
2
O/ns. Even though this is a rough number, it provides a basis for comparison
with experimentally known transport rates of water.
It might be expected that water passes through the lipid nanopores created by the pulse,
and some net water transport undoubtedly does occur this way, perhaps through aquaporin
channels as well. But the water permeability of lipid bilayers and cell membranes (P
f
= 1.3 ×
10
−4
m·s
−1
for POPC [103]; P
f
= 6.7 × 10
−6
m·s
−1
for a more complex bilayer [87]), is high
enough to account for a significant amount of the observed pulse-induced water influx. Lipid
138

nanopores are required for the ion transport and concentration equilibration that occurs after
pulse exposure, and this drives the osmotic flow of water into the cell, but non-permeabilized
membrane leakage pathways appear to be a sufficient mechanism to facilitate post-pulse water
entry.
Molecular simulations reveal how water permeates a lipid bilayer without pore formation
[99, 100]. About once per nanosecond in simulations of 50 nm
2
POPC bilayers at 310 K [187],
an isolated water molecule crosses from one leaflet of an intact bilayer to the other (Figure 5.12).
No pore is involved, and application of a porating electric field has little effect on either the
frequency or the direction of these random, diffusive crossing events. Although detailed
calculations would be premature at the present level of accuracy and understanding of the
correspondence between MD simulations of lipid bilayers and living cell membranes, we can
translate the random and bidirectional 1 H
2
O per nanosecond per 50 nm
2
that we observe in our
POPC simulations to an osmotically driven (initially 100 mOsm) net transport of about 104 H
2
O
per nanosecond per cell, declining exponentially as the ion concentration imbalance is reduced.
Although it is not possible from our data to rule out the possible participation of
membrane protein channels in this process, we have two sets of observations that are consistent
with a major role for lipid nanopores and a minor role for other potential mechanisms. First, we
see no significant inhibition of nanoelectropulse-induced swelling by the ion channel blockers
Gd
3+
and La
3+
[6, 19, 90, 206] at concentrations of 0.1 mM and 1 mM. This is only a non-
specific and imprecise indication, of course, but this result stands in interesting contrast to a
report of Gd
+3
inhibition of osmotic swelling [4]. Pulse characteristics and exposure conditions,
however, are quite different in the two cases. Our pulses are 5 ns, 10 MV/m; the pulses in [4] are
60 ns and 600 ns, 1.2 MV/m. Longer pulses may induce a more extensive membrane
139

restructuring that is Gd
3+
-sensitive. Second, to see whether the opening of aquaporin channels
may be involved in the swelling response (including possibly voltage-sensitive aquaporin
channels [74]) we carried out pulse exposures in the presence of mercuric ion, an inhibitor of
aquaporins [67,135]. The results were similar to those obtained with the lanthanide ions. In these
experiments Hg
2+
had no significant effect on nano-electropulse-induced osmotic swelling of
Jurkat cells.
Our hypothesis is that nano-electropulse-induced cell swelling is caused by the diffusive
loss of intracellular-extracellular ion concentration gradients and osmotically driven water influx.
We have extracted from the swelling data a simple functional expression that describes our
observational data and which corresponds to theoretical mechanisms for this phenomenon
outlined in Appendix A.2:

𝐴 𝑅 ( 𝑡 ) = 𝑏 ∗ 𝑒 𝑥 𝑝 ( − 𝑐 𝑡 ) + 𝑎𝑡 + 𝑑 (5.1)

where A
R
(t) is the time-dependent relative cross-sectional area of the cell, a is the slope of the
linear relative area increase, b is the initial point of exponential relative area increase (at t = 0), d
is the initial point of linear relative area increase, and c is the exponential scaling factor. Note
that at t = 0, b + d = 1 (the area does not begin to increase until just after the pulse exposure).
Figure 5.16 shows how this function can be fit to the swelling response over a wide range of
exposures — 1, 10, and 50 pulses.

140

The nanoelectropulse-induced swelling response can now be analyzed as an overlapping
sequence of phenomena associated with membrane permeabilization, diffusive ion transport, and
the osmotic flow of water. We see now in this functional form the same four steps we proposed
at the beginning of this discussion: 1. lipid nanopore formation (and perhaps other
electropermeabilization-associated restructuring of the cell membrane); 2. equilibration of
intracellular and extracellular concentrations of sodium, potassium, and chloride ions ([Na+]
i
,
[K+]
i
, and [Cl−]
i
] and [Na+]
e
, [K+]
e
, and [Cl−]
e
), at an exponentially decreasing rate
corresponding to the term be
−ct
; 3. net inward H
2
O flux which increases during the period of
small ion equilibration and which creates the osmotic imbalance that drives water flow; and 4. a
continuing inward H
2
O flux at a continually decreasing rate as internal osmolytes are diluted.
Figure 5.16 A simple exponential function describes the time course of nano-electropulse-
induced cell swelling. The comparison between experimental data and the fitting function is
reported for cell swelling after 1, 10 and 50 pulses.
141

The coefficient c has units of inverse time, so we may consider 1/c to represent a time
constant for the overall equilibration of small ions (primarily Na+ and K+) through the
permeabilized membrane. Although we understand that we risk oversimplifying a complex
situation, it may be instructive as a guide to further experimentation and better models to
consider the implications of the value of this constant that we extract from our data. For the
simplest case — single-pulse-induced swelling — we obtain values for τ (1/c) of about 10 s
[132]. Continuing our simple interpretation, we can then say that the ion concentration gradients
are reduced by 1/e over the course of one time constant.

5.5 Conclusions and Summary

Our simulations and analysis show that at the molecular scale electroporation of a
phospholipid bilayer is driven by the restructuring of interfacial water molecules into column-
like structures as their dipole moments align with an external electric field. This electric field-
driven reorganization of the lipid bilayer is associated with an overall drop in the per-molecule
energy for the waters in the growing protrusion. Membrane phospholipids simply follow the
bridging water. This view allows significant reduction of the complexity of the analysis of
electroporation and opens possibilities for applying well-developed analytical and computational
tools to study this problem. Additionally, this insight into the significance of the interfacial water
dynamics can facilitate development of new experimental and technological approaches to better
control and utilize the process of electroporation.
142

Additionally, molecular dynamics simulations of phospholipid bilayers in electric fields
show the sequential, cooperative, electric field-driven rearrangement of water and phospholipid
molecules into a membrane-spanning pore of nanometer dimensions. The time scale of lipid
nanopore formation and the physical properties of the pores are consistent with experiment.
Underlying the poration process we find a more fundamental re-ordering of water molecules at
the boundary of bulk water and a low-permittivity medium (vacuum, hydrocarbon), leading to
energy-minimized water protrusions from the interface that facilitate the escape of water
molecules into a vacuum (electroevaporation) or, in the case of lipid membranes, the
construction of water bridges across the bilayer interior. At the atomic scale the mechanism for
these restructuring processes is stochastic thermal jostling, biased by local electric field
gradients. In phospholipid bilayers the acyl oxygens and charged head groups structurally and
mechanistically contribute to the water-led reorganization of the membrane. Increasingly
available computational power will enable the linking of molecular simulations, continuum
models, and experimental observations with artificial membranes and living cells, strengthening
the theoretical foundation for electropermeabilization-electroporation biotechnology
applications.
We also demonstrated that osmotically driven cell swelling can be used as a sensitive
method for studying plasma membrane electropermeabilization induced by very short and very
intense pulsed electric fields, and that the dynamics of water even after electroporation occurs is
important in determining how cells react to pulsed electric fields. We successfully characterized
the cell swelling response based on simple osmotic pressure models, and we observed in
experiments that swelling is not sensitive to the presence of the lanthanide ions (as ion channel
blockers) or mercuric ions (as aquaporin channel inhibitors) in the external medium. Swelling
143

kinetics can be described with a functional expression that matches a hypothetical mechanism for
permeabilization-induced osmotic swelling — an initial loss of the normal physiological
inorganic ion gradients followed by an influx of water in response to the colloidal ion imbalance
that remains after permeabilization. Although further investigations and more refined models are
needed to clarify the details of the mechanisms underlying electric pulse-induced membrane
permeabilization, the approach adopted in this study allowed us to address the most basic events
of the electroporation phenomenon and to establish a simple but explicit connection between MD
simulations and experimental observations.

5.6 Acknowledgements

We thank Andrei G. Pakhomov, Olga N. Pakhomova, Justin Teissié, D. Peter Tieleman,
James C. Weaver, and Olga Zeni for stimulating discussions, helpful suggestions, and probing
comments. 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. SR was supported in part by the Second University of
Naples.

144

Chapter 6: The Role of Lipids in
Electroporation

6.1 Abstract

Although there is still some disagreement as to which biological molecules are primary
facilitators of electroporation, interface water appears to be the dominant driver for many of the
permeabilization dynamics we observe in molecular dynamics simulations. And yet, in
simulations of the electroporation of phospholipid bilayers, the role of lipids cannot be ignored.
Phosphate atoms for instance, in particular the acyl and phosphoryl oxygens and the amino
hydrogens, respond to an applied electric field independently of water, and because they are
solvated by water, their position and orientation affects all associated water structures. The
electron-dense acyl oxygens on the glycerol backbone of the phospholipid molecule lie deep in
the membrane interface, providing stable local energy minimas for intruding water molecules. To
a much greater extent than the charged head group atoms (phosphate and amine), the acyl
oxygens on the lipid glycerol backbone are associated with the initial, intruding water cone
during water bridge formation. If lipids are not the primary drivers of electroporation they are at
least a close second in determining the characteristics and timescales of cell
electropermeabilization, and so they must also be studied in sufficient detail. Here I will present
research on the effects of temperature on electroporation, specifically how lipid temperature
affects the pore life cycle. Then I will briefly discuss how molecular advection, self-diffusion,
surface tension, edge tension, and steric repulsion affects electroporation by implementing
145

simulations with customized electric pulses and periods to control the pore size and the
expansion and resealing rates of electropores.

6.2 Introduction

Molecular-scale details of the mechanism of electric field-driven pore formation in
phospholipid bilayers are not well understood, in part because experiments cannot access the
nanoscopic domain at which individual pore formation occurs. Analytical and numerical models
can help to fill this void. Previous studies using molecular dynamics (MD) simulations
characterized the life cycle of an electropore as a function of the externally-applied electric field,
from the formation of an initial water defect to the restoration of the intact bilayer [93], but this
is far from a complete description of pore creation and annihilation. Here we continue our
characterization of electropores, modulating the temperature at which electropermeabilization
occurs, allowing us to observe how temperature-related physical properties of lipids affect the
timescales of pore formation and annihilation. These results can be compared to and reconciled
with existing mathematical models of electroporation which incorporate temperature, presenting
a more unified and complete framework for future studies.
The question then that we wish to answer is whether the properties of lipids and lipid
bilayers strongly affects water bridge formation, and are these effects direct or indirect in their
interactions with interface water molecules? Because phospholipid acyl oxygens provide a low-
energy “resting point” which can further increase the likelihood of the formation of an intruding
water cone, it might be expected that they facilitate the formation of lipid nanopores, or water
146

bridges, across a phospholipid bilayer. But the same hydrogen-bond attractions to the electron-
dense acyl oxygens that makes it easier for water to reside deep in the membrane interface also
makes it more difficult for a water molecule to continue into the hydrocarbon interior of the
bilayer, so there is no obvious answer to this question.
Certain lipid properties however do affect pore creation [209] and here we wish to
investigate the roles of temperature and self-diffusion to see if they too affect the detailed
construction of pre-pore pedestals, further impacting the timescales and kinetics of lipid bilayer
electroporation. Lipid tail rigidity for instance may impede the entry of interfacial water
molecules into the membrane interior, subsequently slowing the creation of new water pedestals.
Similarly if a membrane is sufficiently viscous and has high surface tension, the electric pulses
required for permeabilization may need to be modified or have higher amplitudes to ensure that
water bridges can be created and joined together before they are quickly deconstructed by the
hydrophobic forces of the surrounding hydrocarbon interior. Near the phase transition
temperature lipids will also experience changes in trans-gauche isomerization [91] and their lipid
tail order parameters will become more uniform. By using molecular dynamics simulations to
investigate the role of lipid temperature in electroporation we can relate our observations in MD
to well-known experimental results and continuum models which provide further support for
atomistic simulations as a valid and predictive tool in the ongoing study of nanoscopic
biophysical processes.
Other factors which may affect electropore formation are the mechanical properties of
lipids. For instance when dealing with pulses in the low nanosecond and subnanosecond range,
the rise time of the actual waveform delivered to a biological load is a significant percentage of
the entire pulse duration, but experimental and modeling studies of the effects of varying the rise
147

time are rare. If the rise time is long then the membrane may only experience the intended
maximum voltage for fractions of the pulse period, so it becomes necessary to study how
different electric pulses and pulse shapes affect the mechanical properties of lipid electropores.
Molecular dynamics (MD) simulations have contributed to our understanding of lipid bilayer
electropermeabilization, and we show here that when time-varying external fields are introduced
into MD, it becomes possible to modulate geometric electropore properties dynamically by
sampling a large range of frequencies, fields, pulse shapes, and polarities with various charged
and uncharged bilayer systems. This capability also enables the refinement of various models of
nanosecond pore formation by loosening constraints requiring constant external perturbations.
Results are compared to and, to the extent possible at this time, reconciled with existing
mathematical models of electroporation, presenting a more unified and complete framework for
future studies. Our results also allow for the indirect measurement of a handful of mechanical
lipid properties as a basis for more robust molecular models, in addition to providing a novel
way of verifying the validity of current atomistic models.

148

6.3 Lipid Temperature affects Pore Annihilation

6.3.1 Methods

All simulations were performed using the GROMACS 4.0.5 software package on the
University of Southern California High Performance Computing and Communications Linux
cluster. Bilayers consisted of 128 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine
(POPC) lipids parameterized by OPLS/Berger force fields and 35 SPC water molecules per lipid
at an integration time-step of 2 fs. Pore life cycle times were determined using custom Perl
scripts as previously reported. Temperature was modulated using a modified Berendsen
thermostat with corrected velocities known as the velocity rescaling thermostat with a time
constant of 0.1 ps. An NPT ensemble was maintained using the Berendsen barostat which
applied 1 bar in xy:z semi-isotropically with a compressibility of 4.5E-5 bar
-1
. Bond lengths were
constrained using the LINCS algorithm for lipids and SETTLE for water. Short-range
electrostatic and Lennard-Jones interactions were cut off at 1.0 nm. Long-range electrostatics
were calculated by the PME algorithm using fast Fourier transforms and conductive boundary
conditions. Reciprocal-space interactions were evaluated on a 0.12-nm grid with fourth-order B-
spline interpolation. Periodic boundary conditions were employed to mitigate system size effects.
Lipid tail order parameters were calculated for every i
th
atom by using a second-order
Legendre polynomial with an angle theta measured from an axis connecting the i-1 atom and the
i+1 atom to the z-direction which is normal to the plane of the bilayer. Area per lipid was
calculated by dividing the box area (X
box
times Y
box
) by the total number of lipids in one of the
bilayer’s leaflets (64). Heat capacity at constant pressure was calculated by measuring the change
in enthalpy (extracted from the GROMACS tool ‘g_energy’) per unit change in temperature.
149

Lipid lateral diffusion coefficients (restricted in the plane of the membrane) were obtained by
measuring the two-dimensional mean square displacement of lipids as a function of time and
invoking the Einstein diffusion relationship in the linear regime. We set the temperature of lipids
by adjusting the velocity rescaling thermostat to the desired temperature where we then waited
for the temperature, area per lipid, and energies to converge within 5% of their stable final
values. Lipid temperatures varied from 50 C to -5 C. Studies have shown that models of POPC
undergo a phase transition around -2 to -3 C [91] which is consistent with experiment, however
others have found that OPLS/Berger POPC models require much lower temperatures to
transition from the gel to liquid phase. Bilayer thickness was measured by measuring the P-P
distance from mean phosphorus planes calculated only from lipids on opposing leaflets.
Identification of ‘pore waters’ was taken by counting the number of water molecules which lie in
any z-plane containing at most 50% bulk water density (500 kg/m
3
) and within 2 nm in x and y
directions from the axis which traverses normally through the center of the open electropore.

6.3.2 Results and Discussion

To get a feeling for the magnitude of thermal energies present in these systems compared
to the energies delivered from the external electric field, consider for a moment the mean
orientation of a lone dipole caught in a uniform electric field (which was derived in Appendix
A.1): < 𝜇 𝐸 >= 𝜇 �coth �
𝜇 𝐸 𝑘𝑇
� −
𝑘𝑇
𝜇 𝐸 �. This states that the electric dipole orientation will resemble
a hyperbolic tangent at small arguments (rotating randomly around zero orientation as
temperature dominates) and approaching saturation (+ 𝜇 or - 𝜇 ) at large arguments when the field
150

is dominant. In the small electric field limit ( 𝜇𝐸 ≪ 𝑘𝑇 ) we can do a Taylor expansion on the
Langevin function to obtain:

𝐿 ( 𝑥 ) = coth( 𝑥 ) −
1
𝑥 ≈
1
3
𝑥 −
1
4 5
𝑥 3
+ ⋯ (for x << 1) (6.1)
< 𝜇 𝐸 > ≈ 𝜇 �
𝜇 𝐸 3 𝑘𝑇
� = �
𝜇 2
𝐸 3 𝑘𝑇
� (for 𝜇𝐸 << 𝑘𝑇 ) (6.2)

In other words at small electric fields a water molecule (analogous to the electric dipole) orients
itself linearly with the external electric field, but for large electric fields ( 𝜇𝐸 ≫ 𝑘𝑇 ) the dipole
moment is fixed at ± 𝜇 . So how high does the external electric field need to be in order to negate
the contributions from thermal energy? If we take T = 310 K (body temperature), 𝜇 = 2.23 Debye
(the total dipole moment for SPC water with zwitterionic charges 0.82 e and an OH distance of 1
Å), k as Boltzmann’s constant, and E as 400 MV/m (a typical field magnitude from our MD
simulations) we get a dipole moment which is only 22% of its maximum value. It would require
a field as large as 5.76 GV/m to make the argument unity and even then the dipole moment
would still only be at 31.3% of its maximum value. It turns out that even at very high electric
fields thermal effects are important to consider.
The above calculations assume however that we have an ideal gas which follows the law
of equipartition, assumes that bulk water molecules are free to rotate, unhindered by the presence
of nearby hydrogen bonds, and worst of all assumes that water molecules act just as ideal
dipoles. Furthermore we are not interested in the explicit dynamics of water, but rather on the
temperature effects of lipids and lipid bilayers. The purpose of this comparison then is to give the
151

reader an appreciation for the energies and energy magnitudes introduced through small changes
in temperature. In experiments though we are limited to a small range of physiological
temperatures, so we can rarely justify doubling the temperature in a biological system, however
we are always free to double, or even triple the external electric field since that is an uncoupled
variable we have complete control over. Before we observe how temperature explicitly affects
the discrete, transient stages of the electropore life cycle, we will look first at how temperature
affects the static mechanical properties of intact lipid bilayers in molecular dynamics simulations
without an external field present. This will provide us with the proper context to analyze these
same systems when an electric field is later applied and electroporation subsequently occurs.
Although there are many studies which focus on molecular dynamics simulations of
dipalmitoylphosphatidylcholine (DPPC) [113] or dimyristoylphosphatidylcholine (DMPC) lipids
[69] under various temperatures, few studies have looked at the effects of temperature on
simulated POPC. The few studies that have performed these simulations have done so using
replica-exchange molecular dynamics (REMD) [91] or course-grained lipid models [113]
however those force-fields lack the spatial and temporal resolutions we are looking for. Atomic
models on the other hand are difficult to sample for the determination of thermodynamic
quantities because the time-scales required for proper thermalization are usually on the order of
microseconds, however we are left with few options since we wish to describe fully the atomic
interactions of thermalized lipids with interface water molecules to see how the nanoscale
architecture of electroporation changes as a function of temperature.
To do this we performed simulations of POPC at continuously varying temperatures by
repeatedly lowering and raising the temperature from 37 °C (body temperature) to -20 °C,
similar to the simulated annealing technique. To ensure that we disturbed the system as little as
152

possible we changed temperature slowly at a rate of 2 °C per nanosecond which resulted in
minimal energy fluctuations. This rate is even slower than those used in other studies [91]. We
then measured the area per lipid and temperature of the bilayer as a function of time (Figure 6.1).
We observed a sudden drop in the area per lipid when temperature reached 10 °C, a drop which
is consistent with a first-order phase transition, however this drop was not repeatedly observed
when temperature was raised and lowered again across the 10 °C threshold. This unpredictability
emphasizes the weaknesses of current lipid models. Furthermore the experimental phase
transition temperature for POPC is close to -3 °C [91], not 10 °C. DPPC on the other hand has a
well-known freezing temperature around 30 °C in simulations [26], much higher than POPC.
This large difference in behavior between similar lipid types highlights the large disparities
which currently exist between various molecular force fields. Additionally, we also verified that
our solvent, SPC water, remained in liquid form throughout our simulations since the freezing
temperature of SPC water is well below the experimental value of 0 °C [12].

153

In order to address these instabilities, we tried alternatively allowing each system to
equilibrate for at least 40 ns after the temperature was changed to guarantee convergence in the
area per lipid. When we did this we obtained much cleaner and consistent results (Figure 6.2)
which shows first-order phase transition behavior close to -3 °C. Similarly the rate of the area per
lipid change is different below -3 °C compared to the rate of change above -3 °C, and this is
consistent with other studies [91].

Figure 6.1 Area per lipid and lipid temperature in a POPC bilayer simulation as a function of
time while temperature is continuously varied. Although the area per lipid (blue) suddenly
drops around 10000 ps (10 ns) while the temperature (red) is 10 degrees Celsius - a sign of a
phase transition - this does not appear at later times. This highlights some of the
inconsistencies in the OPLS/Berger POPC force field model, especially since the
experimental phase transition temperature (dotted black line) is closer to -3 degrees Celsius.
154

Other changes to the bilayer that we observed at different temperatures was the
hydrocarbon tail rigidity which was measured with the z-directed order parameter S
z
, a second-
order Legendre polynomial related to the experimental carbon-deuterium order parameter S
cd
.
Values for S
z
can be found in Figure 6.3 and 6.4 for the sn-1 and sn-2 lipid tails which
correspond to the palmitoyl and oleoyl groups respectively. The sn-1 palmitoyl tail, containing
no double bonds, appears to become more rigid across carbons 8 to 16 as temperature drops,
however there is little change in the carbon chain order close to the glycerol backbone. For the
sn-2 oleoyl tail which contains a single double bond we see an increase in tail rigidity across all
carbon atoms after the double bond when temperature is lowered. In both the sn-1 and sn-2 tails
we see a non-linear change in the order parameter between 0 and 5 degrees Celsius which also
suggests a phase transition somewhere in between these temperatures. We also measured the heat
Figure 6.2 Area per lipid versus temperature of a POPC bilayer when the system is allowed
to equilibrate for 40 ns at a target temperature. Note the fluctuations present close to -3
degrees Celsius, a sign perhaps of a first-order phase transition from lipid gel phase to liquid.
155

capacity at constant pressure (Figure 6.5) and the lipid lateral diffusion (Figure 6.6) as a function
of temperature to verify that a phase transition had indeed occurred, however these results were
less conclusive. This is likely because we were sampling on sub-microsecond timescales which
does not allow us to obtain robust statistical averages, and thus a phase transition was not evident
in those measurements.

Figure 6.3 Order parameter S
z
for the sn-1 palmitoyl lipid tail at various temperatures. Note
that ‘carbon number’ i in this diagram refers to carbon atom i+3 in a real lipid. This data
confirms that a phase transition may be occurring near -3 °C as in Figure 6.2.
Figure 6.4 Order parameter S
z
for the sn-2 oleoyl lipid tail at various temperatures. Note that
‘carbon number’ i in this diagram refers to carbon atom i+3 in a real lipid. This data confirms
that a phase transition may be occurring near -3 °C as in Figure 6.2.

156

Figure 6.5 Heat capacity at constant pressure was measured by observing the change in
enthalpy per unit temperature. Evidence of a phase transition at -3 °C is less evident here and
longer sampling times may be required for increased accuracy.

Figure 6.6 POPC lateral diffusion decreases as temperature is lowered. Evidence of a phase
transition at -3 °C is unclear.
157

After each system was equilibrated at a given temperature we ran five trials where we
applied an external electric field of 400 MV/m and measured the corresponding average pore
creation time by looking at the total time required for electropore formation and the individual
times required for each constituent process to occur within pore creation (Figure 6.7 for T < 0 °C
and Figure 6.8 for T � 0 °C). For negative temperatures around -3 °C there is great variability in
the time required for the formation of water bridges (pore initiation) from about 19 nanoseconds
to over 100 nanoseconds. Note that each trial was run only up to 100 nanoseconds, so processes
which occur past 100 ns were not studied in detail here. Pore maturation (later renamed to pore
‘expansion’), the time it took for an electropore to expand from a minimal size to a predefined
‘mature’ size, slowed as the temperature decreased. For temperatures above 0 °C, the pore
creation time as a whole did not appear to change very much from body temperature since the
standard deviation was quite large. Pore expansion however was affected by temperature as one
might expect due to changes in lipid viscosity.

158

Figure 6.7 Total and constituent pore creation times for sub-zero lipid temperatures.

Figure 6.8 Total and constituent pore creation times for lipid temperatures which vary from 0
to 50 degrees Celsius. Pore maturation, a loose measurement of the rate of pore expansion,
appears to be the only temperature-dependent pore creation process.

159

Interestingly our results are in contrast with another study performed by Leekumjorn et
al. [91] which shows that temperatures does in fact modify the pore creation time to a large
extent. The methods of Leekumjorn et al. however are questionable and we believe there are
serious errors in their study which invalidate their conclusions. Although they test the effects of
numerous temperatures on electroporation (from 2-74 °C) they do not explicitly state that they
had equilibrated any of their systems before-hand. Let us suppose however that they did
equilibrate properly but failed to mention it in their paper. They find that over a 72 degree
change in temperature there is only a 0.65 ns change in the pore creation time from 0.8 ns to 0.15
ns with an external electric field of 600 MV/m. In fact in one of our previous studies [93] we
showed that at 600 MV/m the pore creation time at body temperature was 0.4 ns � 0.2 ns, almost
precisely the average of the values found by Leekumjorn et al. It is also important to mention
that they ran only a single trial per temperature tested, and found only a 0.05 ns to 0.2 ns change
in the pore creation time over a 72 degree temperature differential, much smaller than the
standard deviations we found in Figures 6.7 and 6.8. As a result we feel we have adequate
reasons to question their findings, and we believe that our observations are in fact valid ones.
As for pore annihilation and resealing, temperature appears to have quite a large effect on
the timescales required for electropore dissipation. At temperatures above body temperature (T >
37 °C) the pore annihilation time (Figure 6.9) shows a linear decrease in the total pore
annihilation time with similar decreases in all of the constituent steps. On the other hand we find
electropores are incredibly long-lived past 200 nanoseconds for temperatures below body
temperature, quite a contrast to the higher temperature results. This does seem logical however
since thermal energies could in fact dominate the local rearrangement of interface water
molecules when the megavolt-per-meter external electric field is removed. During pore creation
160

when an electric field is present, it’s possible that the large electrostatic contributions overpower
any small changes in temperature and that could explain why temperature does not modify pore
creation by a large amount whereas it does strongly affect the timescales of pore annihilation.
Also perhaps the POPC force field is so heavily based on older DPPC force fields that it retains
some of the 30 degree phase transition properties of its predecessor?

In addition to modifying the bilayer thickness, temperature can also affect the speed at
which pore water molecules diffuse back into the bulk compartment. Figure 6.10 shows a typical
measurement of how quickly pore waters exit the pore interior as a function of time and
temperature. For high temperatures the number of pore waters goes to zero very quickly because
the thermal forces present are largely outweighing the local hydrogen bond forces which
maintain the structural integrity of the electropore. When temperature is dropped to values close
to zero, the number of water molecules which remain present in the pore are large and this may
contribute to the longer pore annihilation times we observed. If the amount of thermal energy
Figure 6.9 At high temperatures pores reseal quickly, but at low temperatures pores appear to
be stable out to 200 ns. Pores at every temperature however shrink to a minimum size when
the field is removed.

161

delivered is less than 44 kJ/mol (the energy required to liberate water molecules from bulk) then
it is entirely possible for an electropore to remain ‘frozen in’ as we have observed for low
temperatures. This result also provides us with a unique metric for future studies in live cells by
observing whether ‘cold’ cells experience long-lived permeabilization after pulsing occurs, a
result which could further validate our existing molecular models.

Figure 6.10 ‘Cold’ electropores take longer to minimize than warm ones, and they also
contain more waters in their interior. Bilayer thickness also increases as the temperature
decreases and lipid tails subsequently become more rigid.
162

6.4 Modulating Pore Radius with Alternating Electric Fields

6.4.1 Methods

Simulations were performed using the GROMACS 4.5.4 software package on the
University of Southern California High Performance Computing and Communications Linux
cluster. Bilayers consisted of either 128 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine
(POPC) lipids or 108 POPC lipids mixed with 20 1-palmitoyl-2-oleoyl-sn-glycero-3-
phosphatidylserine (POPS) lipids, all of which were parameterized by OPLS/Berger force fields
and combined with 70 SPC water molecules per lipid at an integration time-step of 2 fs. Pore life
cycle times were determined using custom Perl scripts as previously reported [93]. NPT
ensembles were maintained using the Berendsen barostat which applied 1 bar in xy:z semi-
isotropically with a compressibility of 4.5E-5 bar
-1
. Pore radius measurements were obtained by
identifying a pore axis which passes directly through the pore in the z-dimension (based on
perpendicular x and y water density profiles). Then, bins were assigned between mean lipid
phosphorus planes where, for each bin, maximally-distant water molecules were identified
relative to the pore axis in x and y to obtain local semi-major and semi-minor pore diameters.
Following this the semi-major and semi-minor pore diameters were averaged together in each
bin, and finally all local pore diameters were averaged together to obtain an average pore radius
or volume. Periodic boundary conditions were employed to mitigate system size effects.
Time-dependent electric pulses were applied by modifying the GROMACS source code
and replacing the built-in constant electric field function with a custom routine which delivers
electric pulses of the form:

163

𝐸 𝑧 = 𝐸 0
∗ 𝑐𝑜𝑠 ( 𝜔𝑡 ) ∗ Θ(cos( 𝜔 𝑧 𝑡 )) (6.3)

Here Θ is the Heavyside step function which selects only positive unipolar field amplitudes, t is
time in picoseconds, and 𝜔 𝑧 is the angular frequency for electric fields in the z-direction which
is perpendicular to the plane of the bilayer. In this study we used an angular pulse frequency of
0.00209 radians per second which corresponded to a unipolar pulse width of 1.5 nanoseconds (or
3 nanoseconds bipolar). Measurements of pore radius over time were fit using MATLAB’s
Equations and Systems Solver to obtain optimal coefficients for predefined functional formats.

6.4.2 Results and Discussion

Initial interest in this study came from an early collaboration with a group in Buenos
Aires [44] where lipid bilayer electropores were created in simulations with large external
electric fields, then stabilized using smaller electric fields which stopped pores from growing
(Figure 6.11). These ‘stabilized pores’ were energy minimized structures of constant radius
which depended on the magnitude of the reduced external electric field. The presence of these
pores allowed us to compare theoretical ‘stable’ pore radii in analytical models to our
observations in MD, and often our results were consistent with models, specifically the
asymptotic model of electroporation [86].

164

One problem we faced however was that most continuum models did not focus on the
dynamics of only a single electropore, but rather on the effects of a large distribution of pores,
thus necessitating more observations to closely tie theoretical models to our simulations. Instead
we focused on measuring electropore properties dynamically and we were curious how varying
or time-dependent electric fields would continuously perturb the membrane. Using these fields
would allow us to continuously sample a statistically-significant number of electropores which
we could use to directly compare simulations to models. Many of these continuum models also
characterize electropores by the membrane properties which they occur in such as the membrane
Figure 6.11 The creation of ‘stabilized electropores’ in MD from a collaboration with
Fernández et al. [44]. Electropores were created with an electric field of 750 MV/m, then
stabilized by quickly reducing the field once the pore reached a standard size. Note that there
appears to be a ‘solid-to-liquid-like’ transition between pore-stabilizing fields under 200
MV/m where pore structures are stable, and fields above 200 MV/m where the pore radius
grows indefinitely, analogous to the root-mean-square displacement of particles in a lattice at
increasing temperatures.
165

surface tension, edge tension, lipid viscosity, diffusion coefficient, and steric repulsion force
[86]. All of these are physical properties that are possible to measure in simulations during
electroporation, and the use of time-dependent electric fields also makes it possible to extract the
rate of pore expansion or contraction as a function of these static membrane properties.
Before we settled on the electric pulse function specified in Equation 6.3 we tried many
different functional forms to observe how frequency, amplitude, and pulse shape affect the
electroporation process. Although other studies [146, 150] have thoroughly characterized
mechanical pore formation in vesicles based on continuum models that depend only on water
viscosity or lipid surface tension, no analytical models exist (to the author’s knowledge) which
characterize electropore formation as a function of the external pulse shape or frequency.
We began our studies with bipolar sinusoidal pulses (Figure 6.12) of varying amplitude
and period. We wanted to observe if pores formed more preferentially in high frequency systems.
For all frequencies sampled however, pores formed only when the field was at maximum
amplitude (Figure 6.13) independent of the frequencies used to obtain that amplitude. Thus the
pulse spacing had little effect on the bilayer other than deciding when the maximum field was
delivered. Also because the lipid bilayer is symmetric the field can induce permeabilization
equally at negative or positive polarities. As Figure 6.13 shows, although the polarization of
water is net positive or negative during pore formation, it can rotate by 180 degrees within the
pore if the field requires it, and still continue to provide the necessary aqueous flux in order to
keep the pore open. Amplitude then has the same effect as it did in constant fields where larger
amplitudes created electropores exponentially quicker than smaller amplitudes, consistent with
the stochastic model of electroporation [86]. We also tested triangle pulses (data not shown) but
the results were identical to the data shown for sinusoidal simulations.
166

Figure 6.12 A typical time-dependent bipolar electric field implemented in GROMACS.
Figure 6.13 The net water-dipole response of a POPC bilayer simulation during application
of a time-dependent electric field drawn in Figure 6.12. Electroporation (depicted in this
figure as |polarization| > 1000 debye) occurs at both positive and negative field polarities.
167

We were also concerned that because electropores formed at both positive and negative
polarities, the bilayer was not getting sufficient time to re-polarize in between high-amplitude
pulses. We decided then to fix the period at 3 ns and implement the step function present in
Equation 6.3 in order to select only for positive field polarities. Using unipolar pulses (Figure
6.14) also made it easier to compare our systems to real vesicle or cell experiments in the lab
because all of our experimental electric pulses are unipolar. This introduced adequate time for
the system to re-polarize to pre-pulse levels in between subsequent electropulses (Figure 6.15).

Figure 6.14 A typical time-dependent unipolar electric pulse implemented in GROMACS.
168

As for pulse amplitude, we wanted to use a field which would smoothly open and close
electropores in MD such that our pore radius identification algorithm (described in the methods
section) output smooth datapoints. We tested field amplitudes from 150 MV/m to 500 MV/m,
with bipolar (Figure 6.16A) and unipolar pulses (Figure 6.16B), and found that 400 MV/m,
unipolar pulses resulted in the cleanest pore radius data which also matched up well with the
ideal pulse shape in Figure 6.14. The 400 MV/m pulse also allows electropores in the bilayer to
shrink enough over the 3 ns period (or 1.5 ns ‘field-off’ time) that the average pore radius
remains constant over time. For instance, when a pulse amplitude of 500 MV/m is used instead
the average pore radius increases as time goes on which is undesirable for long-term simulations.

Figure 6.15 The net water-dipole response of a POPC bilayer simulation during application
of a time-dependent electric field drawn in Figure 6.14. Note that the dipole moment profile
follows the electric field profile almost exactly.
169

Figure 6.16 Measurements of the MD time-dependent electropore radius in a POPC
membrane as a function of pulse amplitude and pulse shape (bipolar (A), unipolar (B)).
170

In many ways the functional form of the pore radius for unipolar pulses resembles the
voltage output across an RC circuit with a diode attached to it (Figure 6.17). This analogy of the
membrane being a capacitive electrical device was actually made in Chapter 1 since many
studies [1, 29, 199, 200, 201] have used simple electrical analogs to model the electroporation
process. In the case of the RC circuit if we apply a sinusoidally-varying signal, the positive
voltage will pass through the diode and charge the capacitor over some time-constant τ. When
the sinusoidal source becomes negative the diode stops the transmission of current to the
capacitor, and the capacitor discharges across the resistor, also with a time constant τ. The
resulting voltage across the points on the right-side of Figure 6.17 would be an exponential
increase in voltage with a time constant τ, then an exponential drop in voltage with time constant
τ. This is surprisingly similar to the dynamics observed in Figure 6.16B, and if we model an
analogous ‘time constant’ using Equation 6.4, we find that τ
pore
= 3.33 ns. In other words the
pore radius will increase or decrease by a factor of e over 3.33 ns.

𝑅 = 𝑅 𝑝 𝑒 𝑎𝑘
𝑒 𝑥𝑝 ( −
𝑡 𝜏 𝑝 𝑜 𝑟𝑒
) (6.4)
Figure 6.17 A simple DC circuit with a diode attached which results in output signals similar
to the pore radii observed in Figure 6.16B.
171

Now that we have a standardized pulse period (3 ns), amplitude (400 MV/m), and shape
(unipolar), we can really begin to test how closely electroporation in simulations resembles
continuum models. Figures 6.18 and 6.19 show a side view (without visible lipids) and top view
(without visible water) of a typical electropore sequence to give the reader an idea for how
dynamic electropores created with time-dependent electric pulses can be. Quite visibly, time-
dependent electropores ‘breathe’ with each peak of the electric pulse, enlarging when the
amplitude is large and contracting when the amplitude is small. Water and hydrated lipids alike
move in unison as the pore expands or contracts, and this movement allows us to measure the
rate of pore expansion (or contraction) as a function of time which is directly comparable to
continuum models [86] that provide the steric and diffusive properties of the membrane.

Figure 6.18 A side view of a ‘breathing’ electropore as a time-dependent electric field is
applied. Lipid carbon atoms are not shown for simplicity. Colored atoms are red for both
phosphate and glycerol oxygens, blue for amino nitrogen, and gold for phosphate phosphorus.
172

The stochastic model of electroporation [86] (sometimes called the asymptotic model of
pore formation) tells us that the rate of pore expansion follows a simple functional form:
𝑑𝑟 𝑑𝑡 = 𝜶 𝑉 2
( 𝑟 + 𝜷 )
( 𝑟 + 𝜸 )
+ 𝜹 1
𝑟 5
− 𝜺 + 𝝃 𝑟 (6.5)
where V is voltage, r is pore radius, t is time, and all other variables are constants of the
membrane. The first term represents the local electric force which originates from the constant
external transmembrane potential. In our system V = V
o
*cos(ωt)*Θ(ωt), essentially a voltage
analog to Equation 6.3. The second, third, and fourth terms are attributed to steric repulsion
Figure 6.19 A top-down view of a ‘breathing’ electropore through the lipid bilayer as a time-
dependent electric field is applied. Water is not shown for simplicity. Colored atoms are
phosphate and glycerol oxygen (red), amino nitrogen (dark blue), phosphate phosphorus
(gold) and carbon (cyan).
173

between lipids, line tension present on the pore perimeter, and surface tension of the surrounding
bilayer respectively.
To extract these coefficients from our simulations we need to smooth the pore radius data
since the rate of pore expansion, dr/dt, is a first derivative of the radius, thus it should be a
continuous and well-behaved function. A moving average of pore radius produces a smooth,
well-defined curve (Figure 6.20) which is sufficient for our analyses. From the smoothed pore
radius we take the derivative (dr/dt) and use numerical approximations to extract the best-fit
solution from Equation 6.5. From our data we find very good agreement with theory (Figure
6.21) which allows for the extraction of α, β, γ, δ, ε, and ξ. Our simulations reveal that α =
0.0195 nm/(ns·V
2
), β = 0.31 nm, γ = 0.97 nm, δ = 1.05 nm
6
/ns, ε = 0.011 nm/ns, and ξ = -0.01
ns
-1
. The R-squared value of this fit is 0.9 which signifies a strong correlation between
theoretical models and our observations in MD.

Figure 6.20 A measure of the simulated electropore radius in a POPC bilayer as a function of
time while an external, time-dependent electric field (Figure 6.14) is applied.
174

Theoretical values for these coefficients are, for the most part, similar to the coefficients
obtained from simulations. Using standard input values obtained from Krassowska et al. [86] we
obtain theoretical values of α = 0.0082 nm/(ns·V
2
), β = 0.31 nm, γ = 1.28 nm, δ = 4.42E-4
nm
6
/ns, ε = 0.0013 nm/ns, and ξ = 7.34E-8 ns
-1
. α, β, γ, and ε all appear to be within the same
order of magnitude as the results we obtained through simulations, however δ and ξ are three-
orders of magnitude and six-orders of magnitude off respectively from simulation values. δ
represents the steric lipid repulsion term and it is perhaps not surprising that the OPLS/Berger
force field for lipids in MD gives different repulsion energies than experiment since most
dominant lipid models only guarantee proper lipid density, charge, area per lipid, net dipole
moment, and in some cases electrophoretic mobility [15]. ξ however, the term representing the
surface tension of the bilayer, should be reasonably consistent with simulations. To address this
Figure 6.21 A measure of the rate of simulated electropore expansion in a POPC bilayer as a
function of time while an external, time-dependent electric field is applied. The theoretical
pore advection velocity [86], an analog to the rate of pore expansion, is plotted alongside our
measurements from MD.
175

we found that the calculation in Krassowska et al. for ξ is only valid for electropores which form
on whole cells, thus having an available surface area of A ≈ 4 𝜋 (5 𝐸 − 5 𝑚 )
2
. In our simulations
however we only have a small periodic box with surface area A’ = 49 nm
2
, so when we
substitute A’ for A in the calculation of ξ from Krassowska et al. we obtain ξ = -0.0004 ns
-1
. This
is in much better agreement with the simulation data.

6.5 Conclusions

Small changes in the properties of lipids or lipid bilayers can produce a large change in
the dynamics and self-assembly of phospholipid bilayer electropores. We find that at lower lipid
temperatures which are close to the phase-transition temperature, electropore expansion occurs
more slowly than at higher temperatures, and this dependence appears linear. Temperature did
not have a large effect on the pore initiation time - the time required to construct water bridges
across the membrane, even though many structural changes were observed in the bilayer. The
external electric field on the other hand appears to be the dominant variable which controls pore
initiation, and this is consistent with the stochastic pore hypothesis which states that the rate of
pore formation is proportional to exp(
𝑉 2
𝑘𝑇
). For pore annihilation, temperature had a very large
effect on every stage of pore resealing, and for temperatures of 30 ˚C and lower, electropores
appeared to be stable without an external field present for over 200 nanoseconds. At
temperatures higher than 30 degrees the pore annihilation time appears to decrease linearly as the
temperature increases. Also when the field is removed at lower temperatures after pore creation
176

occurs, the pore radius and the number of water present inside the pore tend to reach a metastable
minimum more slowly than systems at higher lipid temperatures.
Similarly, time-dependent electric fields can be used to control the size and expansion
rate of phospholipid bilayer electropores, and simulations of time-dependent electropores can
facilitate the extraction of a multitude of continuum variables as a basis for comparing molecular
dynamics results to theoretical models. We find good agreement between our simulations and
dominant theoretical models, except for measurements of steric lipid repulsion. Pulse amplitudes
in time-dependent electric fields have a similar effect on pore creation as they did for constant
electric fields, that is they cause electroporation to occur at an exponential rate as the amplitude
grows linearly. Pulse period did not appear to directly affect the pore creation time, however it
does indirectly affect when the maximum electric field is delivered to the bilayer which can
subsequently speed-up or slow-down electroporation. Lipid electropores also appear to oscillate
in time-dependent electric fields in a way which is strikingly similar to RC circuit outputs.
Electrical continuum models were quite accurate at describing the atomic lipid properties that are
present in molecular dynamics simulations during time-dependent phospholipid bilayer
electroporation, a key step in aligning macroscopic theories to discrete computer simulations.

6.6 Acknowledgements

Computation for this work was supported by the University of Southern California Center for
High-Performance Computing and Communications (www.usc.edu/hpcc). Special thanks to
MOSIS for providing additional resources and funding for this project.

177

Conclusions and Outlook

In this dissertation we have shown that performing molecular dynamics simulations of
phospholipid bilayer electroporation is a valid way to study the atomic physical and chemical
interactions which take place during cell membrane electropermeabilization. The use of highly-
paralleled computing clusters and scalable integration techniques allowed for a large number of
statistical studies to be performed at the nanoscale-level which helped bridge the gap between
ideal continuum models and real biological experiments using only ab-initio numerical
simulations. Often these simulations provided unique details about real systems which were
unavailable to experimental studies such as the extraction of molecular free energies or the
summation of electric potentials which gave rise to spatial and temporal perturbations. Together
with robust mathematical modeling, molecular dynamics simulations successfully answered
questions about physical, chemical, and biological phenomena by recreating complex many-body
dynamics from intuitive, first-principle foundations.
These studies found that oxidized lipid bilayers electropermeabilized much more readily
than unoxidized bilayers, a result which could facilitate the development of new
electrochemotherapy techniques for the treatment of cancer through targeted, low-amplitude
electropermeabilization of pretreated, locally-oxidized tumors. These results were confirmed
through experiments of peroxidized live cells which exhibited enhanced electropermeability
compared to unoxidized cells. Molecular dynamics simulations suggested that the cause of this
enhanced permeability was due to the presence of carbonyl oxygen atoms on the oxidized lipid
tail which appeared to mediate additional water transport through the membrane interior.
178

We also investigated the life cycle of an electropore, a standardized temporal sequence of
events which was observed in all electroporation trajectories. This allowed, for the first time,
objective normalization of lipid electropore dynamics so that additional studies could
quantitatively differentiate two physically-distinct pore structures from one another. From this
we found that the pore creation rate in MD depends on the external electric field exponentially, a
result which is consistent with continuum models of electroporation. We also observed that pore
construction, the reorientation of lipids towards the exterior of the electropore, occurs in constant
time independent of the magnitude of the external electric field. Pore annihilation, or resealing,
appears independent of the pore-initiating electric field, pointing instead to stochastic resealing.
Further utilizing electropore life cycle as an objective metric, we found that
palmitoyloleoylphosphatidylserine (a negative-charged lipid) and calcium (a divalent cation)
increases the time required for electroporation to occur and subsequently decreases the time
required for pore resealing. This implies that systems with negatively-charged interfaces, high
surface tensions, and/or thicker membranes exhibit pore creation dynamics which are more time-
consuming while existing electropores in these systems dissipate quickly when the external
electric field is removed. Calcium binding kinetics in our simulations also appeared to agree well
with experimental and theoretical models which use a 1:2 Langmuir binding isotherm to
characterize the binding affinity of ions.
Finally, we studied how individual changes in the physical properties of water and lipid
affect electroporation. Simple water models of electroevaporation provided similar energetics
and structural precursors to water bridge formation when compared to lipid bilayer
electroporation. Observations in MD of atomic water bridge construction suggest that so-called
‘Okuno-Tanioka’ water pedestals facilitate the transition of intact bilayers from one geometric
179

state to another, more energetically-favorable geometric state – the electropore. Lipids on the
other hand, secondary structures in water bridge construction, are also capable of perturbing
electroporation dynamics. Changes in lipid temperature can slow various pore creation processes
such as the expansion of pores once they are created. When the external field is removed, local
thermal fluctuations become the dominant forces which determine the time for subsequent
bilayer resealing. Lipid bilayers can also be manipulated by time-dependent electric fields which
allow for various transient processes to be observed and quantified for comparison to continuum
models. These models validated many of the force fields used in our molecular dynamics
simulations and created a baseline agreement for future simulations to strengthen, thus
supporting the use of further MD studies as reliable and useful tools for biological studies.
Although there is much to improve upon, molecular dynamics approximations provide an
excellent way to link experimental electroporation observations to first-principle physical
models, and it will continue to make great strides as computer hardware and software becomes
quicker and more nimble in the coming decades. Studies of the mechanisms of electroporation
are still in their infancy, and there is much to be discovered by continuing simulations of local
electric interactions with biological systems. As we become more and more dependent on the
utilization of electromagnetic radiation for engineering and medicine, understanding the complex
interactions between biological and physical systems becomes a priority, and I sincerely hope
that future generations of researchers pick up where we (and others) have left off during what is
likely the beginning of a new era in computational biophysics.

180

181

Bibliography

1 I. G. Abidor, V. B. Arakelyan, L. V. Chernomordik, Y. A. Chizmadzhev, V. F.
Pastushenko, and M. R. Tarasevich, "Electric Breakdown of Bilayer Lipid-Membranes
.1. Main Experimental Facts and Their Qualitative Discussion," Bioelectrochemistry and
Bioenergetics 6 (1), 37-52 (1979).

2 J. Alper, "Drug delivery. Breaching the membrane," Science 296 (5569), 838-839 (2002).

3 C. Altenbach and J. Seelig, "Ca-2+ Binding to Phosphatidylcholine Bilayers as Studied
by Deuterium Magnetic-Resonance - Evidence for the Formation of a Ca-2+ Complex
with 2 Phospholipidmolecules," Biochemistry 23 (17), 3913-3920 (1984).

4 F. M. Andre, M. A. Rassokhin, A. M. Bowman, and A. G. Pakhomov, "Gadolinium
blocks membrane permeabilization induced by nanosecond electric pulses and reduces
cell death," Bioelectrochemistry 79 (1), 95-100 (2010).

5 W. Arap, R. Pasqualini, and E. Ruoslahti, "Cancer treatment by targeted drug delivery to
tumor vasculature in a mouse model," Science 279 (5349), 377-380 (1998).

6 C. Aussel, R. Marhaba, C. Pelassy, and J. P. Breittmayer, "Submicromolar La3+
concentrations block the calcium release-activated channel, and impair CD69 and CD25
expression in CD3- or thapsigargin-activated Jurkat cells," Biochemical Journal 313,
909-913 (1996).

7 A. G. Ayuyan and F. S. Cohen, "Lipid peroxides promote large rafts: Effects of excitation
of probes in fluorescence microscopy and electrochemical reactions during vesicle
formation," Biophysical Journal 91 (6), 2172-2183 (2006).

8 V. A. Bakaev and W. A. Steele, "On the computer simulation of a hydrophobic vitreous
silica surface," The Journal of Chemical Physics 111 (21), 9803-9812 (1999).

9 S. J. Beebe, P. M. Fox, L. J. Rec, L. K. Willis, and K. H. Schoenbach, "Nanosecond,
high-intensity pulsed electric fields induce apoptosis in human cells," Faseb Journal 17
(9), 1493-1495 (2003).

10 W. F. Bennett, J. L. MacCallum, and D. P. Tieleman, "Thermodynamic analysis of the
effect of cholesterol on dipalmitoylphosphatidylcholine lipid membranes," J Am Chem
Soc 131 (5), 1972-1978 (2009).

11 R. Benz and U. Zimmermann, "Pulse-Length Dependence of the Electrical Breakdown in
Lipid Bilayer-Membranes," Biochimica Et Biophysica Acta 597 (3), 637-642 (1980).

182

12 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and J. Hermans, "Interaction
models for water in relation to protein hydration," Intermolecular Forces. 331-342
(1981).

13 H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, "The Missing Term in Effective
Pair Potentials," Journal of Physical Chemistry 91 (24), 6269-6271 (1987).

14 H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola, and J. R. Haak,
"Molecular-Dynamics with Coupling to an External Bath," Journal of Chemical Physics
81 (8), 3684-3690 (1984).

15 O. Berger, O. Edholm, and F. Jahnig, "Molecular dynamics simulations of a fluid bilayer
of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant
temperature," Biophysical Journal 72 (5), 2002-2013 (1997).

16 R. A. Bockmann, B. L. de Groot, S. Kakorin, E. Neumann, and H. Grubmuller, "Kinetics,
statistics, and energetics of lipid membrane electroporation studied by molecular
dynamics simulations," Biophysical Journal 95 (4), 1837-1850 (2008).

17 R. A. Bockmann and H. Grubmuller, "Multistep binding of divalent cations to
phospholipid bilayers: A molecular dynamics study," Angewandte Chemie-International
Edition 43 (8), 1021-1024 (2004).

18 J. M. Boettcher, R. L. Davis-Harrison, M. C. Clay, A. J. Nieuwkoop, Y. Z. Ohkubo, E.
Tajkhorshid, J. H. Morrissey, and C. M. Rienstra, "Atomic View of Calcium-Induced
Clustering of Phosphatidylserine in Mixed Lipid Bilayers," Biochemistry 50 (12), 2264-
2273 (2011).

19 S. Boitano, M. J. Sanderson, and E. R. Dirksen, "A Role for Ca2+-Conducting Ion
Channels in Mechanically-Induced Signal-Transduction of Airway Epithelial-Cells,"
Journal of Cell Science 107, 3037-3044 (1994).

20 P. Bonnafous, M. C. Vernhes, J. Teissie, and B. Gabriel, "The generation of reactive-
oxygen species associated with long-lasting pulse-induced electropermeabilisation of
mammalian cells is based on a non-destructive alteration of the plasma membrane,"
Biochimica Et Biophysica Acta-Biomembranes 1461 (1), 123-134 (1999).

21 D. Bostick and M. L. Berkowitz, "The implementation of slab geometry for membrane-
channel molecular dynamics simulations," Biophysical Journal 85 (1), 97-107 (2003).

22 M. Brielmeier, J. M. Bechet, M. H. Falk, M. Pawlita, A. Polack, and G. W. Bornkamm,
"Improving stable transfection efficiency: antioxidants dramatically improve the
outgrowth of clones under dominant marker selection," Nucleic Acids Research 26 (9),
2082-2085 (1998).

183

23 B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M.
Karplus, "Charmm - a Program for Macromolecular Energy, Minimization, and
Dynamics Calculations," Journal of Computational Chemistry 4 (2), 187-217 (1983).

24 G. Bussi, D. Donadio, and M. Parrinello, "Canonical sampling through velocity
rescaling," Journal of Chemical Physics 126 (1) (2007).

25 M. Cemazar, T. Jarm, D. Miklavcic, A. M. Lebar, A. Ihan, N. A. Kopitar, and G. Sersa,
"Effect of electric-field intensity on electropermeabilization and electrosensitivity of
various tumor-cell lines in vitro," Electro- and Magnetobiology 17 (2), 263-272 (1998).

26 R. Chen, D. Poger, and A. E. Mark, "Effect of High Pressure on Fully Hydrated DPPC
and POPC Bilayers," Journal of Physical Chemistry B 115 (5), 1038-1044 (2011).

27 L. V. Chernomordik, S. I. Sukharev, I. G. Abidor, and Y. A. Chizmadzhev, "Breakdown
of Lipid Bilayer-Membranes in an Electric-Field," Biochimica Et Biophysica Acta 736
(2), 203-213 (1983).

28 R. Chiaramonte, E. Bartolini, P. Riso, E. Calzavara, D. Erba, G. Testolin, P. Comi, and
G. V. Sherbet, "Oxidative stress signalling in the apoptosis of jurkat T-lymphocytes,"
Journal of Cellular Biochemistry 82 (3), 437-444 (2001).

29 Y. A. Chizmadzhev and I. G. Abidor, "Bilayer Lipid-Membranes in Strong Electric-
Fields," Bioelectrochemistry and Bioenergetics 7 (1), 83-100 (1980).

30 M. M. A. E. Claessens, F. A. M. Leermakers, F. A. Hoekstra, and M. A. C. Stuart,
"Osmotic shrinkage and reswelling of giant vesicles composed of
dioleoylphosphatidylglycerol and cholesterol," Biochimica Et Biophysica Acta-
Biomembranes 1778 (4), 890-895 (2008).

31 R. J. Clarke, "The dipole potential of phospholipid membranes and methods for its
detection," Adv Colloid Interfac 89, 263-281 (2001).

32 K. S. Cole, "Electric impedance of marine egg membranes.," T Faraday Soc 33 (2), 0966-
0971 (1937).

33 J. A. Cook, D. Gius, D. A. Wink, M. C. Krishna, A. Russo, and J. B. Mitchell, "Oxidative
stress, redox, and the tumor microenvironment," Seminars in Radiation Oncology 14 (3),
259-266 (2004).

34 H. G. L. Coster, "A Quantitative Analysis of Voltage-Current Relationships of Fixed
Charge Membranes and Associated Property of Punch-Through," Biophysical Journal 5
(5), 669-& (1965).

184

35 K. A. DeBruin and W. Krassowska, "Electroporation and shock-induced transmembrane
potential in a cardiac fiber during defibrillation strength shocks," Annals of Biomedical
Engineering 26 (4), 584-596 (1998).

36 K. A. DeBruin and W. Krassowska, "Modeling electroporation in a single cell. I. Effects
of field strength and rest potential," Biophysical Journal 77 (3), 1213-1224 (1999).

37 K. A. DeBruin and W. Krassowska, "Modeling electroporation in a single cell. II. Effects
of ionic concentrations," Biophysical Journal 77 (3), 1225-1233 (1999).

38 W. Dyrka, A. T. Augousti, and M. Kotulska, "Ion flux through membrane channels - An
enhanced algorithm for the Poisson-Nernst-Planck model," Journal of Computational
Chemistry 29 (12), 1876-1888 (2008).

39 E. Egberts, S. J. Marrink, and H. J. C. Berendsen, "Molecular-Dynamics Simulation of a
Phospholipid Membrane," European Biophysics Journal with Biophysics Letters 22 (6),
423-436 (1994).

40 M. Eisenberg, J. E. Hall, and C. A. Mead, "The nature of the voltage-dependent
conductance induced by alamethicin in black lipid membranes," J Membr Biol 14 (2),
143-176 (1973).

41 U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, "A
Smooth Particle Mesh Ewald Method," Journal of Chemical Physics 103 (19), 8577-8593
(1995).

42 D. M. Ferguson, "Parameterization and Evaluation of a Flexible Water Model," Journal
of Computational Chemistry 16 (4), 501-511 (1995).

43 M. L. Fernandez, G. Marshall, F. Sagues, and R. Reigada, "Structural and Kinetic
Molecular Dynamics Study of Electroporation in Cholesterol-Containing Bilayers,"
Journal of Physical Chemistry B 114 (20), 6855-6865 (2010).

44 M. L. Fernandez, M. Risk, R. Reigada, and P. T. Vernier, "Size-controlled nanopores in
lipid membranes with stabilizing electric fields," Biochem Biophys Res Commun 423
(2), 325-330 (2012).

45 W. Frey, J. A. White, R. O. Price, P. F. Blackmore, R. P. Joshi, R. Nuccitelli, S. J. Beebe,
K. H. Schoenbach, and J. F. Kolb, "Plasma membrane voltage changes during
nanosecond pulsed electric field exposure," Biophysical Journal 90 (10), 3608-3615
(2006).

46 B. Gabriel and J. Teissie, "Generation of Reactive-Oxygen Species Induced by
Electropermeabilization of Chinese-Hamster Ovary Cells and Their Consequence on Cell
Viability," European Journal of Biochemistry 223 (1), 25-33 (1994).

185

47 B. Gabriel and J. Teissie, "Direct observation in the millisecond time range of fluorescent
molecule asymmetrical interaction with the electropermeabilized cell membrane,"
Biophysical Journal 73 (5), 2630-2637 (1997).

48 B. Gabriel and J. Teissie, "Fluorescence imaging in the millisecond time range of
membrane electropermeabilisation of single cells using a rapid ultra-low-light
intensifying detection system," European Biophysics Journal with Biophysics Letters 27
(3), 291-298 (1998).

49 E. B. Garon, D. Sawcer, P. T. Vernier, T. Tang, Y. H. Sun, L. Marcu, M. A. Gundersen,
and H. P. Koeffler, "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 (3), 675-682 (2007).

50 J. Gehl, "Electroporation: theory and methods, perspectives for drug delivery, gene
therapy and research," Acta physiologica Scandinavica 177 (4), 437-447 (2003).

51 A. W. Girotti, "Mechanisms of lipid peroxidation," J Free Radic Biol Med 1 (2), 87-95
(1985).

52 R. W. Glaser, S. L. Leikin, L. V. Chernomordik, V. F. Pastushenko, and A. I. Sokirko,
"Reversible Electrical Breakdown of Lipid Bilayers - Formation and Evolution of Pores,"
Biochimica Et Biophysica Acta 940 (2), 275-287 (1988).

53 D. E. Goldman, "Potential, impedance, and rectification in membranes," J Gen Physiol 27
(1), 37-60 (1944).

54 M. Golzio, M. P. Mora, C. Raynaud, C. Delteil, J. Teissie, and M. P. Rols, "Control by
osmotic pressure of voltage-induced permeabilization and gene transfer in mammalian
cells," Biophysical Journal 74 (6), 3015-3022 (1998).

55 M. Golzio, J. Teissie, and M. P. Rols, "Direct visualization at the single-cell level of
electrically mediated gene delivery," Proceedings of the National Academy of Sciences
of the United States of America 99 (3), 1292-1297 (2002).

56 P. F. B. Goncalves and H. Stassen, "Calculation of the free energy of solvation from
molecular dynamics simulations," Pure Appl Chem 76 (1), 231-240 (2004).

57 T. R. Gowrishankar and J. C. Weaver, "An approach to electrical modeling of single and
multiple cells," Proceedings of the National Academy of Sciences of the United States of
America 100 (6), 3203-3208 (2003).

58 A. A. Gurtovenko and I. Vattulainen, "Pore formation coupled to ion transport through
lipid membranes as induced by transmembrane ionic charge imbalance: atomistic
molecular dynamics study," J Am Chem Soc 127 (50), 17570-17571 (2005).

186

59 A. A. Gurtovenko and I. Vattulainen, "Calculation of the electrostatic potential of lipid
bilayers from molecular dynamics simulations: Methodological issues," Journal of
Chemical Physics 130 (21), - (2009).

60 W. A. Hamilton and A. J. H. Sale, "Effects of High Electric Fields on Microorganisms .2.
Mechanism of Action of Lethal Effect," Biochimica Et Biophysica Acta 148 (3), 789-800
(1967).

61 R. L. Harrison, B. J. Byrne, and L. Tung, "Electroporation-mediated gene transfer in
cardiac tissue," FEBS Letters 435 (1), 1-5 (1998).

62 R. Heller, M. J. Jaroszeski, L. F. Glass, J. L. Messina, D. P. Rapaport, R. C. DeConti, N.
A. Fenske, R. A. Gilbert, L. M. Mir, and D. S. Reintgen, "Phase I/II trial for the treatment
of cutaneous and subcutaneous tumors using electrochemotherapy," Cancer 77 (5), 964-
971 (1996).

63 B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije, "LINCS: A linear
constraint solver for molecular simulations," Journal of Computational Chemistry 18
(12), 1463-1472 (1997).

64 B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, "GROMACS 4: Algorithms for
highly efficient, load-balanced, and scalable molecular simulation," Journal of Chemical
Theory and Computation 4 (3), 435-447 (2008).

65 M. Hibino, M. Shigemori, H. Itoh, K. Nagayama, and K. Kinosita, "Membrane
Conductance of an Electroporated Cell Analyzed by Submicrosecond Imaging of
Transmembrane Potential," Biophysical Journal 59 (1), 209-220 (1991).

66 S. E. Hickman, J. Elkhoury, S. Greenberg, I. Schieren, and S. C. Silverstein, "P2z
Adenosine-Triphosphate Receptor Activity in Cultured Human Monocyte-Derived
Macrophages," Blood 84 (8), 2452-2456 (1994).

67 Y. Hirano, N. Okimoto, I. Kadohira, M. Suematsu, K. Yasuoka, and M. Yasui,
"Molecular Mechanisms of How Mercury Inhibits Water Permeation through Aquaporin-
1: Understanding by Molecular Dynamics Simulation," Biophysical Journal 98 (8), 1512-
1519 (2010).

68 J. F. Hoffman, "Physiological Characteristics of Human Red Blood Cell Ghosts," J Gen
Physiol 42 (1), 9-28 (1958).

69 C. J. Hogberg and A. P. Lyubartsev, "A molecular dynamics investigation of the
influence of hydration and temperature on structural and dynamical properties of a
dimyristoylphosphatidylcholine bilayer," Journal of Physical Chemistry B 110 (29),
14326-14336 (2006).

187

70 P. Hojman, H. Gissel, F. M. Andre, C. Cournil-Henrionnet, J. Eriksen, J. Gehl, and L. M.
Mir, "Physiological Effects of High- and Low-Voltage Pulse Combinations for Gene
Electrotransfer in Muscle," Human Gene Therapy 19 (11), 1249-1260 (2008).

71 Q. Hu, S. Viswanadham, R. P. Joshi, K. H. Schoenbach, S. J. Beebe, and P. F.
Blackmore, "Simulations of transient membrane behavior in cells subjected to a high-
intensity ultrashort electric pulse," Phys Rev E Stat Nonlin Soft Matter Phys 71 (3 Pt 1),
031914 (2005).

72 P. Huang and G. H. Loew, "Interaction of an Amphiphilic Peptide with a Phospholipid-
Bilayer Surface by Molecular-Dynamics Simulation Study," Journal of Biomolecular
Structure & Dynamics 12 (5), 937-956 (1995).

73 P. Huang, J. J. Perez, and G. H. Loew, "Molecular-Dynamics Simulations of
Phospholipid-Bilayers," Journal of Biomolecular Structure & Dynamics 11 (5), 927-956
(1994).

74 J. S. Hub, C. Aponte-Santamaria, H. Grubmuller, and B. L. de Groot, "Voltage-Regulated
Water Flux through Aquaporin Channels In Silico," Biophysical Journal 99 (12), L97-
L99 (2010).

75 W. Humphrey, A. Dalke, and K. Schulten, "VMD: Visual molecular dynamics," Journal
of Molecular Graphics 14 (1), 33-38 (1996).

76 T. Idziorek, J. Estaquier, F. Debels, and J. C. Ameisen, "Yopro-1 Permits
Cytofluorometric Analysis of Programmed Cell-Death (Apoptosis) without Interfering
with Cell Viability," Journal of Immunological Methods 185 (2), 249-258 (1995).

77 R. F. Jacob and R. P. Mason, "Lipid peroxidation induces cholesterol domain formation
in model membranes," Journal of Biological Chemistry 280 (47), 39380-39387 (2005).

78 B. Jojart and T. A. Martinek, "Performance of the general amber force field in modeling
aqueous POPC membrane bilayers," Journal of Computational Chemistry 28 (12), 2051-
2058 (2007).

79 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein,
"Comparison of Simple Potential Functions for Simulating Liquid Water," Journal of
Chemical Physics 79 (2), 926-935 (1983).

80 R. P. Joshi, Q. Hu, R. Aly, K. H. Schoenbach, and H. P. Hjalmarson, "Self-consistent
simulations of electroporation dynamics in biological cells subjected to ultrashort
electrical pulses," Physical Review E 64:011913 (1) (2001).

81 K. Kinosita and T. Y. Tsong, "Voltage-Induced Pore Formation and Hemolysis of
Human Erythrocytes," Biochimica Et Biophysica Acta 471 (2), 227-242 (1977).

188

82 M. A. Kol, A. N. C. van Laak, D. T. S. Rijkers, J. A. Killian, A. I. P. M. de Kroon, and B.
de Kruijff, "Phospholipid flop induced by transmembrane peptides in model membranes
is modulated by lipid composition," Biochemistry 42 (1), 231-237 (2003).

83 S. Koronkiewicz, S. Kalinowski, and K. Bryl, "Programmable chronopotentiometry as a
tool for the study of electroporation and resealing of pores in bilayer lipid membranes,"
Biochimica Et Biophysica Acta-Biomembranes 1561 (2), 222-229 (2002).

84 T. Kotnik and D. Miklavcic, "Second-order model of membrane electric field induced by
alternating external electric fields," Ieee Transactions on Biomedical Engineering 47 (8),
1074-1081 (2000).

85 P. Kramar, D. Miklavcic, and A. M. Lebar, "Determination of the lipid bilayer
breakdown voltage by means of linear rising signal," Bioelectrochemistry 70 (1), 23-27
(2007).

86 W. Krassowska and P. D. Filev, "Modeling electroporation in a single cell," Biophysical
Journal 92 (2), 404-417 (2007).

87 A. V. Krylov, P. Pohl, M. L. Zeidel, and W. G. Hill, "Water permeability of asymmetric
planar lipid bilayers: Leaflets of different composition offer independent and additive
resistances to permeation," J Gen Physiol 118 (4), 333-339 (2001).

88 N. Kucerka, D. Marquardt, T. A. Harroun, M. P. Nieh, S. R. Wassall, and J. Katsaras,
"The Functional Significance of Lipid Diversity: Orientation of Cholesterol in Bilayers Is
Determined by Lipid Species," Journal of the American Chemical Society 131 (45),
16358-& (2009).

89 G. Lamoureux, A. D. MacKerell, and B. Roux, "A simple polarizable model of water
based on classical Drude oscillators," Journal of Chemical Physics 119 (10), 5185-5197
(2003).

90 J. Lee, A. Ishihara, G. Oxford, B. Johnson, and K. Jacobson, "Regulation of cell
movement is mediated by stretch-activated calcium channels," Nature 400 (6742), 382-
386 (1999).

91 S. Leekumjorn and A. K. Sum, "Molecular characterization of gel and liquid-crystalline
structures of fully hydrated POPC and POPE bilayers," Journal of Physical Chemistry B
111 (21), 6026-6033 (2007).

92 H. Leontiadou, A. E. Mark, and S. J. Marrink, "Molecular dynamics simulations of
hydrophilic pores in lipid bilayers," Biophysical Journal 86 (4), 2156-2164 (2004).

93 Z. A. Levine and P. T. Vernier, "Life Cycle of an Electropore: Field-Dependent and
Field-Independent Steps in Pore Creation and Annihilation," Journal of Membrane
Biology 236 (1), 27-36 (2010).
189

94 Z. A. Levine and P. T. Vernier, "Calcium and phosphatidylserine inhibit lipid electropore
formation and reduce pore lifetime," J Membr Biol 245 (10), 599-610 (2012).

95 T. J. Lewis, "A model for bilayer membrane electroporation based on resultant
electromechanical stress," Ieee Transactions on Dielectrics and Electrical Insulation 10
(5), 769-777 (2003).

96 J. L. MacCallum, W. F. D. Bennett, and D. P. Tieleman, "Distribution of amino acids in a
lipid bilayer from computer simulations," Biophysical Journal 94 (9), 3393-3404 (2008).

97 M. Maccarrone, N. Rosato, and A. F. Agro, "Electroporation Enhances Cell-Membrane
Peroxidation and Luminescence," Biochemical and Biophysical Research
Communications 206 (1), 238-245 (1995).

98 M. W. Mahoney and W. L. Jorgensen, "A five-site model for liquid water and the
reproduction of the density anomaly by rigid, nonpolarizable potential functions," Journal
of Chemical Physics 112 (20), 8910-8922 (2000).

99 S. J. Marrink and H. J. C. Berendsen, "Simulation of Water Transport through a Lipid-
Membrane," Journal of Physical Chemistry 98 (15), 4155-4168 (1994).

100 S. J. Marrink and H. J. C. Berendsen, "Permeation process of small molecules across
lipid membranes studied by molecular dynamics simulations," Journal of Physical
Chemistry 100 (41), 16729-16738 (1996).

101 S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, and A. H. de Vries, "The
MARTINI force field: coarse grained model for biomolecular simulations," J Phys Chem
B 111 (27), 7812-7824 (2007).

102 M. Marty, G. Sersa, J. R. Garbay, J. Gehl, C. G. Collins, M. Snoj, V. Billard, P. F.
Geertsen, J. O. Larkin, D. Miklavcic, I. Pavlovic, S. M. Paulin-Kosir, M. Cemazar, N.
Morsli, Z. Rudolf, C. Robert, G. C. O'Sullivan, and L. M. Mir, "Electrochemotherapy -
An easy, highly effective and safe treatment of cutaneous and subcutaneous metastases:
Results of ESOPE (European Standard Operating Procedures of Electrochemotherapy)
study," Ejc Supplements 4 (11), 3-13 (2006).

103 J. C. Mathai, S. Tristram-Nagle, J. F. Nagle, and M. L. Zeidel, "Structural determinants
of water permeability through the lipid membrane," J Gen Physiol 131 (1), 69-76 (2008).

104 J. Mattai, H. Hauser, R. A. Demel, and G. G. Shipley, "Interactions of Metal-Ions with
Phosphatidylserine Bilayer-Membranes - Effect of Hydrocarbon Chain Unsaturation,"
Biochemistry 28 (5), 2322-2330 (1989).

190

105 J. P. Mattila, K. Sabatini, and P. K. J. Kinnunen, "Oxidized phospholipids as potential
molecular targets for antimicrobial peptides," Biochimica Et Biophysica Acta-
Biomembranes 1778 (10), 2041-2050 (2008).

106 K. C. Melikov, V. A. Frolov, A. Shcherbakov, A. V. Samsonov, Y. A. Chizmadzhev, and
L. V. Chernomordik, "Voltage-induced nonconductive pre-pores and metastable single
pores in unmodified planar lipid bilayer," Biophysical Journal 80 (4), 1829-1836 (2001).

107 D. Miklavcic, K. Beravs, D. Semrov, M. Cemazar, F. Demsar, and G. Sersa, "The
importance of electric field distribution for effective in vivo electroporation of tissues,"
Biophys J 74 (5), 2152-2158 (1998).

108 L. M. Mir, M. F. Bureau, J. Gehl, R. Rangara, D. Rouy, J. M. Caillaud, P. Delaere, D.
Branellec, B. Schwartz, and D. Scherman, "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 (8), 4262-4267 (1999).

109 S. Miyamoto and P. A. Kollman, "Settle - an Analytical Version of the Shake and Rattle
Algorithm for Rigid Water Models," Journal of Computational Chemistry 13 (8), 952-
962 (1992).

110 Ole G. Mouritsen, Life - as a matter of fat : the emerging science of lipidomics. (Springer,
Berlin, 2005).

111 P. Mukhopadhyay, L. Monticelli, and D. P. Tieleman, "Molecular dynamics simulation
of a palmitoyl-oleoyl phosphatidylserine bilayer with Na+ counterions and NaCl,"
Biophysical Journal 86 (3), 1601-1609 (2004).

112 K. J. Muller, V. L. Sukhorukov, and U. Zimmermann, "Reversible
electropermeabilization of mammalian cells by high-intensity, ultra-short pulses of
submicrosecond duration," Journal of Membrane Biology 184 (2), 161-170 (2001).

113 T. Nagai and Y. Okamoto, "Replica-exchange molecular dynamics simulation of a lipid
bilayer system with a coarse-grained model," Mol Simulat 38 (5), 437-441 (2012).

114 Y. M. Naguib, "A fluorometric method for measurement of peroxyl radical scavenging
activities of lipophilic antioxidants," Anal Biochem 265 (2), 290-298 (1998).

115 M. T. Nelson, W. Humphrey, A. Gursoy, A. Dalke, L. V. Kale, R. D. Skeel, and K.
Schulten, "NAMD: A parallel, object oriented molecular dynamics program," Int J
Supercomput Ap 10 (4), 251-268 (1996).

116 O. M. Nesin, O. N. Pakhomova, S. Xiao, and A. G. Pakhomov, "Manipulation of cell
volume and membrane pore comparison following single cell permeabilization with 60-
and 600-ns electric pulses," Biochimica Et Biophysica Acta-Biomembranes 1808 (3),
792-801 (2011).
191

117 J. C. Neu and W. Krassowska, "Asymptotic model of electroporation," Physical Review
E 59 (3), 3471-3482 (1999).

118 E. Neumann and K. Rosenheck, "Permeability Changes Induced by Electric Impulses in
Vesicular Membranes," Journal of Membrane Biology 10 (3-4), 279-290 (1972).

119 E. Neumann, M. Schaeferridder, Y. Wang, and P. H. Hofschneider, "Gene-Transfer into
Mouse Lyoma Cells by Electroporation in High Electric-Fields," EMBO Journal 1 (7),
841-845 (1982).

120 R. Nuccitelli, X. Chen, A. G. Pakhomov, W. H. Baldwin, S. Sheikh, J. L. Pomicter, W.
Ren, C. Osgood, R. J. Swanson, J. F. Kolb, S. J. Beebe, and K. H. Schoenbach, "A new
pulsed electric field therapy for melanoma disrupts the tumor's blood supply and causes
complete remission without recurrence," International Journal of Cancer 125 (2), 438-445
(2009).

121 R. Nuccitelli, U. Pliquett, X. H. Chen, W. Ford, R. J. Swanson, S. J. Beebe, J. F. Kolb,
and K. H. Schoenbach, "Nanosecond pulsed electric fields cause melanomas to self-
destruct," Biochemical and Biophysical Research Communications 343 (2), 351-360
(2006).

122 Y. Okuno, M. Minagawa, H. Matsumoto, and A. Tanioka, "Simulation study on the
influence of an electric field on water evaporation," Journal of Molecular Structure-
Theochem 904 (1-3), 83-90 (2009).

123 A. G. Pakhomov, A. M. Bowman, B. L. Ibey, F. M. Andre, O. N. Pakhomova, and K. H.
Schoenbach, "Lipid nanopores can form a stable, ion channel-like conduction pathway in
cell membrane," Biochemical and Biophysical Research Communications 385 (2), 181-
186 (2009).

124 A. G. Pakhomov, J. F. Kolb, J. A. White, R. P. Joshi, S. Xiao, and K. H. Schoenbach,
"Long-lasting plasma membrane permeabilization in mammalian cells by nanosecond
pulsed electric field (nsPEF)," Bioelectromagnetics 28 (8), 655-663 (2007).

125 Andrei G. Pakhomov, Damijan Miklavčič, and Marko Markov, Advanced electroporation
techniques in biology and medicine. (CRC Press, Boca Raton, FL, 2010).

126 A. G. Pakhomov, R. Shevin, J. A. White, J. F. Kolb, O. N. Pakhomova, R. P. Joshi, and
K. H. Schoenbach, "Membrane permeabilization and cell damage by ultrashort electric
field shocks," Archives of Biochemistry and Biophysics 465 (1), 109-118 (2007).

127 E. H. Pap, G. P. Drummen, V. J. Winter, T. W. Kooij, P. Rijken, K. W. Wirtz, J. A. Op
den Kamp, W. J. Hage, and J. A. Post, "Ratio-fluorescence microscopy of lipid oxidation
in living cells using C11-BODIPY(581/591)," FEBS Lett 453 (3), 278-282 (1999).

192

128 M. Parrinello and A. Rahman, "Polymorphic Transitions in Single-Crystals - a New
Molecular-Dynamics Method," Journal of Applied Physics 52 (12), 7182-7190 (1981).

129 D. A. Pearlman, D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatham, S. Debolt, D.
Ferguson, G. Seibel, and P. Kollman, "Amber, a Package of Computer-Programs for
Applying Molecular Mechanics, Normal-Mode Analysis, Molecular-Dynamics and Free-
Energy Calculations to Simulate the Structural and Energetic Properties of Molecules,"
Comput Phys Commun 91 (1-3), 1-41 (1995).

130 L. Perera, U. Essmann, and M. L. Berkowitz, "Role of water in the hydration force acting
between lipid bilayers," Langmuir 12 (11), 2625-2629 (1996).

131 Rob Phillips, Jane Kondev, and Julie Theriot, Physical biology of the cell. (Garland
Science, New York, 2009).

132 C. Poignard, A. Silve, F. Campion, L. M. Mir, O. Saut, and L. Schwartz, "Ion fluxes,
transmembrane potential, and osmotic stabilization: a new dynamic electrophysiological
model for eukaryotic cells," European Biophysics Journal with Biophysics Letters 40 (3),
235-246 (2011).

133 D. Popescu, C. Rucareanu, and G. Victor, "A Model for the Appearance of Statistical
Pores in Membranes Due to Selfoscillations," Bioelectrochemistry and Bioenergetics 25
(1), 91-103 (1991).

134 R. D. Porasso and J. J. L. Cascales, "Study of the effect of Na+ and Ca2+ ion
concentration on the structure of an asymmetric DPPC/DPPC plus DPPS lipid bilayer by
molecular dynamics simulation," Colloid Surface B 73 (1), 42-50 (2009).

135 G. M. Preston, J. S. Jung, W. B. Guggino, and P. Agre, "The Mercury-Sensitive Residue
at Cysteine-189 in the Chip28 Water Channel," Journal of Biological Chemistry 268 (1),
17-20 (1993).

136 G. Pucihar, T. Kotnik, D. Miklavcic, and J. Teissie, "Kinetics of transmembrane transport
of small molecules into electropermeabilized cells," Biophysical Journal 95 (6), 2837-
2848 (2008).

137 J. Racine, "gnuplot 4.0: A portable interactive plotting utility," J Appl Econom 21 (1),
133-141 (2006).

138 M. P. Rols, "Electropermeabilization, a physical method for the delivery of therapeutic
molecules into cells," Biochimica Et Biophysica Acta-Biomembranes 1758 (3), 423-428
(2006).

139 M. P. Rols, D. Coulet, and J. Teissie, "Highly Efficient Transfection of Mammalian-Cells
by Electric-Field Pulses - Application to Large Volumes of Cell-Culture by Using a Flow
System," European Journal of Biochemistry 206 (1), 115-121 (1992).
193

140 M. P. Rols, C. Delteil, M. Golzio, P. Dumond, S. Cros, and J. Teissie, "In vivo
electrically mediated protein and gene transfer in murine melanoma," Nature
Biotechnology 16 (2), 168-171 (1998).

141 M. P. Rols and J. Teissie, "Electropermeabilization of Mammalian-Cells - Quantitative-
Analysis of the Phenomenon," Biophysical Journal 58 (5), 1089-1098 (1990).

142 M. P. Rols and J. Teissie, "Experimental-Evidence for the Involvement of the
Cytoskeleton in Mammalian-Cell Electropermeabilization," Biochimica Et Biophysica
Acta 1111 (1), 45-50 (1992).

143 S. Romeo, Y. Wu, Z.A. Levine, M. A. Gundersen, P. T. Vernier, "Water influx and cell
swelling after nanosecond electropermeabilization," Biochimica et Biophysica Acta
(BBA) - Biomembranes 1828 (8), 1715-1722 (2013)
144 B. Rubinsky, "Irreversible electroporation in medicine," Technol Cancer Res Treat 6 (4),
255-260 (2007).

145 J. P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, "Numerical-Integration of Cartesian
Equations of Motion of a System with Constraints - Molecular-Dynamics of N-Alkanes,"
J Comput Phys 23 (3), 327-341 (1977).

146 R. Ryham, I. Berezovik, and F. S. Cohent, "Aqueous Viscosity Is the Primary Source of
Friction in Lipidic Pore Dynamics," Biophysical Journal 101 (12), 2929-2938 (2011).

147 A. J. H. Sale and W. A. Hamilton, "Effects of High Electric Fields on Microorganisms .I.
Killing of Bacteria and Yeasts," Biochimica Et Biophysica Acta 148 (3), 781-& (1967).

148 A. J. H. Sale and W. A. Hamilton, "Effects of High Electric Fields on Micro-Organisms
.3. Lysis of Erythrocytes and Protoplasts," Biochimica Et Biophysica Acta 163 (1), 37-&
(1968).

149 J. M. Sanders, A. Kuthi, Y. H. Wu, P. T. Vernier, and M. A. Gundersen, "A Linear,
Single-stage, Nanosecond Pulse Generator for Delivering Intense Electric Fields to
Biological Loads," Ieee Transactions on Dielectrics and Electrical Insulation 16 (4),
1048-1054 (2009).

150 O. Sandre, L. Moreaux, and F. Brochard-Wyart, "Dynamics of transient pores in
stretched vesicles," Proceedings of the National Academy of Sciences of the United
States of America 96 (19), 10591-10596 (1999).

151 N. Sapay, W. F. D. Bennett, and D. P. Tieleman, "Thermodynamics of flip-flop and
desorption for a systematic series of phosphatidylcholine lipids," Soft Matter 5 (17),
3295-3302 (2009).

194

152 G. Saulis, S. Satkauskas, and R. Praneviciute, "Determination of cell electroporation from
the release of intracellular potassium ions," Analytical Biochemistry 360 (2), 273-281
(2007).

153 P. Schoch, D. F. Sargent, and R. Schwyzer, "Capacitance and Conductance as Tools for
the Measurement of Asymmetric Surface-Potentials and Energy Barriers of Lipid Bilayer
Membranes," Journal of Membrane Biology 46 (1), 71-89 (1979).

154 K. H. Schoenbach, R. P. Joshi, J. F. Kolb, N. Y. Chen, M. Stacey, P. F. Blackmore, E. S.
Buescher, and S. J. Beebe, "Ultrashort electrical pulses open a new gateway into
biological cells," Proceedings of the Ieee 92 (7), 1122-1137 (2004).

155 H. P. Schwan and K. R. Foster, "Microwave Dielectric Properties of Tissue - Some
Comments on Rotational Mobility of Tissue Water," Biophysical Journal 17 (2), 193-197
(1977).

156 K. Schwister and B. Deuticke, "Formation and Properties of Aqueous Leaks Induced in
Human-Erythrocytes by Electrical Breakdown," Biochimica Et Biophysica Acta 816 (2),
332-348 (1985).

157 A. Seelig and J. Seelig, "Dynamic Structure of Fatty Acyl Chains in a Phospholipid
Bilayer Measured by Deuterium Magnetic-Resonance," Biochemistry 13 (23), 4839-4845
(1974).

158 J. Seelig, "Spin Label Studies of Oriented Smectic Liquid Crystals - a Model System for
Bilayer Membranes," Journal of the American Chemical Society 92 (13), 3881-& (1970).

159 D. Sengupta, H. Leontiadou, A. E. Mark, and S. J. Marrink, "Toroidal pores formed by
antimicrobial peptides show significant disorder," Biochimica Et Biophysica Acta-
Biomembranes 1778 (10), 2308-2317 (2008).

160 C. G. Sinn, M. Antonietti, and R. Dimova, "Binding of calcium to phosphatidylcholine-
phosphatidylserine membranes," Colloids and Surfaces A-Physicochemical and
Engineering Aspects 282, 410-419 (2006).

161 K. C. Smith and J. C. Weaver, "Active mechanisms are needed to describe cell responses
to submicrosecond, megavolt-per-meter pulses: Cell models for ultrashort pulses,"
Biophysical Journal 95 (4), 1547-1563 (2008).

162 R. Stampfli and M. Willi, "Membrane Potential of a Ranvier Node Measured after
Electrical Destruction of Its Membrane," Experientia 13 (7), 297-298 (1957).

163 G. Stark, "The Effect of Ionizing-Radiation on Lipid-Membranes," Biochimica Et
Biophysica Acta 1071 (2), 103-122 (1991).

195

164 O. Stern, "The theory of the electrolytic double shift," Z Elktrochem Angew P 30, 508-
516 (1924).

165 D. A. Stewart, I. R. Gowrishankar, and J. C. Weaver, "Transport lattice approach to
describing cell electroporation: Use of a local asymptotic model," Ieee Transactions on
Plasma Science 32 (4), 1696-1708 (2004).

166 I. P. Sugar, "The effects of external fields on the structure of lipid bilayers," J Physiol
(Paris) 77 (9), 1035-1042 (1981).

167 I. P. Sugar and E. Neumann, "Stochastic-Model for Electric Field-Induced Membrane
Pores. Electroporation," Biophysical Chemistry 19 (3), 211-225 (1984).

168 V. L. Sukhorukov, R. Reuss, D. Zimmermann, C. Held, K. J. Muller, M. Kiesel, P.
Gessner, A. Steinbach, W. A. Schenk, E. Bamberg, and U. Zimmermann, "Surviving
high-intensity field pulses: strategies for improving robustness and performance of
electrotransfection and electrofusion," J Membr Biol 206 (3), 187-201 (2005).

169 S. Sun, G. Y. Yin, Y. K. Lee, J. T. Y. Wong, and T. Y. Zhang, "Effects of deformability
and thermal motion of lipid membrane on electroporation: By molecular dynamics
simulations," Biochemical and Biophysical Research Communications 404 (2), 684-688
(2011).

170 T. Tang, F. Wang, A. Kuthi, and M. A. Gundersen, "Diode opening switch based
nanosecond high voltage pulse generators for biological and medical applications," Ieee
Transactions on Dielectrics and Electrical Insulation 14 (4), 878-883 (2007).

171 M. Tarek, "Membrane electroporation: A molecular dynamics simulation," Biophysical
Journal 88 (6), 4045-4053 (2005).

172 B. Tavazzi, D. Di Pierro, A. M. Amorini, G. Fazzina, M. Tuttobene, B. Giardina, and G.
Lazzarino, "Energy metabolism and lipid peroxidation of human erythrocytes as a
function of increased oxidative stress," European Journal of Biochemistry 267 (3), 684-
689 (2000).

173 J. Teissie, M. Golzio, and M. P. Rols, "Mechanisms of cell membrane
electropermeabilization: A minireview of our present (lack of ?) knowledge," Biochimica
Et Biophysica Acta-General Subjects 1724 (3), 270-280 (2005).

174 J. Teissie and T. Y. Tsong, "Electric-Field Induced Transient Pores in Phospholipid-
Bilayer Vesicles," Biochemistry 20 (6), 1548-1554 (1981).

175 D. P. Tieleman, "The molecular basis of electroporation," Biophysical Journal 86 (1),
371a-372a (2004).

196

176 D. P. Tieleman, H. Leontiadou, A. E. Mark, and S. J. Marrink, "Simulation of pore
formation in lipid bilayers by mechanical stress and electric fields," Journal of the
American Chemical Society 125 (21), 6382-6383 (2003).

177 D. P. Tieleman and S. J. Marrink, "Lipids out of equilibrium: energetics of desorption
and pore mediated flip-flop," J Am Chem Soc 128 (38), 12462-12467 (2006).

178 D. P. Tieleman, S. J. Marrink, and H. J. C. Berendsen, "A computer perspective of
membranes: molecular dynamics studies of lipid bilayer systems," Biochimica Et
Biophysica Acta-Reviews on Biomembranes 1331 (3), 235-270 (1997).

179 M. Tokman, J. H. Lee, Z. A. Levine, M. C. Ho, M. E. Colvin, and P. T. Vernier, "Electric
Field-Driven Water Dipoles: Nanoscale Architecture of Electroporation," PLoS One 8 (4)
(2013).

180 R. Vacha, M. L. Berkowitz, and P. Jungwirth, "Molecular Model of a Cell Plasma
Membrane With an Asymmetric Multicomponent Composition: Water Permeation and
Ion Effects," Biophysical Journal 96 (11), 4493-4501 (2009).

181 P. van der Ploeg and H. J. C. Berendsen, "Molecular-Dynamics Simulation of a Bilayer-
Membrane," Journal of Chemical Physics 76 (6), 3271-3276 (1982).

182 D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. Berendsen,
"GROMACS: fast, flexible, and free," J Comput Chem 26 (16), 1701-1718 (2005).

183 Jack Vanderlinde, Classical electromagnetic theory. (Wiley, New York, 1993).

184 Z. Vasilkoski, A. T. Esser, T. R. Gowrishankar, and J. C. Weaver, "Membrane
electroporation: The absolute rate equation and nanosecond time scale pore creation,"
Physical Review E 74 (2) (2006).

185 A. E. Vaughn and M. Deshmukh, "Glucose metabolism inhibits apoptosis in neurons and
cancer cells by redox inactivation of cytochrome c," Nature Cell Biology 10 (12), 1477-
U1228 (2008).

186 R. M. Venable, Y. H. Zhang, B. J. Hardy, and R. W. Pastor, "Molecular-Dynamics
Simulations of a Lipid Bilayer and of Hexadecane - an Investigation of Membrane
Fluidity," Science 262 (5131), 223-226 (1993).

187 P. T. Vernier, Z. A. Levine, and M. A. Gundersen, "Water Bridges in
Electropermeabilized Phospholipid Bilayers," Proceedings of the IEEE 101 (2), 494-504
(2013).

188 P. T. Vernier, Z. A. Levine, Y. H. Wu, V. Joubert, M. J. Ziegler, L. M. Mir, and D. P.
Tieleman, "Electroporating Fields Target Oxidatively Damaged Areas in the Cell
Membrane," PLoS One 4 (11) (2009).
197

189 P. T. Vernier, A. M. Li, L. Marcu, C. M. Craft, and M. A. Gundersen, "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 (5), 795-809 (2003).

190 P. T. Vernier, Y. H. Sun, M. T. Chen, M. A. Gundersen, and G. L. Craviso, "Nanosecond
electric pulse-induced calcium entry into chromaffin cells," Bioelectrochemistry 73 (1),
1-4 (2008).

191 P. T. Vernier, Y. H. Sun, and M. A. Gundersen, "Nanoelectropulse-driven membrane
perturbation and small molecule permeabilization," BMC Cell Biology 7, 37 (2006).

192 P. T. Vernier, Y. H. Sun, L. Marcu, S. Salemi, C. M. Craft, and M. A. Gundersen,
"Calcium bursts induced by nanosecond electric pulses," Biochemical and Biophysical
Research Communications 310 (2), 286-295 (2003).

193 P. T. Vernier and M. J. Ziegler, "Nanosecond field alignment of head group and water
dipoles in electroporating phospholipid bilayers," Journal of Physical Chemistry B 111
(45), 12993-12996 (2007).

194 P. T. Vernier, M. J. Ziegler, and R. Dimova, "Calcium Binding and Head Group Dipole
Angle in Phosphatidylserine-Phosphatidylcholine Bilayers," Langmuir 25 (2), 1020-1027
(2009).

195 P. T. Vernier, M. J. Ziegler, Y. Sun, M. A. Gundersen, and D. P. Tieleman, "Nanopore-
facilitated, voltage-driven phosphatidylserine translocation in lipid bilayers--in cells and
in silico," Phys Biol 3 (4), 233-247 (2006).

196 P. T. Vernier, M. J. Ziegler, Y. H. Sun, W. V. Chang, M. A. Gundersen, and D. P.
Tieleman, "Nanopore formation and phosphatidylserine externalization in a phospholipid
bilayer at high transmembrane potential," Journal of the American Chemical Society 128
(19), 6288-6289 (2006).

197 S. F. Wang, J. X. Chen, M. T. Chen, P. T. Vernier, M. A. Gundersen, and M.
Valderrabano, "Cardiac Myocyte Excitation by Ultrashort High-Field Pulses,"
Biophysical Journal 96 (4), 1640-1648 (2009).

198 A. Warshel, P. K. Sharma, M. Kato, and W. W. Parson, "Modeling electrostatic effects in
proteins," Bba-Proteins Proteom 1764 (11), 1647-1676 (2006).

199 J. C. Weaver, "Electroporation of biological membranes from multicellular to nano
scales," IEEE Transactions on Dielectrics and Electrical Insulation 10 (5), 754-768
(2003).

198

200 J. C. Weaver and Y. A. Chizmadzhev, "Theory of electroporation: A review,"
Bioelectrochemistry and Bioenergetics 41 (2), 135-160 (1996).

201 J. C. Weaver and R. A. Mintzer, "Decreased Bilayer Stability Due to Transmembrane
Potentials," Physics Letters A 86 (1), 57-59 (1981).

202 J. A. White, P. F. Blackmore, K. H. Schoenbach, and S. J. Beebe, "Stimulation of
capacitative calcium entry in HL-60 cells by nanosecond pulsed electric fields," Journal
of Biological Chemistry 279 (22), 22964-22972 (2004).

203 M. Winterhalter, "Black lipid membranes," Current Opinion in Colloid & Interface
Science 5 (3-4), 250-255 (2000).

204 T. K. Wong and E. Neumann, "Electric field mediated gene transfer," Biochem Biophys
Res Commun 107 (2), 584-587 (1982).

205 J. Wong-Ekkabut, Z. T. Xu, W. Triampo, I. M. Tang, D. P. Tieleman, and L. Monticelli,
"Effect of lipid peroxidation on the properties of lipid bilayers: A molecular dynamics
study," Biophysical Journal 93 (12), 4225-4236 (2007).

206 X. C. Yang and F. Sachs, "Block of Stretch-Activated Ion Channels in Xenopus Oocytes
by Gadolinium and Calcium-Ions," Science 243 (4894), 1068-1071 (1989).

207 N. Zamzami, P. Marchetti, M. Castedo, D. Decaudin, A. Macho, T. Hirsch, S. A. Susin,
P. X. Petit, B. Mignotte, and G. Kroemer, "Sequential Reduction of Mitochondrial
Transmembrane Potential and Generation of Reactive Oxygen Species in Early
Programmed Cell-Death," Journal of Experimental Medicine 182 (2), 367-377 (1995).

208 Y. Zhou, C. K. Berry, P. A. Storer, and R. M. Raphael, "Peroxidation of polyunsaturated
phosphatidyl-choline lipids during electroformation," Biomaterials 28 (6), 1298-1306
(2007).

209 M. J. Ziegler and P. T. Vernier, "Interface Water Dynamics and Porating Electric Fields
for Phospholipid Bilayers," Journal of Physical Chemistry B 112 (43), 13588-13596
(2008).

210 U. Zimmermann, "Electric Field-Mediated Fusion and Related Electrical Phenomena,"
Biochimica Et Biophysica Acta 694 (3), 227-277 (1982).

211 U. Zimmermann, G. Pilwat, and F. Riemann, "Dielectric-Breakdown of Cell-
Membranes," Biophysical Journal 14 (11), 881-899 (1974).

199

Appendix A

Analytical Derivations

A.1 A derivation of the Langevin-Debye dipole moment in an electric field

Here we will derive the dipole moment for an ideal electric dipole in the presence of a
+ 𝑧 ̂-directed electric field as a function of field strength. We will make assumptions which are
similar to the case of the three-dimensional paramagnet in an external magnetic field, a typical
problem solved in thermal physics. First we define the interaction energy U of a simple dipole in
an electric field:
𝑈 = − 𝜇 ⃑ ∙ 𝐸 � ⃑
= −| 𝜇 || 𝐸 | 𝑐𝑜𝑠 𝜃
where µ is the maximum dipole moment, E is the magnitude of the electric field, and θ is the
angle between the electric field and the dipole moment. If there exists a large number of electric
dipoles in the system then the number of dipoles with interaction energy U follows a Boltzmann
distribution:
𝑁 ( 𝑈 ) = 𝐴 𝑒 −
𝑈 𝑘𝑇

where k is Boltzmann’s constant, T is temperature, and A is some normalization term. Thus the
average dipole moment in this system can be written as:
200

< 𝜇 > =
∬ 𝜇 𝑁 ( 𝑈 ) 𝑑Ω
∬ 𝑁 ( 𝑈 ) 𝑑Ω

where dΩ = sinθdθdϕ is a differential piece of solid angle. Since we already know that the
electric field is 𝑧 ̂-oriented (E=E
z
) we need only extract µ
z
= µcosθ, thus:
< 𝜇 𝑧 > =
∫ ∫ 𝜇 𝑐𝑜𝑠 𝜃 𝐴 𝑒 +
| 𝜇 || 𝐸 | 𝑐𝑜 𝑠𝜃
𝑘𝑇
𝑠𝑖 𝑛 𝜃 𝑑𝜃 𝑑𝜑 𝜋 0
2 𝜋 0

∫ ∫ 𝐴 𝑒 +
| 𝜇 || 𝐸 | 𝑐𝑜 𝑠𝜃
𝑘𝑇
𝑠𝑖 𝑛 𝜃 𝑑𝜃 𝑑𝜑 𝜋 0
2 𝜋 0

To simplify the above integral we can use the following identity:
∫ d φ
2 𝜋 0
∫ f(cos θ) sin θd θ = 2 𝜋 ∫ 𝑓 ( 𝑐𝑜𝑠 𝜃 ) 𝑑 ( 𝑐𝑜𝑠 𝜃 ) = 2 𝜋 ∫ 𝑓 ( 𝑥 ) 𝑑𝑥 1
− 1

1
− 1
𝜋 0

to arrive at:
< 𝜇 𝑧 > =
2 𝜋𝐴𝜇 ∫ 𝑥 𝑒 𝑐𝑥
𝑑𝑥 1
− 1

2 𝜋𝐴 ∫ 𝑒 𝑐𝑥
𝑑𝑥 1
− 1

where c =
𝜇 𝐸 𝑘𝑇
and x=cosθ. Integrating the numerator by parts gives us:
< 𝜇 𝑧 > = 𝜇 [
𝑥 𝑐 𝑒 𝑐𝑥
|
1
− 1
−
1
𝑐 ∫ 𝑒 𝑐𝑥
𝑑𝑥 1
− 1
]
1
𝑐 𝑒 𝑐𝑥
|
1
− 1

∴ < 𝜇 𝑧 > = 𝜇 �
1
tanh[ 𝑐 ]
−
1
𝑐 �
𝑜𝑟 < 𝜇 𝑧 > = 𝜇 � coth ( 𝑐 ) −
1
𝑐 �
< 𝜇 𝑧 > = 𝜇 � coth (
𝜇𝐸
𝑘𝑇
) −
𝑘𝑇
𝜇𝐸
�
201

< 𝜇 𝑧 > = 𝜇𝐿 (
𝜇𝐸
𝑘𝑇
)
where L(x) is the Langevin function (coth ( 𝑥 ) −
1
𝑥 ) of variable x. Thus the dipole moment in an
external electric field depends only on the relative strength of the electric interaction energy µE
relative to the ambient thermal energy kT.

202

A.2 A derivation of the diffusional flow of water across an intact membrane

Here we will model the diffusional flux of water through a semipermeable barrier where
particle concentrations on each side of the barrier are kept uniform. The primary assumption we
will require throughout this calculation is that we expect a linear change in particle concentration
inside of pores which mediate water flow from one side of the barrier to the other. Starting with
standard Fick diffusion, we have for a divergentless or sourceless flow:
𝜕𝐶
𝜕𝑡
= − ∇ 𝑗 , 𝑗 = − 𝐷 ∇C, thus
𝜕𝐶
𝜕𝑡
= 𝐷 ∇
2
𝐶
where C is particle concentration, t is time, ∇ is the gradient operator, j is particle flow, and D is
diffusion coefficient. Imagining the barrier (or membrane) as a rectangular slab with discrete
cylindrical tubes or pores distributed across its surface, we are free to call the concentration on
one side of the membrane C
1
and the concentration on the other side of the membrane C
2
.
Invoking our linearity assumption about the concentration gradient across the membrane, we
have:
𝜕𝐶
𝜕𝑥
= −
C
1
− C
2
∆ 𝑥 = −
∆C
∆ 𝑥
which when combined with Fick’s equation for a single pore becomes:
𝑗 ( 𝑥 , 𝑡 )
𝑝 𝑜𝑟 𝑒 = − 𝐷 ∂ 𝐶 𝜕𝑥
= 𝐷 ∆C
∆ 𝑥
a description of the flow rate per unit area through a single electropore. However if we consider J
- the total number of particles which traverses through all cylindrical electropores per second (i.e.
J = j(x,t)*𝜋 𝑎 2
) – then we find that:
203

J = 𝑛 ∗ ( 𝜋 𝑎 2
)
𝐷 Δ 𝑥 Δ 𝐶 = 𝓅 Δ 𝐶
where 𝓅 = 𝑛𝜋 𝑎 2
𝐷 Δ 𝑥 is called the membrane ‘permeability’ which depends on the number n and
size of pores which are present, and on the diffusion of particles which pass through these pores.
Now suppose that the membrane partitions two separate compartments of constant
volume V
1
and V
2
which have particle concentrations C
1
and C
2
and total number of particles N
1

and N
2
respectively. If the concentration C
1
is higher than C
2
(i.e. C
1
> C
2
) then:
dN
1
𝑑𝑡 = −
dN
2
𝑑𝑡 = −J ∙ A = −A 𝓅 ( 𝐶 1
− 𝐶 2
)
where 𝑁 1
( 𝑡 ) = 𝐶 1
( 𝑡 ) 𝑉 1
and 𝑁 2
( 𝑡 ) = 𝐶 2
( 𝑡 ) 𝑉 2
. Thus:
dN
1
𝑑𝑡 = V
1
dC
1
𝑑𝑡 = −A 𝓅 ( 𝐶 1
− 𝐶 2
),
dN
2
𝑑𝑡 = V
2
dC
2
𝑑𝑡 = +A 𝓅 ( 𝐶 1
− 𝐶 2
)
Since there is particle conservation it must be required that N
1
+ N
2
= constant ∴
d N
1
𝑑𝑡
+
d N
2
𝑑𝑡
= 0.
This allows us to add the above expressions to zero:
V
1
dC
1
𝑑𝑡 + V
2
dC
2
𝑑𝑡 =
d
𝑑𝑡 (V
1
C
1
+ V
2
C
2
) = 0
It is obvious that V
1
C
1
+ V
2
C
2
is constant in time, so let us call the initial concentrations C
1
(0)
and C
2
(0), and the final uniform concentration C(∞). Then V
1
C
1
(0) + V
2
C
2
(0) = V
1
C
∞
+
V
2
C
∞
, or:
𝐶 ( ∞) =
𝑉 1
𝑉 1
+ 𝑉 2
𝐶 1
(0) +
𝑉 2
𝑉 1
+ 𝑉 2
𝐶 2
(0)
204

which is our first boundary condition. Similarly if we take the difference of the derivatives above
we find that:
dC
1
𝑑𝑡 −
dC
2
𝑑𝑡 =
d(C
1
− C
2
)
𝑑𝑡 = −A 𝓅 𝑉 1
+ 𝑉 2
𝑉 1
𝑉 2
( 𝐶 1
− 𝐶 2
) = −
( 𝐶 1
− 𝐶 2
)
τ
0

which determines the characteristic time constant τ
o
for the equilibrium of concentrations:
τ
0
=
1
A 𝓅 𝑉 1
𝑉 2
𝑉 1
+ 𝑉 2

where:
∆C(t) = ∆C(0) 𝑒 −
𝑡 𝜏 0

(since the concentration derivative returns the concentration times an inverse time constant).
Using our expressions for ∆C = C
1
− C
2
and
d ∆ C
𝑑𝑡
we obtain:
𝐶 1
( 𝑡 ) = 𝐶 ∞
+ [ 𝐶 1
(0) − 𝐶 ∞
] 𝑒 −
𝑡 𝜏 0

𝐶 2
( 𝑡 ) = 𝐶 ∞
+ [ 𝐶 2
(0) − 𝐶 ∞
] 𝑒 −
𝑡 𝜏 0

and because C
1
(0) > C
∞
> C
2
(0), this results in a converging concentration to C
∞
from both initial
concentrations over a sufficient amount of time.

205

In our observations with osmotic cell swelling after electroporation we see a volume
increase that is similar to the theoretical concentration increase of C
2
(0). However in the case of
our cells which have nearly infinite supplies of ‘inner’ and ‘outer’ water molecules, the cell
volume is free to swell linearly until there is no water left in the buffer to continue filling it. To
account for this we can replace C
∞
with a linear function of the form C
∞
= mx+b (a sloping line)
which gives us a consistent model to compare with our experimental observations (Figure 5.16).

Abstract (if available)

Abstract Computer simulations of physical, chemical, and biological systems have improved tremendously over the past five decades. From simple studies of liquid argon in the 1960s to fully atomistic simulations of entire viruses in the past few years, recent advances in high-performance computing have continuously enabled simulations to bridge the gap between scientific theory and experiment. Molecular dynamics simulations in particular have allowed for the direct observation of spatial and temporal events which are at present inaccessible to experiments. For this dissertation I employ all-atom molecular dynamics simulations to study the transient, electric field-induced poration (or electroporation) of phospholipid bilayers at MV/m electric fields. ❧ Phospholipid bilayers are the dominant constituents of cell membranes and act as both a barrier and gatekeeper to the cell interior. This makes their structural integrity and susceptibility to external perturbations an important topic for study, especially as the density of electromagnetic radiation in our environment is increasing steadily. The primary goal of this dissertation is to understand the specific physical and biological mechanisms which facilitate electroporation, and to connect our simulated observations to experiments with live cells and to continuum models which seek to describe the underlying biological processes of electroporation. ❧ In Chapter 1 I begin with a brief introduction to phospholipids and phospholipid bilayers, followed by an extensive overview of electroporation and atomistic molecular dynamics simulations. The following chapters will then focus on peer-reviewed and published work we performed, or on existing projects which are currently being prepared for submission. Chapter 2 looks at how external electric fields affect both oxidized and unoxidized lipid bilayers as a function of oxidation concentration and oxidized lipid type. Oxidative damage to cell membranes represents a physiologically relevant system where lipids can become damaged or severely impacted from interacting with reactive oxygen species, and these events become more frequent with age. The results are then compared to experiments where we show agreement between our simulations, theoretical models, and experiments with peroxidized cells in our lab. In Chapter 3 I outline a set of unique metrics which can be used to quantitatively measure the life cycle of a discrete electropore for the first time, across multiple lipid species, and I compare these results to analytical models where we find good agreement with theory. In Chapter 4 I use the life cycle of an electropore as a tool to measure the effects of electrolyte and lipid headgroup charge on electroporation compared to electrolyte-free and zwitterionic systems, in addition to presenting ion binding isotherms to determine the validity of our simulated electrolyte models. Chapters 5 and 6 focus on the roles of water and lipid respectively on electroporation using simplified water:vacuum systems, osmotic swelling simulations, systems at varying temperature, and systems where we successfully modulated the electropore radius using customized time-dependent electric fields. I conclude this dissertation with a brief summary of these studies followed by a short outlook on the future of electroporation simulations as a whole.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Molecular simulations of water and monovalent ion dynamics in the electroporation of phospholipid bilayers

PDF

Molecular dynamics simulations of lipid bilayers in megavolt per meter electric fields

PDF

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

PDF

The effects of oxidation on bilayer membranes studied using giant unilamellar vesicles

PDF

Biophysical studies of passive transport across synthetic lipid bilayers

PDF

The physics of membrane protein polyhedra

PDF

Hybrid lipid-based nanostructures

PDF

Probing the effects of transmembrane domains on the continuum mechanics of lipid bilayers

PDF

Physical principles of membrane mechanics, membrane domain formation, and cellular signal transduction

PDF

Multi-scale quantum dynamics and machine learning for next generation energy applications

PDF

Out-of-equilibrium dynamics of inhomogeneous quantum systems

PDF

Shock-induced nanobubble collapse and its applications

PDF

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

PDF

Photoexcitation and nonradiative relaxation in molecular systems: methodology, optical properties, and dynamics

PDF

Molecular dynamics study of energetic materials

PDF

Imaging molecular transport across and nanomaterial interaction with lipid membranes

PDF

Neural network for molecular dynamics simulation and design of 2D materials

PDF

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

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

PDF

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

PDF

The effects of divalent counterions on the formation and stabilization of RNA tertiary structure

Asset Metadata

Creator Levine, Zachary Alan (author)

Core Title Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 06/26/2013

Defense Date 06/21/2013

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

Tag biophysics,DNA transfection,electrochemotherapy,electrogenetherapy,electropermeabilization,electroperturbation,electroporation,lipid,lipid bilayer,membrane,membrane biology,membrane biophysics,molecular dynamics,OAI-PMH Harvest,phospholipid,Physics,pulsed power,simulation

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor El-Naggar, Mohamed Y. (committee chair), Haas, Stephan W. (committee member), Malmstadt, Noah (committee member), Nakano, Aiichiro (committee member), Vernier, Paul Thomas (committee member)

Creator Email zachlevine@gmail.com

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

Unique identifier UC11294989

Identifier etd-LevineZach-1707.pdf (filename),usctheses-c3-280283 (legacy record id)

Legacy Identifier etd-LevineZach-1707.pdf

Dmrecord 280283

Document Type Dissertation

Format application/pdf (imt)

Rights Levine, Zachary Alan

Type texts

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

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

Repository Name University of Southern California Digital Library

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