Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Molecular simulations of water and monovalent ion dynamics in the electroporation of phospholipid bilayers
(USC Thesis Other)
Molecular simulations of water and monovalent ion dynamics in the electroporation of phospholipid bilayers
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
MOLECULAR SIMULATIONS OF WATER AND MONOVALENT ION DYNAMICS IN
THE ELECTROPORATION OF PHOSPHOLIPID BILAYERS
by
Ming-Chak Ho
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(PHYSICS)
May 2014
Copyright 2014 Ming-Chak Ho
ii
Dedication
To my family
iii
Acknowledgements
I would like to express my deepest gratitude to those who helped me during the years in
my graduate study.
My researcher advisor Dr. Thomas Vernier provided me an interesting research topic, and
guided me through the world of scientific research with his wisdom and encouragement. His
criticisms improved the quality of my work and shaped me into a better researcher.
My colleague and friend Zach Levine introduced me into the group and the research. His
support in many occasions helped me overcome the obstacles in my study.
Professor Aiichiro Nakano taught me a lot about Molecular Dynamics and scientific
computing. The pursuit of the master’s degree in High Performance Computing and Simulations
under his encouragement was a valuable experience.
My research would not be possible without the USC Center for High Performance
Computing and Communications (HPCC). I am thankful for the staff who work tirelessly to keep
the machines running efficiently, and provide help when I encounter issues.
I also thank Professor Martin Gundersen who allowed me to work in his group and
oversaw my research progress; Professor Stephan Haas who recruited me into the USC Physics
department and provided many years of encouragement and support; Professor Werner Däppen
and Professor Noah Malmstadt for their kindness in participating in my dissertation committee;
and Professor John Shumway from Arizona State University who provided the earliest research
experience in my undergraduate years.
iv
Last but not least, I want to thank my fellow classmates and graduate students: Yung-
Ching Liang, Hsiao-Hsuan Lin, Kien Trinh, Wen Zhang, Yaqi Tao, Yichen Chan, Yung-Hsu
Lin, Yu-Hsuan Sharon Wu, Kok-Wee Song, and Suet-Ying Daisy Mak. You guys endured many
of my frequent visits and engaged in many useful (and useless) discussions on research (and
other totally random topics). We have struggled and laughed together. No matter where you are
and what you decided to do, I hope you will have a successful post-doctorial career.
v
Table of Contents
Dedication ii
Acknowledgements iii
List of Tables viii
List of Figures ix
Abstract xii
Chapter 1: Introduction 1
1.1 Applications of electroporation ........................................................................................ 1
1.2 Electroporation experiments ............................................................................................ 4
1.3 Continuum models of electroporation .............................................................................. 8
Chapter 2: Molecular Dynamics Simulation 14
2.1 Overview ...................................................................................................................... 14
2.2 MD formalism and implementation .............................................................................. 17
2.2.1 MD discretization and integrators........................................................................ 17
2.2.2 Energy function in biological systems ................................................................. 21
2.3 Molecular structures and models ................................................................................... 23
2.3.1 Water models ...................................................................................................... 24
2.3.2 Lipid models ....................................................................................................... 26
2.4 Holonomic constraints .................................................................................................. 28
2.5 Thermodynamics ensembles and coupling algorithms ................................................... 30
2.5.1 Temperature coupling ......................................................................................... 30
2.5.2 Pressure coupling ................................................................................................ 32
2.6 Long-range interactions ................................................................................................ 34
2.7 Previous MD studies on electroporation ........................................................................ 37
2.7.1 Dynamics and life cycle of an electropore ........................................................... 38
2.7.2 Water’s role in electroporation ............................................................................ 41
2.7.3 MD simulation of field-stabilized electropore ..................................................... 43
Chapter 3: Dynamics of water bridge in electroporation 46
3.1 Abstract ........................................................................................................................ 46
3.2 Introduction .................................................................................................................. 48
3.3 Methods........................................................................................................................ 50
3.3.1 Molecular dynamics simulations ......................................................................... 50
3.3.2 Structure ............................................................................................................. 51
vi
3.3.3 Water bridge initiation time, annihilation time .................................................... 51
3.3.4 Measurement of water bridge radius .................................................................... 52
3.4 Results and Discussion ................................................................................................. 53
3.4.1 Internal electric field versus external electric field ............................................... 53
3.4.2 Water bridge formation ....................................................................................... 54
3.4.3 Water bridge stabilization ................................................................................... 58
3.4.4 Water bridge annihilation .................................................................................... 60
3.4.5 Effects of water models in water bridge formation simulation ............................. 61
3.5 Summary ...................................................................................................................... 64
Chapter 4: Ion conductance of nanoscale electropore 66
4.1 Abstract ........................................................................................................................ 66
4.2 Introduction .................................................................................................................. 67
4.3 Methods........................................................................................................................ 69
4.3.1 Molecular dynamics simulations ......................................................................... 69
4.3.2 Structure ............................................................................................................. 70
4.3.3 Conductance simulations ..................................................................................... 70
4.3.4 Pore radius measurement .................................................................................... 71
4.3.5 Current measurement .......................................................................................... 72
4.3.6 Conductance calculation ..................................................................................... 73
4.3.7 Electropore ion concentration measurement ........................................................ 73
4.4 Results and Discussion ................................................................................................. 74
4.4.1 Ion interactions with the phospholipid interface .................................................. 74
4.4.2 Pore formation and stabilization .......................................................................... 75
4.4.3 Na
+
, K
+
, and Cl
-
conductance .............................................................................. 77
4.4.4 Molecular simulations and experimental data ...................................................... 83
4.4.5 Effects of Water model on conductance simulation ............................................. 85
4.5 Summary ...................................................................................................................... 87
Chapter 5: Effects of monovalent ions on phospholipid bilayer and electroporation 88
5.1 Abstract ........................................................................................................................ 88
5.2 Introduction .................................................................................................................. 89
5.3 Methods........................................................................................................................ 92
5.3.1 Molecular dynamics simulations ......................................................................... 92
5.3.2 Structure ............................................................................................................. 93
5.3.3 Evaluation of ion distribution, ion absorption, and bulk ion concentration ........... 93
5.3.4 Electroporation simulations and pore initiation time measurement ...................... 94
5.3.5 Conductance simulations ..................................................................................... 94
5.3.6 Pore radius measurement .................................................................................... 95
5.3.7 Current measurement .......................................................................................... 95
5.3.8 Conductance calculation ..................................................................................... 96
5.4 Results and Discussion ................................................................................................. 97
5.4.1 Ion distribution and ion absorption in lipid bilayer .............................................. 97
5.4.2 Ions’ effect on area per lipid .............................................................................. 101
vii
5.4.3 Water phase separation and pore initiation ......................................................... 103
5.4.4 Ions’ effect on field-stabilized pore radius .......................................................... 106
5.4.5 Electropore conductance in different ion concentration ...................................... 108
5.4.6 Ions’ effect on POPS translocation ..................................................................... 110
5.5 Summary .................................................................................................................... 115
Conclusions and Outlook ...................................................................................................... 117
Bibliography .......................................................................................................................... 122
viii
List of Tables
Table 2.1 Parameters of 3-sites water models............................................................................ 24
Table 3.1 Water bridge initiation time using SPC, SPC/E, TIP3P model ................................... 62
Table 4.1 Electrophoretic mobilities and conductivities in bulk water ....................................... 78
Table 4.2 NaCl and KCl conductance in POPC electropores ..................................................... 80
Table 4.3 KCl concentration inside electropore from MD simulation ........................................ 84
Table 5.1 Bulk ion concentrations in pure POPC bilayer system ............................................. 100
Table 5.2 Bulk ion concentrations in POPC/POPS bilayer system........................................... 101
Table 5.3 Number of translocated POPS toward anode leaflet by external field ...................... 113
Table 5.4 Number of translocated POPS toward cathode leaflet by POPS-cation grouping ..... 114
ix
List of Figures
Figure 1.1 Schematic of electrochemotherapy ............................................................................ 2
Figure 1.2 Microscopic photograph a peel tissue of wine grapes before and after electroporation,
picture of electroporation device in a winery ............................................................................... 3
Figure 1.3 Microscopic photograph of single cell conductance measurement using patch clamp,
electric current versus applied potential obtained by patch clamp ............................................... 4
Figure 1.4 Intakes of YO-PRO-1 and propidium iodide into pulsed cell in fluorescence
experiments ................................................................................................................................. 5
Figure 1.5 Fluorescence images of phosphatidylserine (PS) externalization under intense
nanosecond pulses in experiment................................................................................................. 6
Figure 1.6 Images of Jurkat cells before after electric pulses exposure in cell swelling
experiment .................................................................................................................................. 7
Figure 1.7 Lipid structures of a hydrophobic pore and a hydrophilic pore ................................... 9
Figure 1.8 Energy functions of a transmembrane pore in different phases ................................. 10
Figure 1.9 Energy function (, ) of a pore under different transmembrane potential ......... 11
Figure 2.1 Cross-sectional view of lipid bilayer composed of POPC ......................................... 23
Figure 2.2 Structure of POPC represented by united-atom model ............................................. 27
Figure 2.3 Illustration of modified Lennard-Jones potential functions under different truncation
schemes ... ................................................................................................................................ 35
Figure 2.4 Illustration of the life cycle of an electropore ........................................................... 39
Figure 2.5 Dipole configuration of two water molecules resembles the flat water surface, and
the water protrusion .................................................................................................................. 43
Figure 3.1 Illustration of water bridge radius measurement ....................................................... 52
Figure 3.2 Internal versus external electric field in W-V-W and POPC system ......................... 54
Figure 3.3 Snapshots of water bridge formation in POPC and W-V-W systems ........................ 55
Figure 3.4 Water bridge initiation time versus external electric field strength ........................... 55
Figure 3.5 Total system dipole moment of water in z direction versus time during water bridge
formation .................................................................................................................................. 57
x
Figure 3.6 Time-averaged water bridge radius versus sustaining field ....................................... 58
Figure 3.7 Evolution of water bridge radius during stabilization ............................................... 59
Figure 3.8 Water bridge annihilation time versus initial water bridge radius in POPC system and
W-V-W system ......................................................................................................................... 61
Figure 3.9 Water bridge initiation time using SPC, SPC/E, TIP3P model ................................. 63
Figure 4.1 Illustration of pore radius measurement in a POPC electropore ................................ 72
Figure 4.2 Main ion binding sites of POPC, radial distribution function g(r) between both
glycerol backbone acyl-oxygens and cations ............................................................................. 74
Figure 4.3 Pore radius versus time from conductance simulations, Time-averaged pore radius
versus sustaining field E
s
........................................................................................................... 77
Figure 4.4 Ion conductance versus pore radius .......................................................................... 79
Figure 4.5 Molar pore conductance of K
+
using GROMACS parameters and CHARMM
parameters ................................................................................................................................ 81
Figure 4.6 Snapshots of ion interactions during migration through an electropore ..................... 82
Figure 4.7 Simulated electrophoretic mobility of Na
+
, K
+
and Cl
-
in bulk water versus ion
concentration............................................................................................................................. 82
Figure 4.8 Comparison of simulated KCl conductance with chronopotentiometry data ............. 83
Figure 4.9 Velocity of single Na
+
and K
+
in SPC and SPC/E water under electric field ............. 86
Figure 5.1 Ion distributions in z direction during equilibration .................................................. 98
Figure 5.2 Number of cations absorbed by the lipid bilayer versus total number of cations in the
system ...... ................................................................................................................................ 99
Figure 5.3 Area per lipid under different number of cations in the system ............................... 102
Figure 5.4 Water group separation with different number of cations in the system .................. 103
Figure 5.5 Pore initiation time under E
ext
= 350 MV/m with different number of cations in the
system ...... .............................................................................................................................. 104
Figure 5.6 Water group separation in the system and the corresponding pore initiation time under
E
ext
= 350 MV/m ..................................................................................................................... 105
Figure 5.7 Time-averaged pore radius versus number of Cl
-
of field-stabilized electropore under
E
sus
= 75 MV/m ....................................................................................................................... 107
Figure 5.8 Pore conductance of cations versus number of cations under E
sus
= 75 MV/m ........ 109
Figure 5.9 Pore conductance of Cl
-
versus number of Cl
-
under E
sus
= 75 MV/m .................... 110
xi
Figure 5.10 Snapshots of POPS translocation toward anode leaflet by external field ............... 112
Figure 5.11 Snapshots of POPS translocation toward cathode leaflet by POPS-cations grouping
under external field ................................................................................................................. 112
xii
Abstract
Electroporation provides a controllable method to introduce foreign substances into living
cells. It is widely used by researchers in cell biology and the medical field to manipulate
biological systems at the cellular level. For decades, electroporation has been studied extensively
through experiments and theoretical models, and electroporation-based technologies have been
improved substantially with these efforts. One of the issues in utilizing electroporation is the lack
of understanding in the phenomenon’s molecular mechanism and the microscopic details, mainly
due to the difficulty in the direct experimental observation of the nanosecond-scale electropore
formation process and the nanometer-sized electropore structure. To overcome this issue,
Molecular Dynamics (MD) simulation has become one of the major tools to study
electroporation at the microscopic level.
Recent advancements of high performance computing, such as the increase in processing
power, developments in algorithms and parallelization, have improved the efficiency of MD
simulation substantially. Due to these advancements, MD simulation has become a popular tool
for studying systems that are composed of biomolecules. For nearly a decade of effort, MD
simulation revealed many different aspects of electroporation and it provided a molecular
description of the process. Using MD simulation, we are able to observe the events during the
electropore formation and annihilation, as well as the transport processes of molecules through
the electropore. In addition, MD simulation provides a platform to study the molecular structure
of electropore, and the associated energetic.
xiii
My dissertation is organized as the follows: Chapter 1 provides the motivation of this
research by discussing the applications of electroporation-based technology, electroporation
experiments, and the existing continuum model that describes electroporation. Chapter 2
introduces the MD formalism, models, and various algorithms used in our MD simulation. I will
also discuss some of the previous MD studies of electroporation and their significances at the
end of Chapter 2. In Chapter 3, we will examine the dynamics of water bridge during different
stages of electroporation, and delineate the role of water molecules in electroporation by
comparing the lipid bilayer system with an artificial water-vacuum-water system. In Chapter 4,
we will examine the steady state of an electropore in lipid bilayer and evaluate the pore
conductance of ions. The pore conductance values obtained from the simulations can be
compared with those obtained by experiments. In Chapter 5, we will examine the effects that
monovalent ions impose on lipid bilayer and electropore formation. We will also examine the
pore conductance of ions under various ion concentrations, and the PS translocation process. At
the end, I will summarize the findings in our research and provide a short outlook on MD
simulation in the study of electroporation.
1
Chapter 1: Introduction
1.1 Applications of electroporation
Cell membrane functions as a barrier to separate cell interior from the environment. In
order to govern different biological functions of the cell (such as cell signaling and adhesion),
cell membrane regulates the substances that enter and exit the cell through various mechanisms
(such as osmosis, protein channels, endocytosis and exocytosis). When a cell is exposed to
sufficiently large electric field, structures known as electropores are formed on the cell
membrane, and they allow normally blocked substances to travel into and out of the cell. This
phenomenon is known as Electroporation (or Electropermeabilization). Electroporation provides
a controllable method to manipulate permeability of the cells, and enable many applications for
biological experiments and medical treatments.
The use of electroporation to deliver genetic material into cells can be dated back in the
early 80s, where plasmid DNA was introduced into the mouse lyoma cells using intense
(~kV/cm) and short (~µs) electric pulses. In the experiment, the pulsed, porated cells that
absorbed the plasmid DNA exhibited higher viability in a normally harmful environment [93,
154]. Over the years, electroporation has been utilized to introduce foreign genes into cells in
both in vitro [23, 66, 94, 116] and in vivo [2, 21, 47, 102] experiments. The development of
electroporation-based technology enables researchers to manipulate the structures of cells, and
enhances our understanding in their functionalities.
Electroporation-based technology is currently of great interest in medical research. One
of its most important medical applications is electrochemotherapy. In conventional
2
chemotherapy, chemotherapeutic agents not only kill the cancer cells, but they are also harmful
to other normal cells as well. Patients who undergo conventional chemotherapy often experience
distressful side effects such as decrease in blood cell production, inflammation, and hair loss due
to the casualties from the chemotherapeutic agents. Electrochemotherapy overcomes this issue
by utilizing short and intense pulses to control the doses and localize the absorption of
chemotherapy drugs in the vicinity of the cancer cells [8, 39, 115] (Figure 1.1). In particular,
some of the electrochemotherapy drugs, such as bleomycin, are impermeant to unporated cells
[131]. The use of these electrochemotherapy drugs can prevent the casualties of regular healthy
cells and reduces the undesired side effects during the medical treatment.
Figure 1.1 Schematic of electrochemotherapy [115].
Beside biological and medical research, the usefulness of electroporation is also found in
other applications such as food preservation and food processing. For instance, although the
3
conventional heating method is effective for killing spoilage microorganisms, the thermal effects
often lowers the sensory and nutritional qualities of the food. Electroporation provides an
alternate nonthermal method to inactivate the spoilage microorganisms. Numbers of studies have
shown that juices under electroporation treatment possess higher sensory and nutritional quality
compare to those under conventional heating treatment [31, 83]. Electroporation-based
technology can also be used to extract cell components to enhance food processing. In a
particular study, electric pulses were utilized to extract flavoring substances and pigments from
grapes mash (Figure 1.2) for wine production [113]. The same technique can be applied on other
oil or protein rich plants to extract the nutritional substances. These applications are some of the
examples on how the developments in electroporation-based technology can improve life quality
substantially.
Figure 1.2 Microscopic photographs of a peel tissue of wine grapes a) before and b) after
electroporation. c) Picture of electroporation device in a winery. Images were obtained from ref
[113].
4
1.2 Electroporation experiments
Due to the nanometer size of the electropore, and the nanosecond scale of the pore
formation process, direct observation of electroporation is currently inaccessible by experiment.
Observations of electroporation are based on the changes in the properties of cell membrane.
One of the examples is the change of electric conductance. Although initially considered as the
electric breakdown [157], the rapid increase of electric conductance under sufficiently large
applied voltage was one of the earliest indicators of electropore formations on the cell membrane
[120, 158]. Conductance measurement provides real time information on the kinetics of the pore
opening and closing, and it has been utilized in different experimental settings such as individual
cell [98, 112] (Figure 1.3), cell culture [101], and planar lipid bilayer membrane [65, 70, 71]. By
using ion selective electrode, the selectivity of ions through the electropores can also be
examined by the conductance measurement [99].
Figure 1.3 a) Microscopic photograph of single cell conductance measurement using patch
clamp. b) Plot of electric current versus applied potential obtained by patch clamp, displays rapid
increase of conductance when the applied potential is beyond certain threshold. Both figures
were obtained from ref [112].
5
Another method to observe electroporation is by monitoring the intake of foreign
molecules during the cell’s permeabilized state. For example, fluorescent dyes such as propidium
iodide (PI) and YO-PRO-1 are normally blocked by the cell membrane, but they can travel into
the cell interior after the cell is electro-permeabilized. When these dyes are absorbed by the
porated cell and bind with the nucleic acids, they produce strong fluorescent signals (Figure 1.4)
[143]. Fluorescence intensity in these experiments can be analyzed to extract information such as
the pore density distribution [17, 98]. In addition, voltage-sensitive dyes such as RH292, can be
used to probe the potential profile of the cell, and the decrease in transmembrane potential is also
one of the signatures of a cell undergoing electroporation [53, 54, 67]. Other effects from
electroporation, such as calcium burst, can also be observed by utilizing fluorescent dyes such as
Calcium Green [140].
Figure 1.4 Intakes of YO-PRO-1(left column) and propidium iodide (right column) into pulsed
cell in fluorescence experiments. Images were obtained from ref [143].
6
Structural changes in the cell caused by electroporation can also be used to characterize
electroporation. One of such structural changes is phosphatidylserine (PS) externalization [141].
PS is commonly found in the inner leaflet of cell membrane, and the exposure of PS is one of the
indicators of the cell undergoing apoptosis [139]. Using an intense electric field, PS can be
artificially translocated to the outer leaflet of cell membrane, and the translocated PS can be
observed using fluorescent-tagged PS-binding proteins such as annexin V and FM1-43 [144]
(Figure 1.5).
Figure 1.5 Fluorescence images of phosphatidylserine (PS) externalization under intense
nanosecond pulses. Separate cells were stained with annexin V-FITC (a,c) and FM1-43 (b,d).
Fluorescence intensity exhibits asymmetry of PS externalization between cathode and anode side
of the cell. Images were obtained from ref [144].
7
Cell swelling is another effect that is mediated by electroporation [91, 109]. The
osmolality of a cell is determined by the concentrations of the biomolecules, which are normally
maintained at the physiological values by different protein channels on the cell membrane. When
a cell is electro-permeabilized, the nanoscale electropores allow small molecules such as Na
+
, K
+
and Cl
-
to permeate through the cell membrane, while larger molecules remain trapped in the
cell’s interior. This change of cytoplasm content introduces osmotic imbalance between the cell
interior and exterior. To balance the osmolality, water influxes into the cell in order to dilute the
osmolytes, and causes the cell to swell (Figure 1.6). Cell swelling provides another way to
characterize electroporation, and the change in the cell volume during the swelling provides a
sensitive method to evaluate the electropore sizes [109].
Figure 1.6 Images of Jurkat cells before (A, C) and 60 s after electric pulses exposure (B, D) in
cell swelling experiment. Images were obtained from ref [109].
8
1.3 Continuum models of electroporation
The continuum model of electroporation describes the process as the phase transition
between the nonconductive hydrophobic pore (Figure 1.7a) and the conductive hydrophilic pore
(Figure 1.7b). This model is sometimes referred as the ‘stochastic pore model’ [105, 121]. A
hydrophobic pore is considered as a gap in the lipid membrane, and the energy function ()
associated with a hydrophobic pore is given by
() = 2
ℎ
() (1.3.1)
() =
(∞)
(/)/
(/) (1.3.2)
where is the pore radius, ℎ is the membrane thickness, and
is the interfacial tension between
hydrophobic lipid tails and water, and
() is the n
th
order modified Bessel function with
characteristic length [38, 151]. In the limit of ≪ , the tension term
() can be
approximated by
() =
(∞)/2 and the energy function () can be simplified as the
quadratic function of [92].
In the hydrophilic pore configuration, due to water inside the pore, phospholipid head
groups rotate toward the pore interior to a form pore wall in attempt to shield the hydrophobic
lipid tails from the water. The energy function () associated with a hydrophilic pore is given
by
() = 2
−
+ () (1.3.3)
where is the energy per unit length due to the pore line tension, and is the energy per unit
area of the intact membrane [1, 38, 124, 151]. For a hydrophilic pore with small radius ( ≪ ℎ),
9
the lipid head groups experience steric repulsions from the close packing. The terms () in
equation (1.3.3) represents the energy gain due to the steric repulsions. In the asymptotic model
of Neu and Krassowska [92], the steric repulsion term () is given by () = !"# " ∗
%&
.
The overall energy function () of a pore in the radius space is given by
() = min ((), ()). Figure 1.8 illustrates this energy landscape in the pore radius space.
There are three critical radii in the pore energy function (): The first critical radius
∗
, is
located at the point where () = (), and this is the point where the phase transition occurs
between the hydrophobic pore state and the hydrophilic pore state. The second critical radius
*
is located at the local minima of (), and this being the equilibrium radius of the metastable
hydrophilic pore state. The third critical radius
+
is located at the local maxima of (). If the
radius of a pore is larger than
+
, then the pore will keep expanding and eventually rupture the
membrane (a process known as irreversible breakdown).
Figure 1.7 Lipid structures of a) a hydrophobic pore and b) a hydrophilic pore.
10
Figure 1.8 Energy functions of a transmembrane pore in different phases. The red curve
corresponds to the energy function () of a hydrophobic pore. The blue curve corresponds to
the energy function () of a hydrophilic pore. The combined solid curve represents the overall
energy function () of a pore.
When a cell membrane is placed under an external electric field, charges accumulate on
the membrane surfaces and the membrane is charged as a parallel plate capacitor. If a pore has
formed, part of the membrane is replaced by the aqueous medium and some of the surface charge
is lost. Formation of the pore lowers the capacitance of the membrane, and the corresponding
energy decrease , (, ) is given by
, (, ) = −
-.
(1.3.4)
-.
= /
0
1
0
2
− 14
(1.3.5)
where is the spatially averaged transmembrane voltage, 5
6
is the relative permittivity of water,
5
*
is the relative permittivity of the membrane, and
is the specific capacitance of the pore-
11
free membrane [151]. Typically 5
6
7 80 and 5
*
7 2, so , (, ) decreases with increasing .
Adding this contribution, the pore energy function (, ) is given by
(, ) = () −
-.
(1.3.6)
Figure 1.9 demonstrates the change of energy function (, ) under various voltage [92]. For
sufficiently small , (, )has the same features as (), where the three critical radii
∗
,
*
,
and
+
, can be identified. In this regime, increasing shifts
*
to a larger value and lowers the
energy at
*
, which means larger hydrophilic pores are allowed to form on the membrane and
they are more stable. For sufficiently large ,
*
and
+
merge together and disappear. In this
case, there is no metastable state of the hydrophilic pore. If a hydrophobic pore evolves into a
hydrophilic pore under such voltage, the pore will keep expanding until the membrane
experiences irreversible breakdown.
Figure 1.9 Energy function (, ) of a pore under different transmembrane potential . The
figure was obtained from ref [92].
12
Utilizing the energy function (, ), the pore density function (, ") can be evaluated
using the Smoluchowski equation [151]:
:
:;
− < /
:
:=
/
>
?
@
:A
:=
4 +
:
B
:=
B
4 = C() (1.3.7)
where < is the diffusion coefficient of a pore in the radius space, D
E
is the Boltzmann constant, F
is the absolute temperature, and C() is a source term that represents the pore creation and
annihilation [92, 150]. In the asymptotic model of Neu and Krassowska, the Smoluchowski
equation (1.3.7) is simplified to an ordinal differential equation [92]. Using this formalism,
numerical studies were carried out to investigate various aspects of electroporation, such as the
electropore distribution in the pore radius space [64], the electropore density distribution on a
single cell [28, 29] and the corresponding transmembrane potential profile [107, 118].
Although the continuum model provides a valid theoretical description on
electroporation, there are several drawbacks: One of the main issues is that the equations used in
the model require empirical constants that cannot be measured directly by experiment. These
constants were estimated by theoretical arguments and some of them are only accurate within an
order of magnitude. Secondly, some of the variables in the continuum model, such as the pore
radius and pore number density, cannot be measured directly by experiment, which makes the
validation of the continuum model difficult. Furthermore, the lack of molecular details prohibits
the use of continuum model to study many factors that can affect electroporation, such as the
molecular composition of the cell membrane, the interactions between the electropore and
foreign molecules, and the dynamics of molecules transport through the electropore. These issues
prevent a detail description of electroporation by the continuum model. Molecular simulation, on
the other hand, can be use to resolve some of these issues. By combining the advantage of
13
molecular simulation and the continuum model, we can acquire a better theoretical view on
electroporation, such as a better evaluation of the pore energy [153].
14
Chapter 2: Molecular Dynamics Simulation
2.1 Overview
Molecular dynamics (MD) simulation utilizes numerical integration of Newton’s Second
Law to obtain trajectories of atoms and molecules within a microscopic system. Development of
the MD simulation began in the 1950s [3], with one of the earliest physical systems studied by
the method was the liquid phase of argon, and quantities such as the pair correlation function,
velocity autocorrelation function, and diffusion coefficient were extracted by using system of
864 argon atoms with periodic boundary condition [108]. The main advantage of MD simulation
is that it allows real time observation of events that are not easily observed by experiment. These
events are often occur within the nanosecond and nanometer scale. Beside dynamical study, MD
simulation can also be utilized to study the equilibrium properties of a microscopic system. In
particular, MD simulation provides a clear view on a system’s molecular structure and its effects
on the thermodynamic properties (such as the critical temperature and the order parameters).
Some of these microscopic systems with precise molecular structure (such as a pure lipid bilayer)
are difficult to be attained by experiment.
Over the years, MD simulation has become a popular tool for researchers in different
fields such as material science, physical chemistry, and molecular biology. This popularity is
credited to the astonishing improvement on the method’s efficiency. For instance, computers
nowadays are orders of magnitude faster than the computers decades ago. The developments in
algorithms improved the efficiency of the method by reducing the required computations. Some
15
examples of these algorithms are Particle Mesh Ewald (PME) method [26, 32] and Fast
Multipole Method (FMM) [41] for long-range electrostatic forces evaluation. Another
breakthrough on the MD simulation efficiency is the employment of parallel computing.
Protocols such as Message Passing Interface (MPI) are utilized to decompose the system into
subsystems for parallelization. Developments of network architectures (such as Myrinet and
InfiniBand) also improved the efficiency of MD simulation by reducing the latency between the
processors’ communications during the simulation. With these achievements, researchers
nowadays can easily perform MD simulation for system consists over tens of thousands of
atoms, and in some cases, millions of atoms [110].
There are currently many packages available for MD simulation of biological system.
Some of the popular packages are AMBER [103], CHARMM [18], GROMACS [12, 50], and
NAMD [90]. The primary difference between these packages is the force field parameters used
during their developments (such as all atoms model versus united atoms model, SPC water
model versus TIP3P water model). The main tool we used in our research is the GROMACS
(GROningen MAchine for Chemical Simulations) package. It is a free, open source package
released under the GNU General Public License, and it is used worldwide to study systems that
contain biomolecules such as proteins, lipids, and nucleic acids. The GROMACS package is not
only limited for MD simulation, it can also perform energy minimization, and provides selection
of tools to analyze the molecular structures.
Like any other analytical tool, there are limitations for MD simulation. First, by using
Newton’s Second Law, the nature of MD simulation is purely classical. Quantum effects (such as
tunneling) are not included in the simulation. Luckily, this is usually not a serious concern for
biological systems, since the temperature of interested is usually room temperature (~300 K) or
16
physiological temperature (~310 K), and the behaviors of atoms under these temperatures are
well described classically.
Another limitation of MD method is that the electron densities of molecules are
predefined in the force field parameters, and the charges on the molecules are fixed during the
simulation with the assumption that electrons are in their ground states. This limitation excludes
effects such as electron transfer, electron excitation, and chemical reaction between molecules
during the simulation. Molecules in MD simulation by default also lack polarizability. Although
the effects of polarizability can be introduced by various methods such as adding virtual sites
(Drude particles) or effective pair potential, the treatment of polarizability is often ignored in
order to reduce computing cost.
Beside these methodological limitations, the restriction of computing resources also
introduces additional limitations. For instance, Lennard-Jones interactions are usually being
cutoff beyond certain range and the long-range electrostatic interactions are usually
approximated by Ewald summation in order to reduce computations. In order to facilitate MD
simulation to study the molecular structure that resembles macroscopic system, periodic
boundary condition is often imposed. Artifacts can be introduced by the period images if the
boundary condition is not handled carefully [45]. Challenges to these limitations are currently
topics in scientific computing research.
17
2.2 MD formalism and implementation
2.2.1 MD discretization and integrators
The formalism of MD simulation is based on the Taylor expansions on the particles’
trajectories. For example, consider the following Taylor expansions of a particle’s position:
G(" + H") = G(") + GI(")H" +
GJ(")H"
+
K
G L(")H"
M
+ N(H"
&
) (2.2.1)
G(" − H") = G(") − GI(")H" +
GJ(")H"
−
K
G L(")H"
M
+ N(H"
&
) (2.2.2)
where G is the 3-dimensional position vector, and H" is a small time increment. Adding equation
(2.2.1) and (2.2.2) gives the discretization of position:
G(" + H") = 2G(") − G(" − H") + GJ(")H"
+ N(H"
&
) (2.2.3)
and subtracting equation (2.2.2) from (2.2.1) gives the discretization for velocity:
GI(") =
G(;OP;)%G(;%P;)
QP;
+ N(H"
) (2.2.4)
The acceleration GJ(") in equation (2.2.3) can be obtained from Newton’s Second Law. By
iterating equation (2.2.3) and (2.2.4) in every time step H", the position and velocity of the
particle can be obtained at each time step. The above scheme is known as the Verlet integration.
Notice equation (2.24) is not needed to obtain the trajectory of the particle, and the initial
conditions of Verlet integration are given by G(−H") and G(0), which means the velocity can be
totally omitted in the Verlet integration. This is an issue for simulating a system that is coupled
18
to an external heat bath, since the temperature of a system is provided by the velocities of the
particles based on the Boltzmann distribution.
To incorporate the velocities of particles in the integration, an alternative discretization is
proposed by utilizing the Taylor expansions of position and velocity at half time step H"/2:
G /" +
P;
4 = G(") + GI(")
P;
+
GJ(") /
P;
4
+ N(H"
M
) (2.2.5)
G /" −
P;
4 = G(") − GI(")
P;
+
GJ(") /
P;
4
+ N(H"
M
) (2.2.6)
GI /" +
P;
4 = GI(") + GJ(")
P;
+
G L(") /
P;
4
+ N(H"
M
) (2.2.7)
GI /" −
P;
4 = GI(") − GJ(")
P;
+
G L(") /
P;
4
+ N(H"
M
) (2.2.8)
Subtracting equation (2.2.6) from (2.2.5), we have:
G /" +
P;
4 = G /" −
P;
4 + GI(")H" + N(H"
M
) (2.2.9)
and with the change of variable " ← " +
P;
on equation (2.2.9), the discretization of position is
given by:
G(" + H") = G(") + GI /" +
P;
4 H" + N(H"
M
) (2.2.10)
Subtracting equation (2.2.8) from (2.2.7) gives the discretization of velocity:
GI /" +
P;
4 = GI /" −
P;
4 + GJ(")H" + N(H"
M
) (2.2.11)
19
The algorithm utilizes the discretization (2.2.10) and (2.2.11) is known as the Leap-frog
integration. At each time step, the integrator first calculates the velocity GI /" +
P;
4 using
equation (2.2.11), and then the position G(" + H") is obtained using equation (2.2.10). The name
“Leap-frog” was given because it resembles the children’s game Leap-frog, where one child
(position) leap over the other child (velocity)’s back at each step. Using the Leap-frog
integration, the initial conditions are given by GI /
P;
4 and G(0), and the initial temperature of the
system can be adjusted by assigning the initial velocities of particles at
P;
based on the
Boltzmann distribution. Coupling to external heat bath can also be imposed by scaling the
particles’ velocities at each half time step.
Another scheme that incorporates the velocities of particles in the simulation is the
Velocity Verlet integration. Consider the following Taylor expansions on the position, velocity,
and acceleration:
G(" + H") = G(") + GI(")H" +
GJ(")H"
+ N(H"
M
) (2.2.12)
GI(" + H") = GI(") + GJ(")H" +
G L(")H"
+ N(H"
M
) (2.2.13)
GJ(" + H") = GJ(") + G L(")H" + N(H"
) (2.2.14)
With a little algebra, equation (2.2.14) can be rewritten as:
G L(")H"
=
P;
SGJ(" + H") − GJ(")T + N(H"
M
) (2.2.15)
Substitute (2.2.15) into (2.2.13), we have the following discretization for velocity:
GI(" + H") = GI(") +
P;
SGJ(" + H") + GJ(")T + N(H"
M
) (2.2.16)
20
Define the following term as the velocity at half time step:
GI /" +
P;
4 ≡ GI(") + GJ(")
P;
(2.2.17)
Equation (2.2.12) and (2.2.16) can be rewritten as:
G(" + H") = G(") + GI /" +
P;
4 H" + N(H"
M
) (2.2.18)
GI(" + H") = GI /" +
P;
4 + GJ(" + H")
P;
+ N(H"
M
) (2.2.19)
With the above discretization, the Velocity Verlet integrator first obtains the acceleration GJ(")
from Newton’s Second Law and finds the velocity GI /" +
P;
4 using equation (2.2.17). The new
position G(" + H") is obtained by using equation (2.2.18), and with the new position G(" + H"),
the acceleration GJ(" + H") can be calculated again using Newton’s Second law. Finally, the
velocity GI(" + H") is calculated using equation (2.2.19) to complete the integration cycle. The
advantage of the Velocity Verlet algorithm is that the positions and velocities of particles are
both obtained in every full time step, in contrast to the Leap-frog algorithm, where the positions
are obtained in every full time step and the velocities are obtained in every half time step. The
disadvantage of the Velocity Verlet integration is that it requires extra computing cost for the
calculations of velocities at half time steps.
Both Leap-frog and Velocity Verlet algorithm are implemented in the GROMACS
package. Since we mostly examine the molecular structure of the lipid bilayer and rarely perform
analysis in the velocity space, the Leap-frog integration is more efficient for our research.
21
2.2.2 Energy function in biological systems
In MD simulation, the acceleration of a given coordinate
V
is obtained by Newton’s
Second Law:
V
J = −
*
W
:
:=
X
(2.3.1)
where Y
Z
is the mass of the particle associated with the coordinate
V
and is the potential
energy function of the system. For system that consists of biomolecules, the potential function
has the following general form:
= ∑
&\0
]
W
]
^
=
W^
Z_`
+ ∑
a
W^
=
W^
bB
−
E
W^
=
W^
c Z_`
+ ∑
D
Z`
d
S
Z`
− e
Z`
T
+ ∑
D
Z`>
f
Sg
Z`>
− g
Z`>
T
hijkl d
+l
+ ∑ D
m
/1 + n!S (o − o
)T4
+Zpk+=hjl
(2.3.2)
The interactions between a pair of molecules are provided by the non-bonded (intermolecular)
interactions. The major non-bonded interactions in MD simulation are the Coulomb
(electrostatic) interaction, represented by the first term in equation (2.3.2) and the Lennard-Jones
interaction, represented by the second term in equation (2.3.2). The physical structure of a given
molecule is characterized by the bond lengths, bond angles (defined by triplet of atoms) and the
dihedral angles (defined by the angle between two normal vectors from the planes constructed by
atom triples). The energy cost of altering a molecule’s physical structure is given by the energies
of the bonded (intramolecular) interactions. Namely, the third term in equation (2.3.2) represents
the energy cost of bond stretching, the fourth term represents the energy cost of bond angle
22
bending, and the fifth term represents the energy cost of dihedral angle bending. In general, the
parameters in equation (2.3.2) characterized the molecular structure of a given biological system.
23
2.3 Molecular structures and models
The main structure studied in electroporation simulation is the lipid bilayer. Lipid bilayer
is composed of two leaflets of lipids. Lipids on each leaflet are oriented so that their hydrophobic
tails face the opposite leaflet and shielded from water (see Figure 2.1). In molecular simulation,
water molecules are also included in the system to ensure the hydrations of the lipids. Lipid
bilayer is one of the simplest structures that resemble the cell membrane by providing the basic
barrier function. Pure lipid bilayer such as black lipid membrane (BLM), can also be synthesized
artificially and their properties have been well studied by experiments [73, 88, 20], and
simulations [24, 30, 126].
Figure 2.1 Cross-sectional view of lipid bilayer composed of POPC. The normal vector of the
lipid bilayer is approximately in ẑ. In the figure, the lipid tails of POPC are represented by light
blue stripes. Oxygen, nitrogen and phosphorus atoms of POPC are represented by red, blue and
gold spheres accordingly to indicate the locations of the lipid head groups. Hydrogen and oxygen
of water molecules are represented by white and red dots.
24
2.3.1 Water models
Interactions between biomolecules and water molecules are vital for many biological
functions. In molecular simulation, water molecules also play an important role in maintaining
the proper structure of the biomolecules. There are various water models available for theoretical
study. Some of these models are more complex than others in terms of the number of sites used
to represent a water molecule. 3-sites models are popular for MD simulation due to their
simplicities. In particular, the constraint force of a 3-sites molecule can be solved analytically
and calculated quickly by the SETTLE algorithm [87], the use of 3-sites water model in MD
simulation substantially reduces the computing cost.
Some of the popular 3-site water models are the Simple Point Charge (SPC) model [9],
the TIP3P model [61], and the Extended Simple Point Charge (SPC/E) model [11]. The general
features of these models are the same: Atoms of a water molecule are represented as classical
point particles. Negative charge is assigned on oxygen and positive charge is assigned on
hydrogen, in order to produce the dipole moment of water, and the Lennard-Jones parameters are
assigned on the oxygen atom. Parameters of these models are listed in Table 2.1.
Table 2.1 Parameters of 3-sites water models.
SPC TIP3P SPC/E
r
OH
(Å) 1.0 0.9572 1.0
HOH angle 109.47º 104.52º 109.47º
Lennard-Jones parameter A
(10
-3
kcal Å
12
/mol)
629.4 582.0 629.4
Lennard-Jones parameter B
(10
-3
kcal Å
12
/mol)
625.5 595.0 625.5
Charge on O (e
-
) -0.82 -0.834 -0.8476
Charge on H (e
-
) 0.41 0.417 0.4238
25
Parameterization of the water models were optimized by comparing the results obtained
from simulation to the known properties of water (such as density, diffusion coefficient) and the
radial distribution function obtained from x-ray diffraction [89]. None of the 3-sites water
models mentioned above can reproduce the behavior of real water perfectly, because these
models are lack of polarizability. It is difficult to determine which model is superior to the
others. Each of the models has its own characteristic and it could be more accurate than the
others in a given particular study. For example, the TIP3P model uses the actual bond length and
bond angle of water from experiment, but the difference between its diffusion coefficient and the
known experimental value is larger than the other two models [79]. The SPC/E model has better
diffusion coefficient compares to the other two models, but the difference between its dipole
moment and the experimental value is largest amongst the three.
In our research, we mostly use the SPC water model. It is the default water model for the
GROMACS package and it was used during the modeling of many biomolecular structures,
including the lipid bilayers. We had briefly compared electroporation simulations with the SPC
and the SPC/E model, and the behaviors of water molecules during electropore creation process
are similar regardless of which water model is used. However, quantitative measurements from
the simulation (such as the time scale of electroporation creation, electropore conductance), can
be affected by the difference in the properties of water model.
26
2.3.2 Lipid models
Molecular models of biomolecule can be broadly divided into three categories: all-atom
model, united-atom model, and coarse-grained model. In all-atom model, coordinates of all the
atoms in a molecule are explicitly included to represent the molecule. By using an all-atom
model, full atomistic details of the system are provided. However, for many the studies, some of
these atomistic details are not significant, and they become burdens on the computation. To
improve the efficiency, united-atom models combine the small atoms (such as hydrogen) with
their larger neighboring atoms (such as carbon) into pseudo-atoms in order to reduce the total
number of particles in the system. In general, this is a good approximation if the intramolecular
motions of the system are less important than the intermolecular motions. Coarse-grained model
further reduce the number of particles in the system by combining some of the larger atoms into
large pseudo-atom [80]. It is mainly used for simulation of larger system that require high
efficiency. The choice of molecular model depends on the granularity of the specific study, and
often there is a trade-off between efficiency and resolution.
In our research, we employ the united-atom model in Optimized Potential for Liquid
Simulations (OPLS) force field [62], with the parameters developed by Berger et al. [13]
(optimized under SPC water model). The majority of our simulations are performed on lipid
bilayer composed of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine (POPC, see Figure
2.2). POPC is one of the main constituents of the cell membrane’s outer leaflet and it is one of
the common lipids studied in biophysical experiments. The structure and topologies of the POPC
bilayer can be obtained from the webpage of the Biocomputing Group at University of Calgary
(http://wcm.ucalgary.ca/tieleman/). Electroporation in other types of lipid bilayer were also
studied extensively with MD simulation [104, 123, 128, 156].
27
Figure 2.2 Structure of POPC represented by united-atom model. The constituents of POPC are
nitrogen (blue sphere), phosphorus (gold sphere), oxygen (red sphere) and carbon (light blue
sphere). Hydrogen atoms are combined with the larger neighboring atoms in the united-atom
model.
28
2.4 Holonomic constraints
One of the main control parameters of MD simulation is the integration time step H". In
general, MD simulation is more efficient with larger H". However, H" is limited by the highest
vibrational frequency in the system, given by the terms in equation (2.3.2). If H" is larger than the
characteristic time of a vibrational mode, the system will become unstable and the result
obtained from the simulation will be unphysical. Many of these vibrational modes with high
frequency come from the intramolecular potentials that are used to maintain the bond lengths and
the bond angles of molecules. One way to improve the efficiency is to remove some of these
vibrational modes with high frequencies by using holonomic constraints to maintain the
molecular structure, instead of intramolecular potentials, so that a larger H" can be used [136].
The constraint algorithms used in MD simulation are based on the method of Lagrange
multiplier. A holonomic constraint #(") in MD simulation has the general form
#(") = qG
Z
(") − G
`
(")q
− r
Z`
= 0 (2.4.1)
For system with constraints, the equation of motion for a constrained coordinate
Z
is given by
Y
Z
:
B
:;
B
Z
(") = −
:
:s
W
tuSG
Z
(")T + ∑ v
>
>w
#
>
(")x (2.4.2)
where v
>
is the Lagrange multiplier for the specific constraint #
>
("). In practice, the
unconstrained coordinates of particles are first obtained in each MD step. The goal of the
constraint algorithm is then to obtain the proper v
>
from the unconstrained coordinates, and use
them to find the constrained coordinates that satisfy the equation of motion (2.4.2). This typically
leads to a system of nonlinear equations with unknowns. Different constraint algorithm
29
solves this system of nonlinear equations differently. For example, the SHAKE algorithm [111]
was one of the first constraint algorithm proposed, and it solves the nonlinear equations using
iterative methods such as Newton’s method or Gauss-Seidel method. The LINCS algorithm [49]
formulates the nonlinear equations in the matrix form and approximates the inverse of the matrix
using series expansion. Parallel version of LINCS, known as P-LINCS [51], was developed to
improve the efficiency. For 3-sites molecules (such as water), the SETTLE algorithm [87]
applies the analytical solution for the nonlinear equations to achieve fast calculation of the
constrained coordinates. Since large portion of our systems is composed of water, we typically
use the LINCS algorithm for the constraints of lipids and the SETTLE algorithm for the
constraints of water to optimize the efficiency.
30
2.5 Thermodynamics ensembles and coupling algorithms
By default, MD simulation without any coupling algorithm is in the microcanonical
ensemble. Microcanonical ensemble is a good representation of an isolated system, since the
fixed thermodynamics quantities are the particle number, volume and internal energy (sometimes
the microcanonical ensemble is denoted as NVE). However, system in biological experiment is
usually not isolated. It exchanges heat with the environment and it is influenced by the external
pressure. In order to perform MD simulation that resembles the experiment, coupling algorithms
are introduced to incorporate the thermal effects from the environment.
2.5.1 Temperature coupling
From Maxwell-Boltzmann statistics, the temperature F of a system is related to the
particles’ velocities via the relation:
∑ Y
Z Z
GI Z
=
y
z
D
E
F (2.5.1)
where D
E
is the Boltzmann’s constant and {
+
is the internal degrees of freedom in the system.
The relation (2.5.1) hints that the system’s temperature can be adjusted by scaling the velocities
of the particles. One of the simplest schemes for temperature coupling is the Berendsen
thermostat [10]. In the Berendsen thermostat scheme, temperature of a system is adjusted with
the rate
+@
+;
=
@
|
%@
}
~
(2.5.2)
31
by scaling the velocities with the factor
v = 1 +
~
P;
}
~
/
@
|
@
− 14€
b
B
(2.5.3)
where F
is the target temperature at thermal equilibrium,
@
is the number of time steps to
update the temperature, and �
@
is a parameter that governs the rate of thermalization. Although
the Berendsen thermostat provides an efficient method to adjust the temperature of the system, it
does not produce the canonical ensemble (NVT) correctly.
To resolve this issue, S. Nosé proposed a Hamiltonian with additional variable to
describe the heat bath [95]. The velocities and the time step in Nosé’s proposal are scaled by the
additional variable. W. G. Hoover later modified this approach by introducing a friction term in
the equations of motion to eliminate the scaling of the time step [57]. This temperature coupling
scheme is known as the Nosé-Hoover thermostat.
Another approach to produce the NVT ensemble is the velocity-rescaling temperature
coupling developed by G. Bussi, D. Donadio and M. Parrinello [19]. The scheme is similar to the
Berendsen thermostat, but instead of using a fixed parameter �
@
to control the rate of
thermalization, a stochastic term is included in the calculation of the scaling factor v, in order to
create the correct kinetic energy distribution in the NVT ensemble. Unlike the Berendsen
thermostat and the Nosé-Hoover thermostat, the velocity-rescaling thermostat does not require
any parameter to control the rate of thermalization. Instead, the only parameter it requires is a
seed to generate random numbers for the stochastic term.
Although the Nosé-Hoover thermostat and the velocity-rescaling thermostat yield better
results in the NVT ensemble, the Berendsen thermostat is more efficient due to its simplicity. A
32
mixed of these algorithms is often used to improve the efficiency and accuracy of the simulation.
For example, if the initial configuration of a system is far from thermal equilibrium, Berendsen
thermostat can first be used to thermalize the system quickly near the equilibrium of the desired
temperature, then Nosé-Hoover or velocity-rescaling thermostat can be used to collect the
statistics after the system evolves to the thermal equilibrium.
2.5.2 Pressure coupling
Pressure coupling is based on the rescaling of the particles’ positions (and the simulation
box size). The idea of pressure coupling came from the pressure virial theorem:
‚u = {DF +
M
∑ G
Z`
∙ „
Z`
y
Z,`…Z
(2.5.4)
where ‚ is the pressure, u is the volume, { is the number of the particles in the system. G
Z`
is the
distance between the †
th
and ‡
th
particle, and „
Z`
is the net force between the particles. The term
∑ G
Z`
∙ „
Z`
y
Z,`…Z
is known as the internal virial of the system. Similar to the Berendsen thermostat,
the Berendsen barostat [10] adjusts the pressure of the system with rate
+ˆ
+;
=
ˆ
|
%ˆ
}
‰
(2.5.5)
where the parameter �
Š
takes a similar role as �
@
in the Berendsen thermostat. A major
difference between the temperature and the pressure coupling is that the pressure coupling can be
applied anisotropically to the system to create ensemble with fix membrane tension (sometimes
denoted as NPγT). In the Berendsen barostat, the positions of the particles are scaled by the
tensor:
33
‹
Z`
= H
Z`
−
‰
P;
M}
‰
Œ
Z`
t‚
Z`%
‚
Z`
x (2.5.6)
where H
Z`
is the Kronecker delta,
Š
is the number of time steps between pressure coupling, and
Œ
Z`
is the isothermal compressibility of the system.
Similar to the Berendsen thermostat, the pressure of a system can reach the equilibrium
value using the Berendsen barostat, but the barastat does not yield the correct NPT ensemble.
Additional barostats such as the Parrinello-Rahman barostat [100] and MTTK barostat [78, 132]
were developed to resolve this issue and they are also implemented in the GROMACS package.
34
2.6 Long-range interactions
The size of the system is one of the primary factors that determine the computing cost of
MD simulation. Evaluation of the forces is the most expensive operation in MD simulation due
to the poor scaling with the system size. For a system with { particles, the direct summation of
the electrostatics and Lennard-Jones interactions scales as N({
). With tens of thousands of
particles in the system, direct summation of the forces becomes impractical for most of the
studies.
The simplest treatment to this issue is to truncate the potential function in order to reduce
the number of interactions between particles. In the cutoff scheme, a cutoff radius
�
is defined,
and the modified potential function () is given by the following:
() = Ž
(), ≤
�
0, >
�
‘ (2.6.1)
where
() is the unmodified potential function. Notice in the equation (2.6.1), a discontinuity
is introduced by the scheme at =
�
(See Figure 2.3). To remove this discontinuity, the shift
scheme was proposed. The potential function () in the shift scheme is given by the following:
() = ’
() −
(
�
) − ( −
�
)
‘ :
:=
()“
=
”
, ≤
�
0, >
�
‘ (2.6.2)
The shifted potential function (2.6.2) and its derivative is continuous at =
�
, but the potential
function itself is deviated from the original potential function
() for <
�
(Figure 2.3). To
improve on this issue, the switch scheme utilizes another parameter
and a switch function C()
to control the location of the shifting (see Figure 2.3):
35
() = –
(), ≤
() + C(),
0, >
�
‘
< ≤
�
(2.6.3)
where the switch function C() satisfies the following boundary conditions:
C(
) = 0
‘ :
:=
C()“
=
b
= 0
C(
�
) = −u
(
�
)
‘ :
:=
C()“
=
”
= −
‘ :
:=
u
()“
=
”
(2.6.4)
Figure 2.3 Illustration of the modified Lennard-Jones potential functions under different
truncation schemes. The potential functions are plotted in arbitrary unit with — = 1 and ˜ = 1,
with the cutoffs are set as
= 1.3 and
�
= 2.
Since the Lennard-Jones interaction weakens rapidly with increasing, the above
schemes are sufficient to approximate the potential function within reasonable error. Coulomb
interaction, however, weakens much slower compares to the Lennard-Jones interaction, and the
36
error from using the simple truncation schemes can introduce artifacts in the simulation due to
the lack of long-range interactions [4]. A more sophisticated treatment for the Coulomb
interaction is required.
One of the popular treatments for the long-range Coulomb interaction is the Particle
Mesh Ewald (PME) method [26, 32]. In the PME method, short-range Coulomb interactions
( ≤
�
) are represented by the direct sum in real space just as the cutoff scheme, and the long-
range Coulomb interactions are represented by the Ewald sum (which converges quickly in the
reciprocal space). Utilizing Fast Fourier transform on the Ewald sum, the complexity for
evaluating Coulomb interactions can be reduced from N({
) to N({ ∙ ›œ{). Since the method
uses Fourier transform to change the representation of charge configuration from real space to
reciprocal space, periodic boundary condition is imposed in the simulation by using the PME
method.
37
2.7 Previous MD studies on electroporation
The earliest observation of electropore in MD simulation was reported in 2003.
Electropore formation was observed on a DOPC bilayer under external electric field of 500
MV/m [127]. A more detail MD study on the molecular basis was later presented by D. P.
Tieleman in 2004 [128]. In the study, two systems were examined: DOPC bilayer and a layer of
octane. Both systems were placed under various external field strengths and electropore
formation was observed in both systems when the external field is sufficiently large. The
processes of pore formation in both systems are identical: The formation of electropore begins
with water leaking into the bilayer/octane region and forms a water bridge. The water bridge
continues to expand and the expansion opens up a pore on the membrane.
In addition to D. P. Tieleman’s work, MD study of electroporation was also developed
independently by M. Tarek on DMPC and POPC bilayer [123]. The pore formation process
observed in M. Tarek’s study was similar as those reported previously by D. P. Tieleman. Beside
the electropore formation, M. Tarek also reported pore resealing and the reconstruction of the
lipid bilayer in nanosecond scale after the removal of the external field, which demonstrated that
electroporation is a reversible process for pores with nanometer dimension. A system that
contain peptide nanotube and a system that contains peripheral DNA double strand were also
examined in the study, and an electropore was formed on the lipid bilayer region in all of these
systems.
In later studies, MD simulation was used to investigate the molecular composition of the
lipid bilayers and the effects on electroporation. For example, simulations of PLPC and oxidized
PLPC bilayer revealed that bilayer with oxidized lipids are more permeable for water, and
38
electropore formation is more prone to occur on the oxidized region [146]. Simulations on the
lipid bilayer that contains cholesterols demonstrated that the presences of cholesterols increases
the cohesive force on the membrane, and inhibits the electropore formation process [34].
Electroporation on more complex structures, such as bacterial membranes, were also investigated
using MD simulation [104].
Pore formation on lipid bilayer can be mediated by other mechanisms without the
influence of an external field, and these processes can be observed with MD simulation. For
example, formation of a transient pore can be mediated by ionic charge imbalance, and such
transient pore allows leakage of small ions through the cell membrane [43]. Transmembrane pore
can also be created by applying surface tension on the membrane, and the pore can be stabilized
and controlled by tuning the applied surface tension [74]. Simulation on the surface tension-
stabilized pore allows study of diffusive transport of small molecules (such as water and
monovalent ions) through the transmembrane pore [75].
In the remaining of this section, I shall discuss in detail, some of the previous MD studies
that are closely related to my research.
2.7.1 Dynamics and life cycle of an electropore
In order to perform quantitative study on the electroporation dynamics, Z. A. Levine and
P. T. Vernier proposed a scheme to characterize the stages of an electropore life cycle based on
the events observed in MD simulation [76, 77]. An illustration of the scheme is given in Figure
2.4.
39
Figure 2.4 Illustration of the life cycle of an electropore. In the snapshots, oxygen and hydrogen
of water molecules are represented by red and white dots. Phosphorus atoms of the
phospholipids are represented by gold spheres to indicate the approximate location of the lipid
head groups. The figure was obtained from ref [77].
An electropore’s life cycle can be divided into pore creation and pore annihilation. Pore
creation begins with the pore initiation stage. Pore initiation starts with the application of the
external field. Under the influence of external field, water protrudes into the bilayer region from
both side of the bilayer. The protrusions from top and bottom water group eventually connect
and form a water bridge. At this point, the system enters the pore construction stage.
40
During the pore construction stage, water molecules from the bulk water phases continue
to enter the bilayer region under the influence of the external field, and expand the water bridge.
Meanwhile, lipids rearrange such that the head groups follow the protruding water into the
bilayer center. The movements of lipid rearrangement can be monitored by the positions of the
phosphorus atoms and the nitrogen atoms. Eventually, top and bottom phosphorus groups
connect to each other and the pore wall is formed. At this point, the pore construction is
complete and the system enters the pore expansion stage.
The size of an electropore is determined by the electric field strength. In order to
observed electropore creation within reasonable simulation time (restricted by the computing
resources), the porating field strength is usually a high value (> 350 MV/m). Under such external
field, the electropore will continue to expand after the pore construction until finite size effect
appears and ruptures the lipid bilayer. If the external field is reduced after the pore reaches
certain radius, then the electropore can be stabilized. Discussion on field-stabilized electropore
simulation will be provided in later section.
The dynamics of pore annihilation is analogous to the reverses of the pore creation. Pore
annihilation starts with the removal of the external field. In the absence of the field, lipid head
groups of the pore wall maneuver back to the bilayer and cause the pore to shrink. Eventually
only a single pair of phosphorus connects the top and bottom phosphorus group, and at this point
the pore radius is reduced to zero. The pore wall deconstructs as the last pair of phosphorus
disconnects and the lipids head groups maneuver back to the bilayer. When the lipids reorient
with their tails facing the bilayer center, the once pore region is again hydrophobic.
Subsequently, the water bridge collapse as the water molecules retreat to the water phases, and
41
once all of the water molecules reside in the water phases, the annihilation process is complete
and the system is back to the bilayer configuration.
The elapsed time of each stage in the electropore life cycle was examined by Z. A.
Levine and P. T. Vernier [76]. External fields from 300 to 600 MV/m were used for pore
creation (denoted by
Š
for porating field). Pore initiation is the only stage in the pore creation
process that exhibits dependence of applied field, with its elapsed time decrease rapidly with
Š
(i.e., �
ZZ;
= !"# " ∙ �
%ž(A
‰
)
, with (
Š
) is an increasing function of
Š
). Elapsed times of
other pore creation stages exhibit no dependence on
Š
. For lower fields that were examined (
Š
≤ 400 MV/m), large portion of the pore creation time are occupied by pore initiation. Near the
minimum porating field (
Š
= 300 MV/m), pore initiation time is in general above 10 ns, and
time scales in other pore creation stages are only in order of 0.1 ns.
The value of porating field
Š
imposes no effect on the pore annihilation stages.
However, the annihilation time is determined by the size of the pore at the moment when the
external field is removed. Collapse of the water bridge (last stage of pore annihilation) occupies
the smallest portion of the pore annihilation process, and the average time for the water bridge to
collapse is less than 2 ns, which implies the water bridge is very unstable in the hydrophobic
region of the lipid bilayer interior.
2.7.2 Water’s role in electroporation
One of the objectives of our research is to identify the underlying force responsible for
electroporation. The comparison between DOPC bilayer and octane layer in Tieleman’s study
42
provided some interesting insights: Because octane has no polar head group, the external field
has no effect on its orientation. Later MD studies also demonstrated that the external field
imposes very little effect on the orientations of the lipid head groups before electroporation
occurs [145, 156]. These evidences hint that the orientations of the lipids play little role in the
pore formation, and additional MD simulations suggested that the pore formation process is
promoted by the water defect, generated by interactions between the external field and the water
dipoles at the water/lipid or water/octane interface [156]. The effect of salt was also tested in
Tieleman’s study, and the presence of salt does not change the events in pore formation process
[127].
Protrusion of water into the lipid bilayer region and the formation of water bridge are the
earliest events observed in the pore creation process. To elucidate the role of water in
electroporation, we examined the energetic of the water dynamics during electroporation using
MD simulation [129]. In most of the MD simulations, the correlation coefficient between water
protrusion height and the electrostatic energy is negative, which means the system experiences
reduction in electrostatic energy during water protrusion. In contrast, the Lennard-Jones energy
is relatively steady during the same period. The decrease in electrostatic energy during water
protrusion can be understood with basic electrostatic: Water is a well-known polar molecule, and
the dipole-dipole energy between two water dipoles with distance is given by:
+ZŠ
jk%+ZŠ
jk
=
&\0
|
=
Ÿ
‹ ¡
∙ ‹ ¡
− 3(‹ ¡
∙ ̂
)(‹ ¡
∙ ̂
)£ −
¤ ¡
∙ (‹ ¡
+ ‹ ¡
) (2.7.1)
where ‹ ¡
and ‹ ¡
are the dipole vectors of two water molecules. In the absence of an external
field, there is no ordering of the water dipole moments. Under a strong external field, the water
dipole moments are forced to point along the direction of the field. Consider two configurations
43
of the water dipole vectors: the first represents the water dipoles at the flat surface (Figure 2.5
A). The energy of this configuration is given by:
a
=
¥
B
&\0
|
=
Ÿ
− 2‹ (2.7.2)
The second configuration represents the water dipoles of the water protrusion (Figure 2.5 B), and
its energy is given by:
E
= −
¥
B
\0
|
=
Ÿ
− 2‹ (2.7.3)
Comparing the two configurations, the configuration that resembles the water protrusion (Figure
2.5B) is energetically favorable. The energetic study of the water dynamics also suggested the
dipole-dipole interactions between water molecules are responsible for the construction of water
bridge during the electropore formation process [96, 129].
Figure 2.5 Dipole configuration of two water molecules resembles A) flat water surface, and B)
water protrusion.
2.7.3 MD simulation of field-stabilized electropore
A protocol to simulate field-stabilized electropore was recently developed by Fernández
et al [35]. The protocol is a two-steps process that was briefly used earlier by Böckmann et al to
44
examine the steady state of an electropore [16]. The two-steps protocol is the following: From
the bilayer configuration, a high external field (denoted by
Š
for porating field) is used to create
an electropore in the system. After the pore is formed and reaches certain radius, the external
field is reduced to a lower value (denoted by
l
for sustaining field) to sustain the electropore in
a steady state. Typically,
Š
is larger than 350 MV/m and
l
is around 100 MV/m for our 128
lipids bilayer system. The reason for using this two-steps process is because of the computing
cost of MD simulation. Theoretically, an electropore can be created under the sustaining field
value, but the characteristic time of the electropore formation could be larger than a millisecond,
which is very expensive for the MD simulation in our current resources. The use of a high
porating field at the beginning is intended to stimulate the system over the energy barrier and
speeds up the pore formation process. On the other hand, because the equilibrium pore radius
using the high porating field is larger than the simulation box size, if the external field is retained
in the porating field value after the electropore is formed, the electropore will continue to expand
until finite size defect appears and tears the membrane apart.
Fernández et al. performed detail study on the two-steps protocol for field-stabilized
electropore simulation. The electropore size in the steady stage can be controlled by the
sustaining field
l
, and it does not depend on the porating field strength
Š
nor the configuration
of the system when the field is reduced, as long as the electropore has formed and reaches certain
radius (~ 1 nm). Evolution of pore radius depends on the electropore’s configuration at the
moment when the field is reduced. For the 128 lipids bilayer system, the electropore typically
evolves to the steady state within 20 ns after the field reduction, and the steady state of
electropore can be stabilized using this two-steps protocol for over hundreds of nanoseconds.
45
One of the important properties of a transmembrane pore is its permeability. Diffusive
transport of molecules through a transmembrane pore can be studied with the MD simulation of
a tension-stabilized pore [75]. Similarly, the electrophoretic transport of charged molecules
through a transmembrane pore can be studied with the MD simulation of a field-stabilized
electropore, using the protocol described above. The sustaining field in the conductance
simulation plays both roles of sustaining the pore and driving charged molecules through the
pore. A detail study on the pore conductance of monovalent ions will be presented in Chapter 4.
46
Chapter 3: Dynamics of water bridge in
electroporation
Portions of this chapter are adapted from:
M.C. Ho, Z.A. Levine, and P.T. Vernier, Nanoscale, Electric Field-Driven Water Bridges in
Vacuum Gaps and Lipid Bilayers. J Membr Biol, 2013.
3.1 Abstract
Formation of a water bridge across the lipid bilayer is the first stage of pore formation in
molecular dynamics (MD) simulations of electroporation, suggesting that the intrusion of
individual water molecules into the membrane interior is the initiation event in the sequence that
leads to the formation of a conductive membrane pore. To delineate more clearly the role of
water in membrane permeabilization, we conducted extensive MD simulations of water bridge
formation, stabilization, and collapse in palmitoyloleoylphosphatidylcholine (POPC) bilayers
and in water-vacuum-water (W-V-W) systems, in which two groups of water molecules are
separated by a 2.8 nm vacuum gap, a simple analog of a phospholipid bilayer. Certain features,
such as the exponential decrease in water bridge initiation time with increased external electric
field, are similar in both systems. Other features, such as the relationship between water bridge
lifetime and the diameter of the water bridge, are quite different between the two systems. Data
like these contribute to a better and more quantitative understanding of the relative roles of water
and lipid in membrane electropore creation and annihilation, facilitating a mechanism-driven
47
development of electroporation protocols. These methods can be extended to more complex,
heterogeneous systems that include membrane proteins and intracellular and extracellular
membrane attachments, leading to more accurate models of living cells in electric fields.
48
3.2 Introduction
The cell membrane separates the cell interior from the external environment, enabling
control of substances traveling into and out of the cell. When a sufficiently high electric field is
applied, the cell membrane becomes permeable to charged species and small molecules that
normally cannot cross the membrane barrier, a phenomenon known as electroporation [125,
157]. This nonlethal method for modifying membrane permeability permits the introduction of
foreign materials into living cells, enabling applications like electro-chemotherapy [85] and
electroporation-mediated gene transfer [48]. Electroporation has been studied for decades
experimentally [93] and theoretically [72, 151], but a precise characterization of the structures
responsible for the increased permeability in living cells has not yet been achieved.
Identifying the constituents of the permeabilized membrane is challenging because of the
small dimensions (nanometers) of the putative electropores and the short time scale
(nanoseconds) of pore creation, which makes direct observation of electroporation inaccessible
to current experimental methods. Given this difficulty, molecular dynamics (MD) simulations
have been employed in recent years for the study of electroporation at the molecular level [123,
128], and it has revealed microscopic aspects of the electroporation of lipid bilayers that are
consistent with experiment [16, 34, 77, 104, 145].
Among the mechanistic molecular details of electroporation that remain to be clarified is
the initial step in the reorganization of water and lipid that leads to the formation of a pore in a
lipid bilayer. Previous MD studies have suggested that the intrusion of water molecules into the
hydrophobic, phospholipid tail region of the bilayer is the earliest event that can be identified
49
[128, 145]. After interfacial water molecules enter into and then bridge the bilayer interior, lipids
follow and rearrange to form a hydrophilic wall around the water column.
In an important analytical study, Okuno et al. demonstrated that under the influence of
external electric field, water molecules rearrange themselves (due to their dipole moments) to
form an energy-minimized quasi-conical structure [96]. We observe similar structures in the pore
initiation stage in MD simulations of electroporation [148], and we speculate that this energy-
minimizing configuration of interfacial water in an electric field normal to the interface is the
key initiating structure in the formation of an electropore. In this work, we examine the dynamics
of water bridges in different stages of electroporation in two systems: water-lipid-water (a POPC
bilayer sandwiched by two layers of water), and water-vacuum-water (W-V-W, a vacuum gap
between two layers of water). Although the water-vacuum-water system is somewhat artificial, it
is stable on these time scales, and it provides a simple view of the dynamics of interfacial water
under the influence of an applied electric field.
50
3.3 Methods
3.3.1 Molecular dynamics simulations
All simulations were performed using GROMACS version 4.0.5 [50] on the University of
Southern California High Performance Computing and Communications (HPCC) Linux cluster
(http://www.usc.edu/hpcc/). Lipid topologies were taken from OPLS united-atom parameters
[13] and the Simple Point Charge (SPC) model [9] was used for water. POPC systems were
simulated with constant particle number, pressure, temperature (NPT) and constant particle
number, volume, temperature (NVT) ensembles, as indicated. NVT ensemble was used W-V-W
systems. Note that W-V-W systems cannot be sustained in an NPT configuration since the
pressure coupling would change the volume of the simulation box, push the water layers
together, and collapse the vacuum gap. Temperature was held at 310 K with an external heat bath
using a velocity rescaling algorithm [19]. For the NPT ensemble, pressure was held at 1 bar
using a weak coupling algorithm [10], with a relaxation time of 1 ps and compressibility of 4.5 ×
10
-5
bar
-1
semi-isotropically applied in both normal and in-plane directions relative to the
membrane. Bond lengths were constrained using the LINCS algorithm [49] for lipids and the
SETTLE algorithm [87] for water. Short-range electrostatic and Lennard-Jones interactions were
cut off at 1.0 nm. Long-range electrostatic interactions were calculated with PME algorithm [32]
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 in all directions were employed to mitigate system size effects.
51
3.3.2 Structure
POPC systems contain a lipid bilayer composed of 128 POPC lipids (64 per leaflet) and
approximately 4480 water molecules (35 waters/lipid), with initial system dimensions of
approximately 7 nm × 7 nm × 7 nm. W-V-W systems were created by using a custom Perl script
to remove waters from the central region of a volume of bulk water to produce the desired gap
sizes. Systems were equilibrated for 300 ps so that both kinetic and potential energies reached a
steady value. The planes of the water-lipid interface and the water-vacuum interface are normal
to the z direction of the simulation box. Molecular graphics images were generated with Visual
Molecular Dynamics (VMD) [58].
3.3.3 Water bridge initiation time, annihilation time
Initiation and annihilation times for an electropore in a lipid bilayer system were defined
previously [76], and we establish similar definitions here for water bridge initiation and
annihilation times. For both water-lipid-water and water-vacuum-water systems the water bridge
initiation time is defined as the time for the top and bottom water groups to first merge (without
subsequently separating) after the application of the external electric field. The water bridge
annihilation time is defined as the time for the water to separates into top and bottom groups
(without subsequently reconnecting) after the external field is removed.
52
3.3.4 Measurement of water bridge radius
A Perl script is used to measure the radius of the water bridge. The script first obtains the
water density profile in z and defines the water bridge boundaries as the planes where the water
density is 50% of the bulk value. Then the script subdivides the water bridge region into 10 bins
(See Figure 3.1). Within each bin, the script finds the largest and smallest x coordinate of all the
water molecules within the bin, x
max
and x
min
, and the same is done for the y coordinates. For
each bin, r
x
= (x
max
+ x
min
)/2, r
y
= (y
max
+ y
min
)/2, and r
bin
= (r
x
+ r
y
)/2. The water bridge radius for
a particular time frame is defined as the average of the five smallest values of r
bin
during that
time window.
Figure 3.1 Illustration of water bridge radius measurement. White dots are hydrogen atoms and
red dots are oxygen atoms of the water molecules. Yellow lines are the boundaries of the water
bridge region defined by the radius measurement script. Blue boxes denote the bins defined by
the script and the corresponding x
max
and x
min
values.
53
3.4 Results and Discussion
3.4.1 Internal electric field versus external electric field
We use the magnitude of the electric field in the low-permittivity interior (hydrocarbon
chains or vacuum) as a normalizing reference for comparisons between water-lipid-water and
water-vacuum-water systems. For a POPC system, the distance between the top and bottom
water groups is determined by the properties of the lipid bilayer and does not change more than
15% over the range of field magnitudes used in this work [156], so the internal electric field for a
given POPC system is a function of the applied (external) electric field.
For a W-V-W system the situation is not so simple. The surface charge of the water
molecules at the vacuum-water interface contributes to the net internal electric field, so that for a
given external field, the internal field depends on the vacuum gap thickness, which is defined
when the system is assembled. We choose a separation that result in the same internal electric
field for a given external field in both the POPC and W-V-W systems.
Figure 3.2 shows the linear relationship between external and internal electric fields for
POPC systems and W-V-W systems with different vacuum gap thicknesses. The internal field is
defined as the local field at the center of the system (between the top and bottom water groups),
and it is calculated using the GROMACS function g_potential. As mentioned above, for W-V-W
systems the internal field increases when the vacuum gap decreases. The W-V-W system with a
2.8 nm vacuum gap exhibits an external versus internal field relationship that is very close to that
of a POPC system.
54
Figure 3.2 Internal versus external electric field in W-V-W and POPC system. The internal field
is measured at the midpoint between the upper and lower water groups. Each data point is the
average of 3 trials. Error bars are given as the standard deviation.
3.4.2 Water bridge formation
The processes of water bridge formation are similar in POPC and W-V-W systems
(Figure 3.3). When an external electric field is applied, water dipoles rotate in the direction of the
field, and, after an interval randomly distributed around a time that is exponentially dependent on
the field strength (Figure 3.4), an assembly of water molecules begins to intrude into the
hydrocarbon interior of the lipid bilayer or the vacuum gap of a W-V-W system. For W-V-W
systems the initiation time is also a function of the vacuum gap size. In a very short period of
time the chain of water extends from one interface to the other, forming a bridge. If the field
strength is above a certain value (typically above 150 MV/m in our systems), more and more
water enters the bridging structure, and the radius of the water bridge increases until it reaches
the system boundaries and finite size defects appear due to periodic boundary conditions.
55
Figure 3.3 Water bridge formation in POPC (top and middle row) and W-V-W (bottom row)
systems. Water hydrogen and oxygen atoms are white spheres and red spheres, respectively. In
the top row, nitrogen atoms and phosphorus atoms of POPC head groups are represented as blue
spheres and gold spheres, and lipid tails are represented as light blue lines. The frames in the
middle row are identical to those in the top row, except that only water molecules are displayed.
Figure 3.4 Water bridge initiation time versus external electric field strength. Each data point is
the average of 3 simulations. Error bars are the standard deviation.
56
We examined the initiation time of W-V-W system with two different thicknesses: W-V-
W system with 2.8 nm gap that resemble similar internal vs. external field relationship in POPC
system, and W-V-W system with 4.0 nm gap resemble similar water volume which has
approximately the same distance between top and bottom water groups as in POPC system. The
initiation time in W-V-W system with 2.8 nm gap is much closer to POPC system compares to
W-V-W system with 4.0 nm gap.
No structure that follows the application of the external field and precedes the initial
water intrusion that leads to bridge formation has been identified. Extension of the chain of water
molecules into the lipid/vacuum region is very fast once it begins, typically less than 500 ps for
electric fields equal to or greater than 350 MV/m, while the total initiation time for both lipid and
vacuum systems can be more than 100 ns using a weaker field (< 300 MV/m). Comparing the
two systems, water intrusion time in W-V-W system is much less than POPC system
(approximately 10 ps versus approximately 100 ps), suggesting that the fatty acid chains of the
lipid bilayer behave as a mechanical barrier that slows down the construction of the water bridge.
Water bridge formation can be observed by monitoring the water dipole moment (Figure
3.5). In the absence of an external electric field, water molecules are ordered locally by hydrogen
bonding, and the system has a zero net dipole moment. Application of an external field orients
the water dipoles, producing a finite net dipole moment in the system. The system dipole
moment increases again when the water bridge is formed (indicated by the red lines in Figure 5).
The net dipole moment increases rapidly during the expansion of the water bridge and reaches a
new equilibrium value if the system is stabilized.
57
Figure 3.5: Total system dipole moment of water in z direction versus time during water bridge
formation from selected simulations in W-V-W (upper panel) and POPC (lower panel) systems.
The external field is applied at t = 0. Top and bottom water groups connect at the time indicated
by the red lines.
Recall from basic electrostatics that when two dipoles with are placed parallel to each
other with the same orientation, their dipole-dipole interaction energy is positive. In our
simulations, the water dipoles are forced to align in the z direction by the external field. In order
to minimize the potential energy, the system favors the extension of water intrusions in the z
direction and the lengthening of the growing water bridge, which increases the surface area
normal to the x-y direction and reduces the number of parallel dipole neighbors.
58
3.4.3 Water bridge stabilization
To simplify the discussion, we consider here only the water bridge and subsequent water
column that forms in either POPC or W-V-W systems. A sufficiently large electric field (> 300
MV/m) is required to create a water bridge within a reasonable simulation time (< 25 ns) [156].
If this field continues to be applied after the water bridge forms, the water column expands
radially until anomalies arising from finite size effects appear.
If, however, the external field is reduced after the water bridge forms, the water bridge
can be stabilized with a constant radius, similar to what has been reported previously for
phospholipid bilayers [35]. In our simulations, the water bridge can be stabilized with sustaining
fields ranging from 50 MV/m to 100 MV/m in both POPC and W-V-W systems. Within this
range, the stabilized, time-averaged water bridge radius increases linearly with the sustaining
field similarly for both POPC and W-V-W systems, but without the constraining presence of the
phospholipid pore wall, the water bridge radius in W-V-W systems is larger for a given
sustaining field than the corresponding structure in POPC systems (Figure 3.6).
Figure 3.6 Time-averaged water bridge radius versus sustaining field. Each data point is the
average of 3 simulations. Error bars are the standard deviation.
59
POPC and W-V-W systems differ in the time required for a system to evolve to a new
equilibrium after the external field is reduced. Figure 3.7 shows the evolution of the water bridge
radius under various sustaining fields. The water bridge radius in W-V-W systems reaches an
equilibrium value in a few hundred picoseconds or less after the external field is reduced. POPC
systems take several nanoseconds to stabilize. This difference in the kinetics of water bridge
radius evolution is probably a consequence of the momentum and interactions of the lipids in the
pore wall of the POPC systems, which hinders water bridge expansion and contraction.
Figure 3.7 Evolution of water bridge radius during stabilization from selected simulations. At
" = 0, the external field was reduced to the sustaining field value labeled in the legend.
60
3.4.4 Water bridge annihilation
When the external electric field is removed, water molecules in the system are no longer
preferentially oriented in the z-direction, and the water bridge becomes an energetically
unfavorable configuration. Hydrogen bond associations drive the migration of water molecules
back into the bulk, and the bridge is annihilated. In POPC systems, the deconstruction of the
water bridge involves also the hydrophilic interactions between water and the phospholipid head
groups in the pore wall, as the pore lipids move back into the plane of the bilayer.
Figure 3.8 shows how the water bridge annihilation time (defined as the time between the
field removed and the disconnection of the top and bottom water groups) varies with the radius
of the water column when the field is removed. In W-V-W systems the water bridge collapses
very quickly (< 0.3 ns) when the initial bridge radius is less than 1.2 nm. Water bridge with
radius larger than 1.2 nm are stable for about 10 ns before annihilation, and the annihilation time
is independent of bridge radius for values up to 2.0 nm.
In contrast with W-V-W systems, annihilation time in POPC systems is dependent on
bridge radius. Annihilation of a phospholipid membrane pore involves the rearrangement of both
water and lipids. Since larger pores incorporate more lipids which must move back into the
bilayer, it is not surprising that water bridge annihilation times are much longer in POPC systems
than they are in W-V-W systems. We can conclude that the dynamics of electropore annihilation
in phospholipid bilayers are determined primarily by the properties of the lipids, and that the
lifetime of the water bridge in a POPC system is extended by the hydrophilic interaction between
the bridge water and the pore lipids.
61
Figure 3.8 Water bridge annihilation time versus initial water bridge radius in POPC system (top
panel) and W-V-W system (bottom panel). Each data point is the average of at least three
simulations. Error bars are given by the standard deviation.
3.4.5 Effects of water models in water bridge formation simulation
In most of our MD simulations, water molecules are represented by the SPC water model.
SPC water model is a simple and efficient model that is commonly used for simulation of
biological systems. However, many of its properties deviate from real water molecules due to its
62
simplicity. To examine the effects of water model on electroporation simulation, we simulated
water bridge formation in both POPC and W-V-W system using SPC, SPC/E, and TIP3P model
for water. These are the popular 3-site water models, and the slight differences in their structures
and parameterizations lead to different water properties [79].
The events in water bridge formation process are identical for all three water models, as
described in Section 3.4.2. However, the time scale of water bridge initiation using different
water models can be very different. The water bridge initiation times for simulations using
different water models are presented in Table 3.1.
Table 3.1 Water bridge initiation time using SPC, SPC/E, TIP3P model.
Water bridge initiation time in POPC bilayer (ns)
External field (MV/m) SPC SPC/E TIP3P
300 --- --- 6.49 ± 4.89
350 --- --- 1.84 ± 0.92
450 7.31 ± 4.62 55.18 ± 45.09 0.85 ± 0.23
500 --- --- 0.63 ± 0.33
550 2.03 ± 0.90 7.16 ± 3.29 ---
650 1.12 ± 0.21 2.23 ± 0.90 ---
750 0.68 ± 0.19 0.89 ± 0.30 ---
Water bridge initiation time in W-V-W system (ns)
External field (MV/m) SPC SPC/E TIP3P
300 --- --- 4.60 ± 3.90
350 7.88 ± 6.71 --- 0.66 ± 0.35
400 1.05 ± 1.31 --- 0.069 ±0.068
450 0.19 ± 0.12 6.63 ± 5.82 0.034 ± 0.014
500 0.042 ± 0.016 0.58 ± 0.25 ---
550 --- 0.29 ± 0.09 ---
600 --- 0.078 ± 0.044 ---
Water bridge initiation time of each field is given by the average of 5 trials. Errors are given by
the standard deviation.
63
Regardless of the water model, the initiation time decrease rapidly with increasing
external field, as described in Section 3.4.2 (see Figure 3.9). For a given external field, water
bridge initiation with TIP3P water is fastest amongst the water models, and water bridge
initiation with SPC/E water is the slowest. This order of initiation times (τ
init,TIP3P
< τ
init,SPC
<
τ
init,SPC/E
) is observed in both POPC bilayer and W-V-W system. At lower external field (< 450
MV/m), the initiation times using different water model can be different by orders of magnitude.
One of the differences in the properties of these water models is their self-diffusion
coefficients. The self-diffusion coefficient is 5.6 × 10
-9
m
2
s
-1
for TIP3P water, 4.2 × 10
-9
m
2
s
-1
for
SPC water, and 2.8 × 10
-9
m
2
s
-1
for SPC/E water [79]. This ordering of self-diffusion coefficients
(D
SPC/E
< D
SPC
< D
TIP3P
) is the opposite of the ordering of the initiation times. From Stoke-
Einstein equation, lower value of diffusion coefficient corresponds to higher viscosity of the
fluid. The ordering of diffusion coefficients and the ordering of initiation times suggest that
water bridge formation is more difficult for aqueous medium that possess high viscosity.
Figure 3.9 Water bridge initiation time using SPC, SPC/E, TIP3P model in a) POPC bilayer, and
b) W-V-W system. Each data point is represented by the average 5 trials, and the error is given
by the standard deviation.
64
3.5 Summary
We performed MD simulations to study the water bridge in W-V-W and POPC systems
during different stages of electropore formation in order to delineate the roles of water and lipids.
The processes of water bridge formation in both systems are similar, and both systems exhibit an
exponential decrease in water bridge initiation time with increasing external electric field.
Once a water bridge has formed, the radius of the bridge in both POPC and W-V-W
systems increases linearly with the magnitude of the sustaining electric field over the range of
values we examined. For a given value of the sustaining field, the radius of the water bridge is
larger in a W-V-W system than in a POPC system. Changes in bridge radius in response to
changes in applied electric field are faster in W-V-W systems than in POPC systems.
Annihilation of the water bridge occurs much faster in a W-V-W system than in a POPC
system. Water bridge annihilation time does not increase with the radius when the radius is
above 1.2 nm in W-V-W systems, but it increases significantly with the radius in POPC systems,
suggesting that hydrophilic interactions between the water and the lipid pore wall and among the
pore lipids themselves contribute to the stabilization of a phospholipid electropore, and that the
properties of the lipids are key determinants of the stability of the electropore.
The difference in water model used in the simulation can influence the time scale of
water bridge formation. Pore initiation time can be varied by orders of magnitude using different
water model, and it follows the order τ
init,TIP3P
< τ
init,SPC
< τ
init,SPC/E
for the water models we
examined. The ordering of the self-diffusion coefficient D
SPC/E
< D
SPC
< D
TIP3P
is opposite of the
65
ordering of pore initiation time, indicates that water bridge formation is more difficult for
aqueous medium that possess high viscosity.
66
Chapter 4: Ion conductance of nanoscale
electropore
Portions of this chapter are adapted from:
M.C. Ho, M. Casciola, Z.A. Levine, and P.T. Vernier, Molecular Dynamics Simulations of Ion
Conductance in Field-Stabilized Nanoscale Lipid Electropores. J Phys Chem B, 2013.
4.1 Abstract
Molecular dynamics (MD) simulations of electrophoretic transport of monovalent ions
through field-stabilized electropores in POPC lipid bilayers permit systematic characterization of
the conductive properties of lipid nanopore. The radius of the electropore can be controlled by
the magnitude of the applied sustaining external electric field, which also drives the transport of
ions through the pore. We examined pore conductances for two monovalent salts, NaCl and KCl,
at physiological concentrations. Na
+
conductance is significantly less than K
+
and Cl
-
conductance and is a non-linear function of pore radius over the range of pore radii investigated.
The single pore electrical conductance of KCl obtained from MD simulation is comparable to
experimental values measured by chronopotentiometry.
67
4.2 Introduction
Electroporation, a widely used method for introducing foreign materials into cells [81,
125], has been studied experimentally and with analytical and numerical models for decades [38,
92, 151]. A current goal in improving our understanding of electroporation is the development of
a comprehensive microscopic description of the phenomenon, not an easy task due to the
nanoscale dimensions of the lipid electropore and the short time scale (nanoseconds) of pore
creation, which present challenges to direct experimental observations. For these reasons
molecular dynamics (MD) simulations have become a useful tool for the study of electroporation
in atomic detail. In the past decade, MD simulations have revealed many microscopic aspects of
electroporation including the molecular structure of the electropore [74, 123, 128, 145] and the
dynamics and energetic of pore creation and annihilation [16, 76, 156]. Despite this progress,
MD simulation results are rarely correlated directly and quantitatively with experiments.
A recently developed method for stabilizing electropores in MD simulations [35, 104]
enables the characterization of the properties of these permeabilizing structures under controlled,
steady-state conditions. The pores can be sustained for long periods (hundreds of nanoseconds or
longer) by applying an external electric field with appropriate magnitude along the axis of the
pore. The diameter of the pore can be tuned by varying the external field, enabling a systematic
examination of the ion conductance of lipid electropores.
Historically, an increase in electric conductance was the first indicator of the electric
field-driven breakdown of the membrane barrier function that has come to be known as
electroporation [25, 120], and the measurement of conductance is still a useful method for
68
characterizing membrane permeabilization [98]. In this paper we compare conductance values
from chronopotentiometric studies of planar lipid bilayers, which can measure the conductance
of a single electropore [65, 70, 71], to those obtained from MD simulations of field-stabilized
electropores in 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine (POPC) bilayers in
aqueous solutions containing Na
+
, K
+
, and Cl
-
. The maintenance of intracellular and extracellular
concentrations of these ions ([Na
+
] is much higher outside the cell; [K
+
] is much higher inside
the cell) is critical for the electrophysiology of living cells. Disruption of these intracellular-
extracellular ion gradients, normally maintained by ion channels and pumps at considerable
energy expense [5], by electroporation, which allows diffusive transport of ions and small
molecules through the permeabilized membrane [106, 114], has profound consequences for the
cell, which must restore osmotic and ion balance or die.
Ions affect the properties of lipid bilayers, including area per lipid, head group
orientation, and lateral lipid mobility, as shown both experimentally and in simulations [14, 37,
44, 63, 134, 135]. MD models of ion transport through porated lipid bilayers suggest that Na
+
increases electropore line tension as a result of Na
+
-phospholipid binding, resulting in pore
destabilization and a decrease in pore lifetime [75]. Ion transport through porated lipid bilayers
has been studied to a limited extent with MD methods [43, 75], but here we report the first
molecular measurement of ion electrical conductance through electric field-induced pores, and
we compare these results with those obtained from experimental systems.
69
4.3 Methods
4.3.1 Molecular dynamics simulations
All simulations were performed using GROMACS version 4.0.5 [50] on the University of
Southern California High Performance Computing and Communications (HPCC) Linux cluster
(http://www.usc.edu/hpcc/). Lipid topologies were taken from OPLS united-atom parameters
[13], and the Simple Point Charge (SPC) model
[9] was used for water. For Na
+
, K
+
, Cl
-
, we
employed the set of parameters originally developed by Straatsma and Berendsen [119] that are
supplied with the GROMACS force field. All simulations were performed under the NPT
ensemble. Temperature was held at 310 K with an external heat bath using a velocity rescaling
algorithm [19] with a relaxation time of 0.1 ps. Pressure was held at 1 bar using a weak coupling
algorithm [10], with a relaxation time of 1 ps and compressibility of 4.5 × 10
-5
bar
-1
semi-
isotropically applied in both normal and in-plane directions relative to the membrane. Bond
lengths were constrained using the LINCS algorithm
[49] for lipids and the SETTLE algorithm
[87] for water. Short-range electrostatic and Lennard-Jones interactions were cut off at 1.0 nm.
Long-range electrostatic interactions were calculated with a PME algorithm [32] 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 in all directions were imposed to mitigate system size effects.
70
4.3.2 Structure
All simulations contain a bilayer composed of 128 POPC lipids (64 per leaflet) and
approximately 9000 water molecules (~70 waters/lipid), corresponding to an initial system box
size of approximately 7 nm × 7 nm × 10 nm. The GROMACS function ‘genion’ was used to
replace bulk water molecules with Na
+
, K
+
, and Cl
-
. 40 Na
+
and 40 Cl
-
were added to the NaCl
system and 22 K
+
and 22 Cl
-
were added to the KCl system in order to create systems with free
ion concentrations close to physiological values after equilibration. A custom Perl program was
used to check the binding of ions to the lipid bilayer; an ion within 0.3 nm of a POPC glycerol
backbone acyl oxygen is considered bound. We equilibrate the systems until the numbers of free
and bound ions in the system and the area per lipid reach a steady value. Equilibration time is
typically 100 ns for NaCl systems and 50 ns for KCl systems. Molecular graphics images were
generated with Visual Molecular Dynamics (VMD) [58].
4.3.3 Conductance simulations
Utilizing a protocol similar to the one developed by Fernández et al. [35] to produce
electric field-stabilized pores, we apply an external porating electric field (E
p
= 400 MV/m) to a
POPC bilayer. When the phosphorus groups from the top and bottom leaflets merge
[76], and the
pore radius reaches approximately 1 nm, we reduce the external field to a lower, sustaining field
(E
s
= 50, 75, and 100 MV/m; note that these fields correspond to transmembrane potentials of
roughly 200, 300, and 400 mV. Under these conditions the pore is maintained in a steady state. If
ions are present in the system they are electrophoretically driven through the pore, creating an
ionic current. Conductance simulations were run for 100 ns after the reduction of the external
71
field to E
s
. Three independent trials were carried out for each E
s
by assigning every atom a
randomized velocity (a built-in function of GROMACS) at the time the field was reduced.
In addition, we measured the electrophoretic mobility of the ions in bulk water at several
ion concentrations. A simulation box of bulk water with dimension 7 nm x 7 nm x 10 nm is
created using the GROMACS function ‘genbox’ and ions are inserted with GROMACS function
‘genion’. The system is equilibrated for 300 ps, and simulations of 5 ns with an applied external
field are performed to measure the electrophoretic mobility.
4.3.4 Pore radius measurement
A custom Perl script is used to extract the radius of the electropore. The script first shifts
the configuration to center the pore within the simulation box according to the water density
distribution in all 3 dimensions. Then the upper and lower bounds of the pore region are defined
using the average of the 15 largest and smallest z-coordinates of the POPC phosphorus atoms.
Next the script defines a series of bins of 1.5 nm in the z-direction starting at the low boundary of
the pore region. The analysis proceeds in 0.5 nm increments of the bin center plane, from the i
th
to the i+1
th
bin (see Figure 4.1). For each bin the script examines all the coordinates of the POPC
nitrogen atoms and the POPC phosphorus atoms within the bin, and finds the maximum and
minimum x and y coordinates of those atoms. The x and y coordinates of the pore center of each
bin are defined as the midpoint between the maximum and minimum values in the x and y
directions. The pore radius of each bin is defined as the average distance in the x-y plane from
the pore center to all POPC nitrogen and phosphorus atoms within the bin. At the end the script
72
selects the smallest radius as the pore radius at that time step. When the shape of the pore is not
perfectly cylindrical, the script finds the approximate radius at the minimum point.
Figure 4.1 Pore radius measurement in a POPC electropore. Only phosphorus (gold) and
nitrogen (blue) atoms are displayed for clarity. The red lines represent the boundaries of the pore
region. The solid yellow box represents the i
th
bin and the dashed yellow box represents the i+1
th
bin at a given step of the analysis.
4.3.5 Current measurement
A custom Perl script is used to extract the current from the conductance simulations. The
script first defines the pore region as described above. The pore center plane is defined as the
mid-plane between the upper bound and the lower bound of the pore region. The Perl script
counts the ions of each species that cross through the center plane and divides the count by the
elapsed time to obtain the current. In most cases the current was measured for 50 ns after the
pore radius reached a steady value.
73
4.3.6 Conductance calculation
For each conductance calculation we define the bilayer thickness as the time-averaged
pore region thickness from the current measurement plus 1 nm (the water density is
approximately the bulk water density at a distance of 0.5 nm beyond the pore region). The
potential across the pore is then defined as the bilayer thickness times the value of the external
field. The conductance is obtained by dividing the current by the potential across the pore.
4.3.7 Electropore ion concentration measurement
The ion concentration inside an electropore is evaluated by counting the number of ions
and water molecules inside the pore region, then extracted from the relation
�
W¦§
�
1¨©ª«
=
¬
W¦§
¬
1¨©ª«
,
where the concentration of water n
6h;k=
is assumed to be the bulk value 55.6 M.
74
4.4 Results and Discussion
4.4.1 Ion interactions with the phospholipid interface
In our simulations Na
+
binds more strongly to the phospholipid interface than K
+
or Cl
-
,
with bound Na
+
located deep in the interface at the acyl oxygens of the glycerol backbone and K
+
and Cl
-
showing little affinity for the bilayer, consistent with previous MD studies [44, 63].
Figure 4.2 shows radial distribution functions for both Na
+
and K
+
with respect to the backbone
C=O oxygens (sn-1 and sn-2 positions are merged into one binding data set). With roughly
equivalent ion concentrations in the bulk water (0.11 M Na
+
, 0.12 M K
+
), the Na
+
peak is orders
of magnitude higher than the K
+
peak. Figure 2.2 also demonstrates graphically the validity of
the choice of 0.3 nm as the limiting separation for bound ions and acyl oxygens.
Figure 4.2. Left panel: snapshot of a single POPC molecule. The main ion binding sites (both
glycerol backbone acyl oxygens, C=O) are represented by purple spheres. Right panel: radial
distribution function g(r) between both glycerol backbone acyl oxygens and cations. g(r) is
plotted on a log scale in order to display both curves.
75
This Na
+
-K
+
difference is argued to be exaggerated by the default GROMACS
parameters used in these simulations [44, 69, 75]. Because the van der Waals radius (σ) of Na
+
in
this model is much less than that of K
+
(σ
K
/σ
Na
≈ 2.5), an atom can approach closer to Na
+
than to
K
+
, and the electrostatic attraction between a negatively charged atom or a region of high
electron density and Na
+
is stronger. For other force fields, such as CHARMM [7] and KBFF
[68, 149], σ
K
is much smaller than the GROMACS parameter (σ
K,CHARMM
/σ
K,GROMACS
≈ 0.55 and
σ
K,KBFF
/σ
K,GROMACS
≈ 0.47). With these alternate parameters, K
+
exhibits more significant
attraction to POPC and resides deeper in the membrane [44, 69].
4.4.2 Pore formation and stabilization
In the porating external field E
p
= 400 MV/m, an electropore is created in both NaCl and
KCl POPC systems within 2 ns. In this high electric field the pore will continue to expand until
finite size effects dominate the simulation and the bilayer becomes disorganized. To stabilize the
pore, we reduce the field to a lower sustaining value E
s
when the pore radius reaches
approximately 1 nm. Although the choice of E
p
and the target radius are somewhat arbitrary, we
have shown previously that the initial configuration does not influence the steady-state properties
of the electropore [35].
We investigated sustaining field values E
s
= 30, 40, 50, 75, 100, and 120 MV/m for both
NaCl and KCl systems. For E
s
= 30 MV/m, the external field failed to sustain the pore in the
NaCl system (the pore annihilated within 10 ns in 2 independent trials). For E
s
= 40 MV/m, a
pore is sustained for more than 15 ns in the NaCl system, but we did not observe significant Na
+
transport during this time (probably the pore is too small) and could not obtain a reproducible ion
76
current with this field. In KCl systems the electropore was sustained for over 15 ns at E
s
= 30
MV/m (and also at 40, 50, 75, and 100 MV/m) in 3 independent trials. At E
s
= 120 MV/m, the
bilayer became disorganized within 100 ns in both NaCl and KCl systems, a result of finite size
effects and the periodic boundary conditions.
For E
s
= 50, 75, and 100 MV/m, electropores are stabilized for at least 200 ns in both
NaCl and KCl systems. Figure 4.3 shows the evolution of the pore size after the external field is
reduced to E
s.
The radius (initially about 1 nm in each case) reaches a field-stabilized value
within 30 ns. The stabilized, time-averaged pore radius for each E
s
, calculated from the last 50 ns
of each simulation for both NaCl and KCl systems, increases linearly with E
s
over this range of
sustaining field values (Figure 3). For a given value of E
s
, the pore radius for the NaCl system is
0.14 nm smaller than for the KCl system, probably a result of the more extensive Na
+
binding to
POPC, which decreases the area per lipid and increases the pore line tension [75].
77
Figure 4.3 Left panels: Pore radius versus time from selected conductance simulations. At t = 0
the external field is reduced to the values in the legend. Right panel: Time-averaged pore radius
versus sustaining field E
s
. Results are averaged over three independent trials. Errors are given as
standard deviation of the mean.
4.4.3 Na
+
, K
+
, and Cl
-
conductance
Table 4.1 summarizes the electrophoretic mobilities and molar conductivities in bulk
water extracted from our simulations of Na
+
, K
+
, and Cl
-
. The mobility of Na
+
is less than that of
K
+
and Cl
-
, a consequence of the larger hydrodynamic radius of Na
+
. At infinite dilution, the
ratio between the known values of K
+
molar conductivity and Na
+
molar conductivity [137]
(Λ
0
K
/ Λ
0
Na
) is approximately 1.5. We obtained approximately the same ratio in our simulations,
(Λ
0
K
/ Λ
0
Na
≈ 1.6). At 0.15 M the molar conductivity of ions decreases according to Kohlrausch’s
78
law. As a result of cation-anion interactions, the ratio between the molar conductivities Λ
K
/Λ
Na
in our simulation increases to approximately 1.8.
Table 4.1 Electrophoretic mobilities and conductivities in bulk water.
Ion
Simulated
electrophoretic mobility
(m
2
V
-1
s
-1
)
Simulated molar
conductivity at 0.15 M,
310 K
(S m
2
mol
-1
)
Na
+
0.64 × 10
-7
6.14 × 10
-3
K
+
1.12 × 10
-7
10.78 × 10
-3
Cl
-
(in NaCl) 1.01 × 10
-7
9.73 × 10
-3
Cl
-
(in KCl) 1.06 × 10
-7
10.23 × 10
-3
Ion
Simulated limiting
molar conductivity at
310 K
(S m
2
mol
-1
)
Limiting molar
conductivity at 298 K
(S m
2
mol
-1
)
Na
+
7.7 × 10
-3
5.01 × 10
-3
K
+
12.1 × 10
-3
7.35. × 10
-3
Cl
-
11.1 × 10
-3
7.63 × 10
-3
Electrophoretic mobility is obtained from 5 ns of MD simulation with NVT ensemble. Simulated
molar conductivity is calculated from the electrophoretic mobility.
The conductance of an electroporated lipid bilayer for a given ion species is determined
by the native ion conductance, by ion-pore interactions, and by the radius of the pore, which can
be tuned with the sustaining field E
s
. Ion conductance increases with pore radius, as one would
expect (Table 4.2), but Na
+
conductance is significantly less than K
+
conductance for a given
radius. This is a consequence of stronger Na
+
binding to the interface, which results also in the
nonlinear relation between pore radius and Na
+
conductance at small values for the pore radius
(Figure 4.4). At r
pore
= 0.77 nm (the smallest pore radius measured in the KCl system), the ratio
of the cation conductances, G
K
/G
Na
, is approximately 5.6, a reflection of the stronger interaction
between Na
+
and the phospholipid head groups in the pore wall relative to K
+
. At r
pore
= 1.5 nm
79
(the largest pore radius measured in the NaCl system), G
K
/G
Na
is only 2.3. As the pore size
increases, a larger and larger fraction of the sodium ions traversing the pore sees only bulk-like
water and is not slowed by pore wall binding, thus the Na
+
conductance increases relative to the
K
+
conductance and G
K
/G
Na
decreases. Preliminary data shows that this ratio approaches 1.8, the
value in bulk water simulations, as the pore radius increases further. K
+
and Cl
-
interact very
little with the pore wall; their conductances in lipid electropores are similar to those in bulk
water.
Figure 4.4 Ion conductance versus pore radius. Conductances of K
+
and Cl
-
were obtained from
the same sets of simulations.
80
Table 4.2 NaCl and KCl conductance in POPC electropores.
System
Sustaining
field
(MV/m)
Pore
radius
(nm)
Bilayer
thickness
(nm)
Cation
current
(nA)
Anion
current
(nA)
Cation
molar
conductance
(nS/M)
Anion molar
conductance
(nS/M)
NaCl 50
0.68
± 0.03
4.71
± 0.02
0.03
± 0.02
0.20
± 0.03
0.95 ± 0.47 6.43 ± 0.85
75
1.06
± 0.06
4.87
± 0.06
0.18
± 0.05
1.12
± 0.16
3.76 ± 1.15 23.17 ± 3.39
100
1.50
± 0.03
5.05
± 0.07
0.63
± 0.05
2.48
± 0.14
9.43 ± 0.82 37.31 ± 2.28
KCl 50
0.77
± 0.04
4.62
± 0.01
0.18
± 0.01
0.15
± 0.01
5.28 ± 0.28 4.37 ± 0.56
75
1.20
± 0.02
4.71
± 0.02
0.77
± 0.03
0.81
± 0.02
14.92 ± 0.57 15.71 ± 0.40
100
1.56
± 0.05
5.07
± 0.20
1.67
± 0.07
1.72
± 0.04
22.87 ± 1.00 23.64 ± 0.28
All results are averaged over three independent trials. Errors are the standard deviation.
To evaluate the importance of the ion model, we repeated the conductance simulations
using CHARMM parameters [7] for K
+
(with other atoms in GROMACS parameters). Because
CHARMM K
+
binds to POPC more than GROMACS K
+
, probably because of the smaller van
der Waals radius of the CHARMM K
+
, we added 25 CHARMM K
+
to the system (compared to
22 GROMACS K
+
) in order to achieve a similar free K
+
concentration after equilibration. The
comparison between CHARMM K
+
and GROMACS K
+
conductance is presented in Figure 4.5.
Although the attraction between POPC and K
+
is greater using CHARMM K
+
parameters, the
conductance of K
+
is only slightly lower when the pore radius is small. For a larger pore radius,
the conductance for GROMACS K
+
and CHARMM K
+
is very similar.
81
Figure 4.5 Molar pore conductance of K
+
using GROMACS parameters and CHARMM
parameters.
In addition to the influence of the pore wall, interactions among the ions also affect the
conductance. These interactions become more likely when the ions must pass through nanoscale
pores which can exhibit selectivity and rectification [99]. In our simulations the confinement
increases the likelihood for cation-anion electrostatic attraction, which decreases the net velocity
of the ions. The snapshots in Figure 4.6 capture this behavior. We qualitatively examined the
effect of cation-anion attraction on ion transport by measuring electrophoretic mobilities in a
bulk water system with varying ion concentrations. The results are summarized in Figure 4.7.
For NaCl systems, the ion mobility decreased by nearly half when the concentration increased
from about 0.1 M to 1 M. For KCl systems the decrease in ion mobility is approximately 35%.
82
Figure 4.6 Snapshots of ion interactions during migration through an electropore. External
electric field of 75 MV/m in the +z (upward) direction was applied in the simulation. Na
+
and Cl
-
ions are represented as magenta and green spheres, respectively. Phosphorus atoms and nitrogen
atoms of POPC are shown as brown spheres and blue spheres to delineate the pore. Lipids tails
are represented as blue lines. Water molecules and other atoms are omitted for clarity.
Figure 4.7 Simulated electrophoretic mobility of Na
+
, K
+
and Cl
-
in bulk water versus ion
concentration. Each data point is the average from three independent 5 ns simulations.
83
4.4.4 Molecular simulations and experimental data
One of the goals of our study is to validate our MD simulation results by comparing them
with experimental observations [25, 53, 65, 70, 71, 101, 120]. Here we concentrate on the single-
pore conductance obtained by the method of chronopotentiometry [65, 71]. Figure 4.8 compares
conductance values from our MD simulations with those reported from chronopotentiometry
experiments for POPC electropores in KCl systems. Although the two sets of conductance values
lie within the same order of magnitude, the pore radii in our simulations are smaller than those
reported in the published chronopotentiometry data.
Figure 4.8. Comparison of simulated KCl conductance with chronopotentiometry data of
Kalinowski et al. [65] and Koronkiewicz et al. [71]. Values for pore radii for the
chronopotentiometric data were extracted from the conductance based on the model described in
Ref. [71]. Pore radii for MD data points were determined as described in the Methods.
At this time there is no direct method for measuring electropore radius experimentally.
The pore radius values used for the chronopotentiometric results were extracted from the
84
measured conductance using a simple model where ions are assumed to travel through the pore
as they would through a cylinder of bulk water. In this model the pore diameter [71] is given by
r = 2
®
‰
-
¯\
, where °
Š
is the pore conductance, ± is the length of the pore, and ² is the
conductivity of the electrolyte. Since the accuracy of the extracted pore diameter depends on the
applicability of this simple geometric model, we must place large bands of uncertainty around
the chronopotentiometric results.
An additional assumption embedded in the chronopotentiometric data is that the ion
concentration inside the pore is the same as it is in the electrolyte. To examine this assumption
we evaluated the KCl concentration inside the electropores in our simulations (Table 4.3). These
results show that concentrations of the confined KCl in a conducting electropore are about twice
the bulk value (obtained during equilibration). Although it is not possible to say at this time
which pore model is more correct, the difference in pore ion concentrations between the two
models must also be taken into account when considering the correspondence between our
molecular simulations and this experimental data.
Table 4.3 KCl concentration inside electropore from MD simulation.
Sustaining field (MV/m) K
+
concentration (M) Cl
-
concentration (M)
50 0.36 ± 0.09 0.15 ± 0.02
75 0.23 ± 0.01 0.21 ± 0.01
100 0.25 ± 0.02 0.21 ± 0.02
85
4.4.5 Effects of Water model on conductance simulation
In addition to the ion model, the water model chosen for the simulations also have an
effect on electropore conductance. The SPC model has been used widely in atomistic simulations
of lipid bilayers, but its simplicity carries penalties. For example, the self-diffusion coefficient of
SPC water is about 1.8 times higher than the known experimental value [79], and this implies
lower viscosity (based on Stokes-Einstein equation) and higher ion mobilities. In contrast, the
self-diffusion coefficient of SPC/E water is about 65% of the self-diffusion coefficient of SPC
water.
To examine the importance of the water model for ion conductance in our simulated
electropores, we measured the velocities of single Na
+
and K
+
in bulk SPC and bulk SPC/E water
(Figure 4.9). As expected, the electrophoretic mobility of Na
+
in SPC water is approximately 1.8
times higher than it is in SPC/E water, and the electrophoretic mobility of K
+
in SPC water is
approximately 1.3 times higher than that of K
+
in SPC/E water. In conductance simulations with
pores that are large enough that a significant portion of the electropore consists of bulk water, the
conductance will be determined largely by the mobility of the ions in the bulk water phase. Thus
as the pore radius increases, the choice of water model will have a stronger influence on the
conductance.
86
Figure 4.9 Velocity of single Na
+
and K
+
in SPC and SPC/E water under electric field. Each data
point is the average from three independent 5 ns simulations.
87
4.5 Summary
Electrophoretic ion transport through field-stabilized lipid nanopore can be simulated
using the methods of molecular dynamics. The radius of the electropore can be adjusted by
varying the external field. With a bilayer system of 128 lipids and 9000 water molecules, an
electropore can be stabilized with external field values from 50 MV/m to 100 MV/m,
corresponding to an approximate pore radius range of 0.7 nm to 1.6 nm. Ion conductance is
determined by two factors: pore radius and the extent of ion interaction with the phospholipid
aqueous interface along the pore walls. The stronger binding of Na
+
results in significantly lower
conductance for NaCl systems compared to KCl systems. This difference decreases as the pore
radius is increased and the pore cross section becomes more and more dominated by bulk water,
reducing the probability of pore wall interactions during transport. The electric conductance of a
simulated electropore in a KCl-POPC system is roughly comparable to the single-pore
conductance extracted from chronopotentiometric measurement, but different assumptions
underlying the values used for pore radius and conductance prevent precise comparison of the
two approaches. In addition, the choice of ion and water models also affects the ion conductance
values extracted from simulations. Better models and new experimental methods are needed to
enable quantitative comparisons between atomistic simulations and experiment.
88
Chapter 5: Effect of monovalent ions on
phospholipid bilayer and electroporation
5.1 Abstract
Monovalent ions are common ingredients of the intracellular and extracellular fluid, and
they are important for many biological functions. To examine ions’ effects on electroporation,
we systematically performed molecular dynamics (MD) simulation of electroporation on
phospholipid bilayer with different NaCl and KCl concentrations. We consider two lipid bilayer
systems in our study: pure POPC bilayer and mixed POPC/POPS bilayer. The major difference
between the cations is that Na
+
exhibits strong affinity to the phospholipids and Cl
-
, and the
affinity K
+
exhibits toward other molecules is very weak. Attractions between cations and
phospholipids impose various effects on the lipid bilayer and electroporation, such as decrease in
area per lipid, increase in water group separation and the pore creation time. These effects are
much more noticeable in the systems with NaCl compare the systems with KCl, due to the
stronger cation-lipid attraction of Na
+
. We also examine the pore conductance of the ions under
various ion concentrations. Translocation of POPS via the electropore wall is observed in the
simulation with the POPC/POPS bilayer. Two different modes of POPS translocation are
identified: in one mode, POPS is translocated to the anode-facing leaflet of the bilayer by the
external field. In another mode, POPS forms a group with the cations and the whole group is
translocated toward the cathode-facing leaflet of the bilayer. POPS translocation rates of both
modes are influenced by the ion condition in the system.
89
5.2 Introduction
Electroporation (also known as electropermeabilization) utilizes electric pulse to modify
the permeability of cell membrane and allows absorptions of normally excluded substances into
the cell. The effects of external electric pulse on cell membrane provide the basis for various
applications in biological research (such as gene electrotransfer [36, 84, 93] and cell
electrofusion [22, 82, 159]), and in medical treatment (such as electrochemotherapy [46, 81] and
electrogenetherapy [40, 86]). For decades, electroporation has been studied extensively with
experiments [6, 67, 97, 155] and theoretical models [27, 117, 138, 151, 152]. These efforts
provide the phenomenological understanding of electroporation, but the underlying molecular
mechanism and structure of the permeabilized membrane is still unclear. One of the current goals
in the study is to acquire a comprehensive microscopic description of electroporation. Such
description is currently inaccessible with direct experimental observation, due to the nanometer
size of the electropore and the nanosecond scale of the pore formation process. To overcome this
issue, molecular dynamics (MD) simulation has become one of the major tools to study
electroporation in the microscopic level [42, 59, 123, 128, 130, 122], and molecular model of
electroporation is developed based on the MD studies [76, 129, 148].
In MD simulation, the cell membrane is usually represented by the lipid bilayer. Ions and
other biomolecules are often omitted for simplicity and computing efficiency. However, real
biological membranes are surrounded by aqueous medium that contains ions such as Na
+
, K
+
,
Ca
2+
, Mg
2+
and Cl
-
. Concentrations of these ions are crucial for the electrophysiology of living
cell, and they are maintained at physiological condition by mechanisms such as osmosis and ion
90
channels [52]. The integrity of ion concentrations can be disrupted during the permeabilized state
of the cell, and leads to dramatic effects such as cell swelling [133].
Ions impose various effects on the lipid bilayer, such as decrease in area per lipid, lateral
lipid mobility and head group angles respect to bilayer normal, as shown in experiments and MD
simulations [14, 15, 37, 44, 63, 134, 135, 147]. The effects of ions impose on electroporation
have been studied to a limited extent by MD simulation: Simulation of tension-stabilized pore
showed that the bindings of Na
+
at the water/lipid interface increase the pore line tension and
reduce pore lifetime [75]. Reduction of pore lifetime was also observed in electroporation
simulation with Ca
2+
, and the presents of Ca
2+
inhibit electropore formation process [77]. In this
work, we will extend these analyses by examining the properties of lipid bilayer and the pore
creation process under various concentrations of NaCl and KCl.
Recent development in MD protocol allows simulation of nanometer size field-stabilized
electropore [35]. These electropore can be sustained for over hundreds of nanoseconds by the
external electric field, and the pore conductance can be evaluated using these simulations [55].
The increase of the cell membrane’s electric conductance is one of the indicators of a cell
undergo electropermeabilization [25, 98, 120], and electric conductance of a single pore can be
obtained experimentally through methods such as chronopotentiometry [65, 70, 71]. Previously
we compared the single pore conductance obtained from MD simulation with
chronopotentiometry data, and the results from both methods are comparable within one order of
magnitude [55]. In order to obtain a more comprehensive view on the electrophoretic transport of
ions during electroporation, we examine the conductance of Na
+
, K
+
and Cl
-
in various
concentrations.
91
We examined two lipid bilayer systems in the simulations: A pure 1-palmitoyl-2-oleoyl-
sn-glycero-3-phosphatidylcholine (POPC) bilayer and a bilayer composed of both POPC and 1-
palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylserine (POPS). Normally, phosphatidylserine (PS)
are confined in the inner leaflets of the cell membrane. Externalization of PS is a signal for a cell
undergoes apoptosis and the externalized PS are markers for phagocytosis [33]. Using intense
electric pulse, PS can be artificially translocated from the inner leaflet to the outer leaflet of the
cell membrane within orders of nanoseconds [141, 142]. MD simulation demonstrated that the
electropore formation facilitates PS translocation process by providing the pathway for PS to
maneuver from one leaflet to the other [60, 144]. In this work, we will extract the PS
translocation rate using the field-stabilized electropore simulation and examine the effects of ions
impose on the process.
92
5.3 Methods
5.3.1 Molecular dynamics simulations
All simulations were performed using GROMACS version 4.0.5 [50] on the University of
Southern California High Performance Computing and Communications (HPCC) Linux cluster
(http://www.usc.edu/hpcc/). Lipid topologies were taken from OPLS united-atom parameters
[13], and the Simple Point Charge (SPC) model
[9] was used for water. For Na
+
, K
+
, Cl
-
, we
employed the set of parameters originally developed by Straatsma and Berendsen [119] that are
supplied with the GROMACS force field. We also performed additional simulation with K
+
using the CHARMM parameters [7]. All simulations were performed under constant particle
number, pressure, temperature (NPT) ensemble. Temperature was held at 310 K with an external
heat bath using a velocity-rescaling algorithm [19]. Pressure was held at 1 bar using a weak
coupling algorithm [10], with a relaxation time of 1 ps and compressibility of 4.5 × 10
-5
bar
-1
semi-isotropically applied in both normal and in-plane directions relative to the membrane. Bond
lengths were constrained using the LINCS algorithm
[49] for lipids and the SETTLE algorithm
[87] for water. Short-range electrostatic and Lennard-Jones interactions were cut off at 1.0 nm.
Long-range electrostatic interactions were calculated with PME algorithm [32] 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 in all
directions were imposed to mitigate system size effects.
93
5.3.2 Structure
All simulations contain a bilayer composed of 128 lipids (64 per leaflet) and
approximately 9000 water molecules (~70 waters/lipid), corresponding to an initial system box
size of approximately 7 nm × 7 nm × 10 nm. The mixed POPC-POPS bilayer is composed by
108 POPC and 20 POPS, with the POPS randomly distributed on the +z leaflet. The GROMACS
function ‘genion’ was used to replace bulk water molecules with Na
+
, K
+
, and Cl
-
in the system.
Custom Perl programs were used to check the number of ion absorbed in the bilayer and the
number of ion in the bulk water phase. We equilibrate the systems for at least 150 ns so that the
numbers of free and bound ions in the system and the area per lipid reach a steady value. The
area per lipid is calculated as the lateral area of the simulation box divided by the number of
lipids in a single leaflet (64 lipids per leaflet). Molecular graphics images were generated with
Visual Molecular Dynamics (VMD) [58].
5.3.3 Evaluation of ion distribution, ion absorption, and bulk ion concentration
Ion distributions, ion absorptions and bulk ion concentrations were obtained from the last
100 ns of equilibration. Ion distributions were evaluated using GROMACS function ‘g_density’
with each system being divided into 100 slides in z direction. A custom Perl script was used to
evaluate the absorption of ion to the lipid bilayer and the ion concentrations of the bulk water
phase. The script locates the density peak of the phospholipid nitrogen atoms in z direction for
each bilayer leaflet. Any ion that is located in between the phospholipid nitrogen density peaks is
considered absorbed by the lipid bilayer. Any ion that is more than 0.5 nm outside of the
phospholipid nitrogen density peaks is considered to be in the bulk water phase. The bulk ion
94
concentration is calculated using the relation
�
W¦§
�
1¨©ª«
=
¬
W¦§
¬
1¨©ª«
, with the concentration of water
n
6h;k=
is assumed to be the bulk value 55.6 M.
5.3.4 Electroporation simulations and pore initiation time measurement
Electropore were created by applying an external field with porating field value E
p
= 350
MV/m (approximately the minimum porating field [156]). By assigning every atom a
randomized velocity at the beginning of the simulation (a built-in function of GROMACS), nine
independent trials of electropore formation were performed under each ion condition. Pore
initiation time is defined as the time for the top and bottom water groups to first merge (without
subsequently separating) after the application of the external electric field [76].
5.3.5 Conductance simulations
Conductance simulations were performed utilizing the protocol we developed previously
[55] that was based on the field-stabilized electropore simulation developed by Fernández et al.
[35]. An electropore on the lipid bilayer is first created during the electroporation simulation
described in Section 5.3.4. When the phosphorus groups from top and bottom leaflets merge
[76],
and the pore radius reaches approximately 1.0 nm, we reduced the external field to a lower
sustaining field value of E
s
= 75 MV/m. Under this condition, the pore is maintained in a steady
state. If ions are present in the system, they are electrophoretically driven through the pore,
creating the ionic current. Conductance simulations were run for 300 ns after the reduction of the
95
external field to E
s
, and the ionic currents were sampled during the last 200 ns, when the pore
radius is stabilized in the steady-state value.
5.3.6 Pore radius measurement
A custom Perl script is used to evaluate the radius of the electropore. At each time frame,
the script first shifts the configuration according to the water density distribution in all 3
dimensions, so that the pore is centered in the simulation box. The upper and lower bounds of the
pore region are defined using the average of the 15 largest and smallest z-coordinates of the
POPC phosphorus atoms. Next, the script defines a bin of 1.5 nm in the z-direction center at the
midpoint of the pore region boundaries. The script examines all the coordinates of the POPC
nitrogen atoms and the POPC phosphorus atoms within the bin, and finds the maximum and
minimum x and y coordinates of those atoms. The x and y coordinate of the pore center is
defined as the midpoint between the maximum and minimum values in the x and y directions.
The pore radius of the time frame is defined as the average distance in the x-y plane from the
pore center to all POPC nitrogen and phosphorus atoms within the bin.
5.3.7 Current measurement
A custom Perl script is used to extract the current from the conductance simulations. The
script first defines the pore region as described in the pore radius measurement. The pore center
plane is defined as the mid-plane between the upper bound and the lower bound of the pore
region. The Perl script counts the ions of each species that cross through the pore center plane
96
and divides the count by the elapsed time to obtain the current. The current is evaluated for every
50 ns during the 200 ns sampling window of the conductance simulation.
5.3.8 Conductance calculation
For each conductance calculation, we define the bilayer thickness as the time-averaged
pore region thickness from the current measurement plus 1.0 nm (the water density is
approximately the bulk water density at a distance of 0.5 nm beyond the pore region). The
potential across the pore is then defined as the bilayer thickness times the value of the external
field. The conductance is obtained by dividing the current by the potential across the pore.
97
5.4 Results and Discussion
5.4.1 Ion distribution and ion absorption in lipid bilayer
The relative densities of ions in z direction in each system are illustrated in Figure 5.1.
Ion distribution is influenced by the interactions between the ions and phospholipids. Strong
attraction between Na
+
and glycerol backbone acyl oxygens (C=O) causes the bound Na
+
to
reside in the vicinity of glycerol-carbonyl interface [44, 134], and produce density peaks deep
within the lipid bilayer. K
+
with GROMACS parameter exhibits very little affinity to the
phospholipids due to its large Van der Waals radius [44], and the lack of affinity causes most of
the K
+
to locate in the bulk water phase. This lack of binding between K
+
and phospholipids is
argued to be exaggerated by the GROMACS parameters [44, 69]. In contrast, Van der Waals
radius of K
+
in CHARMM parameters is smaller than GROMACS parameters
(σ
K,CHARMM
/σ
K,GROMACS
≈ 0.55). The smaller Van der Waals radius allows CHARMM K
+
to
exhibit some attraction toward phospholipids and small CHARMM K
+
density peaks are located
inside the lipid bilayer.
Similar to GROMACS K
+
, the attraction between Cl
-
and phospholipid is also very weak,
and majority of the Cl
-
reside outside of the lipid bilayer. Distribution of Cl
-
is influenced by the
cation distribution. If cations are absorbed in the lipid bilayer (such as Na
+
and CHARMM K
+
),
Cl
-
are also accumulated near the bilayer surface to form the double layer structure. In the case of
KCl system with GROMACS K
+
, the double layer structure is absent due to the lack of binding
between the cations and phospholipids.
98
Ion distribution in pure POPC bilayer system is more or less symmetric because of the
symmetry in the system. Because the POPS are embedded in the +z leaflet of the POPC/POPS
bilayer and the negatively charged POPS are more attractive to cations compare to the
zwitterionic POPC, the cation density peak is higher in the +z side of the POPC/POPS bilayer
system.
Figure 5.1 Ion distributions in z direction during equilibration (no external field). Approximately
100 cations and 100 anions are present in each system. The relative density is normalized by the
maximum value (i.e. ρ
relative
(z) = ρ(z) / ρ
max
)
..
Density profile of phospholipid nitrogen (denoted
by N) and phospholipid phosphorus (denoted by P) are displayed to indicate the location of the
lipid bilayer.
99
For those cations that are attractive to phospholipid, their density peaks are located within
the phospholipid nitrogen density peaks (Figure 5.1). Using the nitrogen peaks as criterion to
define the lipid bilayer region, we examined the absorption of cations into the lipid bilayer, and
the results are illustrated in Figure 5.2. Because the attraction of Na
+
toward phospholipids is
much stronger than K
+
, the absorption of Na
+
into the bilayer is much larger than the absorption
of K
+
(following the order of the Hofmeister series [68]). The absorption of K
+
is also larger
while using CHARMM K
+
compare to GROMACS K
+
. Comparing the pure POPC bilayer with
the POPC/POPS bilayer, absorption of cations by the POPC/POPS bilayer is larger due to the
electrostatic attractions from the negatively charged POPS.
Figure 5.2 Number of cations absorbed by the lipid bilayer versus total number of cations in the
system. A cation is considered absorbed if it is located in between the phospholipid nitrogen
density peaks. Each data point is the time-averaged value of 100 ns trajectories during
equilibration. Error is given as the standard deviation amongst the time frames.
100
Bulk ion concentrations of each system are listed in Table 5.1 and Table 5.2. Because of
the low absorption of K
+
and Cl
-
into the bilayer, their concentration in the water phase increase
almost linearly with the number of ions in the system. Physiological concentration of K
+
in
intracellular fluid is about 0.15 M. Our POPC bilayer system with 30 K
+
and POPC/POPS
bilayer system with 50 K
+
are close to this physiological condition. Physiological concentration
of Na
+
in extracellular fluid is about 0.15 M. Our POPC bilayer system with 40 Na
+
and
POPC/POPS bilayer system with 50 Na
+
are close to this physiological condition.
Table 5.1 Bulk ion concentrations in pure POPC bilayer system.
128 POPC
Number of
Na
+
Number of
Cl
-
Bulk Na
+
Concentration (M)
Bulk Cl
-
Concentration (M)
10 10 0.0053 ± 0.0049 0.063 ± 0.009
20 20 0.040 ± 0.009 0.12 ± 0.02
30 30 0.083 ± 0.012 0.17 ± 0.02
40 40 0.14 ±0.02 0.23 ± 0.02
50 50 0.20 ± 0.02 0.29 ± 0.02
70 70 0.31 ± 0.02 0.39 ± 0.03
100 100 0.53 ± 0.02 0.61 ± 0.03
150 150 0.86 ± 0.03 0.92 ± 0.04
Number of K
+
(GROMACS)
Number of
Cl
-
Bulk K
+
Concentration (M)
Bulk Cl
-
Concentration (M)
10 10 0.061 ± 0.012 0.071 ± 0.008
30 30 0.20 ± 0.02 0.21 ± 0.02
50 50 0.33 ± 0.02 0.36 ± 0.02
70 70 0.48 ± 0.02 0.51 ± 0.02
130 130 0.92 ± 0.03 0.95 ± 0.03
Number of K
+
(CHARMM)
Number of
Cl
-
Bulk K
+
Concentration (M)
Bulk Cl
-
Concentration (M)
10 10 0.043 ± 0.010 0.068 ± 0.010
30 30 0.12 ± 0.02 0.20 ± 0.02
50 50 0.29 ± 0.02 0.34 ± 0.02
70 70 0.41 ± 0.02 0.47 ± 0.03
90 90 0.55 ± 0.03 0.61 ± 0.03
140 140 0.90 ± 0.03 0.97 ± 0.03
101
Table 5.2 Bulk ion concentrations in POPC/POPS bilayer system
108 POPS 20 POPS
Number of
Na
+
Number of
Cl
-
Bulk Na
+
Concentration (M)
Bulk Cl
-
Concentration (M)
20 0 0 0
50 30 0.12 ± 0.01 0.18 ± 0.02
80 60 0.28 ± 0.02 0.35 ± 0.03
110 90 0.46 ± 0.02 0.55 ± 0.03
150 130 0.72 ± 0.03 0.78 ± 0.04
Number of K
+
(GROMACS)
Number of
Cl
-
Bulk Na
+
Concentration (M)
Bulk Cl
-
Concentration (M)
20 0 0.05 ± 0.01 0
50 30 0.24 ± 0.02 0.22 ±0.01
80 60 0.45 ± 0.02 0.44 ± 0.02
110 90 0.66 ± 0.03 0.67 ± 0.03
140 120 0.87 ± 0.03 0.89 ± 0.03
Number of K
+
(CHARAMM)
Number of
Cl
-
Bulk Na
+
Concentration (M)
Bulk Cl
-
Concentration (M)
20 0 0.03 ± 0.01 0
50 30 0.18 ± 0.02 0.21 ± 0.01
80 60 0.38 ± 0.02 0.42 ± 0.02
110 90 0.58 ± 0.03 0.62 ± 0.03
140 120 0.79 ± 0.03 0.83 ± 0.03
5.4.2 Ions’ effect on area per lipid
For a cation that exhibit strong attraction toward the phospholipids (such as Na
+
), it can
bind simultaneously with multiple lipids and the cohesive force can draw the lipids closer
together, resulting a reduction in the area per lipid [44, 63, 134]. Figure 5.3 illustrates the area
per lipid with different number of cations in the system, and the decrease in the area per lipid
correlates with the number of cations absorbed by the lipid bilayer (Figure 5.2). Because of the
strong attractions between Na
+
and the phospholipids, the decrease in area per lipid is more
predominant in the NaCl systems compares to the KCl systems. Conversely, the presents of
GROMACS K
+
impose no effect on the area per lipid due to the lack of K
+
-phospholipid
102
binding. Since CHARMM K
+
exhibits some attraction toward the phospholipids, by increasing
the number of CHARMM K
+
in the system, the area per lipid of the bilayer is also decreased
slightly.
Because of the strong attractions between the negatively charged POPS and cations, the
area per lipid in the POPC/POPS bilayer is lower compares to the pure POPC bilayer. The area
per lipid is about 0.2 nm
2
smaller in the POPC/POPS bilayer with a given number of cations in
the system. Because of the steric force between the lipids, the area per lipid cannot be reduced
indefinitely. The value of area per lipid in the POPC/POPS bilayer appears to reach a minimum
of approximately 0.56 nm
2
with 80 Na
+
, and additional Na
+
in the system does not reduce the
area per lipid further.
Figure 5.3 Area per lipid with different number of cations in the system. Each data point is the
time-averaged value of 100 ns trajectories during equilibration. Error is given as the standard
deviation amongst the time frames.
103
5.4.3 Water phase separation and pore initiation
One of the consequences from the reduction in area per lipid is the increase of lipid
bilayer thickness [63, 134], as well as the barrier function. We examined ions’ effect on the
barrier function by measuring the separation between the upper and lower water groups (defined
by the 50% bulk water density points), and the results are illustrated in Figure 5.4. The increase
in the water group separation correlates with the decrease in area per lipid, and similar to area per
lipid, the change of water group separation is also larger in NaCl system, due to the stronger
cation-phospholipid attractions.
Figure 5.4 Water group separations with different number of cations in the system. Water group
separation is defined by the difference in the z coordinates between the upper and lower 50%
bulk water density point. Each data point is the time-averaged value of 100 ns trajectories during
equilibration. Errors are given as the standard deviation amongst the time frames.
104
Figure 5.5 Pore initiation time under E
ext
= 350 MV/m with different number of cations in the
system. Direction of the external field in POPC/POPS system is denoted as “cathode” for -z and
“anode” for +z. Each data point is the time-averaged value of nine independent trials. Error is
given as the standard deviation amongst the trials.
Formation of water bridge across the lipid bilayer is a key event in the electropore
creation process [129, 148]. Previously we examined the water bridge formation between two
water phases using a simple water-vacuum-water system [56]. If the separation between the
water phases is larger, water bridge formation is more difficult, and a longer water bridge
formation time would be required for the process. Pore initiation times (defined by the water
bridge formation [76]) under different ion conditions are illustrated in Figure 5.5. Because the
presences of Na
+
increase the bilayer thickness and the water group separation, pore initiation
time is larger in both lipid bilayers with larger number of Na
+
in the system. In contrast, the
105
presences of K
+
impose almost no effect on the pore initiation time, as they impose very little
effect on the bilayer thickness and water group separation. Water group separation in
POPC/POPS bilayer is typically about 0.1 nm larger than the pure POPC bilayer, and the pore
initiation time is also larger in the POPC/POPS bilayer for any given number of cations in the
system. Figure 5.6 summarizes water group separation and the corresponding pore initiation time
in all of the systems we examined.
Figure 5.6 Water group separation in the system and the corresponding pore initiation time
under E
ext
= 350 MV/m. Data are the identical as Figure 5.4 and Figure 5.5.
By embedding POPS in the +z leaflet of the POPC/POPS bilayer, we introduced an
asymmetry in the system. Because POPS normally reside in the inner leaflet of the cell
membrane, applying external field in the +z direction in our system mimics the anode side of the
106
cell in experiment, and applying external field in the –z direction mimics the cathode side of the
cell. Fluorescence microscopy showed that pore distributions in cell electroporation are
asymmetry, where the pore density is higher on the anode side of the cell compares to the
cathode side [54, 67]. In our simulations of POPC/POPS bilayer, we did not observe any
significant difference in the pore initiation times between the applied external field in +z (anode)
direction and –z (cathode) direction. Theoretical study using continuum model suggested that the
asymmetry of the pore density originated from the ion imbalance between cell interior and
exterior [28, 29]. The lack of asymmetry in the pore initiation time in the simulation could be
due to the lack of physiological ion imbalance in the system. The effect of ion imbalance on
electroporation would be an interesting topic of future study.
5.4.4 Ions’ effect on field-stabilized pore radius
Pore radius of a field-stabilized electropore under different ion condition is illustrated in
Figure 5.7. Under E
sus
= 75 MV/m, the pore radius increase asymptotically with the number of
Cl
-
in the system. Previous MD study of tension-stabilized pore showed that Cl
-
cannot enter the
transmembrane pore through diffusion unless the pore radius is beyond certain value (> 1.1 nm),
and the potential of mean force (PMF) of a Cl
-
inside the pore is higher than the PMF of a Cl
-
in
the water phase [75]. In addition, density of Cl
-
is very low inside or near the bilayer surface
(Figure 5.1). These results suggest that there is repulsion between Cl
-
and the lipid membrane
surface. In the field-stabilized electropore simulation, Cl
-
in the water phase are forcefully pulled
through the pore by the external field. The increase of pore radius could be due to the repulsion
between the Cl
-
inside the pore and the pore wall. Notice the electropores we examined here are
107
about 1.0 nm in radius. The pore radius will be larger if a larger sustaining field is used [35, 55].
Expansion of electropore size due to Cl
-
-lipid repulsion could be less significant for larger pores
because of the larger water volume inside the pore.
Radius of a transmembrane pore is also affected by the cation-lipid bindings. In the
tension-stabilized pore simulation, adding Na
+
in the system increases the pore line tension of the
transmembrane pore and destabilizes the pore [75]. In the field-stabilized electropore simulation,
the pore radius in the steady state is reduced due to the increase in pore line tension from Na
+
-
phospholipid bindings. The asymptotic value of pore radius for electropore in system with NaCl
is about 0.2 nm smaller than the electropore in system with KCl. The steady-state pore radius is
also lower for electropore on the POPC/POPS bilayer compares to the pure POPC bilayer for any
given ion condition.
Figure 5.7 Time-averaged pore radius of field-stabilized electropore versus number of Cl
-
under
E
sus
= 75 MV/m. The cations in the system are a) Na
+
, b) GROMACS K
+
and c) CHARMM K
+
.
Direction of the external field in the POPC/POPS system is denoted as “cathode” for -z and
“anode” for +z. Time-averaged pore radius is evaluated every 50 ns trajectory. Each data point is
the average of four trajectories, and the errors are given by the standard deviation.
108
5.4.5 Electropore conductance in different ion concentration
Two determining factors of electropore conductance are the pore radius and the number
of charger carriers in the system. We have previously examined electropore conductance with
different pore radii under physiological ion concentration [55]. Here we examine the ion
conductance of an electropore under different ion concentrations. Pore conductances of the
cations are illustrated in Figure 5.8 and pore conductances of Cl
-
are illustrated in Figure 5.9.
Pore conductance of Na
+
is much lower than K
+
and C
-
because Na
+
experience strong
pore wall hindering and the cation-anion attractions inside the electropore [55]. Because of the
negatively charged POPS, pore wall hindering is stronger in the POPC/POPS electropore
compares to the POPC electropore, and causes the pore conductance to be lower in the
POPC/POPS system. Pore conductance of Cl
-
is influenced by the cations in the system because
of the cation-anion attractions, namely the pore conductance of Cl
-
is lower in NaCl systems
compares with KCl systems. Intuitively we expect ion conductance to increase with the number
of ions in the system. It is the case for both GROMACS K
+
and Cl
-
because they exhibit little
attraction with other molecules. However, pore conductance of Na
+
does not continue to increase
with number of Na
+
in the system, because the likelihood of a Na
+
encounter Cl
-
inside the
electropore is increased with the number of NaCl in the system, and the electrophoretic motion
of Na
+
and Cl
-
is hindered by the Na
+
-Cl
-
attraction. The pore conductance of Na
+
appears to
reach a maximum of 0.55 nS asymptotically with the number of Na
+
in the POPC system. This
asymptotic behave of pore conductance is also observed in system of POPC/POPS electropore
with CHARMM K
+
.
109
Figure 5.8 Pore conductance of cations versus number of cations under E
sus
= 75 MV/m.
Direction of the external field in POPC/POPS system is denoted as “cathode” for -z and “anode”
for +z. The ion conductance is evaluated every 50 ns trajectory. Each data point is the average of
four trajectories, and the error is given by the standard deviation.
110
Figure 5.9 Pore conductance of Cl
-
versus number of Cl
-
under E
sus
= 75 MV/m. Direction of the
external field in POPC/POPS system is denoted as “cathode” for -z and “anode” for +z. The ion
conductance is evaluated every 50 ns trajectory. Each data point is the average of four
trajectories, and the error is given by the standard deviation.
5.4.6 Ions’ effect on POPS translocation
Electropore behaves as a channel connecting two sides of the cell membrane, and the
pore wall provides a pathway for lipids to maneuver from one leaflet to the other [141, 144].
Because POPS is negatively charged, its motion is influenced by the surrounding charges and the
electric field, and they can be translocated to the opposite leaflet by the external field. The
snapshots in Figure 5.10 capture the motion of POPS translocation along the pore wall in such
event. The motion of POPS translocation is influenced by the cations in the system. Table 5.3
lists the numbers of translocated POPS toward the anode-facing leaflet within 300 ns trajectory
111
under different cation conditions. Because the cation-POPS attraction hinders the motion of
POPS translocation, the number of translocated POPS in systems with Na
+
is much lower than
the number in the systems with K
+
, and translocation rate is also lower in system with higher Na
+
concentration.
Interestingly, POPS translocation toward the cathode-facing leaflet is also observed in the
simulations, and it is counter-intuitive because POPS is negatively charged. This alternate mode
of POPS translocation is caused by the grouping of POPS and the cations. The snapshots in
Figure 5.11 illustrate the translocation process of POPS toward the cathode-facing leaflet. During
the translocation process, the POPS attracts multiple cations simultaneously and they form a
group. The group (POPS and cations) is overall positively charged and the POPS is translocated
toward the cathode-facing leaflet as part of the group. The rate of POPS translocation in this
mode is also influenced by the cation condition. Table 5.4 list the numbers of translocated POPS
toward the cathode-facing leaflet within 300 ns trajectory under different number of cations in
the system. Because of the strong Na
+
-POPS attractions, POPS translocation toward the cathode-
facing leaflet occurs more frequently in system with Na
+
. In contrast, POPS translocation toward
the cathode was not observed in any of the KCl system with GROMACS K
+
. In KCl systems
with CHARMM K
+
, POPS translocation toward cathode was only observed once in the system
with high K
+
concentration.
112
Figure 5.10 Snapshots of POPS translocation toward anode-facing leaflet. External field is
applied in +ẑ (upward). The translocated POPS is represented as the green molecule. Lipid tails
of POPC and POPS are represented by light blue stripes. Nitrogen and phosphorus atoms of
POPC and POPS are represented by blue and gold spheres to indicate the locations of lipid head
groups. Hydrogen and oxygen of water are represented by white and red dots.
Figure 5.11 Snapshots of POPS translocation toward cathode-facing leaflet. External field is
applied in -ẑ (downward). The translocated POPS is represented as the green molecule, and the
bond Na
+
are represented as orange spheres. Nitrogen and phosphorus atoms of POPC and POPS
are represented by blue and gold spheres to indicate the locations of lipid head groups. Lipid tails
of lipids and water are omitted for clarity.
113
Table 5.3 Number of translocated POPS toward anode leaflet by external field
Number of Na
+
Number of translocated POPS
at 100 ns at 200 ns at 300 ns
20 1 2 3
50 3 3 3
80 1 1 1
110 0 1 2
150 0 0 0
Number of K
+
(GROMACS)
Number of translocated POPS
at 100 ns at 200 ns at 300 ns
20 3 7 8
50 4 8 12
80 7 12 12
110 6 10 13
140 6 9 12
Number of K
+
(CHARMM)
Number of translocated POPS
at 100 ns at 200 ns at 300 ns
20 4 6 7
50 5 7 8
80 8 14 18
110 5 7 11
140 2 4 4
114
Table 5.4 Number of translocated POPS toward cathode leaflet by POPS-cation grouping
Number of Na
+
Number of translocated POPS
at 100 ns at 200 ns at 300 ns
20 0 0 0
50 2 0 2
80 0 0 0
110 2 2 2
150 1 3 3
Number of K
+
(GROMACS)
Number of translocated POPS
at 100 ns at 200 ns at 300 ns
20 0 0 0
50 0 0 0
80 0 0 0
110 0 0 0
140 0 0 0
Number of K
+
(CHARMM)
Number of translocated POPS
at 100 ns at 200 ns at 300 ns
20 0 0 0
50 0 0 0
80 0 0 0
110 0 1 1
140 0 0 0
115
5.5 Summary
We performed MD simulations to examine ion’s effects on lipid bilayer and the
electroporation process. Two common monovalent salts, NaCl and KCl, were examined in the
study. The major difference between the cations is that Na
+
exhibits stronger attraction toward
phospholipid, and the attraction K
+
exhibits toward the phospholipid is weak. The absorption of
Na
+
into phospholipid bilayer is larger than the absorption of K
+
. The lack of attraction between
K
+
and phospholipid is more noticeable when the GROMACS parameters are used.
Attractions between a cation and multiple phospholipids draw the phospholipids closer
together, and lead to reduction in the area per lipid. Accompany with the decrease in area per
lipid is the increase in bilayer thickness and the water groups separation. The water bridge
formation process is inhibited by enlarging the water group separation, and subsequently the pore
initiation time is also larger. These effects are influenced by the attractions of the cations with
the phospholipids, and they are more predominant in systems with Na
+
compares to systems with
K
+
.
During field-stabilized electropore simulation, Cl
-
in the water phases are dragged
through the electropore by the external field. Because of the repulsion between Cl
-
and the pore
wall, the size of the electropore is expanded when the Cl
-
enter the electropore. For a small field-
stabilized electropore (r ≈ 1.0 nm), its radius increases asymptotically with the number of Cl
-
in
the system. For cations that exhibit strong attractions with phospholipid (such as Na
+
), their
presences increase the pore line tension and cause the electropore radius to reduced.
116
Pore conductance of ions increase asymptotically with the number of ions in the system.
Because of pore wall hindering and cation-anion attraction, ions that exhibit strong attraction
toward lipid or oppositely charged ion possess much lower pore conductance. For a field-
stabilized electropore with radius about 1.0 nm, pore conductance of Na
+
is much lower than K
+
and Cl
-
. Pore conductance of Cl
-
is lower when paired with Na
+
due to the cation-anion
attractions.
Electropore wall provides a pathway for lipids to maneuver from one leaflet of the bilayer
to the other. Two different modes of POPS translocation are observed in the simulations. In the
first mode, POPS are translocated to the anode-facing leaflet by the external field, and the
motion of POPS translocation is influenced by the cations. For cations that exhibit strong
attractions to the POPS (such as Na
+
), they hinder the motion of POPS translocation and cause
the translocation rate to decrease. The second mode of POPS translocation is mediated by the
grouping of POPS and cations. In this mode, attractions between POPS and cations cause them
to form a group that is overall positively charged, and the POPS is translocated toward the
cathode-facing leaflet as part of the group. POPS translocation in this mode can only occur if the
attractions between the cations and POPS are strong, and the POPS translocation rate in this
mode is larger with higher cation concentration.
117
Conclusions and Outlook
MD simulation has been proven as a valuable tool in enhancing our understanding of
electroporation. Using MD simulation, we can observe the microscopic events during different
stages of electroporation and analyze the electropore structure, which are inaccessible by current
experimental methods. Because of these observations, we are able to develop the molecular
model and microscopic description of the phenomenon.
One of the main goals of our study is to identify the underlying force that is responsible
for electropore formation. Water bridge formation is the earliest event that can be observed in the
pore creation process by MD simulation. Energetic analyses suggest that water bridge formation
is mediated by the dipole-dipole interactions between water molecules under the influence of
external field. These observations imply water protrusion during the construction of the water
bridge plays a major role in puncturing the cell membrane during the electropore formation, and
the membrane behaves as a mechanical barrier against the process. Beside lipid bilayer system,
water bridge formation can also be observed in octane layer system and water-vacuum-water
system. By comparing the dynamics of water bridge in the water-lipid-water system (lipid
bilayer) and the water-vacuum-water system, we delineated the role of water in electroporation.
In both systems, the water bridge initiation time increases exponentially with the magnitude of
the external field, and the overall dipole moment of the system increases rapidly once the water
bridge is formed. The water bridge can be sustained in both lipid bilayer and vacuum system, and
radius of the water bridge can be tuned by adjusting the sustaining field. During the water bridge
annihilation process, lipids play the role of temporary sustaining the water bridge, and the
annihilation process is much slower for water bridge inside an electropore with larger radius.
118
Opening of an electropore allows normally blocked substances to travel through cell
membrane. MD simulation provides the microscopic view of the transport process. Using MD
simulation, we examined the pore conductance of a few common monovalent ions in biological
system: Na
+
, K
+
and Cl
-
. A major determining factor of pore conductance is the pore radius,
which can be controlled by tuning the external field. Pore conductance of a specific type of ion is
determined by the attraction of the ion toward the phospholipids and the attraction of the ion
toward the oppositely charged ions. Na
+
, for instance, exhibit strong attraction toward POPC and
Cl
-
, and its electrophoretic motion is hindered by the pore wall and the cation-anion attraction.
Single pore conductance of KCl through POPC electropore obtained from MD simulation is
comparable within one order of magnitude to those measured in chronopotentiometry, and it
provide a direct, qualitative comparison between simulation and experiment.
Interactions between molecules and phospholipids can change the structural properties of
the cell membrane. Some of these structural changes can be observed in detail using MD
simulation. In our simulation, we observed that a cation can attract multiple phospholipids
simultaneously and the cohesive force draw the lipids closer together. The area per lipid is
decreased and the bilayer thickness is increase by the process. By adding these cations in the
system, formation of water bridge becomes more difficult due to the corresponding increase in
the water group separation, and the pore formation process is also slower. Another structural
change of the cell membrane experiences during electroporation is PS externalization. In our
simulation of the mixed POPC/POPS bilayer, we observed two different modes of POPS
translocation: In the first mode, POPS is translocated toward the anode-facing leaflet by the
external field. In the second mode, attractions between the POPS and cations bundle the
molecules into a group that is overall positively charged, and the POPS is translocated toward the
119
cathode-facing leaflet as part of the group. In both modes of POPS translocation, the electropore
wall serve as a pathway connecting one bilayer leaflet to the other. The rates of POPS
translocation in both modes are influenced by the ion conditions in the system.
MD simulation provided many valuable insights in understanding electroporation, and we
are making steady progress toward a comprehensive microscopic model. However, we are still
distant from obtaining a precise quantitative description of the phenomenon. Because of the scale
of the system, many of the results obtained from MD simulation can only be compared
qualitatively or indirectly with experiment or continuum model. For example, pore density on
permeabilized membrane is one of the variables that have been studied extensively by
experiments and continuum models. Due to the limitation on system size, only one electropore
can be formed on the lipid bilayer in most the MD simulations, and they yield no information on
the pore density. On the other hand, evolution of an electropore can be easily observed using MD
simulation, but it is inaccessible for current experimental method. Our pore conductance
measurement is one of very few cases where quantitative comparison between simulation and
experiment can be achieved. Furthermore, results obtained by MD simulation are influenced by
the modeling and parameterizations of the molecules. In earlier chapters, we saw that pore
initiation time and the pore conductance of ion are influenced by the water model used in the
simulation. Improvements on molecular modeling would be a crucial step in making accurate
quantitative prediction with MD simulation.
Modern scientific research nowadays relies much more on computers. The advancement
in the computing resources expands our capability to explore the physical world. MD simulation
is one of the many examples where the use of computer allows us to explore regime that is
otherwise beyond our reach. Although there remain many challenges in making accurate
120
predictions using MD simulation, the method exhibits great potential in enhancing our
understanding in the microscopic world. Improvements for the method require continuous effort
from experiments, theoretical modeling, and the development in the simulation. Through these
efforts, we can learn more about the physical world, and use the knowledge to improve our lives
for the future.
121
122
Bibliography
1 Abidor, I.G., et al., Electric Breakdown of Bilayer Lipid-Membranes .1. Main
Experimental Facts and Their Qualitative Discussion. Bioelectrochemistry and
Bioenergetics, 1979. 6(1): p. 37-52.
2 Aihara, H. and J. Miyazaki, Gene transfer into muscle by electroporation in vivo. Nature
Biotechnology, 1998. 16(9): p. 867-870.
3 Alder, B.J. and T.E. Wainwright, Studies in Molecular Dynamics .1. General Method.
Journal of Chemical Physics, 1959. 31(2): p. 459-466.
4 Alper, H.E., D. Bassolino, and T.R. Stouch, Computer-Simulation of a Phospholipid
Monolayer-Water System - the Influence of Long-Range Forces on Water-Structure and
Dynamics. Journal of Chemical Physics, 1993. 98(12): p. 9798-9807.
5 Armstrong, C.M. and B. Hille, Voltage-gated ion channels and electrical excitability.
Neuron, 1998. 20(3): p. 371-80.
6 Arnaud-Cormos, D., et al., Microchamber setup characterization for nanosecond pulsed
electric field exposure. IEEE Trans Biomed Eng, 2011. 58(6): p. 1656-62.
7 Beglov, D. and B. Roux, Finite Representation of an Infinite Bulk System - Solvent
Boundary Potential for Computer-Simulations. Journal of Chemical Physics, 1994.
100(12): p. 9050-9063.
8 Belehradek, M., et al., Electrochemotherapy, a new antitumor treatment. First clinical
phase I-II trial. Cancer, 1993. 72(12): p. 3694-700.
9 Berendsen H.J.C., Postma J.P.M., van Gunsteren W.F. and Hermans J. Interaction models
for water in relation to protein hydration. In:Intermolecular Forces (B.Pullman, ed.),
Reidel Publishing Company, Dordrecht , 1981, pp. 331-342.
10 Berendsen, H.J.C., et al., Molecular-Dynamics with Coupling to an External Bath.
Journal of Chemical Physics, 1984. 81(8): p. 3684-3690.
11 Berendsen, H.J.C., J.R. Grigera, and T.P. Straatsma, The Missing Term in Effective Pair
Potentials. Journal of Physical Chemistry, 1987. 91(24): p. 6269-6271.
12 Berendsen, H.J.C., D. Vanderspoel, and R. Vandrunen, Gromacs - a Message-Passing
Parallel Molecular-Dynamics Implementation. Computer Physics Communications,
1995. 91(1-3): p. 43-56.
123
13 Berger, O., O. Edholm, and F. Jahnig, Molecular dynamics simulations of a fluid bilayer
of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant
temperature. Biophysical Journal, 1997. 72(5): p. 2002-2013.
14 Binder, H. and O. Zschornig, The effect of metal cations on the phase behavior and
hydration characteristics of phospholipid membranes. Chem Phys Lipids, 2002. 115(1-
2): p. 39-61.
15 Böckmann, R.A., et al., Effect of sodium chloride on a lipid bilayer. Biophys J, 2003.
85(3): p. 1647-55.
16 Böckmann, R.A., et al., Kinetics, statistics, and energetics of lipid membrane
electroporation studied by molecular dynamics simulations. Biophysical Journal, 2008.
95(4): p. 1837-1850.
17 Bowman, A.M., et al., Analysis of plasma membrane integrity by fluorescent detection of
Tl(+) uptake. J Membr Biol, 2010. 236(1): p. 15-26.
18 Brooks, B.R., et al., CHARMM: The Biomolecular Simulation Program. Journal of
Computational Chemistry, 2009. 30(10): p. 1545-1614.
19 Bussi, G., D. Donadio, and M. Parrinello, Canonical sampling through velocity rescaling.
J Chem Phys, 2007. 126(1): p. 014101.
20 Castellana, E.T. and P.S. Cremer, Solid supported lipid bilayers: From biophysical
studies to sensor design. Surface Science Reports, 2006. 61(10): p. 429-444.
21 Cemazar, M., et al., Effective gene transfer to solid tumors using different nonviral gene
delivery techniques: Electroporation, liposomes, and integrin-targeted vector. Cancer
Gene Therapy, 2002. 9(4): p. 399-406.
22 Chang, D.C., Cell Poration and Cell-Fusion Using an Oscillating Electric-Field.
Biophysical Journal, 1989. 56(4): p. 641-652.
23 Chassy, B.M., A. Mercenier, and J. Flickinger, Transformation of Bacteria by
Electroporation. Trends in Biotechnology, 1988. 6(12): p. 303-309.
24 Chiu, S.W., et al., Combined Monte Carlo and molecular dynamics simulation of fully
hydrated dioleyl and palmitoyl-oleyl phosphatidylcholine lipid bilayers. Biophysical
Journal, 1999. 77(5): p. 2462-2469.
25 Coster, H.G., A quantitative analysis of the voltage-current relationships of fixed charge
membranes and the associated property of "punch-through". Biophys J, 1965. 5(5): p.
669-86.
26 Darden, T., D. York, and L. Pedersen, Particle Mesh Ewald - an N.Log(N) Method for
Ewald Sums in Large Systems. Journal of Chemical Physics, 1993. 98(12): p. 10089-
10092.
124
27 DeBruin, K.A. and W. Krassowska, Electroporation and shock-induced transmembrane
potential in a cardiac fiber during defibrillation strength shocks. Ann Biomed Eng, 1998.
26(4): p. 584-96.
28 DeBruin, K.A. and W. Krassowska, Modeling electroporation in a single cell. I. Effects
of field strength and rest potential. Biophysical Journal, 1999. 77(3): p. 1213-1224.
29 DeBruin, K.A. and W. Krassowska, Modeling electroporation in a single cell. II. Effects
of ionic concentrations. Biophysical Journal, 1999. 77(3): p. 1225-1233.
30 Egberts, E., S.J. Marrink, and H.J.C. Berendsen, Molecular-Dynamics Simulation of a
Phospholipid Membrane. European Biophysics Journal with Biophysics Letters, 1994.
22(6): p. 423-436.
31 El-Hag, A.H., S.H. Jayaram, and M.W. Griffiths, Inactivation of naturally grown
microorganisms in orange juice using pulsed electric fields. Ieee Transactions on Plasma
Science, 2006. 34(4): p. 1412-1415.
32 Essmann, U., et al., A Smooth Particle Mesh Ewald Method. Journal of Chemical
Physics, 1995. 103(19): p. 8577-8593.
33 Fadok, V.A., et al., Loss of phospholipid asymmetry and surface exposure of
phosphatidylserine is required for phagocytosis of apoptotic cells by macrophages and
fibroblasts. J Biol Chem, 2001. 276(2): p. 1071-7.
34 Fernández, M.L., et al., Structural and kinetic molecular dynamics study of
electroporation in cholesterol-containing bilayers. J Phys Chem B, 2010. 114(20): p.
6855-65.
35 Fernández, M.L., et al., Size-controlled nanopores in lipid membranes with stabilizing
electric fields. Biochemical and Biophysical Research Communications, 2012. 423(2): p.
325-330.
36 Ferreira, E., et al., Optimization of a gene electrotransfer method for mesenchymal stem
cell transfection. Gene Ther, 2008. 15(7): p. 537-44.
37 Garcia-Manyes, S., G. Oncins, and F. Sanz, Effect of ion-binding and chemical
phospholipid structure on the nanomechanics of lipid bilayers studied by force
spectroscopy. Biophys J, 2005. 89(3): p. 1812-26.
38 Glaser, R.W., et al., Reversible Electrical Breakdown of Lipid Bilayers - Formation and
Evolution of Pores. Biochimica Et Biophysica Acta, 1988. 940(2): p. 275-287.
39 Gothelf, A., L.M. Mir, and J. Gehl, Electrochemotherapy: results of cancer treatment
using enhanced delivery of bleomycin by electroporation. Cancer Treat Rev, 2003. 29(5):
p. 371-87.
125
40 Goto, T., et al., Highly efficient electro-gene therapy of solid tumor by using an
expression plasmid for the herpes simplex virus thymidine kinase gene. Proceedings of
the National Academy of Sciences of the United States of America, 2000. 97(1): p. 354-
359.
41 Greengard, L.F. and J.F. Huang, A new version of the fast multipole method for screened
Coulomb interactions in three dimensions. Journal of Computational Physics, 2002.
180(2): p. 642-658.
42 Gurtovenko, A.A. and I. Vattulainen, Pore formation coupled to ion transport through
lipid membranes as induced by transmembrane ionic charge imbalance: Atomistic
molecular dynamics study. Journal of the American Chemical Society, 2005. 127(50): p.
17570-17571.
43 Gurtovenko, A.A. and I. Vattulainen, Ion leakage through transient water pores in
protein-free lipid membranes driven by transmembrane ionic charge imbalance.
Biophysical Journal, 2007. 92(6): p. 1878-1890.
44 Gurtovenko, A.A. and I. Vattulainen, Effect of NaCl and KCl on phosphatidylcholine and
phosphatidylethanolamine lipid membranes: insight from atomic-scale simulations for
understanding salt-induced effects in the plasma membrane. J Phys Chem B, 2008.
112(7): p. 1953-62.
45 Gurtovenko, A.A. and I. Vattulainen, Calculation of the electrostatic potential of lipid
bilayers from molecular dynamics simulations: methodological issues. J Chem Phys,
2009. 130(21): p. 215107.
46 Heller, R., et al., Phase I/II trial for the treatment of cutaneous and subcutaneous tumors
using electrochemotherapy. Cancer, 1996. 77(5): p. 964-971.
47 Heller, R., et al., In vivo gene electroinjection and expression in rat liver. Biophysical
Journal, 1997. 72(2): p. Wpm10-Wpm10.
48 Heller, L.C. and R. Heller, In vivo electroporation for gene therapy. Human Gene
Therapy, 2006. 17(9): p. 890-897.
49 Hess, B., et al., LINCS: A linear constraint solver for molecular simulations. Journal of
Computational Chemistry, 1997. 18(12): p. 1463-1472.
50 Hess, B., et al., GROMACS 4: Algorithms for highly efficient, load-balanced, and
scalable molecular simulation. J Chem Theory Comput, 2008. 4(3): p. 435-447.
51 Hess, B., P-LINCS: A parallel linear constraint solver for molecular simulation. Journal
of Chemical Theory and Computation, 2008. 4(1): p. 116-122.
52 Hess, P. and R.W. Tsien, Mechanism of Ion Permeation through Calcium Channels.
Nature, 1984. 309(5967): p. 453-456.
53 Hibino, M., et al., Membrane conductance of an electroporated cell analyzed by
submicrosecond imaging of transmembrane potential. Biophys J, 1991. 59(1): p. 209-20.
126
54 Hibino, M., H. Itoh, and K. Kinosita, Jr., Time courses of cell electroporation as revealed
by submicrosecond imaging of transmembrane potential. Biophys J, 1993. 64(6): p.
1789-800.
55 Ho, M.C., et al., Molecular Dynamics Simulations of Ion Conductance in Field-Stabilized
Nanoscale Lipid Electropores. J Phys Chem B, 2013. 117(39): p. 11633-40.
56 Ho, M.C., Z.A. Levine, and P.T. Vernier, Nanoscale, electric field-driven water bridges
in vacuum gaps and lipid bilayers. J Membr Biol, 2013. 246(11): p. 793-801.
57 Hoover, W.G., Canonical Dynamics - Equilibrium Phase-Space Distributions. Physical
Review A, 1985. 31(3): p. 1695-1697.
58 Humphrey, W., A. Dalke, and K. Schulten, VMD: Visual molecular dynamics. Journal of
Molecular Graphics & Modelling, 1996. 14(1): p. 33-38.
59 Hu, Q., et al., Simulations of transient membrane behavior in cells subjected to a high-
intensity ultrashort electric pulse. Physical Review E, 2005. 71(3).
60 Hu, Q., R.P. Joshi, and K.H. Schoenbach, Simulations of nanopore formation and
phosphatidylserine externalization in lipid membranes subjected to a high-intensity,
ultrashort electric pulse. Phys Rev E Stat Nonlin Soft Matter Phys, 2005. 72(3 Pt 1): p.
031902.
61 Jorgensen, W.L., et al., Comparison of Simple Potential Functions for Simulating Liquid
Water. Journal of Chemical Physics, 1983. 79(2): p. 926-935.
62 Jorgensen, W.L. and J. Tiradorives, The Opls Potential Functions for Proteins - Energy
Minimizations for Crystals of Cyclic-Peptides and Crambin. Journal of the American
Chemical Society, 1988. 110(6): p. 1657-1666.
63 Jurkiewicz, P., et al., Biophysics of lipid bilayers containing oxidatively modified
phospholipids: insights from fluorescence and EPR experiments and from MD
simulations. Biochim Biophys Acta, 2012. 1818(10): p. 2388-402.
64 Joshi, R.P. and K.H. Schoenbach, Electroporation dynamics in biological cells subjected
to ultrafast electrical pulses: A numerical simulation study. Physical Review E, 2000.
62(1): p. 1025-1033.
65 Kalinowski, S., et al., Chronopotentiometric studies of electroporation of bilayer lipid
membranes. Biochimica Et Biophysica Acta-Biomembranes, 1998. 1369(2): p. 204-212.
66 Kanduser, M., D. Miklavcic, and M. Pavlin, Mechanisms involved in gene electrotransfer
using high- and low-voltage pulses--an in vitro study. Bioelectrochemistry, 2009. 74(2):
p. 265-71.
67 Kinosita, K., Jr., et al., Electroporation of cell membrane visualized under a pulsed-laser
fluorescence microscope. Biophys J, 1988. 53(6): p. 1015-9.
127
68 Klasczyk, B. and V. Knecht, Kirkwood-Buff derived force field for alkali chlorides in
simple point charge water. Journal of Chemical Physics, 2010. 132(2).
69 Klasczyk, B. and V. Knecht, Validating Affinities for Ion-Lipid Association from
Simulation against Experiment. Journal of Physical Chemistry A, 2011. 115(38): p.
10587-10595.
70 Koronkiewicz, S., 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, 2002. 1561(2): p. 222-229.
71 Koronkiewicz, S. and S. Kalinowski, Influence of cholesterol on electroporation of
bilayer lipid membranes: chronopotentiometric studies. Biochimica Et Biophysica Acta-
Biomembranes, 2004. 1661(2): p. 196-203.
72 Krassowska, W. and P.D. Filev, Modeling electroporation in a single cell. Biophys J,
2007. 92(2): p. 404-17.
73 Läuger, P., et al., Electrical properties of bimolecular phospholipid membranes. Biochim
Biophys Acta, 1967. 135(1): p. 20-32.
74 Leontiadou, H., A.E. Mark, and S.J. Marrink, Molecular dynamics simulations of
hydrophilic pores in lipid bilayers. Biophysical Journal, 2004. 86(4): p. 2156-2164.
75 Leontiadou, H., A.E. Mark, and S.J. Marrink, Ion transport across transmembrane pores.
Biophysical Journal, 2007. 92(12): p. 4209-4215.
76 Levine, Z.A. and P.T. Vernier, Life Cycle of an Electropore: Field-Dependent and Field-
Independent Steps in Pore Creation and Annihilation. Journal of Membrane Biology,
2010. 236(1): p. 27-36.
77 Levine, Z.A. and P.T. Vernier, Calcium and Phosphatidylserine Inhibit Lipid Electropore
Formation and Reduce Pore Lifetime. Journal of Membrane Biology, 2012. 245(10): p.
599-610.
78 Martyna, G.J., et al., Explicit reversible integrators for extended systems dynamics.
Molecular Physics, 1996. 87(5): p. 1117-1157.
79 Mark, P. and L. Nilsson, Structure and dynamics of the TIP3P, SPC, and SPC/E water
models at 298 K. Journal of Physical Chemistry A, 2001. 105(43): p. 9954-9960.
80 Marrink, S.J., et al., The MARTINI force field: Coarse grained model for biomolecular
simulations. Journal of Physical Chemistry B, 2007. 111(27): p. 7812-7824.
81 Marty, M., et al., 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, 2006. 4(11): p.
3-13.
128
82 Mekid, H. and L.M. Mir, In vivo cell electrofusion. Biochim Biophys Acta, 2000.
1524(2-3): p. 118-30.
83 Min, S., Z.T. Jin, and Q.H. Zhang, Commercial scale pulsed electric field processing of
tomato juice. Journal of Agricultural and Food Chemistry, 2003. 51(11): p. 3338-3344.
84 Mir, L.M., et al., High-efficiency gene transfer into skeletal muscle mediated by electric
pulses. Proceedings of the National Academy of Sciences of the United States of
America, 1999. 96(8): p. 4262-4267.
85 Mir, L.M., et al., Electrochemotherapy: a new treatment of solid tumors. Journal of
Experimental & Clinical Cancer Research, 2003. 22(4): p. 145-148.
86 Mir, L.M., Nucleic acids electrotransfer-based gene therapy (electrogenetherapy): past,
current, and future. Mol Biotechnol, 2009. 43(2): p. 167-76.
87 Miyamoto, S. and P.A. Kollman, Settle - an Analytical Version of the Shake and Rattle
Algorithm for Rigid Water Models. Journal of Computational Chemistry, 1992. 13(8): p.
952-962.
88 Montal, M. and P. Mueller, Formation of bimolecular membranes from lipid monolayers
and a study of their electrical properties. Proc Natl Acad Sci U S A, 1972. 69(12): p.
3561-6.
89 Narten, A.H. and H.A. Levy, Liquid Water - Molecular Correlation Functions from X-
Ray Diffraction. Journal of Chemical Physics, 1971. 55(5): p. 2263-2269.
90 Nelson, M.T., et al., NAMD: A parallel, object oriented molecular dynamics program.
International Journal of Supercomputer Applications and High Performance Computing,
1996. 10(4): p. 251-268.
91 Nesin, O.M., et al., Manipulation of cell volume and membrane pore comparison
following single cell permeabilization with 60- and 600-ns electric pulses. Biochim
Biophys Acta, 2011. 1808(3): p. 792-801.
92 Neu, J.C. and W. Krassowska, Asymptotic model of electroporation. Physical Review E,
1999. 59(3): p. 3471-3482.
93 Neumann, E., et al., Gene-Transfer into Mouse Lyoma Cells by Electroporation in High
Electric-Fields. Embo Journal, 1982. 1(7): p. 841-845.
94 Neumann, E., et al., Calcium-mediated DNA adsorption to yeast cells and kinetics of cell
transformation by electroporation. Biophysical Journal, 1996. 71(2): p. 868-877.
95 Nosé, S., A Molecular-Dynamics Method for Simulations in the Canonical Ensemble.
Molecular Physics, 1984. 52(2): p. 255-268.
96 Okuno, Y., et al., Simulation study on the influence of an electric field on water
evaporation. Journal of Molecular Structure-Theochem, 2009. 904(1-3): p. 83-90.
129
97 Pakhomov, A.G., et al., Long-lasting plasma membrane permeabilization in mammalian
cells by nanosecond pulsed electric field (nsPEF). Bioelectromagnetics, 2007. 28(8): p.
655-63.
98 Pakhomov, A.G., et al., Lipid nanopores can form a stable, ion channel-like conduction
pathway in cell membrane. Biochemical and Biophysical Research Communications,
2009. 385(2): p. 181-186.
99 Pakhomov A.G. and Pakhomova O.N. Nanopores: A distinct transmembrane passageway
in electroporated cells. In:Advanced Electroporation Techniques in Biology in Medicine (A.
G. Pakhomov, D. Miklavcic, and M.S. Markov, eds.), CRC Press, Taylor & Francis Group,
Boca Raton, FL, 2010, pp. 177-194.
100 Parrinello, M. and A. Rahman, Polymorphic Transitions in Single-Crystals - a New
Molecular-Dynamics Method. Journal of Applied Physics, 1981. 52(12): p. 7182-7190.
101 Pavlin, M., et al., Effect of cell electroporation on the conductivity of a cell suspension.
Biophysical Journal, 2005. 88(6): p. 4378-4390.
102 Pavselj, N. and V. Preat, DNA electrotransfer into the skin using a combination of one
high- and one low-voltage pulse. Journal of Controlled Release, 2005. 106(3): p. 407-
415.
103 Pearlman, D.A., et al., 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. Computer Physics
Communications, 1995. 91(1-3): p. 1-41.
104 Piggot, T.J., D.A. Holdbrook, and S. Khalid, Electroporation of the E. coli and S. Aureus
membranes: molecular dynamics simulations of complex bacterial membranes. J Phys
Chem B, 2011. 115(45): p. 13381-8.
105 Popescu, D., C. Rucareanu, and G. Victor, A Model for the Appearance of Statistical
Pores in Membranes Due to Selfoscillations. Bioelectrochemistry and Bioenergetics,
1991. 25(1): p. 91-103.
106 Pucihar, G., et al., Kinetics of transmembrane transport of small molecules into
electropermeabilized cells. Biophys J, 2008. 95(6): p. 2837-48.
107 Pucihar, G., D. Miklavcic, and T. Kotnik, A Time-Dependent Numerical Model of
Transmembrane Voltage Inducement and Electroporation of Irregularly Shaped Cells.
Ieee Transactions on Biomedical Engineering, 2009. 56(5): p. 1491-1501.
108 Rahman, A., Correlations in Motion of Atoms in Liquid Argon. Physical Review a-
General Physics, 1964. 136(2A): p. A405-11.
109 Romeo, S., et al., Water influx and cell swelling after nanosecond
electropermeabilization. Biochim Biophys Acta, 2013. 1828(8): p. 1715-22.
130
110 Rountree, C.L., et al., Atomistic aspects of crack propagation in brittle materials:
Multimillion atom molecular dynamics simulations. Annual Review of Materials
Research, 2002. 32: p. 377-400.
111 Ryckaert, J.P., G. Ciccotti, and H.J.C. Berendsen, Numerical-Integration of Cartesian
Equations of Motion of a System with Constraints - Molecular-Dynamics of N-Alkanes.
Journal of Computational Physics, 1977. 23(3): p. 327-341.
112 Ryttsen, F., et al., Characterization of single-cell electroporation by using patch-clamp
and fluorescence microscopy. Biophys J, 2000. 79(4): p. 1993-2001.
113 Sack, M., et al., Operation of an Electroporation Device for Grape Mash. Ieee
Transactions on Plasma Science, 2010. 38(8): p. 1928-1934.
114 Saulis, G., S. Satkauskas, and R. Praneviciute, Determination of cell electroporation from
the release of intracellular potassium ions. Analytical Biochemistry, 2007. 360(2): p.
273-281.
115 Sersa, G., et al., Electrochemotherapy in treatment of tumours. Eur J Surg Oncol, 2008.
34(2): p. 232-40.
116 Sheng, Y., V. Mancino, and B. Birren, Transformation of Escherichia-Coli with Large
DNA-Molecules by Electroporation. Nucleic Acids Research, 1995. 23(11): p. 1990-
1996.
117 Smith, K.C., J.C. Neu, and W. Krassowska, Model of creation and evolution of stable
electropores for DNA delivery. Biophys J, 2004. 86(5): p. 2813-26.
118 Smith, K.C., et al., The spatially distributed dynamic transmembrane voltage of cells and
organelles due to 10-ns pulses: Meshed transport networks. Ieee Transactions on Plasma
Science, 2006. 34(4): p. 1394-1404.
119 Straatsma, T.P. and H.J.C. Berendsen, Free-Energy of Ionic Hydration - Analysis of a
Thermodynamic Integration Technique to Evaluate Free-Energy Differences by
Molecular-Dynamics Simulations. Journal of Chemical Physics, 1988. 89(9): p. 5876-
5886.
120 Stampfli, R. and M. Willi, Membrane potential of a Ranvier node measured after
electrical destruction of its membrane. Experientia, 1957. 13(7): p. 297-8.
121 Sugar, I.P. and E. Neumann, Stochastic-Model for Electric Field-Induced Membrane
Pores Electroporation. Biophysical Chemistry, 1984. 19(3): p. 211-225.
122 Sun, S., et al., Effects of deformability and thermal motion of lipid membrane on
electroporation: by molecular dynamics simulations. Biochem Biophys Res Commun,
2011. 404(2): p. 684-8.
123 Tarek, M., Membrane electroporation: a molecular dynamics simulation. Biophys J,
2005. 88(6): p. 4045-53.
131
124 Taupin, C., M. Dvolaitzky, and C. Sauterey, Osmotic-Pressure Induced Pores in
Phospholipid Vesicles. Biochemistry, 1975. 14(21): p. 4771-4775.
125 Teissie, J., 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, 2005. 1724(3): p. 270-280.
126 Tieleman, D.P. and H.J.C. Berendsen, Molecular dynamics simulations of a fully
hydrated dipalmitoyl phosphatidylcholine bilayer with different macroscopic boundary
conditions and parameters. Journal of Chemical Physics, 1996. 105(11): p. 4871-4880.
127 Tieleman, D.P., et al., Simulation of pore formation in lipid bilayers by mechanical stress
and electric fields. Journal of the American Chemical Society, 2003. 125(21): p. 6382-
6383.
128 Tieleman, D.P., The molecular basis of electroporation. BMC Biochem, 2004. 5: p. 10.
129 Tokman, M., et al., Electric Field-Driven Water Dipoles: Nanoscale Architecture of
Electroporation. Plos One, 2013. 8(4). p. e61111
130 Tolpekina, T.V., W.K. den Otter, and W.J. Briels, Nucleation free energy of pore
formation in an amphiphilic bilayer studied by molecular dynamics simulations. J Chem
Phys, 2004. 121(23): p. 12060-6.
131 Tounekti, O., et al., Bleomycin, an apoptosis-mimetic drug that induces two types of cell
death depending on the number of molecules internalized. Cancer Res, 1993. 53(22): p.
5462-9.
132 Tuckerman, M.E., et al., A Liouville-operator derived. measure-preserving integrator for
molecular dynamics simulations in the isothermal-isobaric ensemble. Journal of Physics
a-Mathematical and General, 2006. 39(19): p. 5629-5651.
133 Tsong, T.Y., On Electroporation of Cell-Membranes and Some Related Phenomena.
Bioelectrochemistry and Bioenergetics, 1990. 24(3): p. 271-295.
134 Vacha, R., et al., Effects of alkali cations and halide anions on the DOPC lipid
membrane. J Phys Chem A, 2009. 113(26): p. 7235-43.
135 Vacha, R., et al., Mechanism of interaction of monovalent ions with phosphatidylcholine
lipid membranes. J Phys Chem B, 2010. 114(29): p. 9504-9.
136 van Gunsteren, W.F. and M. Karplus, Effect of Constraints on the Dynamics of
Macromolecules. Macromolecules, 1982. 15(6): p. 1528-1544.
137 Vanýsek P. CRC Handbook of Chemistry and Physics 74
th
Edition, ed. Lide D.R., CRC
Press, Boca Raton, U.S.A., 1993, pp 5-90.
138 Vasilkoski, Z., et al., Membrane electroporation: The absolute rate equation and
nanosecond time scale pore creation. Physical Review E, 2006. 74(2).
132
139 Verhoven, B., R.A. Schlegel, and P. Williamson, Mechanisms of phosphatidylserine
exposure, a phagocyte recognition signal, on apoptotic T lymphocytes. J Exp Med, 1995.
182(5): p. 1597-601.
140 Vernier, P.T., et al., Calcium bursts induced by nanosecond electric pulses. Biochem
Biophys Res Commun, 2003. 310(2): p. 286-95.
141 Vernier, P.T., et al., Nanoelectropulse-induced phosphatidylserine translocation. Biophys
J, 2004. 86(6): p. 4040-8.
142 Vernier, P.T., et al., Nanosecond pulsed electric fields perturb membrane phospholipids
in T lymphoblasts. FEBS Lett, 2004. 572(1-3): p. 103-8.
143 Vernier, P.T., et al., Nanoelectropulse intracellular perturbation and
electropermeabilization technology: phospholipid translocation, calcium bursts,
chromatin rearrangement, cardiomyocyte activation, and tumor cell sensitivity. Conf
Proc IEEE Eng Med Biol Soc, 2005. 6: p. 5850-3.
144 Vernier, P.T., et al., Nanopore formation and phosphatidylserine externalization in a
phospholipid bilayer at high transmembrane potential. J Am Chem Soc, 2006. 128(19):
p. 6288-9.
145 Vernier, P.T. and M.J. Ziegler, Nanosecond field alignment of head group and water
dipoles in electroporating phospholipid bilayers. J Phys Chem B, 2007. 111(45): p.
12993-6.
146 Vernier, P.T., et al., Electroporating fields target oxidatively damaged areas in the cell
membrane. PLoS One, 2009. 4(11): p. e7966.
147 Vernier, P.T., M.J. Ziegler, and R. Dimova, Calcium binding and head group dipole
angle in phosphatidylserine-phosphatidylcholine bilayers. Langmuir, 2009. 25(2): p.
1020-7.
148 Vernier, P.T., Z.A. Levine, and M.A. Gundersen, Water Bridges in Electropermeabilized
Phospholipid Bilayers. Proceedings of the Ieee, 2013. 101(2): p. 494-504.
149 Weerasinghe, S. and P.E. Smith, A Kirkwood-Buff derived force field for sodium chloride
in water. Journal of Chemical Physics, 2003. 119(21): p. 11342-11349.
150 Weaver, J.C. and R.A. Mintzer, Decreased Bilayer Stability Due to Transmembrane
Potentials. Physics Letters A, 1981. 86(1): p. 57-59.
151 Weaver, J.C. and Y.A. Chizmadzhev, Theory of electroporation: A review.
Bioelectrochemistry and Bioenergetics, 1996. 41(2): p. 135-160.
152 Weaver, J.C., Electroporation of biological membranes from multicellular to nano
scales. Ieee Transactions on Dielectrics and Electrical Insulation, 2003. 10(5): p. 754-
768.
133
153 Wohlert, J., et al., Free energy of a trans-membrane pore calculated from atomistic
molecular dynamics simulations. Journal of Chemical Physics, 2006. 124(15).
154 Wong, T.K. and E. Neumann, Electric field mediated gene transfer. Biochem Biophys
Res Commun, 1982. 107(2): p. 584-7.
155 Wu, Y.H., et al., Versatile broadband electrode assembly for cell electroporation. Conf
Proc IEEE Eng Med Biol Soc, 2012. 2012: p. 2563-6.
156 Ziegler, M.J. and P.T. Vernier, Interface Water Dynamics and Porating Electric Fields
for Phospholipid Bilayers. Journal of Physical Chemistry B, 2008. 112(43): p. 13588-
13596.
157 Zimmermann, U., G. Pilwat, and F. Riemann, Dielectric breakdown of cell membranes.
Biophys J, 1974. 14(11): p. 881-99.
158 Zimmermann, U., et al., Effects of External Electrical Fields on Cell-Membranes.
Bioelectrochemistry and Bioenergetics, 1976. 3(1): p. 58-83.1. Vernier, P.T., et al.,
Nanosecond pulsed electric fields perturb membrane phospholipids in T lymphoblasts. FEBS Lett,
2004. 572(1-3): p. 103-8.
159 Zimmermann, U., Electric field-mediated fusion and related electrical phenomena.
Biochim Biophys Acta, 1982. 694(3): p. 227-77.
Abstract (if available)
Abstract
Electroporation provides a controllable method to introduce foreign substances into living cells. It is widely used by researchers in cell biology and the medical field to manipulate biological systems at the cellular level. For decades, electroporation has been studied extensively through experiments and theoretical models, and electroporation‐based technologies have been improved substantially with these efforts. One of the issues in utilizing electroporation is the lack of understanding in the phenomenon’s molecular mechanism and the microscopic details, mainly due to the difficulty in the direct experimental observation of the nanosecond‐scale electropore formation process and the nanometer‐sized electropore structure. To overcome this issue, Molecular Dynamics (MD) simulation has become one of the major tools to study electroporation at the microscopic level. ❧ Recent advancements of high performance computing, such as the increase in processing power, developments in algorithms and parallelization, have improved the efficiency of MD simulation substantially. Due to these advancements, MD simulation has become a popular tool for studying systems that are composed of biomolecules. For nearly a decade of effort, MD simulation revealed many different aspects of electroporation and it provided a molecular description of the process. Using MD simulation, we are able to observe the events during the electropore formation and annihilation, as well as the transport processes of molecules through the electropore. In addition, MD simulation provides a platform to study the molecular structure of electropore, and the associated energetic. ❧ My dissertation is organized as the follows: Chapter 1 provides the motivation of this research by discussing the applications of electroporation‐based technology, electroporation experiments, and the existing continuum model that describes electroporation. Chapter 2 introduces the MD formalism, models, and various algorithms used in our MD simulation. I will also discuss some of the previous MD studies of electroporation and their significances at the end of Chapter 2. In Chapter 3, we will examine the dynamics of water bridge during different stages of electroporation, and delineate the role of water molecules in electroporation by comparing the lipid bilayer system with an artificial water‐vacuum‐water system. In Chapter 4, we will examine the steady state of an electropore in lipid bilayer and evaluate the pore conductance of ions. The pore conductance values obtained from the simulations can be compared with those obtained by experiments. In Chapter 5, we will examine the effects that monovalent ions impose on lipid bilayer and electropore formation. We will also examine the pore conductance of ions under various ion concentrations, and the PS translocation process. At the end, I will summarize the findings in our research and provide a short outlook on MD simulation in the study of electroporation.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations
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 physics of membrane protein polyhedra
PDF
Shock-induced nanobubble collapse and its applications
PDF
Molecular dynamics and virial expansions for equilibrium thermodynamics in the solar interior
PDF
Imaging molecular transport across and nanomaterial interaction with lipid membranes
PDF
Electronic correlation effects in multi-band systems
PDF
Converging shocks in water and material effects
PDF
KcsA: mutational analysis of ion selectivity with molecular dynamics
PDF
Towards optimized dynamical error control and algorithms for quantum information processing
PDF
Neural network for molecular dynamics simulation and design of 2D materials
PDF
Modeling the hot-dense plasma of the solar interior in and out of thermal equilibrium
PDF
The effects of divalent counterions on the formation and stabilization of RNA tertiary structure
PDF
Capturing the complexity of ion channels: simulations of long-time dynamics, conformational changes and the effect of the membrane potential
PDF
Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations
PDF
Modeling nanodevices: from semiconductor heterostructures to Josephson junction arrays
PDF
The physics of pulsed streamer discharge in high pressure air and applications to engine technologies
PDF
Photoexcitation and nonradiative relaxation in molecular systems: methodology, optical properties, and dynamics
PDF
Modeling graphene: magnetic, transport and optical properties
Asset Metadata
Creator
Ho, Ming-Chak
(author)
Core Title
Molecular simulations of water and monovalent ion dynamics in the electroporation of phospholipid bilayers
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Physics
Publication Date
04/11/2014
Defense Date
02/20/2014
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
cell membrane,electropermeabilization,electropore,molecular dynamics,monovalent salt,OAI-PMH Harvest
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Gundersen, Martin A. (
committee chair
), Daeppen, Werner (
committee member
), Dappen, Werner (
committee member
), Haas, Stephan W. (
committee member
), Malmstadt, Noah (
committee member
), Vernier, Paul Thomas (
committee member
)
Creator Email
mingcho@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-376394
Unique identifier
UC11295860
Identifier
etd-HoMingChak-2340.pdf (filename),usctheses-c3-376394 (legacy record id)
Legacy Identifier
etd-HoMingChak-2340.pdf
Dmrecord
376394
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Ho, Ming-Chak
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
Tags
cell membrane
electropermeabilization
electropore
molecular dynamics
monovalent salt