SIMULATION STUDY OF LIPID BILAYER/NANOPARTICLE INTERACTIONS
by
Corey Evan Musolff

A DISSERTATION
Submitted
to Michigan State University
in partial fulfillment of the requirements
for the degree of
Physics - Doctor of Philosophy

2013

ABSTRACT
SIMULATION STUDY OF LIPID BILAYER/NANOPARTICLE INTERACTIONS
by
Corey Evan Musolff
A coarse-grained model is used to simulate lipid bilayers (LB) interacting with nanoparticles (NP).
Different equilibrium states are observed as the size of the NP is varied along with the interaction
strength between the NP and LB. Sufficient attraction causes the NP to become wrapped by the
LB. Large NP are able to form pores in the LB. These various outcomes are explained using a
continuum theory.
The same model is used to simulate pore healing in LB. It is observed that the area of the pore
decreases linearly with time as explained by Allen-Cahn theory. Bulk properties of the LB can be
used to predict the closing time of a pore as well as the amount of charge able to flow through the
pore while it is open.

To my wife and advisor for their support and infinite patience.

iii

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vi

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
CHAPTER 1

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 2 REVIEW . . . . . . . . . . .
2.1 Introduction . . . . . . . . . . . . . .
2.2 Experimental Work . . . . . . . . . .
2.3 Atomistic Models . . . . . . . . . . .
2.4 Coarse-Grained Models . . . . . . . .
2.4.1 Implicit Solvent . . . . . . . .
2.4.1.1 Tethered Models . .
2.4.1.2 Non-tethered Models
2.4.1.3 Cooke model . . . .
2.4.2 Explicit Solvent . . . . . . . .
2.4.2.1 MARTINI model . .
2.4.2.2 Klein model . . . .
2.5 Continuum Models . . . . . . . . . .
2.5.1 TGMB Theory . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

1

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

3
3
4
5
6
7
7
8
9
11
13
15
16
19

THERMODYNAMICS OF NANOPARTICLE/LIPID BILAYER INTERACTIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Continuum Model Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Fixed boundary LB with spherical NP . . . . . . . . . . . . . . . . . .
3.1.1.1 Detached liposome (sphere) . . . . . . . . . . . . . . . . . .
3.1.1.2 Embedded liposome (sphere) . . . . . . . . . . . . . . . . .
3.1.1.3 Partially wrapped spherical nanoparticle . . . . . . . . . . . .
3.1.1.4 Single layer micelle (sphere) states . . . . . . . . . . . . . .
3.1.2 Fixed boundary LB with infinite rod . . . . . . . . . . . . . . . . . . .
3.1.2.1 Detached liposome (rod) . . . . . . . . . . . . . . . . . . .
3.1.2.2 Embedded liposome (rod) . . . . . . . . . . . . . . . . . . .
3.1.2.3 Partially wrapped rod-shaped nanoparticle . . . . . . . . . .
3.1.2.4 Single layer micelle (rod) states . . . . . . . . . . . . . . . .
3.1.3 Fixed tension ensemble . . . . . . . . . . . . . . . . . . . . . . . . . .
LB properties from simulation . . . . . . . . . . . . . . . . . . . . . . . . . .
Phase diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Fixed boundaries with spherical NP . . . . . . . . . . . . . . . . . . .
3.3.2 Fixed boundaries with infinite rod NP . . . . . . . . . . . . . . . . . .
3.3.2.1 Small L . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2.2 Large L . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

21
21
22
24
25
26
28
31
32
33
34
35
36
38
41
41
45
48
51

CHAPTER 3
3.1

3.2
3.3

iv

3.4
3.5

Comparison with Ginzburg and Balijepalli
Comparison with MD simulations . . . .
3.5.1 Implementation Details . . . . . .
3.5.2 MD Results . . . . . . . . . . . .
3.5.2.1 Spherical NP . . . . . .
3.5.2.2 Rod-shaped NP . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

51
57
57
60
60
62

CHAPTER 4 LIPID BILAYER PORE DYNAMICS
4.1 Tension required for pore formation . . . . . .
4.1.1 Fracture mechanics . . . . . . . . . . .
4.1.2 Membrane rupturing . . . . . . . . . .
4.2 Dynamics of pore healing . . . . . . . . . . . .
4.2.1 Predicted time dependence . . . . . . .
4.2.2 Simulation results . . . . . . . . . . . .
4.2.2.1 Flat LB . . . . . . . . . . . .
4.2.2.2 Spherical vesicle . . . . . . .
4.2.2.3 Alternative fitting . . . . . . .
4.2.2.4 Ion transfer . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

69
69
69
70
71
72
77
77
80
83
85

CHAPTER 5

.
.
.
.
.
.

.
.
.
.
.
.

CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
Appendix A
Ginzburg and Balijepalli Equations . . . . . . . . . . . . . . . . . . . . 91
Appendix B
Nosé-Hoover Thermostat/Barostat . . . . . . . . . . . . . . . . . . . . 95
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

v

LIST OF TABLES

Table 4.1

Quadratic fit parameters for pore healing simulations . . . . . . . . . . . . . . . 80

Table 4.2

Exponential plus quadratic fit parameters for pore healing simulations . . . . . . 84

Table 4.3

Double exponential fit parameters for pore healing simulations . . . . . . . . . . 85

Table 4.4

Charge transfer through pores . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

vi

LIST OF FIGURES

Figure 2.1

Potential used in the Cooke model. The solid line indicates the attraction
between tail beads. The dashed line indicates the hardcore repulsion for interactions involving head beads. For interpretation of the references to color
in this and all other figures, the reader is referred to the electronic version of
this dissertation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

Figure 3.1

Hydrophilic spherical NP states . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Figure 3.2

Hydrophobic spherical NP states . . . . . . . . . . . . . . . . . . . . . . . . . 30

Figure 3.3

Hydrophilic rod-shaped NP states . . . . . . . . . . . . . . . . . . . . . . . . . 31

Figure 3.4

Hydrophobic rod-shaped NP states . . . . . . . . . . . . . . . . . . . . . . . . 36

Figure 3.5

Surface tension versus area in a flat LB. The green trace represents increasing
area while the black trace represents decreasing area. The red circle indicates
the critical point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

Figure 3.6

Phase diagram of spherical NP states. The x-axis shows the radius of the
nanoparticle. The y-axis shows the attraction between the nanoparticle and
lipid heads. Attraction increases with w. (A) Dissociated reference state, (D)
Embedded liposome without a pore, (E) Embedded liposome with a pore, (I)
Embedded single-layer micelle without a pore, (J) Embedded single-layer
micelle with a pore. Schematics and energies for each phase are given in
figures 3.1 and 3.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

Figure 3.7

Phase diagram of short rod-shaped NP states. The x-axis shows the radius of
the nanoparticle. The y-axis shows the attraction between the nanoparticle
and lipid heads. Attraction increases with w. (A) Dissociated reference state,
(D) Embedded liposome with and without a pore, (E) Embedded liposome
with a spanning pore, (I) Embedded single-layer micelle with and without a
pore, (J) Embedded single-layer micelle with a spanning pore. Schematics
and energies for each phase are given in figures 3.3 and 3.4. . . . . . . . . . . . 50

Figure 3.8

Phase diagram of long rod-shaped NP states. The x-axis shows the radius of
the nanoparticle. The y-axis shows the attraction between the nanoparticle
and lipid heads. Attraction increases with w. (A) Dissociated reference state,
(D) Embedded liposome with and without a pore, (E) Embedded liposome
with a spanning pore, (I) Embedded single-layer micelle with and without a
pore, (J) Embedded single-layer micelle with a spanning pore. Schematics
and energies for each phase are given in figures 3.3 and 3.4. . . . . . . . . . . . 52

vii

Figure 3.9

Phase diagram of spherical NP states. The x-axis shows the radius of the
nanoparticle. The y-axis shows the attraction between the nanoparticle and
lipid heads. (A) Dissociated reference state, (D) Embedded liposome without
a pore, (E) Embedded liposome with a pore, (I) Embedded single-layer micelle without a pore, (J) Embedded single-layer micelle with a pore. Schematics and energies for each phase are given in figures (3.1 and 3.2). Data points
indicate the results of Ginzburg and Balijepalli. (Blue circles): Embedded
micelles with no pores. (Red squares): Embedded liposomes with no pores.
(Green diamonds): Embedded liposomes with pores. . . . . . . . . . . . . . . 55

Figure 3.10 Area per lipid versus temperature from Cooke model simulation. The dark
lines represent running averages of the raw data which is displayed as lighter
dots. In the top trace the temperature is decreasing while in the bottom trace
it is increasing. The gel to liquid transition occurs around kB T = 0.98 ε while
heating. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Figure 3.11 Radial distribution function for a lipid bilayer at various temperatures during
the heating transition. The lack of long range order in kB T = 0.99 ε indicates
the liquid phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
Figure 3.12 Typical starting configuration of my spherical nanoparticle/lipid bilayer simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
Figure 3.13 Cross section of single-layer hybrid micelle. A hydrophobic nanoparticle is
embedded between two single layers of lipids. . . . . . . . . . . . . . . . . . . 61
Figure 3.14 Cross section of bilayer hybrid micelle. A hydrophilic nanoparticle is surrounded by a bilayer and embedded in the membrane. . . . . . . . . . . . . . . 62
Figure 3.15 Strongly hydrophilic nanoparticle is wrapped by a bilayer and causes a pore
in the membrane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
Figure 3.16 Phase diagram of spherical NP states. The x-axis shows the radius of the
nanoparticle. The y-axis shows the attraction between the nanoparticle and
lipid heads. (A) Dissociated reference state, (D) Embedded liposome without
a pore, (E) Embedded liposome with a pore, (I) Embedded single-layer micelle without a pore, (J) Embedded single-layer micelle with a pore. Schematics and energies for each phase are given in figures (3.1 and 3.2). Data points
indicate my molecular dynamics results. (Blue circles): Embedded micelles
with no pores. (Red squares): Embedded liposomes with no pores. (Green
diamonds): Embedded liposomes with pores. . . . . . . . . . . . . . . . . . . 64
Figure 3.17 Typical starting configuration of my rod-shaped nanoparticle/lipid bilayer
simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

viii

Figure 3.18 Phase diagram of rod-shaped NP states. The x-axis shows the radius of the
nanoparticle. The y-axis shows the attraction between the nanoparticle and
lipid heads. (D) Embedded liposome with a small pore, (E) Embedded liposome with a spanning pore, (I) Embedded single-layer micelle with a small
pore or no pore. Schematics and energies for each phase are given in figures
(3.3 and 3.4). Data points indicate my molecular dynamics results. (Blue circles): Embedded micelles with no pores. (Red squares): Embedded micelles
with small pores. (Green diamonds): Embedded liposomes with spanning pores. 66
Figure 3.19 Narrow, hydrophobic, rod-shaped nanoparticle embedded in a bilayer. It is
wrapped by a single layer of lipids and does not form a pore. . . . . . . . . . . 67
Figure 3.20 Thick, hydrophobic, rod-shaped nanoparticle embedded in a bilayer. It is
wrapped by a single layer of lipids and does form a pore. . . . . . . . . . . . . 67
Figure 3.21 Strongly hydrophilic, rod-shaped nanoparticle embedded in a bilayer. It is
wrapped by a double layer of lipids and has caused the membrane to rupture. . . 68
Figure 4.1

Snapshots of pore healing in a flat LB at three different times. . . . . . . . . . . 78

Figure 4.2

Pore area versus time in a flat lipid bilayer under zero tension. This is the
raw data from 20 different runs. . . . . . . . . . . . . . . . . . . . . . . . . . . 78

Figure 4.3

Averaged pore area versus time in a flat lipid bilayer under zero tension. The
data in figure 4.2 was averaged (red circles) and fit to a quadratic function
(solid black line) as well as a linear function (dashed black line). . . . . . . . . 79

Figure 4.4

Snapshots of pore healing in a spherical vesicle at three different times. . . . . . 81

Figure 4.5

Pore area versus time for a spherical vesicle. This is the raw data from 20
different runs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

Figure 4.6

Averaged pore area versus time for a spherical vesicle. The data in figure 4.5
was averaged and fit to both a linear and quadratic function. . . . . . . . . . . . 82

Figure 4.7

Averaged flat lipid bilayer pore data with alternative fits. The data in figure
4.2 was averaged (red circles) and fit to a an exponential plus a linear term
(solid green line) as well as a double exponential (dashed blue line). . . . . . . 84

Figure 4.8

Averaged vesicle pore data with alternative fits. The data in figure 4.5 was
averaged (red circles) and fit to an exponential plus a quadratic function (solid
green line) as well as a double exponential (dashed blue line). . . . . . . . . . . 85

Figure B.1 Histogram of kinetic energy for strong heat bath coupling (Q = 0.648, T damp =
0.005). The scaling variable s fluctuates faster than the system and causes a
bimodal, non-Gaussian distribution. . . . . . . . . . . . . . . . . . . . . . . . . 98

ix

Figure B.2 Temperature versus time for various T damp. . . . . . . . . . . . . . . . . . . . 99

x

Chapter 1
INTRODUCTION

I am interested in how nanomaterials interact with biological structures in the human body, specifically I am focused on interactions between nanoparticles (NP) and lipid bilayers (LB). NP can have
both beneficial and adverse effects on living cells. Nano-sized structures are special because they
are on the same size scale as lipids which are the building blocks of cell walls. Lipid molecules
are roughly a few nanometers long, so that the thickness of a lipid bilayer is about 5 nm. Therefore, certain NP are able to penetrate cells under the right conditions and can disrupt the cell wall
to varying degrees. NP are believed to be ideal drug delivery vehicles for this reason. If a drug
molecule is not able to penetrate a cell on its own it may be able to be attached to a NP which can
then be targeted to the cell. On the other hand certain NP can be toxic and are capable of rupturing
cells. This has encouraged the active study of nanotoxicity. However, even this behavior could be
useful if for instance cancerous cells are targeted.
NP/LB interactions have been looked at from many different angles. In chapter 2 I will give a
review of recent experimental and theoretical work in this field. I have chosen to use molecular dynamics (MD) simulations to study the problem. MD can easily look at a system over a wide range
of conditions without the usual experimental constraints. I would like to identify how NP/LB interactions are affected by fundamental properties such as the size, shape, and charge of NP. A deeper
understanding of these relationships would be useful in explaining effects that have already been
observed as well as predicting effects that have not yet been observed. New NP could be designed
with certain properties to have desirable interactions or to prevent undesirable interactions.
In chapter 3 I look at the thermodynamics of NP/LB interactions. I use a continuum model to
define the energy of various states. I then make phase diagrams showing the lowest energy states
depending on the properties of the system. I have also run MD simulations which agree reasonably
well with the predicted phase diagrams. Unlike the static states studied in chapter 3, chapter 4 looks

1

at the dynamics of pores in LB. I use a simple energy argument to predict the rate at which a pore
will heal under different conditions. These predictions are again compared to MD simulations. My
goal is to understand the general process of pore closing. In chapter 5 I will summarize my work
and draw conclusions.

2

Chapter 2
REVIEW

2.1

Introduction

Understanding the interactions between nanoparticles (NP) and lipid bilayers (LB) is necessary for
studying processes such as viral budding, endo/exocytosis, drug delivery, and harmful effects related to cytotoxicity. I am interested in the conditions under which nanoparticles become wrapped
and possibly rupture lipid bilayers.
This chapter covers the experimental and theoretical work of others in the field, though an
emphasis will be placed on molecular dynamics (MD) simulations. Experiments range from in
vitro studies on small patches of LB to in vivo drug studies on animals. Simulations tend to focus
on individual small scale interactions using a patch of LB or perhaps an entire vesicle or cell.
Simulations can take place over many length scales depending on the level of detail required,
ranging from the atomic scale to that of the entire cell (four orders of magnitude!). A single simulation or model can not hope to bridge that gap. At the small end of the spectrum are atomistic
or all-atom (AA) models. These of course have the advantage of including the most detail. We
are confident that the interactions between individual atoms are understood well enough to create
useful simulations. Unfortunately, these models are limited by computational constraints to looking at rather small systems (thousands of lipids) over short time scales (nanoseconds). Even with
optimistic predictions for future computing power, these models will likely remain inadequate for
looking at large patches of bilayers over long time scales (milliseconds).
At the other end of the spectrum are continuum models which do not have the same kind of size
restrictions. A lipid membrane can be modeled as a simple two dimensional sheet. While these
models can be used to study the large scale structural properties of the membrane, they are forced
to ignore some of the details responsible for complicated interactions.

3

Between these two extremes lie any number of intermediate, mesoscale models. These coarsegrained models are similar to the AA models in that they simulate individual particle interactions,
but instead of each particle representing a single atom these “beads” can represent several atoms.
Typically, one bead will represent a few atoms. Besides a reduction in the number of particles being
simulated, the interactions between them can also be simplified for instance from multi-body to
pairwise forces. Another enormous simplification can be made if the solvent or water molecules
surrounding the membrane are neglected. While water is in a sense responsible for lipid bilayer
formation via the hydrophobic effect, bilayer formation can still be achieved by using effective
forces. These are so called implicit solvent models while those that include water molecules are
explicit solvent models. Since there are many more water molecules surrounding the cell than the
number of lipids, neglecting the water allows much larger and longer simulations.
I will first review experiments in the field in section 2.2. Then I will review various theoretical
models ranging from atomistic to continuum in sections 2.3-2.5. Special emphasis will be given to
a few coarse-grained (CG) models because I have chosen to use a CG model for my simulations.

2.2

Experimental Work

Interactions between NP and LB require a mutual attraction. Hou et al. looked at the accumulation
of functionalized gold NP on LB. They observed that adsorption depended on the specific functionalization as well as the pH of the system [1]. After allowing some time for the NP to adhere,
they measured the total mass of the accumulated NP using inductively coupled plasma- optical
emission spectroscopy (ICP-OES). A wide range of NP sizes (1-100 nm) were tested with larger
NP producing a higher NP mass concentration on the membrane while smaller NP accumulated in
larger numbers producing a higher number concentration.
In addition to a NP simply adsorbing to a LB surface there are various ways in which it can
disrupt the surface, forming pores. Pore formation can be observed either directly by using an
imaging method such as atomic force microscopy (AFM) [2, 3, 4, 5, 6], or indirectly by measuring
the flow of current [7, 8] or dye molecules [5, 6, 7] across the membrane.
4

The method of measuring current flow across a membrane requires an applied voltage. While
the membrane is intact, its natural resistance will prevent current flow. Any spikes in current
indicate a pore. de Planque et al. used this technique to test the dependence of NP size and
functionalization on membrane disruption [8]. They found that small (50 nm) unfunctionalized
silica NP caused more disruption than when when they were amine-functionalized. Conversely,
larger (500 nm) NP caused no disruption unless they were functionalized. This implies a balance
between size and functionalization.
The effect of NP size alone was studied by Roiter et al. [2, 3]. They began by placing various
sized silica NP on a smooth silica substrate. Then a LB was formed over the top. AFM was used
to observe any pores formed in the membrane. NP smaller than 1.2 nm only slightly deformed the
surface and did not produce any pores. Larger NP (1.2-22 nm) did result in pore formation. Smooth
NP larger than 22 nm caused significant deformation of the membrane as it wrapped around the
NP, though they did not produce any pores. Large NP with sharp features or bumps with higher
curvature were also able to form pores. This intermediate range of NP sizes can be explained by
balancing the cost of bending the membrane against the favorable attraction between the NP and
LB. I will discuss this in detail in the next chapter.
Attraction between NP and LB is greatly affected by the charge of the NP. This was demonstrated in an experiment by Hong et al. [6] in which both supported LB and cells in vitro were
exposed to various polymer NP. The main determining factor for pore formation was the charge
concentration on the NP. Positively charged NP are known to be more attracted to LB which tend
to have a slight negative charge.
Medical research often focuses on either the toxicity of NP [9] or their use in drug delivery
[10, 11].

2.3

Atomistic Models

Atomistic models allow simulations on the smallest scale possible under classical mechanics.
Lyubartsev and Rabinovich published a recent review of current simulation techniques and fo5

cused on atomistic models [12]. In these models individual atoms interact through Lennard-Jones
and Coulomb forces. The three most popular models are GROMOS, CHARMM, and AMBER.
The GROMOS model was developed at the University of Groningen in 1978 by van Gunsteren
et al. and is an acronym for GROningen MOlecular Simulation [13, 14]. It makes the simplification of combining hydrogen atoms with the heavy atom they are bound to. The parameters have
been refined over the years. The newest set includes parameters for additional dihedral-angle and
torsional-angle terms [15]. Poger et al. presented a modified set optimized for dipalmitoylphosphatidylcholine (DPPC) lipids which produced lipid area and volume densities in agreement with
experiments [16].
CHARMM was created at Harvard and is an acronym for Chemistry at HARvard Macromolecular Mechanics [17, 18]. In contrast to GROMOS, CHARMM treats hydrogen atoms explicitly
and has its own set of Lennard-Jones, Coulomb, and bond energy parameters. A recent parameter
set was developed specifically for lipids [19].
AMBER was developed at the University of California, San Francisco and is an acronym for
Assisted Model Building with Energy Refinement [20]. It was originally designed to model nucleic
acids and proteins. A newer version called the general Amber force field (GAFF) was created to
simulate drug molecules [21]. It was recently shown to accurately simulate LB [22].
Tieleman used the GROMOS model to simulate electroporation in LB [23]. In the simulation
he applied a voltage across a pure LB surrounded by water. At an electric field strength of 0.5
V/nm he observed pores forming over the span of one nanosecond. Another study by Tarek used
CHARMM to look at similar processes including resealing of the pore after the electric field is
turned off [24].

2.4

Coarse-Grained Models

I am most interested in using coarse-grained (CG) models to study LB. They are able to strike a balance between the detailed realism of atomistic models and the large scale and speed of continuum
models. Several reviews have been written about current CG models [25, 26, 27, 28, 29, 30, 31].
6

The term coarse-grained is somewhat broad and can include a wide range of detail. Different CG
models have varying levels of complexity depending on the properties being studied. It is possible for very simple models to accurately simulate large scale structural properties. They are near
the continuum model end of the spectrum. Many levels of complexity can be added including
the number of beads, the types of interactions between beads, and the inclusion of solvent beads.
These details bring us closer to the atomistic end of the spectrum and are necessary for capturing
nanoscale processes.
Certain interesting NP/LB interactions take place on patches of LB which are hundreds of
nanometers wide (hundreds of thousands of lipids) and occur over millisecond time scales. These
scales are completely beyond the current capabilities of AA models. I will discuss the current CG
models available which I classify as either implicit or explicit solvent models. I will describe three
of the most popular models in more detail.

2.4.1

Implicit Solvent

The simplest and fastest models available employ an implicit solvent, meaning that the solvent
molecules, usually water, are not included in the simulation. These are a good place to start when
looking at LB systems if you are primarily interested in the large scale structural properties of
the LB itself. At first glance it seems impossible to model LB without including water. Water is
responsible for creating and keeping LB together through the hydrophobic effect. Lipids are not
held together by their own attraction to each other. Rather the hydrophobic tails are repelled by the
water and form micelles and bilayers in order to wall themselves off from it. To solve this problem,
the LB structure must be maintained by other means. The typical methods are either to tether the
head beads or to include an effective attractive force between tail beads.

2.4.1.1

Tethered Models

In tethered models the head beads are either tethered to their neighbors or else confined to move
along a two dimensional sheet with certain given properties [32]. Similar to continuum models,
7

tethered models are tuned to reproduce structural surface properties. While the level of detail is
greater than in continuum models and one is able to look at lipid diffusion and surface dynamics,
it is still fairly limited. LB formation and surface pore formation are not observable with such
models. Brannigan et al. discuss recent tethered model applications in their review [26].

2.4.1.2

Non-tethered Models

Without tethering or solvent molecules, effective forces are needed to hold the LB together. If an
attractive force is used between tail beads, lipids will naturally aggregate into micelles and bilayers.
These models can be used to study such formations. Effective force models can capture much more
dynamic behavior, but it can be tricky to tune the attractive forces just right. Early such models
could not produce a liquid phase. The Lennard-Jones type pair forces used proved to be too short
ranged. Forces strong enough to hold the lipids together tend to form solid structures, but fall apart
at higher temperatures. LB in their natural state should have a robust liquid phase. It was shown
by Drouffe in the early 1990’s that multibody forces could be used to produce a liquid phase [33].
In the early 2000’s a few groups proposed pair potential models also capable of producing a liquid
phase. Pair potentials are much more desirable because they are computationally simpler and can
be easily implemented by popular molecular dynamics software packages.
The Farago model uses rigid lipids composed of three beads and a different Lennard-Jones
potential for each combination of bead type [34]. Farago was able to produce a liquid-gel phase
transition and to look at membrane elasticity and pore formation. The bending rigidity found using
this model was significantly higher than was expected for a single component membrane. Other
downsides of the model are that it can not spontaneously form bilayers and that the number of
potential combinations is somewhat cumbersome.
Brannigan and Brown proposed a model using spherocylinders to represent each lipid [35]. The
pair potential between lipids depended on their distances from each other as well as their relative
angles. Their model could reproduce the same behavior as the Farago model as well as being able
to self assemble.

8

2.4.1.3

Cooke model

Around the same time, Cooke et al. [36, 37] proposed their own model. Each lipid in the Cooke
model consists of three beads, one hydrophilic head bead and two hydrophobic tail beads. Interactions are limited to simple pair potentials which I will describe presently. All non-bonded beads
experience a hardcore repulsion via a Weeks-Chandler-Andersen potential

Vrep (r; b) =


 4ε


b
r

12

6
− b +1 ,
r
4




r ≤ rc
(2.1)
r > rc

0,

which is simply a shifted Lennard-Jones potential cut off at the minimum rc = 21/6 b. The
length scale in the model is set by the diameter of a tail bead σ = btail,tail

0.9 nm. The authors

claim that the head bead should be slightly smaller in order to ensure the correct lipid shape.
Therefore they use bhead,head = bhead,tail = 0.95σ

0.86 nm. I have run most of my simulations

at kB T = 0.9ε. Assuming this is near room temperature, a rough energy scale is found to be ε ≈
2-3 kJ/mol.
Individual lipids are held together by strong FENE (finitely extensible nonlinear elastic) bonds
of the form
r 2
1
2
Vbond (r) = − kbond r∞ log 1 −
2
r∞

(2.2)

where kbond = 30ε/σ 2 is the stiffness and r∞ = 1.5σ is the divergence length. Lipids are kept
straight by a harmonic potential between the two end beads
1
Vbend (r) = kbend (r − 4σ )2
2

(2.3)

which has a rest length of 4σ . kbend = 10ε/σ 2 is the bending stiffness of the lipid. In addition
to the aforementioned hardcore repulsion experienced by all non-bonded beads, tail beads feel a
broad attraction to each other of the form

9

Vattr (r) =













−ε,
−ε cos2

π(r−rc )
2wc

r < rc
,

rc ≤ r ≤ rc + wc

(2.4)

r > rc + wc

0,

This is a square well potential of depth ε below rc which then tapers to zero between rc and
rc + wc . The total inter-lipid bead potentials are shown in figure 2.1. The key to the Cooke model
is the width of the potential, wc , which determines the properties of the bilayer such as the melting
temperature. It is the only tunable parameter. The reason why previous pair potential models failed
to produce a liquid phase was that the potential wells were too narrow. Cooke et al. argued that
the form of the potential is unimportant as long as it is wide [37]. To demonstrate this they tried
two different forms of the potential. One simply inserted a flat extension from the minimum of the
normal LJ potential while the other slowly tapered the energy from the minimum to zero following
the smooth cosine squared shape shown in the figure. They showed that both potentials behaved
similarly and were both able to produce a liquid phase.
An estimated timescale can be calculated based on the other parameters σ

m/ε ≈ 10 ps.

However, the effective timescale can be much longer for coarse-grained models with relatively few
degrees of freedom [38, 39]. Cooke and Deserno estimated the effective timescale for the Cooke
model to be τ ≈ 10 ns [37, 40] by comparing their simulated diffusion constant with experiment
[41].
The Cooke model was originally used to study several LB properties including lipid density,
diffusion, and bending energy. The authors looked at the effects of stress on LB observing pore
formation and buckling which was previously impossible using tethered models [36, 37]. Deserno
has further used the model to simulate interactions with proteins [40]. Ruiz-Herrero et al. have used
the Cooke model to simulate the wrapping of NP by flat LB under zero tension. The Cooke model
appears to contain enough detail to simulate realistic nanoparticle interactions. At the same time
its simplicity and coarseness allow simulations over large size (100 nm) and time (milliseconds)
scales. That is why I have chosen to use this model in my own simulations. I will discuss my

10

Figure 2.1: Potential used in the Cooke model. The solid line indicates the attraction between tail
beads. The dashed line indicates the hardcore repulsion for interactions involving head beads. For
interpretation of the references to color in this and all other figures, the reader is referred to the
electronic version of this dissertation.

implementation of the Cooke model in chapter 3.

2.4.2

Explicit Solvent

As the name suggests, explicit solvent models differ in that they include explicit solvent beads,
usually water. The size of the solvent beads are generally chosen to be similar to that of the lipid
beads. For that reason water beads tend to represent multiple water molecules. I will describe
some of the current models below.
One of the first explicit solvent models was developed by Goetz and Lipowsky [42]. They used
three types of beads. Lipids consisted of one or two tail chains containing four hydrophobic beads
each. Head groups were made up of either one or three hydrophilic beads depending on whether
11

the lipid was single or double tailed. The solvent beads were treated the same as hydrophilic head
beads. All three beads had the same mass. Hydrophilic and hydrophobic beads interacted through
a strictly repulsive soft-core (1/r9 ) potential. Hydrophilic beads, both water and head beads, were
attracted by a Lennard-Jones potential. Lipids were held together by harmonic potentials. They
looked at both rigid and floppy lipids by including a bond bending potential. Their model was able
to self assemble fluid bilayers and micelles and they were able to measure the properties of those
structures. A very similar model was used by Stevens et al. to study liposome fusion [43].
Loison et al. used a similar model to study pore formation in LB [44, 45]. The Loison model
represents each lipid with two hydrophilic head beads and two hydrophobic tail beads. Solvent
beads consist of a single bead which represents approximately three water molecules. Individual beads are bound together by strong FENE bonds. All beads are repelled at short range by a
Lennard-Jones like potential. Beyond the Lennard-Jones minimum, beads from different lipids
and solvent beads are attracted by the potential

U(r) =

φ
cos αr2 + β − 1
2

(2.5)

α and β determine the cutoff radius at which the potential becomes zero and above which the
potential remains zero. The potential depth, φ , determines the strength of the potential. It is the
same for all pairs of beads of the same type as well as between head and solvent beads. Tail beads
are not attracted to head or solvent beads, making the corresponding φ value zero. Apart from the
addition of water molecules, this looks very similar to the Cooke model [36, 37]. The shape of
the attractive part of the potential is slightly different. The more gradual increase approaching the
cutoff radius does affect the stability of the membrane. Loison et al. found that the bilayer state
was unstable in the presence of large amounts of solvent [45].
The explicit solvent models listed thus far share the fact that they use simple effective forces
which are are tuned to recreate large scale behavior. They also tend to use very few types of
beads. These characteristics make them similar to the implicit solvent models. Also, the models
mentioned thus far have not used any electrostatic forces. The next group of explicit solvent

12

models differ in that they try to simulate atomic scale interactions more accurately. They use
a larger variety of beads as well as a variety of potentials to simulate the different interactions
between species. Electrostatic forces are also included. The added details make these models more
complicated and less general, but arguably more realistic depending on the system. The two most
popular models in this regime are the MARTINI [46, 47] and Klein [48, 49, 25] models. Each uses
beads consisting of four heavy atoms that interact through Lennard-Jones and electrostatic forces.
I will explain both of these models in more detail below.
Dissipative particle dynamics (DPD) was developed as an alternative to traditional molecular
dynamics (MD) by Hoogerbrugge and Koelman [50, 51]. Groot and Rabone adapted this method
to simulate LB [52]. Their model uses beads consisting of three heavy atoms which is comparable
to the MARTINI and Klein models. Each bead experiences a conservative force, a random force,
and a drag force. The conservative force defines a cutoff radius, below which beads feel a repulsion
that increases linearly as they get closer. This extremely soft repulsion greatly increases diffusion.
The time step used in MD is limited by the usual hard-core repulsion between beads. Trajectories
must be updated frequently in MD to make sure beads don’t wander too close together. Without
having to worry about the infinite energy of overlapping beads, DPD can use a time step 1000 times
larger than in MD. The coarse-graining gives an additional effective speedup, allowing Groot and
Rabone to use a time step of about 5 ps.

2.4.2.1

MARTINI model

Marrink et al. developed a popular coarse grained model called MARTINI [46, 47]. MARTINI
classifies 18 different building blocks for lipids and solvents. Each bead represents about four
heavy atoms. A typical lipid is therefore composed of about twelve beads. The main distinction between beads is polarity. Beads are classified as being polar, nonpolar, apolar, or charged.
Nonpolar and charged beads are subdivided according to their affinity for hydrogen bonding. The
subtypes are donor, acceptor, both, or none. Polar and apolar beads are assigned a number from
one to five indicating its degree of polarity.

13

Lipids are held together by a weak harmonic potential between nearest neighbors
1
Vbond (R) = Kbond (R − Rbond )2
2

(2.6)

where Kbond = 1250 kJ mol−1 nm−2 and Rbond = 0.47 nm for every combination of beads.
Next nearest neighbors experience a bending potential
1
Vangle (θ ) = Kangle [cos(θ ) − cos(θ0 )]2
2

(2.7)

where Kangle = 25 kJ mol−1 and θ0 = 180◦ for most types of bonds. Besides nearest neighbors,
every combination of beads experiences one of ten different Lennard-Jones potentials.

ULJ (r) = 4εi j

σi j 6
σi j 12
−
r
r

(2.8)

The interaction strengths, εi j , vary from 2.0-5.6 kJ/mol. The closest approach distance, σi j ,
equals 0.47 nm for all but a couple interactions. The final contribution is the Coulomb potential
for charged beads.

Uel (r) =

qi q j
4πε0 εr r

(2.9)

Screening is enforced by setting εr = 15. The original model had the problem of water beads
freezing near room temperature. The newest version inhibits nucleation by replacing 10% of the
water beads with “antifreeze” beads that have an effective radius of 0.57 nm instead of the typical
0.47 nm.
Model parameters have been tuned to reproduce experimental results. The use of 18 different
building blocks allows the MARTINI model to capture a wide variety of interactions while keeping
a manageable number of parameters. The effective reduction in the number of atoms along with
the use of simple pair interactions makes it possible to simulate larger systems over longer times
compared to all-atom models. Simulations can be run with a time step of about 40 fs. The effective
timescale is about four times longer than that due to the lack of friction.

14

The MARTINI model has been used to simulate NP/LB interactions. A molecular dynamics
study by Lee and Larson found membrane pore formation caused by small dendrimers (R ≈ 2 nm)
[53]. Nanoparticle shape is another important factor. Nangia and Sureshkumar performed MD
simulations of NP with various shapes and charge densities interacting with LB [54]. They found
that shapes with high aspect ratios like rice were able to translocate more easily than spheres.

2.4.2.2

Klein model

Another earlier model was developed by Shelley et al. to specifically simulate dimyristoylphosphatidylcholine (DMPC) lipids [48, 49, 25]. Water beads represent four water molecules and
interact with each other via a Lennard-Jones 6-4 potential
σαβ 6

15
V (ri j ) = εαβ
4

ri j

−

σαβ 4
ri j

(2.10)

The parameters were calibrated to reproduce the correct melting temperature of DMPC. The
lipid tails consist of two different beads, (CH2 )3 and (CH2 )2 − CH3 . The heads are made up of
choline, phosphate, glycerol backbone, and ester group beads. Adjacent beads are bound together
by a harmonic potential

V (ri j ) =

kb;αβ
2

(ri j − r0;αβ )2

(2.11)

The bond angles among tail beads are governed by a cosine potential

V (θi jk ) = kθ ;αβ γ 1 − cos(π − θi jk )

(2.12)

These parameters are tuned by comparing with atomistic simulations. Bond angles involving
head beads are governed by a quartic bond angle potential

VQ (θi jk ) = AB2 /4
A = Kθ ;αβ γ /(2d0 )
15

(2.13)

B = d0 − (θi jk − π)2
d0 = (θ0;αβ γ − π)2
All non-bonded interactions other than water-water are given by a Lennard-Jones 9-6 potential
15
V (ri j ) = εαβ
4

σαβ 9
ri j

−

σαβ 6
ri j

(2.14)

These parameters were determined by matching densities from experiments and radial distribution functions from atomistic simulations. Finally, the choline and phosphate beads were assigned
respective charges of +e and -e and a dielectric constant of 78.
The Klein model is a niche model designed to study a very specific type of lipid (DMPC).
Its parameters have been finely tuned to accurately represent that system. The coarseness of the
model is about the same the MARTINI model with an entire lipid consisting of thirteen beads. The
bond length and bond angle potentials must be calculated every femtosecond while the non-bonded
interactions less than 11 Å are updated every 2 f s and longer non-bonded interactions are updated
every 40 f s.

2.5

Continuum Models

Continuum models are useful tools for studying the shape and large scale dynamics of LB [55, 56,
26]. LB lend themselves to continuum models due to their relatively simple morphology and the
large number of lipids usually involved. Bulk parameters such as bending energy and surface tension act as inputs to the models. Some recent reviews of continuum models are given by Brannigan
et al. [26] and Brown [57].
Canham was one of the first to use such a model to understand the shapes of red blood cells
[55]. He was able to produce a variety of shapes by minimizing bending energy alone. Shortly
after, Helfrich developed a similar model which has become the standard [56]. He considered three
types of strain: stretching, tilt, and curvature. He found that the curvature modulus is much smaller
than the stretching or tilt moduli. That means that the surface area and tilt of the LB tend to remain
16

fixed and that the shape of the LB is determined mostly by the curvature. He also argued that tilt
is less important than stretching. The Canham-Helfrich (C-H) Hamiltonian of the LB, neglecting
tilt, is given by [27, 55, 56, 58]

H =

d2A γ +

κ
¯
(C +C2 −C0 )2 + κC1C2
2 1

(2.15)

¯
where γ is the lateral membrane tension. κ and κ are the bending rigidity and saddle-splay
modulus. C1 and C2 are the local curvatures. For a sphere of radius R, C1 = C2 = 1/R. C0 is the
spontaneous curvature which is zero as long as the LB composition is the same on both sides. I will
make that assumption in later calculations. For a LB of fixed topology the last term in the integral
will contribute only a constant and can therefore be neglected. Yoo et al. calculated the energy
of bending that occurs around the edge of a pore to find the equilibrium shape and compared with
MD simulations [59].
In addition to bending energy there is a stretching energy associated with stretching or compressing bonds in the membrane. To lowest order it is assumed that the stretching energy is proportional to the square of the change in area. If the membrane is interacting with another surface
such as a NP, there will also be an adhesion energy proportional to the contact area.
Deserno has used continuum models to find equilibrium configurations involving NP and LB
[60, 61, 62, 40]. A paper by Deserno and Gelbart [60] looked at the interactions between spherical
colloids and spherical vesicles. They considered the following factors: the sizes of the colloid and
vesicle, the adhesion energy between them, and also the surface tension and bending energy of the
vesicle. The colloid was treated as a hard sphere while the vesicle was able to deform while binding
with the colloid. The energy contributions they considered were the adhesion energy between the
colloid and vesicle (Ead ), the tension energy of the stretched vesicle (Eten ), and the curvature
energy for bending the vesicle (Ec ).

Etotal = Ead + Eten + Ec

(2.16)

Ead = −kad Aad

(2.17)

17

Eten =

kten [max(A; A1 ) − A1 ]2
2
A1

(2.18)

1
kc (κ1 + κ2 )2 dA
2

(2.19)

Ec =

where kad is the adhesion constant and Aad is the area of contact. kten is the lateral compressibility, A1 is the unstretched vesicle surface area, and A is equilibrium stretched area. The max
function in Eten ensures that the vesicle will only stretch and can not be compressed. Ec is simply the Helfrich energy defined above, neglecting the Gaussian and spontaneous curvatures. The
energy gain through adhesion must be balanced against the deformation costs.
For a given set of vesicle and colloid sizes and membrane parameters they were able to minimize the energy to find how far the colloid penetrated as well as the curvature of the vesicle near
separation. They found that vesicles greater than a certain minimum (R ∼ 300 nm) were able to
completely envelop colloids within a certain range of sizes. Colloids that were too large were
only partially wrapped. Vesicles below the critical size could only achieve partial wrapping regardless of the size of the colloid. Very small colloids were not able to penetrate the vesicle at all
and remained unbound. In that case the surface tension cost overwhelmed the potential gain from
adhesion.
Other studies have looked at the interactions of NP with flat membranes and found similar
conditions for partial or complete wrapping [62, 61, 63]. In a more recent paper Deserno calculated
the tension required to form a pore in a membrane in the absence of a NP which agreed with coarsegrained MD simulations [40].
Van Lehn and Alexander-Katz also used the C-H Hamiltonian to model NP interacting with
a LB under zero tension [64]. Their NP consisted of two different phases which were either hydrophilic or hydrophobic and were free to mix. This caused different adhesion energies depending
on which part of the NP was in contact with either the lipid head or tail region. In addition to the
bending and adhesion terms, they included the mixing energy for the NP surface. They predicted
varying degrees of penetration depending on the adhesion strengths and ratio or phases on the NP
which they compared to their own MD simulations.

18

2.5.1

TGMB Theory

Although different from the continuum models above, Ginzburg and Balijepalli used a combination
of self-consistent field theory and density functional theory to model the interaction between NP
and LB [65]. Their work is significant because it proposes a phase diagram of NP/LB states. They
predict outcomes based on NP properties which can be compared to other theories and experiments.
I have run simulations closely following the conditions used by Ginzburg and Balijepalli. I will
present a comparison in chapter 3.
The technique was developed earlier by Thompson, Ginzburg, Matsen, and Balazs (TGMB) to
model the interaction between NP and block polymers [66, 67]. In the Ginzburg and Balijepalli
adaptation, diblock polymers are replaced by lipid diblocks (D) which have hydrophilic (H) and
lipophilic (L) parts. The other two components are the water solvent (W) and the nanoparticle (P).
The total free energy can be written as
FND
= f1 + f2 + f3 + f4
kB T ρ0V
ϕ
QP
f1 = P ln
α
V

f2 =

1
V

− ϕD ln

QD
V

− ND ϕW

(2.20)
QW
V

(2.21)


dr  ∑ χαβ ND φα (r)φβ (r) − ξ (r) 1 − ∑ φα (r) 
α,β


f3 =

1
V

(2.22)

α

dr −



∑

wα (r)φα (r) − wP (r)ρP (r)

(2.23)

α(=P)

f4 =

1
V

¯
dr ρP (r)ΨCS φP (r)

(2.24)

The fist term, f1 (2.21), contains the entropic free energies. The second term, f2 (2.22), contains the local interactions determined by Flory-Huggins parameters [68, 69]. I will explain each
of the contributions in detail in appendix 5.
Ginzburg and Balijepalli solved the above field equations on a two dimensional lattice 120x120.
The size of each lattice site was 0.36 nm on a side making the total lattice roughly 43x43 nm2 .
They used nanoparticles of radii (R p = 1.6, 2.0, 2.4, 2.8, and 3.2 nm). The nanoparticles were
19

made stationary and their densities were defined by sharp Gaussian functions. The Flory-Huggins
parameters for the interactions between nanoparticle (P), water (W), hydrophilic (H), and lipophilic
(L) blocks were as follows: χHL = χW L = 1.0 meaning that the lipophilic blocks will separate from
the water and hydrophilic blocks. χW H = χPL = 0.0 meaning that water and hydrophilic blocks as
well as nanoparticle and lipophilic blocks are as likely to mix as not. χPW = 1.0 while χPH varies
from +1.0 (repulsive) to -3.0 (strongly attractive). They argue that χPH = +1.0 corresponds to uncharged or hydrophobic nanoparticles while χPH = −3.0 corresponds to highly charged nanoparticles. It is known that nanoparticles are attracted more strongly to bilayers when charged [70, 71].
The values of χPL and χPW make it more favorable for the nanoparticle to become absorbed into
the lipophilic region than to remain separated from the membrane. The other possible configuration is for the nanoparticle to become surrounded by the bilayer in contact with the hydrophilic
blocks. The most favorable configuration is determined by the size of the nanoparticle and the
interaction between nanoparticle and hydrophilic block.
They found that electrically neutral nanoparticles were absorbed into the center of the bilayer
for all particle sizes tested. As the charge density was increased (χPH decreased) the nanoparticlehydrophilic block interactions became more favorable. Charged nanoparticles were coated by a
bilayer of lipids instead of the monolayer needed for the uncharged nanoparticles. That caused
twice as many lipids to be pulled away from the membrane surface resulting in thinning and even
rupture. Small charged nanoparticles formed micelles while remaining embedded in the membrane. Charged nanoparticles above a certain critical size, in addition to forming micelles, were
also found to form pores in the membrane. While the study was done over a narrow range of
nanoparticle sizes they assume that the trend of large, highly charged nanoparticles causing pores
would continue. I will present the details of their results in chapter 3 as they pertain to my results.
Ginzburg and Balijepalli acknowledged that their model only deals with the thermodynamics of
the system and that a dynamical model is needed to understand the details of the pore formation
such as energy barriers and time scales.

20

Chapter 3
THERMODYNAMICS OF NANOPARTICLE/LIPID BILAYER INTERACTIONS

Interactions between nanoparticles (NP) and lipid bilayers (LB) can be understood in terms of the
membrane properties and the attraction between the NP and LB. In this chapter I will outline a
simple continuum model and define the free energies for various NP/LB states. I consider both
spherical and rod-shaped NP. By comparing those energies I will make a phase diagram to identify
when NP become wrapped and when pore formation becomes favorable. I will then compare my
phase diagram with the DFT results of Ginzburg and Balijepalli [65]. As I mentioned in section
2.5.1, I have chosen to compare with their work because they have also proposed a phase diagram.
Finally, I will present MD simulations I have run which I will compare with my predicted phase
diagram.

3.1

Continuum Model Theory

The energy of a membrane can be characterized in terms of large scale deformations like bending
and stretching. In this section I will present the energies for various NP/LB states. All energies are
defined with respect to an initial reference state in which the NP and LB are separated. I will first
define the energies and then compare them in order to make a phase diagram.
The choice of ensemble is an important factor in determining the wrapping behavior. In the case
of an isolated cell or vesicle interacting with small NP, a fixed area ensemble may be appropriate.
While there can be significant deformation of the membrane near the NP, the large scale dimensions
of the LB remain relatively unchanged. This assumes that there are insufficient free lipids to be
added to the LB or that the process occurs much more slowly than the NP wrapping. This means
that wrapping will increase tension in the LB and that exceeding a critical tension could rupture the
LB. I calculate free energies for the various states in this ensemble and later present corresponding
fixed boundary MD simulations.
21

Alternatively, it may be more appropriate to fix the lateral tension of the LB. Cells are able
to regulate a small but nonzero tension by adding or removing lipids. While wrapping is still
possible, stable pore formation is impossible under low tension. One could imagine a very large
or highly charged NP causing severe stretching and inhibiting the cell’s tension regulation. That
would cause a spike in tension and temporary pore formation. Ruiz-Herrero et al. observed such
transient pores in their simulations [63]. One must be very careful in defining the barostat in
MD simulations. The barostat controls how quickly the pressure can relax and can therefore be
responsible for pressure/tension spikes.
In this section I will consider a fixed area and fixed tension LB interacting with a spherical NP.
I will also discuss the infinite rod NP case as it appears to be similar to a 2D system.

3.1.1

Fixed boundary LB with spherical NP

Here I list the possible states for a spherical NP interacting with a LB with fixed boundaries. I
first consider hydrophilic NP meaning that the NP will adhere to the lipid heads, opposed to the
tails. I will separately consider hydrophobic NP which prefer to be in contact with the lipid tails.
For each state I define the free energy with respect to the reference state in which the NP and LB
are dissociated. Figures 3.1 and 3.2 show schematics for each of the states along with their free
energies.

22

Figure 3.1:
Hydrophilic spherical NP states

A Reference state defined as zero energy

B Detached liposome without a pore

C Detached liposome with a pore

F1a =

K
2

2
4πR2 −πR2
n
p
A0

+ 8πκ + 2πR p γ − 4πR2 w
n

23

Figure 3.1 (cont’d)

D Embedded liposome without a pore

E Embedded liposome with a pore

F2a =

K
2

2
3πR2 −πR2
n
p
A0

+ 8πκ + 2πR p γ − 4πR2 w
n

F Partially wrapped liposome

(4πR2 f 2 )2
n
F3a = K
+ 8πκ f − 4πR2 f w
n
2
A
0

3.1.1.1

Detached liposome (sphere)

The first nonzero energy state involves the nanoparticle being completely wrapped and detached
from the bilayer. I call this state “detached”. Its energy is
2
2
K 4πRn − πR p
F1a =
2
A0

2

+ 8πκ + 2πR p γ − 4πR2 w
n

(3.1)

The first term is the stretching energy of the bilayer. This is an harmonic approximation appropriate for relatively small changes to the membrane area [34, 72, 37, 40]. K is the stretching
24

modulus and Rn is the radius of the NP. A0 is the starting unstretched area of the bilayer and is a
constant in this ensemble. The squared term in parentheses is the change in area, in this case the
surface area of the NP minus the area of a possible pore of radius R p .
The second term is the bending energy of the bilayer around the nanoparticle and is independent
of the nanoparticle size in the case of total spherical coverage [56]. κ is the bending modulus of
the LB. The third term is the line energy around the edge of the pore. The cost comes from lipid
heads having to curve around the edge of the pore to protect the hydrophobic tails. The constant,
γ, is a material property and is directly proportional to K. In section 3.2 I will calculate values for
K, A0 , κ, and γ based on my own MD simulations.
The last term is the adhesion energy proportional to the nanoparticle’s surface area and is the
only favorable term. The constant, w, is a parameter in my simulations on the order of 20 ε/σ 2 .
While I treat w as a constant, it should depend on the lipid concentration in contact with the NP.
A weakly attractive NP will allow thinning of the LB and result in less adhesion energy. On the
other hand, a strongly attractive NP could cause a tight aggregation of lipids in contact with it and
result in a higher adhesion energy. Assuming a constant w is reasonable for small changes in lipid
concentration. This issue will come up again later when analyzing my phase diagram.
This energy applies to states B and C in figure 3.1 depending on whether or not there is a pore.

3.1.1.2

Embedded liposome (sphere)

The second state consists of a completely covered nanoparticle which is embedded in the bilayer.
The energy is
2
2
K 3πRn − πR p
F2a =
2
A0

2

+ 8πκ + 2πR p γ − 4πR2 w
n

(3.2)

and corresponds to states D and E in figure 3.1 depending on whether or not there is a pore.
The change in area in the first term is the surface area of the NP minus the cross sectional area
of the NP displacing the membrane as well as the area of the pore. The only difference from the

25

detached state is a reduction in stretching energy. For that reason, the embedded states will always
be more favorable.
Equation (3.2) assumes that the pore is physically separated from the embedded liposome. It
is possible and in some cases preferable for the pore to be adjacent to the liposome. By sharing
an edge, the line energy around the pore can be reduced. This decrease in energy is somewhat
mitigated by a decrease in pore area which increases the stretching energy. After calculating the
changes to the stretching and line energy terms, I found that the effect on the total energy is negligible.
An additional energy contribution I have neglected is the LB curvature around the interface of
the embedded liposome. For small NP, whose size is comparable to the LB thickness, the interface
curvature will be negligible. For larger NP, the curvature will increase until the radius of curvature
around the interface is similar to that around the edge of a pore. In this large NP limit, an additional
line energy term should be included in which γ increases with the NP radius. I believe equation
(3.2) is appropriate for describing my simulations which involve small NP.

3.1.1.3

Partially wrapped spherical nanoparticle

Another possible state is the partially wrapped NP given by state F in figure 3.1. If f represents
the fraction of the NP covered, its energy is

F3a =

K (4πR2 f 2 )2
n
+ 8πκ f − 4πR2 f w
n
2
A0

(3.3)

The bending and adhesion terms are simply proportional to f . This equation assumes that
the NP is less than half wrapped (0 ≤ f ≤ 1/2). Further wrapping requires a more complicated
form. The stretching term contains the surface area of a spherical cap minus the excluded circular
intersection area. The derivative with respect to f must be taken to find the optimal wrapping
fraction.

26

dF3a
= 4π
df

8πR4 f 3 K
n
+ 2κ − R2 w
n
A0

(3.4)

Initial wrapping depends on the derivative as f approaches zero. Setting the limit of the derivative as f approaches zero equal to zero gives the following condition.

dF3a
=0
d f f →0
⇒ 4πR2 w = 8πκ
n

(3.5)

The two remaining terms are the adhesion and bending energies. If the adhesion is stronger,
then at least partial wrapping will occur. For a given value of w a critical nanoparticle radius can
be found.

Rn1 =

2κ
w

(3.6)

The condition for at least a metastable half wrapped state can be found by taking the limit of
the derivative as f approaches 1/2.

dF3a
=0
d f f →1/2
(2πR2 )2
n
⇒ 4πR2 w = 8πκ + K
n
A0

(3.7)

In this case the adhesion not only has to overcome the cost of bending but also some stretching.
Again, assuming a constant w, the critical minimum radius would be

Rn2 =

A0
2πK

w−

w2 −

8πκK
A0

(3.8)

In order to see if total wrapping is favorable, let’s compare the half wrapped state and the
embedded state with no pore

27

F2a |R p →0 − F3a | f →1/2 = 2π

2πR4 K
n
+ 2κ − R2 w
n
A0

(3.9)

The condition for complete wrapping can be determined by finding when equation (3.9) is negative. There are two positive roots indicating a range over which complete wrapping is favorable.
1
√
πK

1−

1
1 − ξ < Rn < √
πK

1+

1−ξ

(3.10)

where

ξ=

16πκK
A0 w2

(3.11)

Very small NP can not gain enough through adhesion to pay the cost of bending. Very large
NP simply require too much stretching.

3.1.1.4

Single layer micelle (sphere) states

The previous sections assumed a hydrophilic NP in contact with lipid heads and wrapped with a
bilayer. If instead, the NP is hydrophobic or attracted more to the lipid tails, then it will prefer
to be wrapped by a single layer and the energies will be slightly different. Again I consider the
embedded and detached states with and without pores. First, the detached states, corresponding to
states G and H in figure 3.2 will have an energy similar to F1a ,
2
2
K 2πRn − πR p
F4a =
2
A0

2

+ 4πκ + 2πR p γ − 4πR2 v
n

(3.12)

The bending and stretching terms are smaller because only a single layer of lipids is being
deformed. The adhesion term has a new constant, v, corresponding to the NP/tail interaction.
There is an increase in energy from the separation of monolayers around the nanoparticle. Since
that is also proportional to the surface area of the nanoparticle I have combined that effect with the
nanoparticle adhesion by lowering v appropriately. Although v and the NP/head attraction strength,

28

w, can be independent, what matters is their relative magnitudes. In later sections I will treat v as
a constant while allowing w to vary.
The embedded states corresponding to states I and J in figure 3.2 will be similar to F2a ,
2
2
K πRn − πR p
F5a =
2
A0

2

+ 4πκ + 2πR p γ − 4πR2 v
n

29

(3.13)

Figure 3.2:
Hydrophobic spherical NP states

G Detached micelle without a pore

H Detached micelle with a pore

F4a =

K
2

2
2πR2 −πR2
n
p
A0

+ 4πκ + 2πR p γ − 4πR2 v
n

I Embedded micelle without a pore

J Embedded micelle with a pore

F5a =

K
2

2
πR2 −πR2
n
p
A0

+ 4πκ + 2πR p γ − 4πR2 v
n

30

3.1.2

Fixed boundary LB with infinite rod

Now I will list the possible states of the same fixed boundary LB, instead interacting with an infinite
rod. Schematics of the side views of these states, looking down the length of the rod, are the same
as those for spherical NP. The energies are listed in figures 3.3 and 3.4. The subscripts, “b” and
“c”, denote the rod states opposed to the “a” previously used for the spherical NP. I will still use Rn
to now represent the radius of the rod. Additionally, I define the length of the rod as L. By infinite
rod, I mean that it spans the length of my simulation box which has periodic boundaries. The area
of the unstretched LB is therefore L2 .

Figure 3.3:
Hydrophilic rod-shaped NP states

A Reference state defined as zero energy

B Detached liposome with small or no pore

K
F1b =
2

πR2
p
2πRn −
L

2

+

C Detached liposome with a spanning pore

L
πκ
Rn

L
F1c = Rn πκ + 2Lγ − 2πRn Lw

+ 2πR p γ − 2πRn Lw

31

Figure 3.3 (cont’d)
D Embedded liposome with small or no pore

K
F2b =
2

πR2
p
2(π − 1)Rn −
2L

− 2πRn Lw + πR p γ +

E Embedded liposome with a spanning pore

2
L
F2c = Rn πκ + Lγ − 2πRn Lw

L
πκ
Rn

F Partially wrapped liposome

F3b = K 4R2
n
2

3.1.2.1

π f −2

f (1 − f )

2

L
+ Rn πκ f − 2πRn Lw f

Detached liposome (rod)

The detached states (B and C in figure 3.3) again consist of a completely wrapped NP (this time
rod-shaped) and a LB which may have a pore.
In the spherical NP cases, we assumed that the pore area which was proportional to the NP
area was much smaller than the area of the LB. In the infinite rod case, the length of the rod is
comparable (in our case equal) to the length of the LB. It is therefore possible for a pore to span
the entire length of the membrane. I separately consider the cases of small pores and pores that
span the entire length of the LB. The energy for a detached liposome state with a small round pore
is

32

K
F1b =
2

πR2
p
2πRn −
L

2

L
πκ + 2πR p γ − 2πRn Lw
Rn

+

(3.14)

The area term now contains the surface area of the rod as well as the pore area. Both terms in
parentheses are divided by L because of the division by the LB area. The bending term is now size
dependent because the curvature is still 1/Rn while the rod’s surface area is 2πRn L. Integrating
the curvature then gives

κ
2

Eb =

dA

=

1 2
Rn

=

κ
2

2πRn L
R2
n

L
πκ
Rn

(3.15)

The line energy term remains the same while the adhesion term is again proportional to the
rod’s surface area. Equation (3.14) also describes the state without a pore if R p is equal to zero.
In the case of a large spanning pore, the LB is free to retract and relieve all of the stretching
energy. The free energy of the spanning pore state is

F1c =

L
πκ + 2Lγ − 2πRn Lw
Rn

(3.16)

Besides the absence of stretching energy, the only other difference is that the line energy is
proportional to 2L because there are two edges equal to the length of the LB.

3.1.2.2

Embedded liposome (rod)

The embedded states (D and E in figure 3.3) can also have either a small or spanning pore. The
energy of the small pore state is
K
F2b =
2

πR2
p
2(π − 1)Rn −
2L

2

− 2πRn Lw + πR p γ +

L
πκ
Rn

(3.17)

The change in area is equal to the surface area of the rod (2πRn L) minus the displaced area
(2Rn L) and the pore area. The most energy efficient pore shape is a half circle adjacent to the
33

liposome. It allows a smaller line energy for a given pore area. That is why the line energy is
proportional to πR p . The bending and adhesion terms remain the same.
The energy for the spanning pore case is

F2c =

L
πκ + Lγ − 2πRn Lw
Rn

(3.18)

The line energy is half of that in equation (3.16) because there is only one edge.

3.1.2.3

Partially wrapped rod-shaped nanoparticle

The partially wrapped rod-shaped nanoparticle has an energy of

F3b =

K
4R2
n
2

π f −2

f (1 − f )

2

+

L
πκ f − 2πRn Lw f
Rn

(3.19)

The area term contains the surface area of the LB in contact with the NP minus the displaced
area. The bending and adhesion terms are the same as in states B-E except that they are proportional to the wrapping fraction, f .
As in the spherical NP case we can find the condition for wrapping by taking the derivative
with respect to f

dF3b
1−2f
= 4K π −
(fπ −2
df
f (1 − f )

f (1 − f ))R2 + L
n

πκ
− 2πRn w
Rn

(3.20)

Since I am assuming an infinite rod (L >> Rn ), only the terms proportional to L are significant.
Therefore the derivative is negative and wrapping is favorable if

Rn >

κ
2w

(3.21)

Notice that this condition is independent of f meaning that total wrapping will occur for rods
above the critical radius. For a finite rod, there is a stable partially wrapped state which I’ll describe
in section 3.3.2.

34

3.1.2.4

Single layer micelle (rod) states

Similar to the spherical hydrophobic NP states, I will now list the states in which the rod prefers to
be in contact with lipid tails and is therefore wrapped by a single layer of lipids, forming a micelle.
I begin with the detached rod-shaped micelle state with a small pore (G in figure 3.4). The energy
is given by
K
F4b =
2

πR2
p
πRn −
L

2

+ πκ

L
+ 2πR p γ − 2πRn Lv
2Rn

(3.22)

The stretching and bending terms are smaller than in state B (F1b ) because it is easier to bend
or stretch a single layer of lipids. The energy for a detached micelle with a spanning pore (H in
figure 3.4) is

F4c = πκ

L
+ 2Lγ − 2πRn Lv
2Rn

(3.23)

The bending energy is half of that for state C (F1c ). An embedded rod-shaped micelle with a
small pore (I in figure 3.4) has an energy of
K
F5b =
2

πR2
p
(π − 2)Rn −
2L

2

+ πκ

L
+ πR p γ − 2πRn Lv
2Rn

(3.24)

The stretching and bending terms are similarly smaller than in state D (F2b ) because of the
single layer wrapping. An embedded rod-shaped micelle with a spanning pore (J in figure 3.4) has
an energy of

F5c = πκ

L
+ Lγ − 2πRn Lv
2Rn

35

(3.25)

Figure 3.4:
Hydrophobic rod-shaped NP states
G Detached micelle with small or no pore

K
F4b =
2

πR2
p
πRn −
L

2

+ πκ

H Detached micelle with a spanning pore

L
2Rn

L
F4c = πκ 2Rn + 2Lγ − 2πRn Lv

+ 2πR p γ − 2πRn Lv
I Embedded micelle with small or no pore

K
F5b =
2

πR2
p
(π − 2)Rn −
2L

+ πκ

3.1.3

J Embedded micelle with a spanning pore

2
L
F5c = πκ 2Rn + Lγ − 2πRn Lv

L
+ πR p γ − 2πRn Lv
2Rn

Fixed tension ensemble

If the tension of the LB is held constant then it becomes impossible for a stable pore to exist. The
only benefit of a pore is to relieve tension. If the LB is kept at zero tension then any pore will heal

36

itself to avoid the unfavorable line energy around the edge of the pore. If the LB is under positive
tension, then there will be competition between tension relief and line energy.
As an example, let us consider the detached spherical liposome state. The only change from
F1a (equation 3.1) is in the stretching term. Under constant tension this term is simply equal to the
tension, Σ, times the change in area with respect to the reference state.

F1d = Σ(4πR2 − πR2 ) + 8πκ + 2πR p γ − 4πR2 w
n
p
n

(3.26)

First let us determine when wrapping is favorable before a pore forms. For that we need to find
when the energy is negative.

F1d (R p = 0) = π(4R2 (Σ − w) + 8κ) ≤ 0
n

⇒ Rnc ≥

2κ
w−Σ

(3.27)

A minimum requirement for wrapping is that the benefit from adhesion outweighs the cost of
stretching (w > Σ). Additionally, the nanoparticle must be larger than the critical size in equation
(3.27) in order to overcome the cost of bending which dominates for small nanoparticles. RuizHerrero et al. performed simulations at zero tension and found the same condition for wrapping
[63].
To see the dependence of pore size on the energy we can differentiate with respect to R p
dF1d
= 2π(γ − ΣR p )
dR p

(3.28)

The slope is initially positive meaning that there is a barrier to pore formation. The derivative
becomes negative for pores larger than

R pc = γ/Σ

37

(3.29)

Any pore larger than this will continue to grow. The size of the energy barrier can be found by
comparing the energy at the critical pore radius and the energy for no pore.

Ebarrier = F1d (R p = γ/Σ) − F1d (R p = 0)
= πγ 2 /Σ

(3.30)

While the critical radii and size of the energy barriers will be different for other states, the
general behavior remains the same. Very small NP will not become wrapped and pores will be
unstable for any value of Σ. Under positive tension small pores will heal themselves while pores
larger than a critical size will grow until the LB tears apart.

3.2

LB properties from simulation

Several LB properties can be determined from simulation. If the LB is slowly stretched biaxially,
the surface tension is expected to increase linearly with the area. The critical tension at which a
pore forms is related to the line tension around the pore. Cooke and Deserno performed such tests
on their LB model [37, 40]. I have run my own simulations using their model to measure these
properties. As discussed in section 3.1, the energy of a LB with a pore is

E=

K (A − A0 − πR2 )2
p
+ 2πγR p
2
A0

(3.31)

Tolpekina et al. [72] used the following relation between the membrane free energy and surface
tension in the canonical ensemble
dE
dL

= L⊥ L (2pzz − pxx − pyy ) = 2L Σ

(3.32)

N,V,T

where L and L⊥ are the simulation box dimensions parallelel and perpendicular to the membrane. pi j represent the pressure tensor components. The right-hand side serves as the definition
of the surface tension, Σ. The area of the LB, A0 , is equal to (L )2 . Differentiating equation (3.31)
gives
38

dE
dL

= 2L
N,V,T

⇒Σ=

K
(A − A0 − A p )
A0

K
(A − A0 − A p )
A0

(3.33)

The optimal pore size occurs when the derivative with respect to R p is zero.
KR p
dE
= 2π γ −
A − A0 − πR2
p
dR p
A0

=0

(3.34)

Equation (3.34) only has a positive root if its discriminant is positive.

Discriminant

4π 5 K 2
dE
=
4K 2 (A − A0 )3 − 27πγ 2 A2
0
4
dR p
A0

(3.35)

At the critical point (A = Ac ), the discriminant is equal to zero and can be solved for the line
energy constant, γ.

γ=

2K
3A0

(Ac − A0 )3
3π

(3.36)

Figure 3.5 shows surface tension versus area for my simulations. I measured the tension for
both slowly increasing (green) and decreasing (black) areas. Tension was calculated from the
pressure components which are easily accessible from the simulation. Although equation (3.33)
assumes a constant volume, it is still appropriate if the box is stretched very slowly.
Under low tension, the relation is linear as indicated by equation (3.33). By fitting my data to
that equation I was able to determine K ≈ 22 ε/σ 2 and A0 ≈ 1820 σ 2 . The hysteresis between the
green and black traces is due to the energy barrier for pore formation and is rate dependent. As the
rate of area change is decreased, the two traces should converge. The true critical point should lie
between the two peaks and depend only on the membrane properties indicated by equation (3.36).
By running constant area simulations over a range of areas, I observed the true critical point to
be around 1970 σ 2 , indicated by the red point in figure 3.5. Along with the other previously
determined parameters, equation (3.36) gives a value of γ ≈ 4.0 ε/σ . The values for K and γ are
39

Figure 3.5: Surface tension versus area in a flat LB. The green trace represents increasing area
while the black trace represents decreasing area. The red circle indicates the critical point.

in good agreement with those found by Deserno [40]. My values are slightly larger because my
simulations were run at a lower temperature.
The bending modulus, κ, can be estimated from thin plate theory as

κ=

Kh2
48(1 − ν 2 )

(3.37)

where h is the thickness of the membrane and ν is the Poisson ratio measuring compressibility.
Assuming a perfectly incompressible membrane, ν = 0.5. A LB thickness of 5 nm ≈ 5.6 σ , along
with my value for K gives a bending modulus of κ ≈ 19 ε. This value is also in good agreement
with the bending modulus found by Cooke and Deserno using fluctuation measurements [37]. I

40

will use these values when calculating phase diagrams in the next section.

3.3

Phase diagram

I will now compare the various states described in section 3.1 to see which are the most favorable
for different parameters. I will then create a phase diagram to compare with my simulation results
and with those of Ginzburg and Balijepalli [65].

3.3.1

Fixed boundaries with spherical NP

The embedded states will always have a lower energy because the stretching energy is smaller.
This is a result of my neglecting the curvature around the embedded liposome or micelle interface.
As previously mentioned, a large embedded NP will cause significant bending. In the limiting case
where the line energy around the interface is equal to that around the edge of a pore, the detached
states clearly become more favorable. In all of my small NP simulations, only embedded states
were observed.
My phase diagrams will depend on the properties of the NP, namely the size and strength of
attraction. By comparing the energies I will determine the transitions. First I will find the critical
NP size required for pore formation. A pore will form in the embedded liposome case when the
energies of states D and E are equal. It is useful to define the difference between the energies as

∆F2a = F2a (R p ) − F2a (R p = 0)
= πR p

πR p K
(−6R2 + R2 ) + 2γ
n
p
2A0

(3.38)
(3.39)

This function is by definition equal to zero at (R p = 0) and increases as R p increases. The pore
state is favorable if the function has a positive root, meaning that it eventually becomes negative
for a positive pore radius. The factor of R p corresponds to the zero at the origin. The zero we
are interested in must lie in the cubic function in parentheses. Cubic equations have multiple roots

41

only if their discriminant is positive. Therefore finding the discriminant of the term in parentheses
gives another constraint.

Discriminant[∆F2a ] =

27π 2 K 2
(2π 2 K 2 R6 − A2 γ 2 )
n
0
4
A0

(3.40)

Setting the discriminant equal to zero and solving for Rn determines the critical NP radius
which transitions from D to E. Nanoparticles larger than this will form an embedded liposome and
pore.

REa =

1/3
γA
√ 0
2πK

(3.41)

The same procedure can be used to compare the embedded micelle states with and without
pores. Again, I start with the difference between energies

∆F3a = F3a (R p ) − F3a (R p = 0)
=

πR p
πR p K(−2R2 + R2 ) + 4A0 γ
n
p
2A0

(3.42)
(3.43)

Taking the discriminant of the cubic factor gives

Discriminant[∆F5a ] = 256π 4 R6 K 4 − 432π 2 A2 γ 2 K 2
n
0

(3.44)

The critical NP radius at which an embedded micelle forms a pore, transitioning from I to J, is
therefore

RJa =

1/3
√
γA
3 √ 0
2πK

(3.45)

RE a and RJ a give clear transitions for micelle and liposome pore formation. Now I will find
when total wrapping becomes more favorable than the dissociated reference state. Since I am using
a single adhesion strength between the NP and tail beads, v, it is easier to find the transition from
the dissociated state (A) to the embedded micelle state (I). Setting F5a with no pore equal to zero
and solving for Rn gives
42

F5a (R p = 0) = π

⇒ RIa = 2

πR4 K
n
− 2R2 v + 4κ
n
2A0

A0 v
πK

1−

1−

πκK
2A0 v2

=0

(3.46)

The sizes of the three transition radii are as follows

RIa < REa < RJa

(3.47)

This defines four regions in each of which there are two possible states. For NP smaller than
RIa the possible states are A and D. For NP between RIa and REa the possible states are D and I.
Between REa and RJa the possible states are E and I. For NP larger than RJa the possible states
are E and J. To find the adhesion strength (w) dependence of the transitions we must compare the
energies in each region.
The transition from the dissociated reference state (A) to the embedded liposome (D) can be
found by setting F2a with no pore equal to zero and solving for w

wADa =

9πR2 K 2κ
n
+ 2
8A0
Rn

(3.48)

There is also a partially wrapped state (F). In section 3.1 I found the minimum NP radius
required for partial wrapping. If equation(3.5) is solved instead for w, we find the transition from
detached (A) to partially wrapped (F).
2κ
wAFa = 2
Rn

(3.49)

which is the second term in equation (3.48). Once the adhesion strength reaches wADa , the
fully wrapped, embedded state is favorable. The range of w over which a NP becomes partially
wrapped is given by the first term in equation (3.48). The partially wrapped phase becomes larger
as the surface area of the NP increases with respect to the LB area. For small NP (R2 << A0 ), the
n
partially wrapped phase is negligible.
43

The transition from state D to I in the second region can be found by setting F2a equal to F5a
and solving for w

wDIa =

κ
πR2 K
n
+ 2 +v
A0
Rn

(3.50)

To find the remaining two transitions we must first solve for the optimal pore sizes for states E
and J. This is done by setting the derivative of F2a with respect to R p equal to zero and solving for
R p.
dF2a
= 2π
dR p

πR p K 2
(R p − 3R2 ) + γ
n
A0

=0

2
( αa − 1 − αa )2/3 + 1
2
( αa − 1 − αa )1/3

⇒ R p2a = Rn

(3.51)

where

αa =

A0 γ
2πR3 K
n

(3.52)

The optimal micelle pore radius can be found similarly.
πR p K 2
(R p − R2 ) + γ
n
A0

dF5a
= 2π
dR p
Rn
⇒ R p5a = √
3

=0

2
( βa − 1 − βa )2/3 + 1
2
( βa − 1 − βa )1/3

(3.53)

where
√
3 3A0 γ
βa =
2πR3 K
n

(3.54)

The transition from state E to I in the third region can be found by setting F2a equal to F5a with
no pore and solving for w.

F2a (R p = R p2a ) = F5a (R p = 0)
44

⇒ wEIa =

πK
4A0

4R2 − 3R2 +
n
p2a

R4
p2a
2R2
n

1
+v+ 2
Rn

R p2a
γ +κ
2

(3.55)

where R p2a is given by equation (3.51). The transition from state E to J in the fourth region can
be found by setting F2a equal to F5a with the appropriate pore size substitutions (3.51) and (3.53)
and solving for w.

F2a (R p = R p2a ) = F5a (R p = R p5a )

πK
⇒ wEJa =
4A0

4R2 +
n

1
+v + 2
Rn

R4 − R4
p2a
p5a
2R2
n

+ R2 − 3R2
p2a
p5a

(3.56)

R p2a − R p5a
γ +κ
2

Figure 3.6 shows the phase diagram of w versus Rn with wADa , wAFa , wDIa , wEIa , and wEJa
plotted over their respective regions. The partially wrapped state is present but is narrower than
the line thickness of the plot. As mentioned before, the favorable region for partial wrapping is
negligible for NP small compared to the LB.

3.3.2

Fixed boundaries with infinite rod NP

Similar to the previous section I will now explain the phase transitions for the case of an infinite
rod NP interacting with a LB with fixed boundaries. The energies are given in figures 3.3 and 3.4.
The phase diagram looks similar to that of the spherical NP case with a couple differences. The
relative sizes of the critical radii depend on the size of the membrane. For that reason I will make
one phase diagram for large L and one for smaller L.
Embedded states are always preferable to detached states because there will be less stretching.
Again, this is because I have neglected the curvature around the interface of the embedded NP
which I believe is appropriate for small NP. We are left with states A, D, E, F, I, and J as possibili-

45

Figure 3.6: Phase diagram of spherical NP states. The x-axis shows the radius of the nanoparticle.
The y-axis shows the attraction between the nanoparticle and lipid heads. Attraction increases
with w. (A) Dissociated reference state, (D) Embedded liposome without a pore, (E) Embedded
liposome with a pore, (I) Embedded single-layer micelle without a pore, (J) Embedded singlelayer micelle with a pore. Schematics and energies for each phase are given in figures 3.1 and
3.2.

ties. To see when a small pore becomes favorable with an embedded liposome we find where F2b
is smaller than its value at R p = 0.

∆F2b = F2b − F2b (R p = 0)
= πR p

K
L

πR3
p
+ (π − 1)Rn R p + γ
8L

(3.57)
(3.58)

This has a positive root if the discriminant of the cubic expression in brackets is positive.
πK 2
Discriminant[∆F2b ] =
32(π − 1)3 R3 K 2 − 27πLγ 2
n
5
64L
46

(3.59)

πLγ 2
4K 2

3
⇒ RDb =
2(π − 1)

1/3

(3.60)

The same procedure can be used to find the critical NP size required to form a small pore in
the embedded micelle case.

∆F5b = F5b − F5b (R p = 0)
= πR p γ −

Discriminant[∆F5b ] =

R pK
2L

(3.61)

(π − 2)Rn −

πR2
p
4L

πK 2
4(π − 2)3 R3 K 2 − 27πLγ 2
n
5
64L

3
⇒ RIb1 =
π −2

πLγ 2
4K 2

(3.62)

(3.63)

1/3

(3.64)

RDb is always smaller than RIb1 . To find the minimum NP radius required to form an embedded
micelle with no pore we can set F5b equal to zero (the reference energy) and solve for Rn .
πLκ
1
=0
F5b (R p = 0) = π 2 R2 K − 2πLRn v +
n
2
2Rn

⇒ RIb2 =

1
4φ

κ
v

+3 sin

√
√
1
3 − 3 cos cos−1 1 − 2φ 2
3
1
cos−1 1 − 2φ 2
3

(3.65)

(3.66)
(3.67)

where

φ=

3(π − 2)2 K
16πL

3κ
v3

(3.68)

These three critical points, RDb , RIb1 , and RIb2 , are again independent of w and will appear as
vertical segments in the phase diagram. RIb2 can be either larger or smaller than RIb1 depending
on the size of L.
47

The optimal pore radii in states D and I can be found by seeing when the derivatives of F2b and
F5b with respect to R p are equal to zero. By substituting the optimal pore radius back into F2b and
comparing with the energy of the spanning pore state (F2c ), we can find the NP transition radius
from state D to E, REb .
The NP transition radius from state I to J, RJb , can similarly be found by substituting the
optimal pore radius into F5b and setting the energy equal to F5c . RJb is always larger than both
RIb1 and REb , but the order of RIb1 and REb depends on the size of L.
Now that I have calculated all of the transitions which are strictly radius dependent, I will
calculate the transitions which depend on the adhesion constant, w.

3.3.2.1

Small L

For relatively small L, the order of the radial transitions is

RDb < RIb2 < REb < RIb1 < RJb

(3.69)

which defines six regions in which w dependent transitions are required. For NP smaller than
RDb , the possible states are A, D, and F. The transition from state A to D can be found by setting
F2b with no pore equal to zero and solving for w.

F2b (R p = 0) = 2(π − 1)2 R2 K + πL
n

κ
− 2Rn w = 0
Rn

(π − 1)2 Rn K
κ
⇒ wADb1 =
+ 2
πL
2Rn

(3.70)

As in the spherical case, there is also a stable, partially wrapped state which can be found by
finding when the derivative of F3b with respect to f as f approaches zero, is negative.
dF3b
df

= 8R2 K + πL
n
f →0

⇒ wAFb =

κ
− 2Rn w = 0
Rn

4Rn K
κ
+ 2
πL
2Rn
48

(3.71)

For NP between RDb and RIb2 , the possible states are A and D with a small pore. The transition
can be found by substituting the optimal pore radius into F2b and finding where it becomes zero.

F2b R p = R p2b = 0
⇒ wADb2 =

2
K
1
Lκ
πR2 − 4(π − 1)Rn L + R p2b γ +
p2b
2
2Rn L 8πL
Rn

(3.72)

where R p2b is the optimal small pore radius. For NP between RIb2 and REb , the possible states
are I without a pore and D with a small pore. The transition can be found by substituting the
optimal pore radius into F2b and setting the energy equal to F5b without a pore. Solving for w
gives

F2b R p = R p2b = F5b R p = 0
⇒ wD2I1b =

K
1
2Rn L − R2
p2b
2
4Rn L 4L
κ
+2R p2b γ + L 4Rn v +
Rn

2(3π − 4)Rn L − πR2
p2b
(3.73)

For NP between REb and RIb1 , the possible states are I without a pore and E. The transition
can be found by setting F5b without a pore equal to F2c and solving for w.

F2c = F5b R p = 0
⇒ wEI1b =

(π − 2)2 R3 K
1
n
2Rn (2πRn v + γ) + πκ −
L
4πR2
n

(3.74)

For NP between RIb1 and RJb , the possible states are I with a small pore and E. The transition
can be found by setting F5b with its optimal pore radius equal to F2c and solving for w.

F2c = F5b R p = R p5b
⇒ wEI2b =

1
πLκ
(2πRn v + γ)L − πR p5b γ +
2πRn L
2Rn
K
− 2 πR2 − 2(π − 2)Rn L
p5b
8L
49

(3.75)

Figure 3.7: Phase diagram of short rod-shaped NP states. The x-axis shows the radius of the
nanoparticle. The y-axis shows the attraction between the nanoparticle and lipid heads. Attraction
increases with w. (A) Dissociated reference state, (D) Embedded liposome with and without a
pore, (E) Embedded liposome with a spanning pore, (I) Embedded single-layer micelle with and
without a pore, (J) Embedded single-layer micelle with a spanning pore. Schematics and energies
for each phase are given in figures 3.3 and 3.4.

Finally, for NP larger than RJb , the possible states are E and J. The transition can be found by
setting F2c equal to F5c and solving for w.

F2c = F5c
⇒ wEJb =

κ
+v
4R2
n

(3.76)

Figure 3.7 shows the w vs. Rn phase diagram for L = 40σ . The partially wrapped state between
A and D is narrower than the line width.

50

3.3.2.2

Large L

For larger L, the order of the radial transitions becomes

RIb2 < RDb < RIb1 < REb < RJb

(3.77)

This introduces two new transitions. For NP between RIb2 and RDb , the possible states are D
without a pore and I without a pore. The transition can be found by setting F2b equal to F5b and
solving for w.

F2b R p = 0 = F5b R p = 0
⇒ wD1I1b =

κ
(3π − 4)Rn K
+ 2 +v
4L
4Rn

(3.78)

For NP between RIb1 and REb , the possible states are D and I, each with a small pore. The
transition can be found by substituting the optimal pore radii into F2b and F5b and finding when
they are equal.

F2b R p = R p2b = F5b R p = R p5b
⇒ wD2I2b =

K
1
R2 − R2 − 2Rn L
p2b
p5b
2
4Rn L L
+2 R p2b − R p5b γ + L 4Rn v +

8Rn L + π R2 + R2 − 6Rn L
p2b
p5b
κ
Rn

(3.79)

where R p2b and R p5b are the optimal small pore radii for the embedded doubly and singly
wrapped states respectively. Figure 3.8 shows the w vs. Rn phase diagram for a much longer rod
(L = 1000σ ).

3.4

Comparison with Ginzburg and Balijepalli

I will now see how well my phase diagram agrees with the DFT results of Ginzburg and Balijepalli
[65]. As explained in section 2.5.1 they calculated the favorable states for NP/LB interactions while

51

Figure 3.8: Phase diagram of long rod-shaped NP states. The x-axis shows the radius of the
nanoparticle. The y-axis shows the attraction between the nanoparticle and lipid heads. Attraction
increases with w. (A) Dissociated reference state, (D) Embedded liposome with and without a
pore, (E) Embedded liposome with a spanning pore, (I) Embedded single-layer micelle with and
without a pore, (J) Embedded single-layer micelle with a spanning pore. Schematics and energies
for each phase are given in figures 3.3 and 3.4.

varying the size and charge density of the NP. In order to compare with Ginzburg and Balijepalli’s
results I first need to relate the Florry-Huggins parameter, χ, to the adhesion energy density, w. Let
us start with the definition of the Flory-Huggins parameter between two types of material (say A
and B).

χAB =

z
1
wAB − (wAA + wBB )
τ
2

(3.80)

where z is the coordination of a particular bead, τ is equal to the Boltzmann constant times
the temperature, and w is the interaction energy between two beads. The Flory-Huggins parameter

52

indicates the energy difference between mixed and unmixed states. The three types of beads I use
below are hydrophilic (H), lipophilic (L), and nanoparticle (P), corresponding to the head, tail, and
NP beads respectively. The attraction between lipid beads is fixed by the Cooke model, so I will
use those constraints to set the energy scale. The hydrophilic-lipophilic Flory-Huggins parameter
is given by

χHL =

1
z
wHL − (wHH + wLL ) = 1.0
τ
2

(3.81)

The value on the right hand side of this equation as well as those used in the following equations
are the values used by Ginzburg and Balijepalli. The Cooke model sets wHL = wHH = 0 and
wLL = −ε. This gives us the relation between the energy scales.
z 2
=
τ
ε

(3.82)

Now let us look at the particle-tail interaction.

χPL =

z
1
wPL − (wPP + wLL ) = 0.0
τ
2

(3.83)

While wLL is fixed, we have the freedom to choose wPL = wPP = −ε to satisfy the relation.
Finally, we can determine the nanoparticle-head interaction.
z
1
wPH − (wPP + wHH )
τ
2

(3.84)

ε
wPH = (χPH − 1)
2

χPH =

(3.85)

I can now appropriately map Flory-Huggins parameters in Ginzburg and Balijepalli’s model to
interaction strengths in the Cooke model. Since the Cooke model does not have an explicit solvent,
we can safely think of it as having solvent beads with effectively zero interaction strength with all
other beads. For the sake of completeness let us look at the Flory-Huggins parameters involving
the solvent.

53

χLW =

z
1
wLW − (wLL + wWW ) = 1.0
τ
2

(3.86)

The new subscript, W , denotes the solvent beads. wLW and wWW are not uniquely determined,
but we can choose wLW = wWW = 0 to satisfy the relation. Similarly,

χHW =

1
z
wHW − (wHH + wWW ) = 0.0
τ
2

(3.87)

χPW =

z
1
wPW − (wPP + wWW ) = 1.0
τ
2

(3.88)

give wHW = wPW = 0.
The only interaction varied by Ginzburg and Balijepalli is χPH which is related to the NP
bead/head bead attraction, wPH . The adhesion energy density referred to in the phase diagrams of
the previous section can be calculated by estimating the number of nearest neighbors to be

w = 1.89 − 9.51wPH

(3.89)

Substituting equation (3.85) gives a linear relation between the Florry-Huggins parameter and
adhesion energy density.

w = 6.6 − 4.8χPH

(3.90)

Ginzburg and Balijepalli varied χPH from -3 to +1 corresponding to w ranging from 2 to 20.
Their NP radii ranged from 1.6 nm to 3.2 nm. Figure 3.9 shows my predicted phase diagram along
with Ginzburg and Balijepalli’s results. The blue circles indicate an embedded micelle with no
pore which corresponds to state I. The red squares indicate an embedded liposome with no pore,
corresponding to state D. The green diamonds indicate an embedded liposome state with a pore,
corresponding to state E.
One discrepancy between my phase diagram and Ginzburg and Balijepalli’s results is the three
points in the lower left of figure 3.9. My phase diagram appears to overestimate the wrapping

54

Figure 3.9: Phase diagram of spherical NP states. The x-axis shows the radius of the nanoparticle.
The y-axis shows the attraction between the nanoparticle and lipid heads. (A) Dissociated reference
state, (D) Embedded liposome without a pore, (E) Embedded liposome with a pore, (I) Embedded
single-layer micelle without a pore, (J) Embedded single-layer micelle with a pore. Schematics
and energies for each phase are given in figures (3.1 and 3.2).
Data points indicate the results of Ginzburg and Balijepalli. (Blue circles): Embedded micelles
with no pores. (Red squares): Embedded liposomes with no pores. (Green diamonds): Embedded liposomes with pores.

transitions to states D and I. I believe this comes from an overestimate of the bending energy
which dominates for small NP wrapping. I assumed bending around a closed spherical surface.
However, there is a band around the equator of the NP where it intersects with the LB that should
not be counted. This region can be ignored if the NP is much thicker than the LB, but it makes up
a significant fraction of the surface area for small NP. Appropriate reduction of the bending energy
would lower the wrapping transition.
Another discrepancy in my phase diagram is related to pore formation. My prediction of a

55

larger critical NP size required for pore formation goes back to my neglect of the curvature around
the NP/LB interface. As the NP gets larger that neglected curvature becomes comparable to the
curvature around the edge of a pore. By including that curvature cost, the transition from state D
to E would occur for slightly smaller NP.
The final discrepancy is that I have predicted the pore transition to be independent of w while
Ginzburg and Balijepalli only observed pore formation for strongly attractive NP. I assumed that
stretching occurred uniformly throughout the LB with a small uniform decrease in lipid concentration and that the adhesion was independent of stretching. It’s true that a weakly attractive NP will
cause uniform stretching throughout the LB regardless of whether or not it is in contact with the
NP. However, as w increases, lipids will pack more tightly around the NP, increasing the stretching away from the NP. In that way, a strongly attractive NP could increase the stretching above a
critical tension, causing a pore. Combining the stretching and adhesion energy terms may fix the
discrepancy.
One could argue that a more appropriate comparison with Ginzburg and Balijepalli’s results
would be with the infinite rod phase diagram shown in figure 3.7. While the predicted transition
from singly wrapped micelles to doubly wrapped liposomes is in agreement, there are significant
discrepancies with the predicted pore transitions. The spanning pore transition for hydrophilic NP
should occur around Rn = 1.1 nm according to my predictions compared to the 2.4 nm observed
by Ginzburg and Balijepalli. I would also expect a spanning pore transition for hydrophobic NP to
occur around Rn = 2.1 nm while such a state would presumably exist above 3.2 nm in the Ginzburg
and Balijepalli model.
While qualitative comparisons can be made, differences in the models prevent quantitative
agreement regarding transition points. Ginzburg and Balijepalli stressed that their results represented only a qualitative view of the dependency of these transitions on NP size and surface treatment. Hydrophobic NP will become singly wrapped while hydrophilic NP will become doubly
wrapped. There is also a critical NP size in the nanometer range which causes pore formation. On
these general trends, our models are in agreement. A more careful comparison would require de-

56

termining the the bending and stretching moduli which correspond to the Ginzburg and Balijepalli
model.

3.5

Comparison with MD simulations

I have used molecular dynamics to simulate the NP/LB interactions described in the previous
sections for both spherical and rod-shaped NP. Interactions were simulated using the Cooke model
[36, 37] which was explained in chapter 2. As was previously stated, the Cooke model was chosen
for it’s ability to handle large systems. For the purpose of comparison, I tried to mimic the same
system properties used by Ginzburg and Balijepalli [65]. In this section I will first go over the
details of my particular implementation and then I will present my results.

3.5.1

Implementation Details

Simulations were run using LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator), which is an open source molecular dynamics package distributed by Sandia National Laboratories. One benefit of using this package is that it can easily be run on parallel systems. My
simulations were run on the High Performance Computer Cluster (HPCC) at Michigan State University. The HPCC consists of multiple clusters with both AMD and Intel processors, both dual
and quad cores, ranging from 2.2 GHz to 2.4 GHz per core. I used generic submissions to the
cluster which used a combination of any available nodes.
Nanoparticles were chosen to have the same radii as those used by Ginzburg and Balijepalli
(1.6-3.2 nm) [65]. A square patch of LB containing 3200 lipids was used. The simulation box
had fixed periodic boundary conditions with dimensions roughly 40x40x27 nm3 . The horizontal
dimensions were initially equilibrated at zero tension before the NP was introduced. The temperature was held constant at kB T = 1.0 ε using a Nosé-Hoover thermostat [73, 74]. That temperature
was shown to be in the liquid phase by Cooke and Deserno [37]. I confirmed this by measuring the
projected area as a function of temperature. The phase diagram is shown in figure 3.10. The lighter

57

Figure 3.10: Area per lipid versus temperature from Cooke model simulation. The dark lines
represent running averages of the raw data which is displayed as lighter dots. In the top trace the
temperature is decreasing while in the bottom trace it is increasing. The gel to liquid transition
occurs around kB T = 0.98 ε while heating.

dots indicate the raw data while the dark lines are running averages to clearly show the transition.
The top trace represents the cooling transition while the bottom trace represents the heating transition. Cooke and Deserno observed a similar hysteresis. The gel-fluid transition occurs around
kB T = 0.98 ε in the lower trace.
I also looked at the radial distribution function of the bilayer at various temperatures during the
heating transition as shown in figure 3.11. It is clear from the disappearance of long range order
that the liquid transition occurs between kB T = 0.97 ε and 0.99 ε.
The main tunable parameter in the Cooke model is the potential width, wc , which I set to 1.6 σ .

58

Figure 3.11: Radial distribution function for a lipid bilayer at various temperatures during the
heating transition. The lack of long range order in kB T = 0.99 ε indicates the liquid phase.

That set the overall cutoff distance to 2.6 σ ≈ 2.3 nm. Nanoparticles consisted of beads which were
the same size as the lipid beads (σ ), arranged in a spherical cutout of an FCC lattice. The relative
positions of the NP beads were held rigid throughout the simulation.
Interactions between NP and lipid beads were also governed by the Cooke potential. Different strengths of the potential were used to correspond to the Flory-Huggins parameters used by
Ginzburg and Balijepalli. The relation between the Florry-Huggins parameter and the interaction
between individual beads was explained in the previous section.
Simulations were run in the canonical ensemble, keeping the temperature, volume, and number
of particles constant, similar to Ginzburg and Balijepalli’s work. I used a Nosé-Hoover thermostat

59

[73, 74] which is explained in appendix 5. A single damping parameter, Q, is used to determine
how strongly the system is coupled to the heat bath. It is the effective “mass” associated with the
virtual variable s. Great care must be taken to select appropriate Q values [75]. I explain my choice
of Q in appendix 5.
I used a time step of 0.01 τ which was used by Cooke and Deserno [37]. The bilayer was
allowed to equilibrate for 104 τ (≈ 93 ns) before introducing the NP. The NP was initially placed
at rest about 1 nm from the bilayer. Simulations were run for another 104 τ after the NP was
introduced. Each simulation ran in parallel on nine cores of the HPCC for approximately one hour.
LAMMPS divided the simulation box into a 3x3 grid and calculated interactions in each region on
a separate processor.

3.5.2

MD Results

In this section I will describe my simulation results. As explained in the previous section I tried to
reproduce the same conditions used by Ginzburg and Balijepalli [65]. Again, the bilayer contained
3200 lipids and had an area of approximately 40x40 nm2 . The temperature was set to kT = ε and,
like the area, was held constant. I used five different sizes of spherical nanoparticles with radii
ranging from 1.6 nm to 3.2 nm. I also ran simulations on infinite rods with the same radii. The
Flory-Huggins parameter governing the attraction between NP and lipid head beads ranged from
-3 to +1.

3.5.2.1

Spherical NP

Figure 3.12 shows the initial configuration for the spherical NP system. For every combination of
parameters the NP became embedded in the bilayer. For each configuration, I observed whether
the NP became wrapped by a single layer of lipids which Ginzburg and Balijepalli called a singlelayer hybrid micelle or by a bilayer which they called a bilayer hybrid micelle. Figures 3.13 and
3.14 show cross sections of those two outcomes.
In addition to the type of micelle formed I also observed the formation of pores by certain NP.
60

Figure 3.12: Typical starting configuration of my spherical nanoparticle/lipid bilayer simulations.

Figure 3.13: Cross section of single-layer hybrid micelle. A hydrophobic nanoparticle is embedded
between two single layers of lipids.

In each of those cases the NP was wrapped by a bilayer. Figure 3.15 shows the top view of a NP
embedded in a bilayer and the accompanying pore that is located adjacent to the NP as explained
in section 3.1.1.2. The pore is mostly round while sharing an edge with the liposome. That allows
it to maximize its area while minimizing the line energy around the edge of the pore.
Figure 3.16 shows my predicted phase diagram along with my simulation results. The symbols
have the same meaning as in figure 3.9. The same discrepancies are present in the lower left and
upper right. The same explanation given at the end of section 3.4 applies.

61

Figure 3.14: Cross section of bilayer hybrid micelle. A hydrophilic nanoparticle is surrounded by
a bilayer and embedded in the membrane.

3.5.2.2

Rod-shaped NP

I performed the same simulations with “infinite” rod nanoparticles which spanned the simulation box with periodic boundaries. The radii of the rods were the same as those of the spherical
nanoparticles. Figure 3.17 shows the starting configuration for the 3.2 nm rod.
In all cases the NP became embedded as either a singly wrapped micelle or a doubly wrapped
liposome. Figure 3.18 includes the same phase diagram as figure 3.7 except scaled for the sizes I
have used. The data points represent my simulation results.
The blue circles in the lower left represent narrow hydrophobic rods which form embedded
micelles without a pore. These agree with state I predicted by my phase diagram. Figure 3.19
shows the resulting equilibrium state.
The red squares in the lower right represent thicker hydrophobic rods which form embedded
micelles but also form small pores. These also agree with state I and are above the critical NP
radius required for pore formation. Figure 3.20 shows one of these states. Notice the semicircular
pore located adjacent to the micelle which is the most efficient shape.
The green diamonds in the upper region represent hydrophilic rods. In all cases the rod formed

62

Figure 3.15: Strongly hydrophilic nanoparticle is wrapped by a bilayer and causes a pore in the
membrane.

an embedded liposome and caused a pore which spanned the entire length of the LB. This is in
agreement with state E predicted by my phase diagram. Figure 3.21 shows one such final state.
The left side of the figure shows the severed edge that continues straight across the membrane. The
right side of the liposome is smooth where it detached from the membrane.

63

Figure 3.16: Phase diagram of spherical NP states. The x-axis shows the radius of the nanoparticle.
The y-axis shows the attraction between the nanoparticle and lipid heads. (A) Dissociated reference
state, (D) Embedded liposome without a pore, (E) Embedded liposome with a pore, (I) Embedded
single-layer micelle without a pore, (J) Embedded single-layer micelle with a pore. Schematics
and energies for each phase are given in figures (3.1 and 3.2).
Data points indicate my molecular dynamics results. (Blue circles): Embedded micelles with
no pores. (Red squares): Embedded liposomes with no pores. (Green diamonds): Embedded
liposomes with pores.

64

Figure 3.17: Typical starting configuration of my rod-shaped nanoparticle/lipid bilayer simulations.

65

Figure 3.18: Phase diagram of rod-shaped NP states. The x-axis shows the radius of the nanoparticle. The y-axis shows the attraction between the nanoparticle and lipid heads. (D) Embedded
liposome with a small pore, (E) Embedded liposome with a spanning pore, (I) Embedded singlelayer micelle with a small pore or no pore. Schematics and energies for each phase are given in
figures (3.3 and 3.4).
Data points indicate my molecular dynamics results. (Blue circles): Embedded micelles with no
pores. (Red squares): Embedded micelles with small pores. (Green diamonds): Embedded
liposomes with spanning pores.

66

Figure 3.19: Narrow, hydrophobic, rod-shaped nanoparticle embedded in a bilayer. It is wrapped
by a single layer of lipids and does not form a pore.

Figure 3.20: Thick, hydrophobic, rod-shaped nanoparticle embedded in a bilayer. It is wrapped by
a single layer of lipids and does form a pore.

67

Figure 3.21: Strongly hydrophilic, rod-shaped nanoparticle embedded in a bilayer. It is wrapped
by a double layer of lipids and has caused the membrane to rupture.

68

Chapter 4
LIPID BILAYER PORE DYNAMICS

4.1

Tension required for pore formation

There is an energy cost associated with membrane pores due to the curvature around the edge
of the pore. This energy increases linearly with the radius of the pore. Therefore a pore is not
energetically favorable unless a surface tension is applied to the membrane. Under low tension a
pore will not form, or if formed will heal over time. For sufficiently large tension, a pore can lower
the total energy by relieving this tension. In that case the energy will continue to decrease as the
pore grows until all of the tension is released. In this section I will calculate the critical relation
between surface tension and pore size.

4.1.1

Fracture mechanics

This problem is similar to the fracturing of materials. Griffith studied the effect of cracks on
stressed materials [76]. Consider an object with an applied stress, σ , which experiences a corresponding strain, ε. The Young’s modulus, E, is equal to the ratio σ /ε and is an intrinsic property
of the material. The energy stored in the stretched material per unit volume is

U∗ =
U∗ =

σ dε

Eε 2 σ 2
=
2
2E

(4.1)

Now consider a crack of length a perpendicular to the applied stress which spans the entire
height of the material (t). The crack will lower the stress energy for part of the object while adding
energy due to the broken bonds along the inner surface of the crack. If the crack starts at the outer
edge of the object, the volume of the region of relieved stress will be a triangular prism of height
t. The area of the triangular base turns out to be πa2 . The relieved stress energy is then
69

U =−

σ2
(πa2t)
2E

(4.2)

The additional surface energy is equal to the inner surface area of the crack, 2at, times a
material specific surface energy constant, γ.

S = 2γat

(4.3)

The critical crack size can be found by setting the derivative of the total energy equal to zero
πσ 2t
d(U + S)
=−
a + 2γt = 0
da
E
2γE
⇒ ac =
πσ 2

(4.4)
(4.5)

For a given stress, ac is the critical crack size. Cracks larger than this will continue to propagate.
Solving equation (4.4) instead for σ gives the critical stress required for a starting crack of length
a to grow.
2γE
πa

σc =

(4.6)

Notice that these relations do not depend on the thickness of the object which means a similar
approach can be used to analyze the continuum membrane model.

4.1.2

Membrane rupturing

In section 3.1.3 I used a similar energy balancing to look at the criteria for pore formation in a
membrane under tension. I will briefly rederive the critical pore radius here. Considering only the
pore dependent terms from equation (3.26) the free energy consists of a stretching term and a line
energy term

F = −πR2 Σ + 2πR p γ
p

70

(4.7)

where Σ is the surface tension and γ is a constant related to the bending modulus. The line
energy is due to the fact that lipid heads will curve around the edge of the pore to cover the
hydrophobic tails and is proportional to the circumference of the pore. Differentiating this energy
with respect to R p gives
dF
= −2πΣR p + 2πγ = 0
dR p
γ
⇒ R pc =
Σ

(4.8)
(4.9)

For a given applied tension, R pc gives the critical pore size. Small pores will heal over time.
Pores larger than the critical size will grow indefinitely. The metastable pore can be observed if
the membrane area is held constant. If a membrane is stretched to sufficient tension a pore will
form. Figure 3.5 from section 3.2 shows tension versus membrane area for a simulation in which
the area was slowly increased. The sharp drop indicates pore formation. Beyond that transition the
tension decreases like A−1/2 as shown by Cooke and Deserno [37, 40]. While the entire membrane
does stretch, most of the increased area is added to the pore. Therefore the critical pore radius and
critical tension are in fact inversely proportional. For any area above the critical value in figure
3.5 there is a pore and a corresponding tension given by equation (4.9). If the simulation is run
instead at constant tension, then the pore will be unstable. If either the pore radius or tension are
decreased, the pore will heal. If either are instead increased, then the pore will grow until it spans
the entire membrane.

4.2

Dynamics of pore healing

In this section I will determine theoretically how the size of a pore changes over time and compare
with MD simulations.

71

4.2.1

Predicted time dependence

I will now calculate how the time required for pore healing scales with the size of the pore. In order
for the pore to heal in the presence of tension, the pore radius must not exceed the critical value
in equation (4.9). In the absence of tension, the pore energy will always be positive and it will
therefore always heal. A flat lipid bilayer with a pore can be viewed as a two dimensional system
with two phases. The main phase is the homogeneous region containing lipids. The second phase
is the pore which contains no lipids. We can use lipid concentration (c) as the order parameter to
distinguish the two phases. In the case of zero tension the concentration is not conserved. Filling
in the pore (increasing the concentration in one phase) does not change the concentration in the
surrounding area. It is known that the interface between two phases in a non-conserved system
will move with a speed proportional to the mean curvature of the interface [77]. I will rederive that
result below.
Cahn and Hilliard showed that the free energy of a phase segregated system contains a homogeneous term representing the energy within each region as well as gradient term representing the
energy of the interface [78]. The total free energy is then

F = Nv

V

[ f0 (c) + κ(∇c)2 + . . . ]dV

(4.10)

where Nv is the number of particles per unit volume, f0 (c) is the homogeneous volume energy
which depends only on the concentration, c. Only even powers of the gradient are present because
the scalar energy must be invariant with respect to the direction of the gradient. While the energy
can depend on higher powers of the gradient, only the lowest order term is kept if we assume the
gradient to be much smaller than the reciprocal of the distance between particles. The variational
derivative of the free energy with respect to the concentration is

δF
δc

=

∂ f0
∂
−∇
κ(∇c)2
∂c
∂ (∇c)

= f0 (c) − 2κ(∇2 c)

72

(4.11)
(4.12)

The variational derivative is a local thermodynamic force which is assumed to be proportional
to the rate of change of the order parameter [77].
δF
∂c
= −α
∂t
δc

(4.13)

This is a phenomenological assertion. It contains only the first time derivative because the
system is in an overdamped liquid state. Combining equations (4.11) and (4.13) leads to
∂c
= α(2κ(∇2 c) − f0 (c))
∂t

(4.14)

f (c)
1 ∂c
= ∇2 c − 0
2ακ ∂t
2κ

(4.15)

⇒

By absorbing 2ακ into the time scale and defining
f (c)
V (c) = 0
2κ

(4.16)

we arrive at the time-dependent Ginzburg-Landau equation
∂c
= ∇2 c −V (c)
∂t

(4.17)

Now let us apply this equation to the case of a pore centered at the origin in a two dimensional
membrane. I will follow the same procedure used by Krapivsky et al. to look at the shrinking
of a spherical droplet [79]. The system can be described in polar coordinates with no angular
dependence. We can assume the concentration has the following form

c(r,t) = m(r − R(t))

(4.18)

where r is the radial coordinate and R(t) is the time dependent pore radius. If we define
x = r − R(t), then m(x) is defined over the domain x ≥ −R(t) with a sharp transition at x = 0. The
concentration is zero for negative x and becomes a positive constant for positive x. The left-hand
side of equation (4.17) becomes
73

∂m ∂m ∂x
∂R
=
= −m
∂t
∂ x ∂t
∂t

(4.19)

Using polar coordinates, the Laplacian term becomes
∂m
1 ∂
r
r ∂r
∂r

=m +

m
m
=m +
r
R(t) + x

(4.20)

Equation (4.17) is then

m +

1
∂R
+
m −V (m) = 0
R(t) + x ∂t

(4.21)

Multiplying equation (4.21) by m and integrating over x from −R(t) to ∞ gives
∞
∞
∞
1
1
∂R
2
2 dx −V (x)
(m )
+
+
(m )
=0
2
∂t
−R(t) R(t) + x
−R(t)
−R(t)

(4.22)

The first is zero because m is constant away from the interface. The third term is proportional
to the difference in volume energy in the membrane compared to in the pore. It is equal to zero
in the pore and equal to a constant in the membrane which I will call −V . It is negative because
of the negative potential energy from lipid attraction. Since (m )2 is large near the interface and
zero away from it, the second term is equal to the expression in square brackets evaluated at the
interface (x = 0). Equation (4.22) becomes
1
∂R
+
+V = 0
R(t) + x ∂t x→0

(4.23)

which finally gives the equation of motion for the pore radius
dR −1
=
−V
dt
R

(4.24)

Remember that multiple material dependent constants were absorbed into the time scale. To
make this clear I will add constants to equation (4.24).
dR −A
=
−B
dt
R
74

(4.25)

where A is proportional to the line energy around the pore and B is proportional to the volume
energy inside the membrane. Since both A and B are positive, both terms are driving the pore shut.
Separating variables and integrating, gives
t
−R
dR =
dt
0
R0 BR + A
R

(4.26)

where R0 is the initial pore radius. The resulting equation is
1
A + BR
B(R0 − R) − Log
A + BR0
B2

=T

(4.27)

While this does not have an analytic solution, it can be solved numerically. To understand its
general behavior I will consider the two limiting cases where A or B are zero. B = 0 means that the
volume energy is the same both inside the membrane and the pore. This corresponds to a magnetic
system in which a region of spins with a certain orientation is surrounded by spins of the opposite
orientation. The energy within each region is the same and the dynamics are governed completely
by the curvature at the interface. Equation (4.25) reduces to the well known Allen-Cahn equation
[77] in two dimensions. In a general non-conserved system, a non-equilibrium interface will move
with a speed proportional to the mean curvature of the interface. In this limit, the equation of
motion for the pore radius is

R=

R2 − 2At
0

(4.28)

It is known that the size of domains in non-conserved systems scale as t 1/2 . The area of the
pore is

A p = A p0 − 2πAt

(4.29)

where A p0 is the initial pore area. The area decreases linearly with time.
If on the other hand we consider A = 0, we are completely neglecting the line energy around
the pore. The system is governed by the difference in volume energies. It behaves like a pressure
difference across the interface, closing the pore. The equation of motion becomes
75

R = R0 − Bt

(4.30)

In this case the radius decreases linearly and the area is
2
√
A p0 − πBt

Ap =

(4.31)

We can expect the actual solution to equation (4.25) to lie somewhere between the two extremes
of equations (4.29) and (4.31). In the next section I will show simulation results for pore healing.
The area does appear to decrease quadratically. I will now show how the fit parameters from a
second order polynomial fit can be related to the constants A and B in equation (4.25). If the initial
pore radius at time zero is R0 , then equation (4.25) can be used to find the initial first and second
derivatives.
−A
˙
R(0) =
−B
R0
˙
AR(0) −A A
¨
= 2
R(0) =
+B
R2
R0 R0
0
Now assume that the pore area can be approximately fit by

A p = C0t 2 +C1t +C2

(4.32)

(4.33)

(4.34)

C2 is clearly the starting pore area (πR2 ). The pore radius can be expressed as
0

R=

C0 2 C1
t + t + R2
0
π
π

(4.35)

The initial first and second derivatives of equation 4.35 are
C
˙
R(0) = 1
2πR0
1
¨
R(0) =
2π

C2
2C0
− 13
R0
2πR0
76

(4.36)

(4.37)

By setting equation (4.32) equal to (4.36) and equation (4.33) equal to (4.37), C0 and C1 can
be solved for in terms of A and B.

C0 = πB

A
+B
R0

C1 = −2π(A + BR0 )

(4.38)

(4.39)

A quadratic approximation of the solution to equation (4.25) is therefore

Ap = π B

A
+ B t 2 − 2(A + BR0 )t + R2
0
R0

(4.40)

In order to determine the relative contributions to healing, one can fit data to see if it is mostly
linear, as predicted by Allen-Cahn, or if there is a significant quadratic component.

4.2.2

Simulation results

I have run simulations on both a flat LB with periodic boundaries and a spherical LB vesicle and
observed how a pore heals over time. The Cooke lipid model [36, 37] was again used.

4.2.2.1

Flat LB

I began with an equilibrated 80x80 LB with periodic boundaries under zero tension. Unlike the
fixed area simulations in chapter 3 these simulations required free tensionless boundaries to allow
closing of the pore. A pore was made in the center by removing lipids within a radius of 8 σ (≈ 9
nm). That left tail beads initially exposed around the edge of the pore. Figure 4.1 shows three
different states during healing.
I calculated the area of the pore by dividing the LB into a 2D grid of bins approximately 1 σ 2
each. Any bin which doesn’t contain a bead contributes to the pore area. Figure 4.2 shows the pore
area calculated every 1 τ for 20 different runs. There is initially a rapid decrease in area as lipids

77

Figure 4.1: Snapshots of pore healing in a flat LB at three different times.

Figure 4.2: Pore area versus time in a flat lipid bilayer under zero tension. This is the raw data
from 20 different runs.

78

Figure 4.3: Averaged pore area versus time in a flat lipid bilayer under zero tension. The data in
figure 4.2 was averaged (red circles) and fit to a quadratic function (solid black line) as well as a
linear function (dashed black line).

curve around to coat the interior of the pore with lipid heads. This occurs within a few τ. There is
a large variation in closing times, which I define as when the pore area first becomes zero.
Figure 4.3 shows an average of the 20 runs. Each run is a piecewise function which is decreasing at early times and becomes zero beyond the closing time. I believe the averaged data contains
three distinct regions. There is a rapid decrease during the first 10 τ in which lipids are bending
around the edge. After this initial relaxation there is a more gradual decrease in area lasting about
50 τ. This middle region is most likely to correspond to the dynamics discussed in the previous
section. At late times there is a tail which I believe is mostly due to averaging data with different
closing times.

79

Flat LB
Vesicle

C0
0.017 ± 0.002
0.007 ± 0.002

C1
−3.3 ± 0.2
−3.5 ± 0.3

C2
162 ± 9
342 ± 5

Table 4.1: Parameters for the pore area versus time data, fit to quadratic equation (4.34).

In order to isolate the middle region, I removed the very early time data as well as the zeros
after the pore is closed before fitting each run separately to equation (4.34). The averages of those
fit parameters are shown in table 4.1. The averaged quadratic fit is represented by the solid black
line in figure 4.3.
While the coefficient for the quadratic term is not zero, the linear term is dominant throughout
the closing process. The average closing time, equal to the first root of equation (4.34), is 92 τ. The
time at which the quadratic and linear terms are equal is 198 τ. Therefore a linear approximation
should be reasonable. Figure 4.3 also shows the average linear fit to the data represented by the
dashed black line. The closing time predicted by the linear fit is 73 τ which is an underestimate of
about 20% compared to the quadratic fit.

4.2.2.2

Spherical vesicle

I also ran simulations on a spherical vesicle with an initial radius of about 20 σ (≈18 nm). I
created a pore by removing 10% of the lipids within a certain polar angle. Figure 4.4 shows the
initial configuration as well as an intermediate state and the fully healed vesicle.
I calculated the pore area by dividing the surface of the spherical vesicle into equal area bins
identified by their angular coordinates, θ and φ . Bins which did not contain any beads were
counted as part of the pore. Figure 4.5 shows the pore area versus time for 20 different runs.
Figure 4.6 shows the average of the data in figure 4.5. There are again three distinct regions:
an initial sharp drop, then a long gradual decline, and finally a rounded tail. I once more removed
the early time data and late time zeros from each run before fitting them to a quadratic function.
The averaged quadratic fit parameters are listed in table 4.1. This time the quadratic coefficient is

80

Figure 4.4: Snapshots of pore healing in a spherical vesicle at three different times.

Figure 4.5: Pore area versus time for a spherical vesicle. This is the raw data from 20 different
runs.

81

Figure 4.6: Averaged pore area versus time for a spherical vesicle. The data in figure 4.5 was
averaged and fit to both a linear and quadratic function.

even smaller. The estimated closing times from the quadratic and linear fits are 131 τ and 122 τ
respectively. The linear approximation still underestimates the closing time, but only by 7% in this
case.
It appears that the area of pores decrease linearly at early times and that the process can be
explained by a simple Allen-Cahn theory. The pore closing is driven primarily by the membrane
curvature around the pore. I believe there is a different process at work later on. Lipids which are
initially oriented with their heads pointing toward the center of the pore must reorient themselves.
The large variation in closing times among my simulation runs indicates that this is an activated
process. It is known that Allen-Cahn theory is only valid when the thickness of the interface, which

82

in this case is on the order of the LB thickness, is much smaller than the radius of curvature of the
pore. In my flat LB simulation, the radius of curvature is only about twice as large as the interface
thickness. In the vesicle simulation the ratio is closer to four, which may explain why the vesicle
data is more linear.
While Allen-Cahn theory can explain the healing dynamics of pores which are much wider
than the LB thickness, a different theory is needed to understand the final closing process.

4.2.2.3

Alternative fitting

While the quadratic and linear fits of the previous two sections are meant to describe the middle
region of decreasing area, I have also fit the entire averaged data sets shown in figures 4.3 and
4.6. In an attempt to capture the nonlinear early and late time behavior, I fit each data set to the
following two equations

A p = A0 · exp (−A1 · t) + A2 · t 2 + A3 · t

(4.41)

A p = B0 · exp (−B1 · t) + B2 · exp (−B3 · t)

(4.42)

which I refer to as “exponential plus quadratic” and “double exponential”. Figure 4.7 shows
the flat lipid bilayer data fit to both the exponential plus quadratic represented by the solid green
line and the double exponential represented by the dashed blue line. Tables 4.2 and 4.3 list the corresponding fit parameters. After finding an insignificant quadratic coefficient I performed another
fit excluding that term. Since equations (4.41) and (4.42) each contain four fit parameters they do
a comparable job fitting the data. The two fit lines are nearly overlapping throughout the entire
range.
Figure 4.8 shows the averaged vesicle data again fit to equations (4.41) and (4.42). While the
quadratic coefficient was rather small, including the term significantly improved the fit in this case.
Both new functions do a reasonable job fitting the middle and late time behavior, but still miss
the early time relaxation. Perhaps such models could be used to represent the entire evolution of
83

Figure 4.7: Averaged flat lipid bilayer pore data with alternative fits. The data in figure 4.2 was
averaged (red circles) and fit to a an exponential plus a linear term (solid green line) as well as a
double exponential (dashed blue line).

Flat LB
Vesicle

A0
175
366

A1
0.0257
0.00560

A2
N/A
0.00934

A3
-0.182
-2.43

Table 4.2: Parameters for the pore area versus time data, fit to an exponential plus quadratic equation (4.41).

a healing pore. However, I believe it would be much more useful to understand and model the
behavior in each of the three regions separately as part of a comprehensive model.

84

Flat LB
Vesicle

B0
9400
-1730

B1
0.0145
0.0322

B2
-9230
2090

B3
0.0142
0.0280

Table 4.3: Parameters for the pore area versus time data, fit to a double exponential equation (4.42).

Figure 4.8: Averaged vesicle pore data with alternative fits. The data in figure 4.5 was averaged
(red circles) and fit to an exponential plus a quadratic function (solid green line) as well as a double
exponential (dashed blue line).

4.2.2.4

Ion transfer

When considering the possibly harmful effects of pores, it is not only important to know how long
a pore remains open, but also what can happen while it is open. Cells maintain a difference in
ion concentration inside compared to outside. Knowing how the pore area changes with time,

85

we can calculate how much charge can pass through the pore. If we assume the difference in
ion concentration remains fairly constant during transfer, we can also assume a constant current
density through the pore, J. The total charge that passes through the pore is then
T

Q=
0

J · A p (t)dt

(4.43)

where T is closing time and A p (t) is the pore area for which we can use the linear or quadratic
approximations considered in the previous sections. I have calculated the total charge transfer
for both the flat LB and vesicle cases. Table 4.4 shows the results for the linear and quadratic
approximations.
Total charge transferred (J · σ 2 /τ)
Flat LB
Vesicle
Quadratic 5.3E3
19.8E3
Linear
5.2E3
19.6E3
Table 4.4: Total calculated charge transferred through pores using linear and quadratic approximations.

The differences in charge transfer between the linear and quadratic approximations are much
smaller compared to differences in closing times. One reason is that the underestimates of area
that take place in the linear approximations at early and late times, are compensated for by overestimates in the middle (see figures 4.2 and 4.6). Also, most of the charge transfer occurs at early
times when the decrease in area is mostly linear. The differences between the linear and quadratic
approximations are only 1.7% for the flat LB and 1.1% for the vesicle cases. It appears that the
linear approximation predicted by Allen-Cahn works quite well for calculating charge transfer.

86

Chapter 5
CONCLUSION

I have looked at the interactions between nanoparticles (NP) and lipid bilayers (LB). Simple continuum models can be used to determine both the thermodynamics of NP/LB states and the dynamics
of pore healing in LB. I used a coarse-grained lipid model to simulate the behavior in both cases.
In chapter 3 I looked at various equilibrium states resulting from NP/LB interactions. Using
a continuum theory, I defined a free energy for each possible state (figures 3.1-3.4). The theory
depends on the size and attraction strength of the NP and the macroscopic properties of the membrane such as the bending and stretching moduli. By comparing the free energies, I identified the
equilibrium states under different conditions and constructed a phase diagram (figures 3.6-3.8). I
varied the properties of the NP while fixing the membrane parameters, This was done for both
spherical and rod-shaped nanoparticles.
The driving force behind these interactions is the attraction between the NP and the lipid heads
and tails. This attraction favors wrapping of the NP by the LB. It can be wrapped by a double or
single layer of lipids depending on whether it is attracted more to the lipid heads or tails. Total
wrapping involves a bending energy which is independent of the NP size and is therefore the
dominant contribution for small NP. In that case the adhesion, which is proportional to the surface
area of the NP, can not overcome the cost of bending.
The stretching energy is proportional to the square of the NP surface area and therefore dominates for large NP. The stretching energy can be relieved by forming a pore in the LB at the cost
of line energy around the edge of the pore. Very small NP will not become wrapped due to the
overwhelming cost of bending. Slightly larger NP will become embedded in the LB. NP larger
than a critical size cause pore formation in the LB.
Ginzburg and Balijepalli studied these small NP interactions using density functional theory
[65]. In section 3.4 I compared their results with my continuum model predictions. The gen-

87

eral range of transitions is in decent agreement with my phase diagram (figure 3.9). The slight
discrepancies in the transition points are due to simplifications in my model which I explained.
In section 3.5 I used a coarse-grained particle model to simulate the interactions under the
same conditions. The NP sizes and attraction strengths I used were the same as those used by
Ginzburg and Balijepalli. The resulting states were reasonably well predicted by the continuum
model (figures 3.16 and 3.18).
I have compared three different techniques for studying NP/LB interactions which are all in
general agreement about the range of transitions. Mutual agreement lends credibility to all three
methods. Such a multiscale approach is a good way to study these kinds of systems. Simulations
which can quickly and easily be run over a wide range of conditions are an excellent supplement
to experiments which can be difficult and expensive. Continuum theory can identify the interesting
regions of phase space and help direct the conditions for simulation and experiment.
I have found a particular set of NP properties which result in absorption and pore formation.
Given the membrane properties of another system, this method could be used to predict the effects
of NP interaction.
In chapter 4 I looked at the dynamics of pore healing. The evolution of a pore is determined by
the line energy around the edge of the pore as well as the volume energy of the LB away from the
pore. The volume energy is related to the stretching or compression of the LB. A compressed LB
will want to spread out and will accelerate the closing of the pore. A stretched lipid bilayer will
instead work to increase the size of the pore. The line energy is always a positive energy related
to the bending modulus of the LB and is proportional to the circumference of the pore. The line
energy acts to close the pore.
I showed that the volume energy contribution will change the area of the pore quadratically with
time while the line energy contribution produces a linear decrease in pore area. If both energies are
present, the evolution will have both a linear and a quadratic component. Using the same coarsegrained lipid model, I simulated the closing of a pore in a flat LB and also in a spherical vesicle. By
fitting the pore area data, I was able to compare the relative importance of the linear and quadratic

88

contributions. In both the flat LB and the spherical vesicle cases, the pore area decreased quite
linearly. The quadratic fit determined the pore closing time more accurately than the linear fit.
The linear fit underestimated the closing time by 20% in the flat LB case and by only 7% in the
spherical vesicle case.
Charge flow through a pore is proportional to the time integral of the pore area curve. The
total charge flow predicted by the linear and quadratic fits agreed even better. The predicted charge
flows differed by less than 2% between the flat LB and spherical vesicle cases (table 4.4).

89

APPENDICES

90

APPENDIX A:

GINZBURG AND BALIJEPALLI EQUATIONS

Here I will elaborate on the TGMB equations used in section 2.5.1. I am mostly rehashing the
thorough explanation given in the paper by Ginzburg and Balijepalli [65]. Once again, the total
free energy can be written as
FND
= f1 + f2 + f3 + f4
kB T ρ0V

(1)

ϕ
QP
QD
QW
f1 = P ln
− ϕD ln
− ND ϕW ln
α
V
V
V


1
dr  ∑ χαβ ND φα (r)φβ (r) − ξ (r) 1 − ∑ φα (r) 
f2 =
V
α
α,β

f3 =

1
V

dr −

(2)

(3)



∑

wα (r)φα (r) − wP (r)ρP (r)

(4)

α(=P)

f4 =

1
V

¯
dr ρP (r)ΨCS φP (r)

(5)

In the first equation (1) ND is the degree of polymerization of the lipid diblocks or the number
of segments in each lipid, either hydrophilic (H) or lipophilic (L). Each segment has a Khun length
a =0.4 nm and therefore a volume of (ρ0 )−1 =0.064 nm3 . V is simply the system volume.
The fist term, f1 (2), contains the entropic free energies. ϕP , ϕD , and ϕW are the volume
fractions of the nanoparticle, lipids, and water respectively which add to one. QP , QD , and QW are
the individual partition functions for each component. α is the particle to diblock volume ratio
4πR3 ρ0
p
α=
3ND

(6)

where R p is the radius of the nanoparticle.
The second term, f2 (3), contains the local interactions determined by Flory-Huggins parameters [68, 69] as discussed in section 3.5.1. φα (r) represents the local volume fraction of each

91

component. The last part of (3) also enforces incompressibility. ξ (r) is a pressure field which I
will define later in this section.
The third term, f3 (4), introduces the chemical potential fields, w, which depend on the local
interactions. ρP (r) in (4) is the particle center density which is related to the particle’s local volume
¯
fraction, φP (r), and “smoothed density” function, φP (r), as follows
r−r
(RP /RgD )

(7)

r−r
(2RP /RgD )

(8)

1
ND

dr ρP (r )Θ 1 −

1
¯
φP (r) = d
2 ND

dr ρP (r )Θ 1 −

φP (r) =

d is the dimensionality of the simulation box. Θ is the Heaviside step function. RgD = a(ND /6)1/2
is the radius of gyration of the diblock.
Finally, the fourth term, f4 (5), contains the non-ideal hard-sphere interactions. It involves a
density functional theory approximation used by Tarazona [80]. The hard-sphere fluid free energy
is given by the Carnahan-Starling [81] equation of state
4x − 3x2
ΨCS (x) =
(1 − x)2

(9)

The partition functions used in f1 (2) are given by

QW =

exp {−wW (r)} dr

(10)

QP =

exp {−wP (r)} dr

(11)

QD =

q (r, 1) dr

(12)

The propagator q (r, s) and its counterpart q † (r, s) are defined by the following diffusion equations
∂ q (r, s)
= ∇2 − wt(s) q (r, s)
∂s
∂ q † (r, s)
= ∇2 − wt(s) q † (r, s)
∂s

92

(13)
(14)

s is an index denoting the position along the diblock which ranges from 0 to 1. t(s) = L if s < f ,
and H otherwise. These equations are subject to the boundary conditions q (r, 0) = q † (r, 1) = 1.
The chemical potential fields, w are given by
wH (r) = χHL ND φL (r) + χHP ND φP (r) + ξ (r)

(15)

wL (r) = χHL ND φH (r) + χW L ND φW (r) + ξ (r)

(16)

wW (r) = χW L ND φL (r) + χW P ND φP (r) + ξ (r)

(17)

1
¯
wP (r) = ΨCS φP (r) +
ND

dr Θ 1 −

r−r
(RP /RgD )

×

χHP ND φH (r ) + χW P ND φW (r ) + ξ (r ) +
1
2d ND

dr Θ 1 −

r−r
(2RP /RgD )

¯
× ρP (r )ΨCS φP (r )

The pressure field, ξ (r), is given by
ξ (r) = P(r) + ε (φP (r) + φW (r) + φH (r) + φL (r) − 1)

(18)

H (r) + HH (r) + HL (r)
P(r) = W
3

(19)

where HH (r), HL (r), and HW (r) are the right hand sides of equations (15-17) except for the ξ
terms. ε is a heuristic parameter which Ginzburg and Balijepalli found to be convergent when set
equal to 60.
The algorithm used by Ginzburg and Balijepalli was adapted from the one used by Drolet and
Fredrickson [82] and in the original TGMB paper [66]. The first step is to initialize the chemical
potential and pressure fields. Those can then be used to calculate the density fields with the help
of equations (7) and (8) as well as these relations
ϕ V
ρP (r) = P
exp [−wP (r)]
α QP

(20)

V
exp [−wW (r)]
QP

(21)

f
V
q(r, s)q † (r, s)ds
QD 0

(22)

φW (r) = ϕW
φL (r) = ϕD

93

φH (r) = ϕD

1
V
>
q(r, s)q † (r, s)ds
QD
f

φH (r) + φL (r) + φP (r) + φW (r) = 1

(23)
(24)

Then the chemical potentials are updated according to
wt+1 (r) = (1 − λi )wt (r) + λi µi (r)
i
i

(25)

where t is the iteration number. The index i refers to H, L, W , or P. µi are the new chemical
potentials calculated from equations (15-18). λi acts as an effective time step. Ginzburg and
Balijepalli set λP = 0 in order to maintain a constant particle field while setting all other λi = 0.025.
Finally, the pressure field is updated by equations (18) and (19). These steps are iterated until the
fields converge.

94

APPENDIX B:

NOSÉ-HOOVER THERMOSTAT/BAROSTAT

My simulations were run in the isothermal-isobaric ensemble fixing the number of particles, temperature, and pressure (NPT). The natural ensemble for molecular dynamics is the microcanonical
ensemble in which the number of particles, volume, and total energy (NVE) are conserved. The
isothermal-isobaric ensemble is more desirable because it can be directly compared to experiments
in which temperature and pressure are controlled. The temperature and pressure are fixed using a
Nosé-Hoover thermostat and barostat [73, 74] which I will explain now.
The standard Hamiltonian consists of kinetic energy including coordinates q i and masses mi
as well as potential energy φ (q ). In this system the total energy will be conserved while the
kinetic energy (temperature) and pressure are free to fluctuate. Nosé and Hoover introduced a
different Hamiltonian which allows the total energy to fluctuate while conserving the temperature
and pressure. The temperature and pressure can be included separately and I will discuss the
thermostat first.

Thermostat
The Nosé-Hoover formulation transforms the real variables (q ,p ,t ) where p are the momenta, to
virtual variables (q,p,t) through the following transformations

qi = qi

(26)

pi = pi /s

(27)

t

t =

dt
s

(28)

Nosé interpreted these transformations as a rescaling of time by dt = dt/s. The new Hamiltonian in terms of the virtual variables is

95

p2
pi 2
+ φ (q) + s + gkT ln s
(29)
2Q
2mi s2
i
where g is the number of degrees of freedom, k is the Boltzmann constant, and T is the fixed
H =∑

external temperature. ps is the conjugate momentum for s. Q behaves like an effective mass for the
variable s. It determines the strength of the coupling between the system and the heat bath and sets
a time scale for temperature fluctuations. A small value for Q corresponds to strong coupling and
causes frequent velocity rescaling. Extremely small values can make the simulation unstable. A
large value for Q corresponds to weak coupling and can cause very slow temperature equilibration.
It has been shown that the most efficient equilibration occurs when the heat bath fluctuations are in
resonance with the natural system fluctuations. Nosé derived an estimate of the heat bath frequency
[83, 84] which I will now show. The Hamiltonian equations of motion area
p
dqi ∂ H
=
= i2
dt
∂ pi mi s
∂H
dpi
∂φ
=−
=−
dt
∂ qi
∂ qi

(30)
(31)

ps
ds ∂ H
=
=
dt
∂ ps
Q
d ps
∂H
=−
=
dt
∂s

p2
∑ miis2 − gkT
i

(32)
/s

(33)

Setting the time derivative of (32) equal to (33) a second order differential equation for s is
obtained.

Q

p2

d2s
∂H
=−
=
2
∂s
dt

∑ miis2 − gkT

/s

(34)

i

When the system is in equilibrium s will fluctuate around its average value < s >. We can then
replace s with < s > +δ s where δ s contains the time dependence. A linear expansion of (34) gives

Q

d2δ s
1
=
<s>
dt 2

p2
∑ mi < is >2 − gkT
i

96

−

δs
< s >2

p2
i
− gkT
mi < s >2
i

3∑

(35)

For small values of Q or strong coupling to the heat bath the variable s will fluctuate more
quickly than the particles in the system. The temperature can then be defined as

gkT = ∑

p2
i
2
i mi < s >

(36)

−2gkT
d2δ s
=
δs
2
dt
Q < s >2

(37)

and (35) reduces to

This is simply a harmonic oscillator equation whose angular frequency is

ω=

2gkT
Q < s >2

(38)

The period of fluctuations in s is therefore

τ=

2π
=π <s>
ω

2Q
gkT

(39)

Nosé found through simulations that < s > was on the order of one [83]. Multiple studies have
been done to determine the effect of the parameter Q on simulations [83, 84, 75]. Fluctuations in
s show up in the temperature and kinetic energy and can be on a much different time scale than
the natural system fluctuations. If the frequency of s fluctuations are significantly larger than that
of the system, it can result in an isolated mode and cause a non-Gaussian distribution of kinetic
energy.
Figure B.1 shows a histogram of the kinetic energy of a simulated lipid bilayer for a small value
of Q. Notice the slight dip at the center. This bimodal distribution was predicted by Nosé [84].
Extremely small values of Q can even make the simulations unstable, similar to using too large of
a time step. Large values for Q correspond to weak coupling to the heat bath and result in long
equilibration times. The optimal value for Q is one which places the fluctuations of s in resonance
with those of the system. It will maintain a Gaussian distribution while being able to equilibrate in
the least amount of time.

97

Figure B.1: Histogram of kinetic energy for strong heat bath coupling (Q = 0.648, T damp =
0.005). The scaling variable s fluctuates faster than the system and causes a bimodal, non-Gaussian
distribution.

The characteristic timescale of the system can be estimated from the potential. For small
oscillations about the equilibrium positions the motion will appear harmonic. The spring constant
for an equivalent harmonic potential can be found by taking the second derivative of the original
potential and evaluating it at the equilibrium position. The potential between tail beads in the
Cooke model is given by the sum of equations (2.1) and (2.4) and is represented in figure 2.1. The
second derivative of that potential evaluated at r = 21/6 ∗ b is approximately 57 ε/σ 2 . The period
then follows as

98

Figure B.2: Temperature versus time for various T damp.

τ = 2π

m
≈ 0.83
k

(40)

The molecular dynamics package, LAMMPS, which I’ve used for my simulations uses a parameter called T damp which is related to Q by

T damp =

Q
gkt

(41)

They define T damp as a rough timescale for temperature equilibration. It is directly proportional to the period derived in (39). Using the period in (40) implies a T damp value of around
0.2.

99

I ran simulations of a 40x40 lipid bilayer at a temperature of kT = ε. I varied T damp from
0.02 to 2 to observe the effects of the damping parameter. The units for T damp are the standard
Lennard-Jones time units (

mσ 2 ).
ε

There were 9600 particles making g = 28800 degrees of free-

dom. Figure B.2 shows temperature versus time for three different values of T damp. All three
simulations were run at the same temperature and are simply displaced for clarity.
The small features in the top two traces are due to the normal fluctuations of the system. The
larger features in the top two traces are due to fluctuations of the heat bath. The bottom trace
(T damp = 0.02) contains fluctuations which are faster than the normal system fluctuations and
could prevent equilibration. In the middle trace (T damp = 0.2) the heat bath fluctuations are on
the same order though slightly longer than the system fluctuations which is ideal. The system will
be able to equilibrate and the temperature and kinetic energy will have the proper distribution. The
heat bath coupling used in the top trace (T damp = 2) would also work but the system would take
ten times longer to equilibrate. The LAMMPS documentation recommends trying a T damp value
equal to 100 time steps. Since my time step was 0.01, that would suggest a T damp equal to 1
which I believe is larger than necessary in this case. Again, it would produce reasonable results,
but simulations would have to be run for much longer. All simulations shown in other chapters
were run with T damp = 0.2.

Barostat
Pressure can also be held constant through use of a barostat. The method was originally used by
Andersen [85] and involves a rescaling of the coordinates by V 1/3 , where V is the volume of the
simulation box. Similar to equations (26-28)

qi = qiV 1/3
pi =

pi

sV 1/3
t dt
t =
s
100

(42)
(43)
(44)

The Hamiltonian for the isothermal/isobaric ensemble is
p2
p2
(45)
+ φ (qV 1/3 ) + s + gkT ln s + V + PexV
2Q
2W
2mi s2V 1/3
i
where Pex is the externally applied pressure and pV is the conjugate momentum of V . The fifth
H =∑

pi 2

term on the right looks like a kinetic energy making W act as an effective mass similar to Q. It
determines the strength of the pressure coupling and hence the frequency of volume fluctuations.
W is sometimes referred to as the mass of an imagined piston that is controlling the pressure. While
a normal piston compresses the volume along one dimension, these volume changes occur in all
three dimensions.
Care must be taken in choosing a value for W . A value which is too small can again cause
rapid, independent oscillation that produce a non-Gaussian ensemble distribution. Fluctuations
that are short relative to the time step of the simulation can also cause instability. Similar to the
derivation of temperature fluctuations above, Nosé and Klein [86] derived a relation between the
period of volume fluctuations and W as follows

τ0 = 2π

W
3LB

(46)

where L is the length of the simulation box and B is the bulk modulus. In the same paper Nosé
and Klein offer another approximation for the fluctuation period

τ0 = 2π

L
v1

(47)

where v1 is the speed of sound in the simulated material. If equations (46) and (47) are set
equal, W can be solved for.

W=

3BL3
v2
1

(48)

This can be further simplified because the bulk modulus and speed of sound are related as
follows

101

B=ρ

∂P
∂ρ

(49)

v2 =
1

∂P
∂ρ

(50)

where P is pressure and ρ is the mass density. Therefore

B = ρv2
1

(51)

Equation 48 can then be simplified to

W = 3L3 ρ =

3L3 mN
V

(52)

where m is the mass per particle which is equal to one. L is the length of the simulation box
along which pressure waves propagate. In my simulations I’m interested in fluctuations parallel to
the lipid bilayer surface or the xy plane. I also use square bilayers so that

L = Lx = Ly

(53)

The volume of the bilayer is equal to the area times the height or thickness of the bilayer (h)

V = L2 h

(54)

3NL
h

(55)

W can once again be reduced to

W=

LAMMPS sets W through another parameter called Pdamp defined as

Pdamp =

W
NkT

(56)

They describe it as the timescale for pressure equilibration. Substituting (55) into (56) gives

102

Pdamp =

3L
kT h

The units for Pdamp are again the standard Lennard-Jones time units (

(57)
mσ 2 ).
ε

Using the

approximate dimensions of my bilayer (40x40x5 σ 3 ) and my temperature of kT =0.9 ε implies
that Pdamp should be around 5. That corresponds to 1000 time steps which is exactly what the
LAMMPS documentation recommends. Since all of my simulations use the same size bilayer, I
have used this Pdamp value for all of them.

103

BIBLIOGRAPHY

104

BIBLIOGRAPHY

[1] Wen-Che Hou, Babak Yaghoubi Moghadam, Charlie Corredor, Paul Westerhoff, and
Jonathan D Posner. Distribution of functionalized gold nanoparticles between water and lipid
bilayers as model cell membranes. Environmental science & technology, 46(3):1869–76,
February 2012.
[2] Yuri Roiter, Maryna Ornatska, Aravind R Rammohan, Jitendra Balakrishnan, David R Heine,
and Sergiy Minko. Interaction of lipid membrane with nanostructured surfaces. Langmuir :
the ACS journal of surfaces and colloids, 25(11):6287–99, June 2009.
[3] Yuri Roiter, Maryna Ornatska, Aravind R Rammohan, Jitendra Balakrishnan, David R Heine,
and Sergiy Minko. Interaction of nanoparticles with lipid membrane. Nano letters, 8(3):941–
4, March 2008.
[4] Pascale R Leroueil, Stephanie a Berry, Kristen Duthie, Gang Han, Vincent M Rotello,
Daniel Q McNerny, James R Baker, Bradford G Orr, and Mark M Banaszak Holl. Wide
varieties of cationic nanoparticles induce defects in supported lipid bilayers. Nano letters,
8(2):420–4, February 2008.
[5] Pascale R Leroueil, Seungpyo Hong, Almut Mecke, James R Baker, Bradford G Orr, and
Mark M Banaszak Holl. Nanoparticle interaction with biological membranes: does nanotechnology present a Janus face? Accounts of chemical research, 40(5):335–42, May 2007.
[6] Seungpyo Hong, Pascale R Leroueil, Elizabeth K Janus, Jennifer L Peters, Mary-Margaret
Kober, Mohammad T Islam, Bradford G Orr, James R Baker, and Mark M Banaszak Holl.
Interaction of polycationic polymers with supported lipid bilayers and cells: nanoscale hole
formation and enhanced membrane permeability. Bioconjugate chemistry, 17(3):728–34,
2006.
[7] S a Klein, S J Wilk, T J Thornton, and J D Posner. Formation of nanopores in suspended lipid
bilayers using quantum dots. Journal of Physics: Conference Series, 109:012022, March
2008.
[8] Maurits R R de Planque, Sara Aghdaei, Tiina Roose, and Hywel Morgan. Electrophysiological Characterization of Membrane Disruption by Nanoparticles - Supporting Information.
ACS nano, pages 8–10, April 2011.
[9] Nastassja Lewinski, Vicki Colvin, and Rebekah Drezek. Cytotoxicity of nanoparticles. Small
(Weinheim an der Bergstrasse, Germany), 4(1):26–49, January 2008.
[10] Wim H De Jong and Paul J a Borm. Drug delivery and nanoparticles:applications and hazards.
International journal of nanomedicine, 3(2):133–49, January 2008.
[11] M Wang and M Thanou. Targeting nanoparticles to cancer. Pharmacological research : the
official journal of the Italian Pharmacological Society, 62(2):90–9, August 2010.

105

[12] Alexander P. Lyubartsev and Alexander L. Rabinovich. Recent development in computer
simulations of lipid bilayers. Soft Matter, 7(1):25, 2011.
[13] Jan Hermans, Herman J. C. Berendsen, Wilfred F. Van Gunsteren, and Johan P. M. Postma. A
consistent empirical potential for water-protein interactions. Biopolymers, 23(8):1513–1518,
August 1984.
[14] W F van Gunsteren, SR Billeter, AA Eising, PH Hunenberger, P Kruger, AE Mark, WRP
Scott, and IG Tironi. Biomolecular Simulation: The GROMOS96 Manual and User Guide.
Hochschulverlag AG, ETH Zurich, 1996.
[15] Nathan Schmid, Andreas P Eichenberger, Alexandra Choutko, Sereina Riniker, Moritz
Winger, Alan E Mark, and Wilfred F van Gunsteren. Definition and testing of the GROMOS
force-field versions 54A7 and 54B7. European biophysics journal : EBJ, 40(7):843–56, July
2011.
[16] David Poger, Wilfred F Van Gunsteren, and Alan E Mark. A new force field for simulating phosphatidylcholine bilayers. Journal of computational chemistry, 31(6):1117–25, April
2010.
[17] Alexander D. MacKerell, Joanna Wiorkiewicz-Kuczera, and Martin Karplus. An all-atom
empirical energy function for the simulation of nucleic acids. Journal of the American Chemical Society, 117(48):11946–11975, December 1995.
[18] A. D. MacKerell, D Bashford, R. L. Dunbrack, J D Evanseck, M J Field, S Fischer, J Gao,
H Guo, S Ha, D. Joseph-McCarthy, L Kuchnir, K Kuczera, F T K Lau, C Mattos, S Michnick,
T Ngo, D T Nguyen, B Prodhom, W E Reiher, B Roux, M Schlenkrich, J C Smith, R Stote,
J Straub, M Watanabe, J. Wiórkiewicz-Kuczera, D Yin, and M Karplus. All-Atom Empirical
Potential for Molecular Modeling and Dynamics Studies of Proteins. The Journal of Physical
Chemistry B, 102(18):3586–3616, April 1998.
[19] Jeffery B Klauda, Richard M Venable, J Alfredo Freites, Joseph W O’Connor, Douglas J
Tobias, Carlos Mondragon-Ramirez, Igor Vorobyov, Alexander D MacKerell, and Richard W
Pastor. Update of the CHARMM all-atom additive force field for lipids: validation on six lipid
types. The journal of physical chemistry. B, 114(23):7830–43, June 2010.
[20] Scott J Weiner, Peter A Kollman, David A Case, U Chandra Singh, Caterina Ghio, Guliano
Alagona, Salvatore Profeta, and Paul Weiner. A new force field for molecular mechanical simulation of nucleic acids and proteins. Journal of the American Chemical Society,
106(3):765–784, February 1984.
[21] Junmei Wang, Romain M Wolf, James W Caldwell, Peter a Kollman, and David a Case.
Development and testing of a general amber force field. Journal of computational chemistry,
25(9):1157–74, July 2004.
[22] Balázs Jójárt and Tamás A Martinek. Performance of the general amber force field in modeling aqueous POPC membrane bilayers. Journal of computational chemistry, 28(12):2051–8,
September 2007.

106

[23] D Peter Tieleman. The molecular basis of electroporation. BMC biochemistry, 5:10, July
2004.
[24] Mounir Tarek. Membrane electroporation: a molecular dynamics simulation. Biophysical
journal, 88(6):4045–53, June 2005.
[25] Steve O Nielsen, Carlos F Lopez, Goundla Srinivas, and Michael L Klein. Coarse grain
models and the computer simulation of soft materials. Journal of Physics: Condensed Matter,
16(15):R481–R512, April 2004.
[26] Grace Brannigan, Lawrence C-L Lin, and Frank L H Brown. Implicit solvent simulation
models for biomembranes. European biophysics journal : EBJ, 35(2):104–24, January 2006.
[27] M Muller, K Katsov, and M Schick. Biological and synthetic membranes: What can be
learned from a coarse-grained description? Physics Reports, 434(5-6):113–176, November
2006.
[28] M Venturoli, M Maddalenasperotto, M Kranenburg, and B Smit. Mesoscopic models of
biological membranes. Physics Reports, 437(1-2):1–54, December 2006.
[29] Gary S Ayton, Will G Noid, and Gregory A Voth. Systematic Coarse Graining of Biomolecular and Soft-Matter Systems. MRS Bulletin, 32(November):929–934, 2007.
[30] Sandra V Bennun, Matthew I Hoopes, Chenyue Xing, and Roland Faller. Coarse-grained
modeling of lipids. Chemistry and physics of lipids, 159(2):59–66, June 2009.
[31] T. Murtola, A. Bunker, I. Vattulainen, M. Deserno, and M. Karttunen. Multiscale modeling
of emergent materials: biological and soft matter. Physical Chemistry Chemical Physics,
11(12):1869–1892, 2009.
[32] Gary Ayton and Gregory a Voth. Bridging microscopic and mesoscopic simulations of lipid
bilayers. Biophysical journal, 83(6):3357–70, December 2002.
[33] J M Drouffe, a C Maggs, and S Leibler. Computer simulations of self-assembled membranes.
Science (New York, N.Y.), 254(5036):1353–6, November 1991.
[34] Oded Farago. “Water-free” computer model for fluid bilayer membranes. The Journal of
Chemical Physics, 119(1):596, 2003.
[35] Grace Brannigan and Frank L H Brown. Solvent-free simulations of fluid membrane bilayers.
The Journal of chemical physics, 120(2):1059–71, January 2004.
[36] Ira Cooke, Kurt Kremer, and Markus Deserno. Tunable generic model for fluid bilayer membranes. Physical Review E, 72(1):011506, July 2005.
[37] Ira R Cooke and Markus Deserno. Solvent-free model for self-assembling fluid bilayer membranes: stabilization of the fluid phase based on broad attractive tail potentials. The Journal
of chemical physics, 123(22):224710, December 2005.

107

[38] Berk Hess, S León, N Van Der Vegt, and Kurt Kremer. Long time atomistic polymer trajectories from coarse grained simulations: bisphenol-A polycarbonate. Soft Matter, 2(5):409,
2006.
[39] V. A. Harmandaris, N. P. Adhikari, N. F. A. van der Vegt, K. Kremer, B. A. Mann, R. Voelkel,
H. Weiss, and CheeChin Liew. Ethylbenzene Diffusion in Polystyrene: United Atom Atomistic/Coarse Grained Simulations and Experiments. Macromolecules, 40(19):7026–7035,
September 2007.
[40] Markus Deserno. Mesoscopic Membrane Physics: Concepts, Simulations, and Selected Applications. Macromolecular Rapid Communications, 30(9-10):752–771, May 2009.
[41] P F Fahey and W W Webb. Lateral diffusion in phospholipid bilayer membranes and multilamellar liquid crystals. Biochemistry, 17(15):3046–53, July 1978.
[42] Rudiger Goetz and Reinhard Lipowsky. Computer simulations of bilayer membranes: Selfassembly and interfacial tension. The Journal of Chemical Physics, 108(17):7397, 1998.
[43] Mark Stevens, Jan Hoh, and Thomas Woolf. Insights into the Molecular Mechanism of
Membrane Fusion from Simulation: Evidence for the Association of Splayed Tails. Physical
Review Letters, 91(18):188102, October 2003.
[44] C. Loison, M. Mareschal, K. Kremer, and F. Schmid. Thermal fluctuations in a lamellar
phase of a binary amphiphile–solvent mixture: A molecular-dynamics study. The Journal of
Chemical Physics, 119(24):13138, 2003.
[45] C Loison, M Mareschal, and F Schmid. Pores in bilayer membranes of amphiphilic
molecules: coarse-grained molecular dynamics simulations compared with simple mesoscopic models. The Journal of chemical physics, 121(4):1890–900, July 2004.
[46] Siewert J Marrink, H Jelger Risselada, Serge Yefimov, D Peter Tieleman, and Alex H
de Vries. The MARTINI force field: coarse grained model for biomolecular simulations.
The journal of physical chemistry. B, 111(27):7812–24, July 2007.
[47] Siewert J. Marrink, Alex H. de Vries, and Alan E. Mark. Coarse Grained Model for Semiquantitative Lipid Simulations. The Journal of Physical Chemistry B, 108(2):750–760, January 2004.
[48] John C. Shelley, Mee Y. Shelley, Robert C. Reeder, Sanjoy Bandyopadhyay, and Michael L.
Klein. A Coarse Grain Model for Phospholipid Simulations. The Journal of Physical Chemistry B, 105(19):4464–4470, May 2001.
[49] John C. Shelley, Mee Y. Shelley, Robert C. Reeder, Sanjoy Bandyopadhyay, Preston B.
Moore, and Michael L. Klein. Simulations of Phospholipids Using a Coarse Grain Model.
The Journal of Physical Chemistry B, 105(40):9785–9792, October 2001.
[50] P. J Hoogerbrugge and J. M. V. A Koelman. Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics. Europhysics Letters (EPL), 19(3):155–160, June
1992.
108

[51] J. M. V. A Koelman and P. J Hoogerbrugge. Dynamic Simulations of Hard-Sphere Suspensions Under Steady Shear. Europhysics Letters (EPL), 21(3):363–368, January 1993.
[52] R D Groot and K L Rabone. Mesoscopic simulation of cell membrane damage, morphology
change and rupture by nonionic surfactants. Biophysical journal, 81(2):725–36, August 2001.
[53] Hwankyu Lee and R.G. Larson. Membrane Pore Formation Induced by Acetylated and
Polyethylene Glycol-Conjugated Polyamidoamine Dendrimers. The Journal of Physical
Chemistry C, pages 5316–5322, 2011.
[54] Shikha Nangia and Radhakrishna Sureshkumar. Effects of Nanoparticle Charge and Shape
Anisotropy on Translocation through Cell Membranes. Langmuir : the ACS journal of surfaces and colloids, 28(51):17666–71, December 2012.
[55] P B Canham. The minimum energy of bending as a possible explanation of the biconcave
shape of the human red blood cell. Journal of theoretical biology, 26(1):61–81, January 1970.
[56] W. Helfrich. Elastic properties of lipid bilayers: theory and possible experiments. Z. Naturforsch, 28(11):693–703, 1973.
[57] Frank L H Brown. Elastic modeling of biomembranes and lipid bilayers. Annual review of
physical chemistry, 59:685–712, January 2008.
[58] E Evans. Bending Resistance and Chemically Induced Moments in Membrane Bilayers.
Biophysical Journal, 14(12):923–931, December 1974.
[59] Jejoong Yoo, Meyer B Jackson, and Qiang Cui. A comparison of coarse-grained and continuum models for membrane bending in lipid bilayer fusion pores. Biophysical journal,
104(4):841–52, February 2013.
[60] Markus Deserno and William M. Gelbart. Adhesion and Wrapping in ColloidVesicle Complexes. The Journal of Physical Chemistry B, 106(21):5543–5552, May 2002.
[61] Markus Deserno. Elastic deformation of a fluid membrane upon colloid binding. Physical
Review E, 69(3):031903, March 2004.
[62] Markus Deserno. When do fluid membranes engulf sticky colloids? Journal of Physics:
Condensed Matter, 16(22):S2061–S2070, June 2004.
[63] Teresa Ruiz-Herrero, Enrique Velasco, and Michael F Hagan. Mechanisms of budding of nanoscale particles through lipid bilayers. The journal of physical chemistry. B,
116(32):9595–603, August 2012.
[64] Reid C. Van Lehn and Alfredo Alexander-Katz. Penetration of lipid bilayers by nanoparticles
with environmentally-responsive surfaces: simulations and theory. Soft Matter, 7(24):11392,
2011.
[65] Valeriy V Ginzburg and Sudhakar Balijepalli. Modeling the thermodynamics of the interaction of nanoparticles with cell membranes. Nano letters, 7(12):3716–22, December 2007.

109

[66] R B Thompson, V V Ginzburg, M W Matsen, and a C Balazs. Predicting the mesophases
of copolymer-nanoparticle composites. Science (New York, N.Y.), 292(5526):2469–72, June
2001.
[67] Russell B. Thompson, Valeriy V. Ginzburg, Mark W. Matsen, and Anna C. Balazs. Block
Copolymer-Directed Assembly of Nanoparticles: Forming Mesoscopically Ordered Hybrid
Materials. Macromolecules, 35(3):1060–1071, January 2002.
[68] Maurice L. Huggins. Solutions of Long Chain Compounds. The Journal of Chemical Physics,
9(5):440, 1941.
[69] Paul J. Flory. Thermodynamics of High Polymer Solutions. The Journal of Chemical Physics,
9(8):660, March 1941.
[70] Yang Li and Ning Gu. Thermodynamics of charged nanoparticle adsorption on charge-neutral
membranes: a simulation study. The journal of physical chemistry. B, 114(8):2749–54, March
2010.
[71] Jiaqi Lin, Hongwu Zhang, Zhen Chen, and Yonggang Zheng. Penetration of lipid membranes
by gold nanoparticles: insights into cellular uptake, cytotoxicity, and their relationship. ACS
nano, 4(9):5421–9, September 2010.
[72] T V Tolpekina, W K den Otter, and W J Briels. Simulations of stable pores in membranes:
system size dependence and line tension. The Journal of chemical physics, 121(16):8014–20,
October 2004.
[73] Shuichi Nose. A unified formulation of the constant temperature molecular dynamics methods. The Journal of Chemical Physics, 81(1):511, 1984.
[74] WG Hoover. Canonical dynamics: Equilibrium phase-space distributions. Physical Review
A, 31(3):1695–1697, 1985.
[75] K. Cho and JD Joannopoulos. Ergodicity and dynamical properties of constant-temperature
molecular dynamics. Physical Review A, 45(10):7089, 1992.
[76] AA Griffith. The phenomena of rupture and flow in solids. Philosophical Transactions of the
Royal Society of London, A, 221:163–198, 1921.
[77] SM Allen and JW Cahn. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. Acta Metallurgica, 27:1085–1095, 1979.
[78] John W. Cahn and John E. Hilliard. Free Energy of a Nonuniform System. I. Interfacial Free
Energy. The Journal of Chemical Physics, 28(2):258, 1958.
[79] Pavel L. Krapivsky, Sidney Redner, and Eli Ben-Naim. A Kinetic View of Statistical Physics.
2010.
[80] P. Tarazona. A density functional theory of melting. Molecular physics, 52(1):81–96, August
1984.

110

[81] Norman F Carnahan and Kenneth E Starling. Equation of state for nonattracting rigid spheres.
The Journal of Chemical Physics, 51(2):635, 1969.
[82] F Drolet and Glenn H Fredrickson. Combinatorial screening of complex block copolymer
assembly with self-consistent field theory. Physical review letters, 83(1):4317–4320, 1999.
[83] Shichi Nosé. A molecular dynamics method for simulations in the canonical ensemble.
Molecular Physics, 52(2):255–268, June 1984.
[84] Shuichi Nosé. Constant temperature molecular dynamics methods. Progress of theoretical
physics. Supplement, (103):1–46, 1991.
[85] Hans C. Andersen. Molecular dynamics simulations at constant pressure and/or temperature.
The Journal of Chemical Physics, 72(4):2384, 1980.
[86] Shuichi Nosé and M.L. Klein. Constant pressure molecular dynamics for molecular systems.
Molecular Physics, 50(5):1055–1076, December 1983.

111