1 |
mmeineke |
971 |
|
2 |
|
|
|
3 |
|
|
\chapter{\label{chapt:lipid}Phospholipid Simulations} |
4 |
|
|
|
5 |
|
|
\section{\label{lipidSec:Intro}Introduction} |
6 |
|
|
|
7 |
mmeineke |
1001 |
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 |
mmeineke |
1061 |
dipalmitoylphosphatidylcholine (DPPC),\cite{lindahl00} to the |
11 |
mmeineke |
1001 |
spontaneous aggregation of DPPC molecules into fluid phase |
12 |
mmeineke |
1061 |
($L_{\alpha}$) bilayers.\cite{marrink01} With the exception of a few |
13 |
|
|
ambitious |
14 |
|
|
simulations,\cite{marrink01:undulation,marrink:2002,lindahl00} most |
15 |
mmeineke |
1001 |
investigations are limited to 64 to 256 |
16 |
mmeineke |
1061 |
phospholipids.\cite{lindahl00,sum:2003,venable00,gomez:2003,smondyrev:1999,marrink01} |
17 |
mmeineke |
1001 |
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 |
mmeineke |
1061 |
molecules. Added to the difficulty is the electrostatic nature of the |
22 |
mmeineke |
1001 |
phospholipid head groups and water, requiring the computationally |
23 |
|
|
expensive Ewald sum or its slightly faster derivative particle mesh |
24 |
mmeineke |
1061 |
Ewald sum.\cite{nina:2002,norberg:2000,patra:2003} These factors all |
25 |
|
|
limit the potential size and time lengths of bilayer simulations. |
26 |
mmeineke |
1001 |
|
27 |
|
|
Unfortunately, much of biological interest happens on time and length |
28 |
mmeineke |
1061 |
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 |
mmeineke |
1001 |
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 |
35 |
|
|
bilayer. A simulation of this size is impractical with current |
36 |
|
|
atomistic models. |
37 |
|
|
|
38 |
|
|
Another class of simulations to consider, are those dealing with the |
39 |
|
|
diffusion of molecules through a bilayer. Due to the fluid-like |
40 |
|
|
properties of a lipid membrane, not all diffusion across the membrane |
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 |
mmeineke |
1061 |
nanoseconds) within the bilayer. Such long simulation times are |
45 |
mmeineke |
1001 |
difficulty to obtain when integrating the system with atomistic |
46 |
|
|
detail. |
47 |
|
|
|
48 |
|
|
Addressing these issues, several schemes have been proposed. One |
49 |
mmeineke |
1061 |
approach by Goetz and Liposky\cite{goetz98} is to model the entire |
50 |
mmeineke |
1001 |
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 |
mmeineke |
1061 |
model proposed by Marrinck \emph{et. al}.\cite{marrink04}~uses a |
55 |
mmeineke |
1001 |
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 |
58 |
|
|
bilayer. While the solvent spheres are kept charge-neutral and |
59 |
|
|
interact with the bilayer solely through an attractive Lennard-Jones |
60 |
|
|
potential. |
61 |
|
|
|
62 |
|
|
The model used in this investigation adds more information to the |
63 |
mmeineke |
1026 |
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 |
mmeineke |
1061 |
phospholipids, allowing for the ability to scale the parameters to |
67 |
mmeineke |
1026 |
reflect various sized chain configurations while keeping the number of |
68 |
|
|
interactions small. What sets this model apart, however, is the use |
69 |
mmeineke |
1061 |
of dipoles to represent the electrostatic nature of the |
70 |
mmeineke |
1026 |
phospholipids. The dipole electrostatic interaction is shorter range |
71 |
mmeineke |
1061 |
than Coulombic ($\frac{1}{r^3}$ versus $\frac{1}{r}$), eliminating the |
72 |
mmeineke |
1026 |
need for a costly Ewald sum. |
73 |
mmeineke |
1001 |
|
74 |
mmeineke |
1026 |
Another key feature of this model, is the use of a dipolar water model |
75 |
mmeineke |
1061 |
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 |
mmeineke |
1026 |
|
82 |
|
|
The following is an outline of this chapter. |
83 |
mmeineke |
1061 |
Sec.~\ref{lipidSec:Methods} is an introduction to the lipid model |
84 |
mmeineke |
1026 |
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 |
87 |
|
|
Sec.~\ref{lipidSec:ExpSetup}. Sec.~\ref{lipidSec:Results} and |
88 |
|
|
Sec.~\ref{lipidSec:Discussion} give a summary of the results and |
89 |
|
|
interpretation of those results respectively. Finally, the |
90 |
|
|
conclusions of this chapter are presented in |
91 |
|
|
Sec.~\ref{lipidSec:Conclusion}. |
92 |
|
|
|
93 |
mmeineke |
971 |
\section{\label{lipidSec:Methods}Methods} |
94 |
|
|
|
95 |
mmeineke |
1026 |
\subsection{\label{lipidSec:lipidModel}The Lipid Model} |
96 |
|
|
|
97 |
mmeineke |
971 |
\begin{figure} |
98 |
mmeineke |
1061 |
\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 |
mmeineke |
971 |
\label{lipidFig:lipidModel} |
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 |
mmeineke |
1061 |
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 |
mmeineke |
971 |
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 |
mmeineke |
1061 |
scaled by $\sigma_{\text{tail}}$ and $\epsilon_{\text{tail}}$. |
116 |
mmeineke |
971 |
|
117 |
|
|
The long range interactions between lipids are given by the following: |
118 |
|
|
\begin{equation} |
119 |
mmeineke |
1061 |
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 |
mmeineke |
971 |
\label{lipidEq:LJpot} |
125 |
|
|
\end{equation} |
126 |
|
|
and |
127 |
|
|
\begin{equation} |
128 |
mmeineke |
1061 |
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 |
mmeineke |
971 |
\label{lipidEq:dipolePot} |
136 |
|
|
\end{equation} |
137 |
|
|
Where $V_{\text{LJ}}$ is the Lennard-Jones potential and |
138 |
|
|
$V_{\text{dipole}}$ is the dipole-dipole potential. As previously |
139 |
|
|
stated, $\sigma_{ij}$ and $\epsilon_{ij}$ are the Lennard-Jones |
140 |
|
|
parameters which scale the length and depth of the interaction |
141 |
|
|
respectively, and $r_{ij}$ is the distance between beads $i$ and $j$. |
142 |
|
|
In $V_{\text{dipole}}$, $\mathbf{r}_{ij}$ is the vector starting at |
143 |
|
|
bead$i$ and pointing towards bead $j$. Vectors $\mathbf{\Omega}_i$ |
144 |
|
|
and $\mathbf{\Omega}_j$ are the orientational degrees of freedom for |
145 |
|
|
beads $i$ and $j$. $|\mu_i|$ is the magnitude of the dipole moment of |
146 |
|
|
$i$, and $\boldsymbol{\hat{u}}_i$ is the standard unit orientation |
147 |
|
|
vector of $\boldsymbol{\Omega}_i$. |
148 |
|
|
|
149 |
mmeineke |
1061 |
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 |
mmeineke |
971 |
The bends are subject to a harmonic potential: |
153 |
|
|
\begin{equation} |
154 |
mmeineke |
1061 |
V_{\text{bend}}(\theta) = k_{\theta}( \theta - \theta_0 )^2 |
155 |
mmeineke |
971 |
\label{lipidEq:bendPot} |
156 |
|
|
\end{equation} |
157 |
mmeineke |
1061 |
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 |
mmeineke |
971 |
\begin{equation} |
166 |
mmeineke |
1061 |
V_{\text{torsion}}(\phi) = |
167 |
|
|
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0 |
168 |
mmeineke |
971 |
\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 |
mmeineke |
1061 |
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 |
mmeineke |
971 |
directly involved in a bonded pair are calculated. |
177 |
|
|
|
178 |
mmeineke |
1026 |
All simulations presented here use a two chained lipid as pictured in |
179 |
mmeineke |
1061 |
Fig.~\ref{lipidFig:lipidModel}. The chains are both eight beads long, |
180 |
mmeineke |
1026 |
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 |
mmeineke |
1061 |
are summarized in Table~\ref{lipidTable:tcBendParams} and |
184 |
|
|
\ref{lipidTable:tcTorsionParams}. |
185 |
mmeineke |
971 |
|
186 |
mmeineke |
1061 |
\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 |
mmeineke |
1026 |
|
202 |
mmeineke |
1061 |
\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 |
mmeineke |
1026 |
As mentioned previously, the water model used throughout these |
233 |
mmeineke |
1061 |
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 |
mmeineke |
1026 |
|
242 |
|
|
|
243 |
|
|
\subsection{\label{lipidSec:ExpSetup}Experimental Setup} |
244 |
|
|
|
245 |
|
|
Two main starting configuration classes were used in this research: |
246 |
|
|
random and ordered bilayers. The ordered bilayer starting |
247 |
|
|
configurations were all started from an equilibrated bilayer at |
248 |
|
|
300~K. The original configuration for the first 300~K run was |
249 |
|
|
assembled by placing the phospholipids centers of mass on a planar |
250 |
|
|
hexagonal lattice. The lipids were oriented with their long axis |
251 |
|
|
perpendicular to the plane. The second leaf simply mirrored the first |
252 |
|
|
leaf, and the appropriate number of waters were then added above and |
253 |
|
|
below the bilayer. |
254 |
|
|
|
255 |
|
|
The random configurations took more work to generate. To begin, a |
256 |
|
|
test lipid was placed in a simulation box already containing water at |
257 |
|
|
the intended density. The waters were then tested for overlap with |
258 |
|
|
the lipid using a 5.0~$\mbox{\AA}$ buffer distance. This gave an |
259 |
|
|
estimate for the number of waters each lipid would displace in a |
260 |
|
|
simulation box. A target number of waters was then defined which |
261 |
|
|
included the number of waters each lipid would displace, the number of |
262 |
|
|
waters desired to solvate each lipid, and a fudge factor to pad the |
263 |
|
|
initialization. |
264 |
|
|
|
265 |
|
|
Next, a cubic simulation box was created that contained at least the |
266 |
|
|
target number of waters in an FCC lattice (the lattice was for ease of |
267 |
|
|
placement). What followed was a RSA simulation similar to those of |
268 |
|
|
Chapt.~\ref{chapt:RSA}. The lipids were sequentially given a random |
269 |
|
|
position and orientation within the box. If a lipid's position caused |
270 |
|
|
atomic overlap with any previously adsorbed lipid, its position and |
271 |
|
|
orientation were rejected, and a new random adsorption site was |
272 |
|
|
attempted. The RSA simulation proceeded until all phospholipids had |
273 |
|
|
been adsorbed. After adsorption, all water molecules with locations |
274 |
|
|
that overlapped with the atomic coordinates of the lipids were |
275 |
|
|
removed. |
276 |
|
|
|
277 |
|
|
Finally, water molecules were removed one by one at random until the |
278 |
|
|
desired number of waters per lipid was reached. The typical low final |
279 |
|
|
density for these initial configurations was not a problem, as the box |
280 |
|
|
would shrink to an appropriate size within the first 50~ps of a |
281 |
|
|
simulation in the $\text{NPT}_{xyz}$ ensemble. |
282 |
|
|
|
283 |
|
|
\subsection{\label{lipidSec:Configs}The simulation configurations} |
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 |
mmeineke |
1061 |
in an ordered bilayer and observed over a period of 10~ns. Simulation |
288 |
mmeineke |
1026 |
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 |
mmeineke |
1061 |
spontaneous aggregation of an inverted hexagonal lamellar phase. |