7 |
|
In the past 10 years, computer speeds have allowed for the atomistic |
8 |
|
simulation of phospholipid bilayers. These simulations have ranged |
9 |
|
from simulation of the gel phase ($L_{\beta}$) of |
10 |
< |
dipalmitoylphosphatidylcholine (DPPC), \cite{Lindahl:2000} to the |
10 |
> |
dipalmitoylphosphatidylcholine (DPPC),\cite{lindahl00} to the |
11 |
|
spontaneous aggregation of DPPC molecules into fluid phase |
12 |
< |
($L_{\alpha}$ bilayers. \cite{Marrinck:2001} With the exception of a |
13 |
< |
few ambitious |
14 |
< |
simulations,\cite{Marrinch:2001b,Marrinck:2002,Lindahl:2000} most |
12 |
> |
($L_{\alpha}$) bilayers.\cite{marrink01} With the exception of a few |
13 |
> |
ambitious |
14 |
> |
simulations,\cite{marrink01:undulation,marrink:2002,lindahl00} most |
15 |
|
investigations are limited to 64 to 256 |
16 |
< |
phospholipids.\cite{Lindal:2000,Sum:2003,Venable:2000,Gomez:2003,Smondyrev:1999,Marrinck:2001a} |
16 |
> |
phospholipids.\cite{lindahl00,sum:2003,venable00,gomez:2003,smondyrev:1999,marrink01} |
17 |
|
This is due to the expense of the computer calculations involved when |
18 |
|
performing these simulations. To properly hydrate a bilayer, one |
19 |
|
typically needs 25 water molecules for every lipid, bringing the total |
20 |
|
number of atoms simulated to roughly 8,000 for a system of 64 DPPC |
21 |
< |
molecules. Added to the difficluty is the electrostatic nature of the |
21 |
> |
molecules. Added to the difficulty is the electrostatic nature of the |
22 |
|
phospholipid head groups and water, requiring the computationally |
23 |
|
expensive Ewald sum or its slightly faster derivative particle mesh |
24 |
< |
Ewald sum.\cite{Nina:2002,Norberg:2000,Patra:2003} These factors all |
25 |
< |
limit the potential size and time lenghts of bilayer simulations. |
24 |
> |
Ewald sum.\cite{nina:2002,norberg:2000,patra:2003} These factors all |
25 |
> |
limit the potential size and time lengths of bilayer simulations. |
26 |
|
|
27 |
|
Unfortunately, much of biological interest happens on time and length |
28 |
< |
scales unfeasible with current simulation. One such example is the |
29 |
< |
observance of a ripple phase ($P_{\beta'}$) between the $L_{\beta}$ |
30 |
< |
and $L_{\alpha}$ phases of certain phospholipid |
31 |
< |
bilayers.\cite{Katsaras:2000,Sengupta:2000} These ripples are shown to |
28 |
> |
scales infeasible with current simulation. One such example is the |
29 |
> |
observance of a ripple phase ($P_{\beta^{\prime}}$) between the |
30 |
> |
$L_{\beta}$ and $L_{\alpha}$ phases of certain phospholipid |
31 |
> |
bilayers.\cite{katsaras00,sengupta00} These ripples are shown to |
32 |
|
have periodicity on the order of 100-200~$\mbox{\AA}$. A simulation on |
33 |
|
this length scale would have approximately 1,300 lipid molecules with |
34 |
|
an additional 25 water molecules per lipid to fully solvate the |
41 |
|
happens at pores. Some molecules of interest may incorporate |
42 |
|
themselves directly into the membrane. Once here, they may possess an |
43 |
|
appreciable waiting time (on the order of 10's to 100's of |
44 |
< |
nanoseconds) within the bilayer. Such long simulation times are |
44 |
> |
nanoseconds) within the bilayer. Such long simulation times are |
45 |
|
difficulty to obtain when integrating the system with atomistic |
46 |
|
detail. |
47 |
|
|
48 |
|
Addressing these issues, several schemes have been proposed. One |
49 |
< |
approach by Goetz and Liposky\cite{Goetz:1998} is to model the entire |
49 |
> |
approach by Goetz and Liposky\cite{goetz98} is to model the entire |
50 |
|
system as Lennard-Jones spheres. Phospholipids are represented by |
51 |
|
chains of beads with the top most beads identified as the head |
52 |
|
atoms. Polar and non-polar interactions are mimicked through |
53 |
|
attractive and soft-repulsive potentials respectively. A similar |
54 |
< |
model proposed by Marrinck \emph{et. al.}\cite{Marrinck:2004}~ uses a |
54 |
> |
model proposed by Marrinck \emph{et. al}.\cite{marrink04}~uses a |
55 |
|
similar technique for modeling polar and non-polar interactions with |
56 |
|
Lennard-Jones spheres. However, they also include charges on the head |
57 |
|
group spheres to mimic the electrostatic interactions of the |
63 |
|
interactions than the previous two models, while still balancing the |
64 |
|
need for simplifications over atomistic detail. The model uses |
65 |
|
Lennard-Jones spheres for the head and tail groups of the |
66 |
< |
phopholipids, allowing for the ability to scale the parameters to |
66 |
> |
phospholipids, allowing for the ability to scale the parameters to |
67 |
|
reflect various sized chain configurations while keeping the number of |
68 |
|
interactions small. What sets this model apart, however, is the use |
69 |
< |
of dipoles to represent the electrosttaic nature of the |
69 |
> |
of dipoles to represent the electrostatic nature of the |
70 |
|
phospholipids. The dipole electrostatic interaction is shorter range |
71 |
< |
than coulombic ($\frac{1}{r^3}$ versus $\frac{1}{r}$), eliminating the |
71 |
> |
than Coulombic ($\frac{1}{r^3}$ versus $\frac{1}{r}$), eliminating the |
72 |
|
need for a costly Ewald sum. |
73 |
|
|
74 |
|
Another key feature of this model, is the use of a dipolar water model |
75 |
< |
to represent the solvent. The soft sticky dipole ({\scssd}) |
76 |
< |
water \cite{Liu:1996a} relies on the dipole for long range |
77 |
< |
electrostatic effects, butalso contains a short range correction for |
78 |
< |
hydrogen bonding. In this way the systems in this research mimic the |
79 |
< |
entropic contribution to the hydrophobic effect due to hydrogen-bond |
80 |
< |
network deformation around a non-polar entity, \emph{i.e.}~ the |
81 |
< |
phospholipid. |
75 |
> |
to represent the solvent. The soft sticky dipole ({\sc ssd}) water |
76 |
> |
\cite{liu96:new_model} relies on the dipole for long range electrostatic |
77 |
> |
effects, but also contains a short range correction for hydrogen |
78 |
> |
bonding. In this way the systems in this research mimic the entropic |
79 |
> |
contribution to the hydrophobic effect due to hydrogen-bond network |
80 |
> |
deformation around a non-polar entity, \emph{i.e.}~the phospholipid. |
81 |
|
|
82 |
|
The following is an outline of this chapter. |
83 |
< |
Sec.~\ref{lipoidSec:Methods} is an introduction to the lipid model |
83 |
> |
Sec.~\ref{lipidSec:Methods} is an introduction to the lipid model |
84 |
|
used in these simulations. As well as clarification about the water |
85 |
|
model and integration techniques. The various simulation setups |
86 |
|
explored in this research are outlined in |
92 |
|
|
93 |
|
\section{\label{lipidSec:Methods}Methods} |
94 |
|
|
96 |
– |
|
97 |
– |
|
95 |
|
\subsection{\label{lipidSec:lipidModel}The Lipid Model} |
96 |
|
|
97 |
|
\begin{figure} |
98 |
< |
|
99 |
< |
\caption{Schematic diagram of the single chain phospholipid model} |
100 |
< |
|
98 |
> |
\centering |
99 |
> |
\includegraphics[width=\linewidth]{twoChainFig.eps} |
100 |
> |
\caption[The two chained lipid model]{Schematic diagram of the double chain phospholipid model. The head group (in red) has a point dipole, $\boldsymbol{\mu}$, located at its center of mass. The two chains are eight methylene groups in length.} |
101 |
|
\label{lipidFig:lipidModel} |
105 |
– |
|
102 |
|
\end{figure} |
103 |
|
|
104 |
|
The phospholipid model used in these simulations is based on the |
105 |
|
design illustrated in Fig.~\ref{lipidFig:lipidModel}. The head group |
106 |
|
of the phospholipid is replaced by a single Lennard-Jones sphere of |
107 |
< |
diameter $fix$, with $fix$ scaling the well depth of its van der Walls |
108 |
< |
interaction. This sphere also contains a single dipole of magnitude |
109 |
< |
$fix$, where $fix$ can be varied to mimic the charge separation of a |
107 |
> |
diameter $\sigma_{\text{head}}$, with $\epsilon_{\text{head}}$ scaling |
108 |
> |
the well depth of its van der Walls interaction. This sphere also |
109 |
> |
contains a single dipole of magnitude $|\boldsymbol{\mu}|$, where |
110 |
> |
$|\boldsymbol{\mu}|$ can be varied to mimic the charge separation of a |
111 |
|
given phospholipid head group. The atoms of the tail region are |
112 |
|
modeled by unified atom beads. They are free of partial charges or |
113 |
|
dipoles, containing only Lennard-Jones interaction sites at their |
114 |
|
centers of mass. As with the head groups, their potentials can be |
115 |
< |
scaled by $fix$ and $fix$. |
115 |
> |
scaled by $\sigma_{\text{tail}}$ and $\epsilon_{\text{tail}}$. |
116 |
|
|
117 |
|
The long range interactions between lipids are given by the following: |
118 |
|
\begin{equation} |
119 |
< |
EQ Here |
119 |
> |
V_{\text{LJ}}(r_{ij}) = |
120 |
> |
4\epsilon_{ij} \biggl[ |
121 |
> |
\biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12} |
122 |
> |
- \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6} |
123 |
> |
\biggr] |
124 |
|
\label{lipidEq:LJpot} |
125 |
|
\end{equation} |
126 |
|
and |
127 |
|
\begin{equation} |
128 |
< |
EQ Here |
128 |
> |
V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i}, |
129 |
> |
\boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[ |
130 |
> |
\boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j} |
131 |
> |
- |
132 |
> |
\frac{3(\boldsymbol{\hat{u}}_i \cdot \mathbf{r}_{ij}) % |
133 |
> |
(\boldsymbol{\hat{u}}_j \cdot \mathbf{r}_{ij}) } |
134 |
> |
{r^{2}_{ij}} \biggr] |
135 |
|
\label{lipidEq:dipolePot} |
136 |
|
\end{equation} |
137 |
|
Where $V_{\text{LJ}}$ is the Lennard-Jones potential and |
146 |
|
$i$, and $\boldsymbol{\hat{u}}_i$ is the standard unit orientation |
147 |
|
vector of $\boldsymbol{\Omega}_i$. |
148 |
|
|
149 |
< |
The model also allows for the bonded interactions of bonds, bends, and |
150 |
< |
torsions. The bonds between two beads on a chain are of fixed length, |
151 |
< |
and are maintained according to the {\sc rattle} algorithm. \cite{fix} |
149 |
> |
The model also allows for the bonded interactions bends, and torsions. |
150 |
> |
The bond between two beads on a chain is of fixed length, and is |
151 |
> |
maintained according to the {\sc rattle} algorithm.\cite{andersen83} |
152 |
|
The bends are subject to a harmonic potential: |
153 |
|
\begin{equation} |
154 |
< |
eq here |
154 |
> |
V_{\text{bend}}(\theta) = k_{\theta}( \theta - \theta_0 )^2 |
155 |
|
\label{lipidEq:bendPot} |
156 |
|
\end{equation} |
157 |
< |
where $fix$ scales the strength of the harmonic well, and $fix$ is the |
158 |
< |
angle between bond vectors $fix$ and $fix$. The torsion potential is |
159 |
< |
given by: |
157 |
> |
where $k_{\theta}$ scales the strength of the harmonic well, and |
158 |
> |
$\theta$ is the angle between bond vectors |
159 |
> |
(Fig.~\ref{lipidFig:lipidModel}). In addition, we have placed a |
160 |
> |
``ghost'' bend on the phospholipid head. The ghost bend adds a |
161 |
> |
potential to keep the dipole pointed along the bilayer surface, where |
162 |
> |
$theta$ is now the angle the dipole makes with respect to the {\sc |
163 |
> |
head}-$\text{{\sc ch}}_2$ bond vector. The torsion potential is given |
164 |
> |
by: |
165 |
|
\begin{equation} |
166 |
< |
eq here |
166 |
> |
V_{\text{torsion}}(\phi) = |
167 |
> |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0 |
168 |
|
\label{lipidEq:torsionPot} |
169 |
|
\end{equation} |
170 |
|
Here, the parameters $k_0$, $k_1$, $k_2$, and $k_3$ fit the cosine |
171 |
|
power series to the desired torsion potential surface, and $\phi$ is |
172 |
< |
the angle between bondvectors $fix$ and $fix$ along the vector $fix$ |
173 |
< |
(see Fig.:\ref{lipidFig:lipidModel}). Long range interactions such as |
174 |
< |
the Lennard-Jones potential are excluded for bead pairs involved in |
175 |
< |
the same bond, bend, or torsion. However, internal interactions not |
172 |
> |
the angle the two end atoms have rotated about the middle bond |
173 |
> |
(Fig.:\ref{lipidFig:lipidModel}). Long range interactions such as the |
174 |
> |
Lennard-Jones potential are excluded for atom pairs involved in the |
175 |
> |
same bond, bend, or torsion. However, internal interactions not |
176 |
|
directly involved in a bonded pair are calculated. |
177 |
|
|
178 |
|
All simulations presented here use a two chained lipid as pictured in |
179 |
< |
Fig.~\ref{lipidFig:twochain}. The chains are both eight beads long, |
179 |
> |
Fig.~\ref{lipidFig:lipidModel}. The chains are both eight beads long, |
180 |
|
and their mass and Lennard Jones parameters are summarized in |
181 |
|
Table~\ref{lipidTable:tcLJParams}. The magnitude of the dipole moment |
182 |
|
for the head bead is 10.6~Debye, and the bend and torsion parameters |
183 |
< |
are summarized in Table~\ref{lipidTable:teBTParams}. |
183 |
> |
are summarized in Table~\ref{lipidTable:tcBendParams} and |
184 |
> |
\ref{lipidTable:tcTorsionParams}. |
185 |
|
|
186 |
< |
\section{label{lipidSec:furtherMethod}Further Methodology} |
186 |
> |
\begin{table} |
187 |
> |
\caption{The Lennard Jones Parameters for the two chain phospholipids.} |
188 |
> |
\label{lipidTable:tcLJParams} |
189 |
> |
\begin{center} |
190 |
> |
\begin{tabular}{|l|c|c|c|} |
191 |
> |
\hline |
192 |
> |
& mass (amu) & $\sigma$($\mbox{\AA}$) & $\epsilon$ (kcal/mol) \\ \hline |
193 |
> |
{\sc head} & 72 & 4.0 & 0.185 \\ \hline |
194 |
> |
{\sc ch}\cite{Siepmann1998} & 13.02 & 4.0 & 0.0189 \\ \hline |
195 |
> |
$\text{{\sc ch}}_2$\cite{Siepmann1998} & 14.03 & 3.95 & 0.18 \\ \hline |
196 |
> |
$\text{{\sc ch}}_3$\cite{Siepmann1998} & 15.04 & 3.75 & 0.25 \\ \hline |
197 |
> |
{\sc ssd}\cite{liu96:new_model} & 18.03 & 3.051 & 0.152 \\ \hline |
198 |
> |
\end{tabular} |
199 |
> |
\end{center} |
200 |
> |
\end{table} |
201 |
|
|
202 |
+ |
\begin{table} |
203 |
+ |
\caption[Bend Parameters for the two chain phospholipids]{Bend Parameters for the two chain phospholipids. All alkane parameters are based off of those from TraPPE.\cite{Siepmann1998}} |
204 |
+ |
\label{lipidTable:tcBendParams} |
205 |
+ |
\begin{center} |
206 |
+ |
\begin{tabular}{|l|c|c|} |
207 |
+ |
\hline |
208 |
+ |
& $k_{\theta}$ ( kcal/($\text{mol deg}^2$) ) & $\theta_0$ ( deg ) \\ \hline |
209 |
+ |
{\sc ghost}-{\sc head}-$\text{{\sc ch}}_2$ & 0.00177 & 129.78 \\ \hline |
210 |
+ |
$x$-{\sc ch}-$y$ & 58.84 & 112.0 \\ \hline |
211 |
+ |
$x$-$\text{{\sc ch}}_2$-$y$ & 58.84 & 114.0 \\ \hline |
212 |
+ |
\end{tabular} |
213 |
+ |
\end{center} |
214 |
+ |
\end{table} |
215 |
+ |
|
216 |
+ |
\begin{table} |
217 |
+ |
\caption[Torsion Parameters for the two chain phospholipids]{Torsion Parameters for the two chain phospholipids. Alkane parameters based on TraPPE.\cite{Siepmann1998}} |
218 |
+ |
\label{lipidTable:tcTorsionParams} |
219 |
+ |
\begin{center} |
220 |
+ |
\begin{tabular}{|l|c|c|c|c|} |
221 |
+ |
\hline |
222 |
+ |
All are in kcal/mol $\rightarrow$ & $k_3$ & $k_2$ & $k_1$ & $k_0$ \\ \hline |
223 |
+ |
$x$-{\sc ch}-$y$-$z$ & 3.3254 & -0.4215 & -1.686 & 1.1661 \\ \hline |
224 |
+ |
$x$-$\text{{\sc ch}}_2$-$\text{{\sc ch}}_2$-$y$ & 5.9602 & -0.568 & -3.802 & 2.1586 \\ \hline |
225 |
+ |
\end{tabular} |
226 |
+ |
\end{center} |
227 |
+ |
\end{table} |
228 |
+ |
|
229 |
+ |
|
230 |
+ |
\section{\label{lipidSec:furtherMethod}Further Methodology} |
231 |
+ |
|
232 |
|
As mentioned previously, the water model used throughout these |
233 |
< |
simulations was the {\scssd} model of |
234 |
< |
Ichiye.\cite{liu:1996a,Liu:1996b,Chandra:1999} A discussion of the |
235 |
< |
model can be found in Sec.~\ref{oopseSec:SSD}. As for the integration |
236 |
< |
of the equations of motion, all simulations were performed in an |
237 |
< |
orthorhombic periodic box with a thermostat on velocities, and an |
238 |
< |
independent barostat on each cartesian axis $x$, $y$, and $z$. This |
239 |
< |
is the $\text{NPT}_{xyz}$. ensemble described in Sec.~\ref{oopseSec:Ensembles}. |
233 |
> |
simulations was the {\sc ssd} model of |
234 |
> |
Ichiye.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md} A |
235 |
> |
discussion of the model can be found in Sec.~\ref{oopseSec:SSD}. As |
236 |
> |
for the integration of the equations of motion, all simulations were |
237 |
> |
performed in an orthorhombic periodic box with a thermostat on |
238 |
> |
velocities, and an independent barostat on each Cartesian axis $x$, |
239 |
> |
$y$, and $z$. This is the $\text{NPT}_{xyz}$. ensemble described in |
240 |
> |
Sec.~\ref{oopseSec:Ensembles}. |
241 |
|
|
242 |
|
|
243 |
|
\subsection{\label{lipidSec:ExpSetup}Experimental Setup} |
284 |
|
|
285 |
|
Table ~\ref{lipidTable:simNames} summarizes the names and important |
286 |
|
details of the simulations. The B set of simulations were all started |
287 |
< |
in an ordered bilayer and observed over a period of 10~ns. Simulution |
287 |
> |
in an ordered bilayer and observed over a period of 10~ns. Simulation |
288 |
|
RL was integrated for approximately 20~ns starting from a random |
289 |
|
configuration as an example of spontaneous bilayer aggregation. |
290 |
|
Lastly, simulation RH was also started from a random configuration, |
291 |
|
but with a lesser water content and higher temperature to show the |
292 |
< |
spontaneous aggregation of an inverted hexagonal lamellar phase. |
292 |
> |
spontaneous aggregation of an inverted hexagonal lamellar phase. |