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 |
|
|
dipalmitoylphosphatidylcholine (DPPC), \cite{Lindahl:2000} 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 |
15 |
|
|
investigations are limited to 64 to 256 |
16 |
|
|
phospholipids.\cite{Lindal:2000,Sum:2003,Venable:2000,Gomez:2003,Smondyrev:1999,Marrinck:2001a} |
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 |
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. |
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 |
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 |
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 |
|
|
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 |
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 |
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 |
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 |
|
|
phopholipids, 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 |
70 |
|
|
phospholipids. The dipole electrostatic interaction is shorter range |
71 |
|
|
than coulombic ($\frac{1}{r^3}$ versus $\frac{1}{r}$), eliminating the |
72 |
|
|
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 |
|
|
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. |
82 |
|
|
|
83 |
|
|
The following is an outline of this chapter. |
84 |
|
|
Sec.~\ref{lipoidSec:Methods} is an introduction to the lipid model |
85 |
|
|
used in these simulations. As well as clarification about the water |
86 |
|
|
model and integration techniques. The various simulation setups |
87 |
|
|
explored in this research are outlined in |
88 |
|
|
Sec.~\ref{lipidSec:ExpSetup}. Sec.~\ref{lipidSec:Results} and |
89 |
|
|
Sec.~\ref{lipidSec:Discussion} give a summary of the results and |
90 |
|
|
interpretation of those results respectively. Finally, the |
91 |
|
|
conclusions of this chapter are presented in |
92 |
|
|
Sec.~\ref{lipidSec:Conclusion}. |
93 |
|
|
|
94 |
mmeineke |
971 |
\section{\label{lipidSec:Methods}Methods} |
95 |
|
|
|
96 |
|
|
|
97 |
mmeineke |
1026 |
|
98 |
|
|
\subsection{\label{lipidSec:lipidModel}The Lipid Model} |
99 |
|
|
|
100 |
mmeineke |
971 |
\begin{figure} |
101 |
|
|
|
102 |
|
|
\caption{Schematic diagram of the single chain phospholipid model} |
103 |
|
|
|
104 |
|
|
\label{lipidFig:lipidModel} |
105 |
|
|
|
106 |
|
|
\end{figure} |
107 |
|
|
|
108 |
|
|
The phospholipid model used in these simulations is based on the |
109 |
|
|
design illustrated in Fig.~\ref{lipidFig:lipidModel}. The head group |
110 |
|
|
of the phospholipid is replaced by a single Lennard-Jones sphere of |
111 |
|
|
diameter $fix$, with $fix$ scaling the well depth of its van der Walls |
112 |
|
|
interaction. This sphere also contains a single dipole of magnitude |
113 |
|
|
$fix$, where $fix$ can be varied to mimic the charge separation of a |
114 |
|
|
given phospholipid head group. The atoms of the tail region are |
115 |
|
|
modeled by unified atom beads. They are free of partial charges or |
116 |
|
|
dipoles, containing only Lennard-Jones interaction sites at their |
117 |
|
|
centers of mass. As with the head groups, their potentials can be |
118 |
|
|
scaled by $fix$ and $fix$. |
119 |
|
|
|
120 |
|
|
The long range interactions between lipids are given by the following: |
121 |
|
|
\begin{equation} |
122 |
|
|
EQ Here |
123 |
|
|
\label{lipidEq:LJpot} |
124 |
|
|
\end{equation} |
125 |
|
|
and |
126 |
|
|
\begin{equation} |
127 |
|
|
EQ Here |
128 |
|
|
\label{lipidEq:dipolePot} |
129 |
|
|
\end{equation} |
130 |
|
|
Where $V_{\text{LJ}}$ is the Lennard-Jones potential and |
131 |
|
|
$V_{\text{dipole}}$ is the dipole-dipole potential. As previously |
132 |
|
|
stated, $\sigma_{ij}$ and $\epsilon_{ij}$ are the Lennard-Jones |
133 |
|
|
parameters which scale the length and depth of the interaction |
134 |
|
|
respectively, and $r_{ij}$ is the distance between beads $i$ and $j$. |
135 |
|
|
In $V_{\text{dipole}}$, $\mathbf{r}_{ij}$ is the vector starting at |
136 |
|
|
bead$i$ and pointing towards bead $j$. Vectors $\mathbf{\Omega}_i$ |
137 |
|
|
and $\mathbf{\Omega}_j$ are the orientational degrees of freedom for |
138 |
|
|
beads $i$ and $j$. $|\mu_i|$ is the magnitude of the dipole moment of |
139 |
|
|
$i$, and $\boldsymbol{\hat{u}}_i$ is the standard unit orientation |
140 |
|
|
vector of $\boldsymbol{\Omega}_i$. |
141 |
|
|
|
142 |
|
|
The model also allows for the bonded interactions of bonds, bends, and |
143 |
|
|
torsions. The bonds between two beads on a chain are of fixed length, |
144 |
|
|
and are maintained according to the {\sc rattle} algorithm. \cite{fix} |
145 |
|
|
The bends are subject to a harmonic potential: |
146 |
|
|
\begin{equation} |
147 |
|
|
eq here |
148 |
|
|
\label{lipidEq:bendPot} |
149 |
|
|
\end{equation} |
150 |
|
|
where $fix$ scales the strength of the harmonic well, and $fix$ is the |
151 |
|
|
angle between bond vectors $fix$ and $fix$. The torsion potential is |
152 |
|
|
given by: |
153 |
|
|
\begin{equation} |
154 |
|
|
eq here |
155 |
|
|
\label{lipidEq:torsionPot} |
156 |
|
|
\end{equation} |
157 |
|
|
Here, the parameters $k_0$, $k_1$, $k_2$, and $k_3$ fit the cosine |
158 |
|
|
power series to the desired torsion potential surface, and $\phi$ is |
159 |
|
|
the angle between bondvectors $fix$ and $fix$ along the vector $fix$ |
160 |
|
|
(see Fig.:\ref{lipidFig:lipidModel}). Long range interactions such as |
161 |
|
|
the Lennard-Jones potential are excluded for bead pairs involved in |
162 |
|
|
the same bond, bend, or torsion. However, internal interactions not |
163 |
|
|
directly involved in a bonded pair are calculated. |
164 |
|
|
|
165 |
mmeineke |
1026 |
All simulations presented here use a two chained lipid as pictured in |
166 |
|
|
Fig.~\ref{lipidFig:twochain}. The chains are both eight beads long, |
167 |
|
|
and their mass and Lennard Jones parameters are summarized in |
168 |
|
|
Table~\ref{lipidTable:tcLJParams}. The magnitude of the dipole moment |
169 |
|
|
for the head bead is 10.6~Debye, and the bend and torsion parameters |
170 |
|
|
are summarized in Table~\ref{lipidTable:teBTParams}. |
171 |
mmeineke |
971 |
|
172 |
mmeineke |
1026 |
\section{label{lipidSec:furtherMethod}Further Methodology} |
173 |
|
|
|
174 |
|
|
As mentioned previously, the water model used throughout these |
175 |
|
|
simulations was the {\scssd} model of |
176 |
|
|
Ichiye.\cite{liu:1996a,Liu:1996b,Chandra:1999} A discussion of the |
177 |
|
|
model can be found in Sec.~\ref{oopseSec:SSD}. As for the integration |
178 |
|
|
of the equations of motion, all simulations were performed in an |
179 |
|
|
orthorhombic periodic box with a thermostat on velocities, and an |
180 |
|
|
independent barostat on each cartesian axis $x$, $y$, and $z$. This |
181 |
|
|
is the $\text{NPT}_{xyz}$. ensemble described in Sec.~\ref{oopseSec:Ensembles}. |
182 |
|
|
|
183 |
|
|
|
184 |
|
|
\subsection{\label{lipidSec:ExpSetup}Experimental Setup} |
185 |
|
|
|
186 |
|
|
Two main starting configuration classes were used in this research: |
187 |
|
|
random and ordered bilayers. The ordered bilayer starting |
188 |
|
|
configurations were all started from an equilibrated bilayer at |
189 |
|
|
300~K. The original configuration for the first 300~K run was |
190 |
|
|
assembled by placing the phospholipids centers of mass on a planar |
191 |
|
|
hexagonal lattice. The lipids were oriented with their long axis |
192 |
|
|
perpendicular to the plane. The second leaf simply mirrored the first |
193 |
|
|
leaf, and the appropriate number of waters were then added above and |
194 |
|
|
below the bilayer. |
195 |
|
|
|
196 |
|
|
The random configurations took more work to generate. To begin, a |
197 |
|
|
test lipid was placed in a simulation box already containing water at |
198 |
|
|
the intended density. The waters were then tested for overlap with |
199 |
|
|
the lipid using a 5.0~$\mbox{\AA}$ buffer distance. This gave an |
200 |
|
|
estimate for the number of waters each lipid would displace in a |
201 |
|
|
simulation box. A target number of waters was then defined which |
202 |
|
|
included the number of waters each lipid would displace, the number of |
203 |
|
|
waters desired to solvate each lipid, and a fudge factor to pad the |
204 |
|
|
initialization. |
205 |
|
|
|
206 |
|
|
Next, a cubic simulation box was created that contained at least the |
207 |
|
|
target number of waters in an FCC lattice (the lattice was for ease of |
208 |
|
|
placement). What followed was a RSA simulation similar to those of |
209 |
|
|
Chapt.~\ref{chapt:RSA}. The lipids were sequentially given a random |
210 |
|
|
position and orientation within the box. If a lipid's position caused |
211 |
|
|
atomic overlap with any previously adsorbed lipid, its position and |
212 |
|
|
orientation were rejected, and a new random adsorption site was |
213 |
|
|
attempted. The RSA simulation proceeded until all phospholipids had |
214 |
|
|
been adsorbed. After adsorption, all water molecules with locations |
215 |
|
|
that overlapped with the atomic coordinates of the lipids were |
216 |
|
|
removed. |
217 |
|
|
|
218 |
|
|
Finally, water molecules were removed one by one at random until the |
219 |
|
|
desired number of waters per lipid was reached. The typical low final |
220 |
|
|
density for these initial configurations was not a problem, as the box |
221 |
|
|
would shrink to an appropriate size within the first 50~ps of a |
222 |
|
|
simulation in the $\text{NPT}_{xyz}$ ensemble. |
223 |
|
|
|
224 |
|
|
\subsection{\label{lipidSec:Configs}The simulation configurations} |
225 |
|
|
|
226 |
|
|
Table ~\ref{lipidTable:simNames} summarizes the names and important |
227 |
|
|
details of the simulations. The B set of simulations were all started |
228 |
|
|
in an ordered bilayer and observed over a period of 10~ns. Simulution |
229 |
|
|
RL was integrated for approximately 20~ns starting from a random |
230 |
|
|
configuration as an example of spontaneous bilayer aggregation. |
231 |
|
|
Lastly, simulation RH was also started from a random configuration, |
232 |
|
|
but with a lesser water content and higher temperature to show the |
233 |
|
|
spontaneous aggregation of an inverted hexagonal lamellar phase. |