4 |
|
|
5 |
|
\section{\label{lipidSec:Intro}Introduction} |
6 |
|
|
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 |
7 |
> |
In the past 10 years, increasing computer speeds have allowed for the |
8 |
> |
atomistic simulation of phospholipid bilayers for increasingly |
9 |
> |
relevant lengths of time. These simulations have ranged from |
10 |
> |
simulation of the gel phase ($L_{\beta}$) of |
11 |
|
dipalmitoylphosphatidylcholine (DPPC),\cite{lindahl00} to the |
12 |
|
spontaneous aggregation of DPPC molecules into fluid phase |
13 |
|
($L_{\alpha}$) bilayers.\cite{marrink01} With the exception of a few |
14 |
|
ambitious |
15 |
|
simulations,\cite{marrink01:undulation,marrink:2002,lindahl00} most |
16 |
< |
investigations are limited to 64 to 256 |
16 |
> |
investigations are limited to a range of 64 to 256 |
17 |
|
phospholipids.\cite{lindahl00,sum:2003,venable00,gomez:2003,smondyrev:1999,marrink01} |
18 |
< |
This is due to the expense of the computer calculations involved when |
19 |
< |
performing these simulations. To properly hydrate a bilayer, one |
18 |
> |
The expense of the force calculations involved when performing these |
19 |
> |
simulations limits the system size. To properly hydrate a bilayer, one |
20 |
|
typically needs 25 water molecules for every lipid, bringing the total |
21 |
|
number of atoms simulated to roughly 8,000 for a system of 64 DPPC |
22 |
|
molecules. Added to the difficulty is the electrostatic nature of the |
23 |
< |
phospholipid head groups and water, requiring the computationally |
24 |
< |
expensive Ewald sum or its slightly faster derivative particle mesh |
25 |
< |
Ewald sum.\cite{nina:2002,norberg:2000,patra:2003} These factors all |
26 |
< |
limit the potential size and time lengths of bilayer simulations. |
23 |
> |
phospholipid head groups and water, requiring either the |
24 |
> |
computationally expensive Ewald sum or the faster, particle mesh Ewald |
25 |
> |
sum.\cite{nina:2002,norberg:2000,patra:2003} These factors all limit |
26 |
> |
the system size and time scales of bilayer simulations. |
27 |
|
|
28 |
|
Unfortunately, much of biological interest happens on time and length |
29 |
< |
scales infeasible with current simulation. One such example is the |
30 |
< |
observance of a ripple phase ($P_{\beta^{\prime}}$) between the |
31 |
< |
$L_{\beta}$ and $L_{\alpha}$ phases of certain phospholipid |
32 |
< |
bilayers.\cite{katsaras00,sengupta00} These ripples are shown to |
33 |
< |
have periodicity on the order of 100-200~$\mbox{\AA}$. A simulation on |
34 |
< |
this length scale would have approximately 1,300 lipid molecules with |
35 |
< |
an additional 25 water molecules per lipid to fully solvate the |
36 |
< |
bilayer. A simulation of this size is impractical with current |
37 |
< |
atomistic models. |
29 |
> |
scales well beyond the range of current simulation technology. One |
30 |
> |
such example is the observance of a ripple phase |
31 |
> |
($P_{\beta^{\prime}}$) between the $L_{\beta}$ and $L_{\alpha}$ phases |
32 |
> |
of certain phospholipid bilayers.\cite{katsaras00,sengupta00} These |
33 |
> |
ripples are shown to have periodicity on the order of |
34 |
> |
100-200~$\mbox{\AA}$. A simulation on this length scale would have |
35 |
> |
approximately 1,300 lipid molecules with an additional 25 water |
36 |
> |
molecules per lipid to fully solvate the bilayer. A simulation of this |
37 |
> |
size is impractical with current atomistic models. |
38 |
|
|
39 |
< |
Another class of simulations to consider, are those dealing with the |
40 |
< |
diffusion of molecules through a bilayer. Due to the fluid-like |
41 |
< |
properties of a lipid membrane, not all diffusion across the membrane |
42 |
< |
happens at pores. Some molecules of interest may incorporate |
43 |
< |
themselves directly into the membrane. Once here, they may possess an |
44 |
< |
appreciable waiting time (on the order of 10's to 100's of |
45 |
< |
nanoseconds) within the bilayer. Such long simulation times are |
46 |
< |
difficulty to obtain when integrating the system with atomistic |
46 |
< |
detail. |
39 |
> |
The time and length scale limitations are most striking in transport |
40 |
> |
phenomena. Due to the fluid-like properties of a lipid membrane, not |
41 |
> |
all diffusion across the membrane happens at pores. Some molecules of |
42 |
> |
interest may incorporate themselves directly into the membrane. Once |
43 |
> |
here, they may possess an appreciable waiting time (on the order of |
44 |
> |
10's to 100's of nanoseconds) within the bilayer. Such long simulation |
45 |
> |
times are nearly impossible to obtain when integrating the system with |
46 |
> |
atomistic detail. |
47 |
|
|
48 |
< |
Addressing these issues, several schemes have been proposed. One |
48 |
> |
To address these issues, several schemes have been proposed. One |
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 |
68 |
|
interactions small. What sets this model apart, however, is the use |
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 |
72 |
< |
need for a costly Ewald sum. |
71 |
> |
than Coulombic ($\frac{1}{r^3}$ versus $\frac{1}{r}$), and eliminates |
72 |
> |
the 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 ({\sc ssd}) water |
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. |
80 |
> |
deformation around a non-polar entity, \emph{i.e.}~the phospholipid |
81 |
> |
molecules. |
82 |
|
|
83 |
|
The following is an outline of this chapter. |
84 |
< |
Sec.~\ref{lipidSec: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 |
90 |
< |
conclusions of this chapter are presented in |
91 |
< |
Sec.~\ref{lipidSec:Conclusion}. |
84 |
> |
Sec.~\ref{lipidSec:Methods} is an introduction to the lipid model used |
85 |
> |
in these simulations, as well as clarification about the water model |
86 |
> |
and integration techniques. The various simulations explored in this |
87 |
> |
research are outlined in |
88 |
> |
Sec.~\ref{lipidSec:ExpSetup}. Sec.~\ref{lipidSec:resultsDis} gives a |
89 |
> |
summary and interpretation of the results. Finally, the conclusions |
90 |
> |
of this chapter are presented in Sec.~\ref{lipidSec:Conclusion}. |
91 |
|
|
92 |
|
\section{\label{lipidSec:Methods}Methods} |
93 |
|
|
108 |
|
contains a single dipole of magnitude $|\boldsymbol{\mu}|$, where |
109 |
|
$|\boldsymbol{\mu}|$ can be varied to mimic the charge separation of a |
110 |
|
given phospholipid head group. The atoms of the tail region are |
111 |
< |
modeled by unified atom beads. They are free of partial charges or |
112 |
< |
dipoles, containing only Lennard-Jones interaction sites at their |
113 |
< |
centers of mass. As with the head groups, their potentials can be |
114 |
< |
scaled by $\sigma_{\text{tail}}$ and $\epsilon_{\text{tail}}$. |
111 |
> |
modeled by beads representing multiple methyl groups. They are free |
112 |
> |
of partial charges or dipoles, and contain only Lennard-Jones |
113 |
> |
interaction sites at their centers of mass. As with the head groups, |
114 |
> |
their potentials can be scaled by $\sigma_{\text{tail}}$ and |
115 |
> |
$\epsilon_{\text{tail}}$. |
116 |
|
|
117 |
|
The long range interactions between lipids are given by the following: |
118 |
|
\begin{equation} |
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] |
132 |
> |
3(\boldsymbol{\hat{u}}_i \cdot \mathbf{\hat{r}}_{ij}) % |
133 |
> |
(\boldsymbol{\hat{u}}_j \cdot \mathbf{\hat{r}}_{ij} \biggr] |
134 |
|
\label{lipidEq:dipolePot} |
135 |
|
\end{equation} |
136 |
|
Where $V_{\text{LJ}}$ is the Lennard-Jones potential and |
139 |
|
parameters which scale the length and depth of the interaction |
140 |
|
respectively, and $r_{ij}$ is the distance between beads $i$ and $j$. |
141 |
|
In $V_{\text{dipole}}$, $\mathbf{r}_{ij}$ is the vector starting at |
142 |
< |
bead$i$ and pointing towards bead $j$. Vectors $\mathbf{\Omega}_i$ |
142 |
> |
bead $i$ and pointing towards bead $j$. Vectors $\mathbf{\Omega}_i$ |
143 |
|
and $\mathbf{\Omega}_j$ are the orientational degrees of freedom for |
144 |
|
beads $i$ and $j$. $|\mu_i|$ is the magnitude of the dipole moment of |
145 |
|
$i$, and $\boldsymbol{\hat{u}}_i$ is the standard unit orientation |
146 |
< |
vector of $\boldsymbol{\Omega}_i$. |
146 |
> |
vector rotated with Euler angles: $\boldsymbol{\Omega}_i$. |
147 |
|
|
148 |
|
The model also allows for the bonded interactions bends, and torsions. |
149 |
|
The bond between two beads on a chain is of fixed length, and is |
158 |
|
(Fig.~\ref{lipidFig:lipidModel}). In addition, we have placed a |
159 |
|
``ghost'' bend on the phospholipid head. The ghost bend adds a |
160 |
|
potential to keep the dipole pointed along the bilayer surface, where |
161 |
< |
$theta$ is now the angle the dipole makes with respect to the {\sc |
162 |
< |
head}-$\text{{\sc ch}}_2$ bond vector. The torsion potential is given |
163 |
< |
by: |
161 |
> |
$\theta$ is now the angle the dipole makes with respect to the {\sc |
162 |
> |
head}-$\text{{\sc ch}}_2$ bond vector |
163 |
> |
(Fig.~\ref{lipidFig:ghostBend}). The torsion potential is given by: |
164 |
|
\begin{equation} |
165 |
|
V_{\text{torsion}}(\phi) = |
166 |
|
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0 |
174 |
|
same bond, bend, or torsion. However, internal interactions not |
175 |
|
directly involved in a bonded pair are calculated. |
176 |
|
|
177 |
+ |
\begin{figure} |
178 |
+ |
\centering |
179 |
+ |
\includegraphics[width=\linewidth]{ghostBendFig.eps} |
180 |
+ |
\caption[Depiction of the ``ghost'' bend]{The ``ghost'' bend is a bend potential added to constrain the motion of the dipole on the {\sc head} group. The potential follows Eq.~\ref{lipidEq:bendPot} where $\theta$ is now the angle that the dipole makes with the {\sc head}-$\text{{\sc ch}}_2$ bond vector.} |
181 |
+ |
\label{lipidFig:ghostBend} |
182 |
+ |
\end{figure} |
183 |
+ |
|
184 |
|
All simulations presented here use a two chained lipid as pictured in |
185 |
|
Fig.~\ref{lipidFig:lipidModel}. The chains are both eight beads long, |
186 |
|
and their mass and Lennard Jones parameters are summarized in |
243 |
|
performed in an orthorhombic periodic box with a thermostat on |
244 |
|
velocities, and an independent barostat on each Cartesian axis $x$, |
245 |
|
$y$, and $z$. This is the $\text{NPT}_{xyz}$. ensemble described in |
246 |
< |
Sec.~\ref{oopseSec:Ensembles}. |
246 |
> |
Sec.~\ref{oopseSec:integrate}. |
247 |
|
|
248 |
|
|
249 |
|
\subsection{\label{lipidSec:ExpSetup}Experimental Setup} |
265 |
|
estimate for the number of waters each lipid would displace in a |
266 |
|
simulation box. A target number of waters was then defined which |
267 |
|
included the number of waters each lipid would displace, the number of |
268 |
< |
waters desired to solvate each lipid, and a fudge factor to pad the |
269 |
< |
initialization. |
268 |
> |
waters desired to solvate each lipid, and a factor to pad the |
269 |
> |
initial box with a few extra water molecules. |
270 |
|
|
271 |
|
Next, a cubic simulation box was created that contained at least the |
272 |
|
target number of waters in an FCC lattice (the lattice was for ease of |
273 |
|
placement). What followed was a RSA simulation similar to those of |
274 |
|
Chapt.~\ref{chapt:RSA}. The lipids were sequentially given a random |
275 |
|
position and orientation within the box. If a lipid's position caused |
276 |
< |
atomic overlap with any previously adsorbed lipid, its position and |
277 |
< |
orientation were rejected, and a new random adsorption site was |
276 |
> |
atomic overlap with any previously placed lipid, its position and |
277 |
> |
orientation were rejected, and a new random placement site was |
278 |
|
attempted. The RSA simulation proceeded until all phospholipids had |
279 |
< |
been adsorbed. After adsorption, all water molecules with locations |
280 |
< |
that overlapped with the atomic coordinates of the lipids were |
281 |
< |
removed. |
279 |
> |
been adsorbed. After placement of all lipid molecules, water |
280 |
> |
molecules with locations that overlapped with the atomic coordinates |
281 |
> |
of the lipids were removed. |
282 |
|
|
283 |
< |
Finally, water molecules were removed one by one at random until the |
284 |
< |
desired number of waters per lipid was reached. The typical low final |
285 |
< |
density for these initial configurations was not a problem, as the box |
286 |
< |
would shrink to an appropriate size within the first 50~ps of a |
287 |
< |
simulation in the $\text{NPT}_{xyz}$ ensemble. |
283 |
> |
Finally, water molecules were removed at random until the desired |
284 |
> |
number of waters per lipid was reached. The typical low final density |
285 |
> |
for these initial configurations was not a problem, as the box shrinks |
286 |
> |
to an appropriate size within the first 50~ps of a simulation in the |
287 |
> |
$\text{NPT}_{xyz}$ ensemble. |
288 |
|
|
289 |
< |
\subsection{\label{lipidSec:Configs}The simulation configurations} |
289 |
> |
\subsection{\label{lipidSec:Configs}Configurations} |
290 |
|
|
291 |
< |
Table ~\ref{lipidTable:simNames} summarizes the names and important |
292 |
< |
details of the simulations. The B set of simulations were all started |
293 |
< |
in an ordered bilayer and observed over a period of 10~ns. Simulation |
294 |
< |
RL was integrated for approximately 20~ns starting from a random |
295 |
< |
configuration as an example of spontaneous bilayer aggregation. |
296 |
< |
Lastly, simulation RH was also started from a random configuration, |
297 |
< |
but with a lesser water content and higher temperature to show the |
298 |
< |
spontaneous aggregation of an inverted hexagonal lamellar phase. |
291 |
> |
The first class of simulations were started from ordered |
292 |
> |
bilayers. They were all configurations consisting of 60 lipid |
293 |
> |
molecules with 30 lipids on each leaf, and were hydrated with 1620 |
294 |
> |
{\sc ssd} molecules. The original configuration was assembled |
295 |
> |
according to Sec.~\ref{lipidSec:ExpSetup} and simulated for a length |
296 |
> |
of 10~ns at 300~K. The other temperature runs were started from a |
297 |
> |
frame 7~ns into the 300~K simulation. Their temperatures were reset |
298 |
> |
with the thermostating algorithm in the $\text{NPT}_{xyz}$ |
299 |
> |
integrator. All of the temperature variants were also run for 10~ns, |
300 |
> |
with only the last 5~ns being used for accumulation of statistics. |
301 |
> |
|
302 |
> |
The second class of simulations were two configurations started from |
303 |
> |
randomly dispersed lipids in a ``gas'' of water. The first |
304 |
> |
($\text{R}_{\text{I}}$) was a simulation containing 72 lipids with |
305 |
> |
1800 {\sc ssd} molecules simulated at 300~K. The second |
306 |
> |
($\text{R}_{\text{II}}$) was 90 lipids with 1350 {\sc ssd} molecules |
307 |
> |
simulated at 350~K. Both simulations were integrated for more than |
308 |
> |
20~ns, and illustrate the spontaneous aggregation of the lipid model |
309 |
> |
into phospholipid macro-structures: $\text{R}_{\text{I}}$ into a |
310 |
> |
bilayer, and $\text{R}_{\text{II}}$ into a inverted rod. |
311 |
> |
|
312 |
> |
\section{\label{lipidSec:resultsDis}Results and Discussion} |
313 |
> |
|
314 |
> |
\subsection{\label{lipidSec:diffusion}Lateral Diffusion Constants} |
315 |
> |
|
316 |
> |
The lateral diffusion constant, $D_L$, is the constant characterizing |
317 |
> |
the diffusive motion of the lipid within the plane of the bilayer. It |
318 |
> |
is given by the following Einstein relation valid at long |
319 |
> |
times:\cite{allen87:csl} |
320 |
> |
\begin{equation} |
321 |
> |
2tD_L = \frac{1}{2}\langle |\mathbf{r}(t) - \mathbf{r}(0)|^2\rangle |
322 |
> |
\end{equation} |
323 |
> |
Where $\mathbf{r}(t)$ is the position of the lipid at time $t$, and is |
324 |
> |
constrained to lie within a plane. For the bilayer simulations the |
325 |
> |
plane of constrained motion was that perpendicular to the bilayer |
326 |
> |
normal, namely the $xy$-plane. |
327 |
> |
|
328 |
> |
Fig.~\ref{lipidFig:diffusionFig} shows the lateral diffusion constants |
329 |
> |
as a function of temperature. There is a definite increase in the |
330 |
> |
lateral diffusion with higher temperatures, which is exactly what one |
331 |
> |
would expect with greater fluidity of the chains. However, the |
332 |
> |
diffusion constants are all two orders of magnitude smaller than those |
333 |
> |
typical of DPPC.\cite{Cevc87} This is counter-intuitive as the DPPC |
334 |
> |
molecule is sterically larger and heavier than our model. This could |
335 |
> |
be an indication that our model's chains are too interwoven and hinder |
336 |
> |
the motion of the lipid, or that a simplification in the model's |
337 |
> |
forces has led to a slowing of diffusive behavior within the |
338 |
> |
bilayer. In contrast, the diffusion constant of the {\sc ssd} water, |
339 |
> |
$9.84\times 10^{-6}\,\text{cm}^2/\text{s}$, compares favorably with |
340 |
> |
that of bulk water ($2.2999\times |
341 |
> |
10^{-5}\,\text{cm}^2/\text{s}$\cite{Holz00}). |
342 |
> |
|
343 |
> |
\begin{figure} |
344 |
> |
\centering |
345 |
> |
\includegraphics[width=\linewidth]{diffusionFig.eps} |
346 |
> |
\caption[The lateral diffusion constants versus temperature]{The lateral diffusion constants for the bilayers as a function of temperature.} |
347 |
> |
\label{lipidFig:diffusionFig} |
348 |
> |
\end{figure} |
349 |
> |
|
350 |
> |
\subsection{\label{lipidSec:densProf}Density Profile} |
351 |
> |
|
352 |
> |
Fig.~\ref{lipidFig:densityProfile} illustrates the densities of the |
353 |
> |
atoms in the bilayer systems normalized by the bulk density as a |
354 |
> |
function of distance from the center of the box. The profile is taken |
355 |
> |
along the bilayer normal, in this case the $z$ axis. The profile at |
356 |
> |
270~K shows several structural features that are largely smoothed out |
357 |
> |
by 300~K. The left peak for the {\sc head} atoms is split at 270~K, |
358 |
> |
implying that some freezing of the structure might already be occurring |
359 |
> |
at this temperature. From the dynamics, the tails at this temperature |
360 |
> |
are very much fluid, but the profile could indicate that a phase |
361 |
> |
transition may simply be beyond the length scale of the current |
362 |
> |
simulation. In all profiles, the water penetrates almost |
363 |
> |
5~$\mbox{\AA}$ into the bilayer, completely solvating the {\sc head} |
364 |
> |
atoms. The $\text{{\sc ch}}_3$ atoms although mainly centered at the |
365 |
> |
middle of the bilayer, show appreciable penetration into the head |
366 |
> |
group region. This indicates that the chains have enough mobility to |
367 |
> |
bend back upward to allow the ends to explore areas around the {\sc |
368 |
> |
head} atoms. It is unlikely that this is penetration from a lipid of |
369 |
> |
the opposite face, as the lipids are only 12~$\mbox{\AA}$ in length, |
370 |
> |
and the typical leaf spacing as measured from the {\sc head-head} |
371 |
> |
spacing in the profile is 17.5~$\mbox{\AA}$. |
372 |
> |
|
373 |
> |
\begin{figure} |
374 |
> |
\centering |
375 |
> |
\includegraphics[width=\linewidth]{densityProfile.eps} |
376 |
> |
\caption[The density profile of the lipid bilayers]{The density profile of the lipid bilayers along the bilayer normal. The black lines are the {\sc head} atoms, red lines are the {\sc ch} atoms, green lines are the $\text{{\sc ch}}_2$ atoms, blue lines are the $\text{{\sc ch}}_3$ atoms, and the magenta lines are the {\sc ssd} atoms.} |
377 |
> |
\label{lipidFig:densityProfile} |
378 |
> |
\end{figure} |
379 |
> |
|
380 |
> |
|
381 |
> |
\subsection{\label{lipidSec:scd}$\text{S}_{\text{{\sc cd}}}$ Order Parameters} |
382 |
> |
|
383 |
> |
The $\text{S}_{\text{{\sc cd}}}$ order parameter is often reported in |
384 |
> |
the experimental characterizations of phospholipids. It is obtained |
385 |
> |
through deuterium NMR, and measures the ordering of the carbon |
386 |
> |
deuterium bond in relation to the bilayer normal at various points |
387 |
> |
along the chains. In our model, there are no explicit hydrogens, but |
388 |
> |
the order parameter can be written in terms of the carbon ordering at |
389 |
> |
each point in the chain:\cite{egberts88} |
390 |
> |
\begin{equation} |
391 |
> |
S_{\text{{\sc cd}}} = \frac{2}{3}S_{xx} + \frac{1}{3}S_{yy} |
392 |
> |
\label{lipidEq:scd1} |
393 |
> |
\end{equation} |
394 |
> |
Where $S_{ij}$ is given by: |
395 |
> |
\begin{equation} |
396 |
> |
S_{ij} = \frac{1}{2}\Bigl\langle(3\cos\Theta_i\cos\Theta_j |
397 |
> |
- \delta_{ij})\Bigr\rangle |
398 |
> |
\label{lipidEq:scd2} |
399 |
> |
\end{equation} |
400 |
> |
Here, $\Theta_i$ is the angle the $i$th axis in the reference frame of |
401 |
> |
the carbon atom makes with the bilayer normal. The brackets denote an |
402 |
> |
average over time and molecules. The carbon atom axes are defined: |
403 |
> |
$\mathbf{\hat{z}}\rightarrow$ vector from $C_{n-1}$ to $C_{n+1}$; |
404 |
> |
$\mathbf{\hat{y}}\rightarrow$ vector that is perpendicular to $z$ and |
405 |
> |
in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$; |
406 |
> |
$\mathbf{\hat{x}}\rightarrow$ vector perpendicular to |
407 |
> |
$\mathbf{\hat{y}}$ and $\mathbf{\hat{z}}$. |
408 |
> |
|
409 |
> |
The order parameter has a range of $[1,-\frac{1}{2}]$. A value of 1 |
410 |
> |
implies full order aligned to the bilayer axis, 0 implies full |
411 |
> |
disorder, and $-\frac{1}{2}$ implies full order perpendicular to the |
412 |
> |
bilayer axis. The {\sc cd} bond vector for carbons near the head group |
413 |
> |
are usually ordered perpendicular to the bilayer normal, with tails |
414 |
> |
farther away tending toward disorder. This makes the order parameter |
415 |
> |
negative for most carbons, and as such $|S_{\text{{\sc cd}}}|$ is more |
416 |
> |
commonly reported than $S_{\text{{\sc cd}}}$. |
417 |
> |
|
418 |
> |
Fig.~\ref{lipidFig:scdFig} shows the $S_{\text{{\sc cd}}}$ order |
419 |
> |
parameters for the bilayer system at 300~K. There is no appreciable |
420 |
> |
difference in the plots for the various temperatures, however, there |
421 |
> |
is a larger difference between our models ordering, and that of |
422 |
> |
DMPC. As our values are closer to $-\frac{1}{2}$, this implies more |
423 |
> |
ordering perpendicular to the normal than in a real system. This is |
424 |
> |
due to the model having only one carbon group separating the chains |
425 |
> |
from the top of the lipid. In DMPC, with the flexibility inherent in a |
426 |
> |
multiple atom head group, as well as a glycerol linkage between the |
427 |
> |
head group and the acyl chains, there is more loss of ordering by the |
428 |
> |
point when the chains start. |
429 |
> |
|
430 |
> |
\begin{figure} |
431 |
> |
\centering |
432 |
> |
\includegraphics[width=\linewidth]{scdFig.eps} |
433 |
> |
\caption[$\text{S}_{\text{{\sc cd}}}$ order parameter for our model]{A comparison of $|\text{S}_{\text{{\sc cd}}}|$ between our model (blue) and DMPC\cite{petrache00} (black) near 300~K.} |
434 |
> |
\label{lipidFig:scdFig} |
435 |
> |
\end{figure} |
436 |
> |
|
437 |
> |
\subsection{\label{lipidSec:p2Order}$P_2$ Order Parameter} |
438 |
> |
|
439 |
> |
The $P_2$ order parameter allows us to measure the amount of |
440 |
> |
directional ordering that exists in the bilayer. Each lipid molecule |
441 |
> |
can be thought of as a cylindrical tube with the head group at the |
442 |
> |
top. If all of the cylinders are perfectly aligned, the $P_2$ order |
443 |
> |
parameter will be $1.0$. If the cylinders are completely dispersed, |
444 |
> |
the $P_2$ order parameter will be 0. For a collection of unit vectors, |
445 |
> |
the $P_2$ order parameter can be solved via the following |
446 |
> |
method.\cite{zannoni94} |
447 |
> |
|
448 |
> |
Define an ordering matrix $\mathbf{Q}$, such that, |
449 |
> |
\begin{equation} |
450 |
> |
\mathbf{Q} = \frac{1}{N}\sum_i^N % |
451 |
> |
\begin{pmatrix} % |
452 |
> |
u_{ix}u_{ix}-\frac{1}{3} & u_{ix}u_{iy} & u_{ix}u_{iz} \\ |
453 |
> |
u_{iy}u_{ix} & u_{iy}u_{iy}-\frac{1}{3} & u_{iy}u_{iz} \\ |
454 |
> |
u_{iz}u_{ix} & u_{iz}u_{iy} & u_{iz}u_{iz}-\frac{1}{3} % |
455 |
> |
\end{pmatrix} |
456 |
> |
\label{lipidEq:po1} |
457 |
> |
\end{equation} |
458 |
> |
Where the $u_{i\alpha}$ is the $\alpha$ element of the unit vector |
459 |
> |
$\mathbf{\hat{u}}_i$, and the sum over $i$ averages over the whole |
460 |
> |
collection of unit vectors. This allows the matrix element |
461 |
> |
$Q_{\alpha\beta}$ to be written: |
462 |
> |
\begin{equation} |
463 |
> |
Q_{\alpha\beta} = \langle u_{\alpha}u_{\beta} - |
464 |
> |
\frac{1}{3}\delta_{\alpha\beta} \rangle |
465 |
> |
\label{lipidEq:po2} |
466 |
> |
\end{equation} |
467 |
> |
|
468 |
> |
Having constructed the matrix, diagonalizing $\mathbf{Q}$ yields three |
469 |
> |
eigenvalues and eigenvectors. The eigenvector associated with the |
470 |
> |
largest eigenvalue, $\lambda_{\text{max}}$, is the director for the |
471 |
> |
system of unit vectors. The director is the average direction all of |
472 |
> |
the unit vectors are pointing. The $P_2$ order parameter is then |
473 |
> |
simply |
474 |
> |
\begin{equation} |
475 |
> |
\langle P_2 \rangle = \frac{3}{2}\lambda_{\text{max}} |
476 |
> |
\label{lipidEq:po3} |
477 |
> |
\end{equation} |
478 |
> |
|
479 |
> |
Table~\ref{lipidTab:blSummary} summarizes the $P_2$ values for the |
480 |
> |
bilayers, as well as the dipole orientations. The unit vector for the |
481 |
> |
lipid molecules was defined by finding the moment of inertia for each |
482 |
> |
lipid, then setting $\mathbf{\hat{u}}$ to point along the axis of |
483 |
> |
minimum inertia. For the {\sc head} atoms, the unit vector simply |
484 |
> |
pointed in the same direction as the dipole moment. For the lipid |
485 |
> |
molecules, the ordering was consistent across all temperatures, with |
486 |
> |
the director pointed along the $z$ axis of the box. More |
487 |
> |
interestingly, is the high degree of ordering the dipoles impose on |
488 |
> |
the {\sc head} atoms. The directors for the dipoles consistently |
489 |
> |
pointed along the plane of the bilayer, with the directors |
490 |
> |
anti-aligned on the top and bottom leaf. |
491 |
> |
|
492 |
> |
\begin{table} |
493 |
> |
\caption[Structural properties of the bilayers]{Bilayer Structural properties as a function of temperature.} |
494 |
> |
\label{lipidTab:blSummary} |
495 |
> |
\begin{center} |
496 |
> |
\begin{tabular}{|c|c|c|c|c|} |
497 |
> |
\hline |
498 |
> |
Temperature (K) & $\langle L_{\perp}\rangle$ ($\mbox{\AA}$) & % |
499 |
> |
$\langle A_{\parallel}\rangle$ ($\mbox{\AA}^2$) & % |
500 |
> |
$\langle P_2\rangle_{\text{Lipid}}$ & % |
501 |
> |
$\langle P_2\rangle_{\text{{\sc head}}}$ \\ \hline |
502 |
> |
270 & 18.1 & 58.1 & 0.253 & 0.494 \\ \hline |
503 |
> |
275 & 17.2 & 56.7 & 0.295 & 0.514 \\ \hline |
504 |
> |
277 & 16.9 & 58.0 & 0.301 & 0.541 \\ \hline |
505 |
> |
280 & 17.4 & 58.0 & 0.274 & 0.488 \\ \hline |
506 |
> |
285 & 16.9 & 57.6 & 0.270 & 0.616 \\ \hline |
507 |
> |
290 & 17.0 & 57.6 & 0.263 & 0.534 \\ \hline |
508 |
> |
293 & 17.5 & 58.0 & 0.227 & 0.643 \\ \hline |
509 |
> |
300 & 16.9 & 57.6 & 0.315 & 0.536 \\ \hline |
510 |
> |
\end{tabular} |
511 |
> |
\end{center} |
512 |
> |
\end{table} |
513 |
> |
|
514 |
> |
\subsection{\label{lipidSec:miscData}Further Bilayer Data} |
515 |
> |
|
516 |
> |
Also summarized in Table~\ref{lipidTab:blSummary}, are the bilayer |
517 |
> |
thickness and area per lipid. The bilayer thickness was measured from |
518 |
> |
the peak to peak {\sc head} atom distance in the density profiles. The |
519 |
> |
area per lipid data compares favorably with values typically seen for |
520 |
> |
DMPC (60.0~$\mbox{\AA}^2$ at 303~K)\cite{petrache00}. Although are |
521 |
> |
values are lower this is most likely due to the shorter chain length |
522 |
> |
of our model (8 versus 14 for DMPC). |
523 |
> |
|
524 |
> |
\subsection{\label{lipidSec:randBilayer}Bilayer Aggregation} |
525 |
> |
|
526 |
> |
A very important accomplishment for our model is its ability to |
527 |
> |
spontaneously form bilayers from a randomly dispersed starting |
528 |
> |
configuration. Fig.~\ref{lipidFig:blImage} shows an image sequence for |
529 |
> |
the bilayer aggregation. After 3.0~ns, the basic form of the bilayer |
530 |
> |
can already be seen. By 7.0~ns, the bilayer has a lipid bridge |
531 |
> |
stretched across the simulation box to itself that will turn out to be |
532 |
> |
very long lived ($\sim$20~ns), as well as a water pore, that will |
533 |
> |
persist for the length of the current simulation. At 24~ns, the lipid |
534 |
> |
bridge is dispersed, and the bilayer is still integrating the lipid |
535 |
> |
molecules from the bridge into itself, and has still been unable to |
536 |
> |
disperse the water pore. |
537 |
> |
|
538 |
> |
\begin{figure} |
539 |
> |
\centering |
540 |
> |
\includegraphics[width=\linewidth]{bLayerImage.eps} |
541 |
> |
\caption[Image sequence of the bilayer aggregation]{Image sequence of the bilayer aggregation. The blue beads are the {\sc head} atoms the grey beads are the chains, and the red and white bead are the water molecules. A box has been drawn around the periodic image.} |
542 |
> |
\label{lipidFig:blImage} |
543 |
> |
\end{figure} |
544 |
> |
|
545 |
> |
\subsection{\label{lipidSec:randIrod}Inverted Rod Aggregation} |
546 |
> |
|
547 |
> |
Fig.~\ref{lipidFig:iRimage} shows a second aggregation sequence |
548 |
> |
simulated in this research. Here the fraction of water had been |
549 |
> |
significantly decreased to observe how the model would respond. After |
550 |
> |
1.5~ns, The main body of water in the system has already collected |
551 |
> |
into a central water channel. By 10.0~ns, the channel has widened |
552 |
> |
slightly, but there are still many sub channels permeating the lipid |
553 |
> |
macro-structure. At 23.0~ns, the central water channel has stabilized |
554 |
> |
and several smaller water channels have been absorbed to main |
555 |
> |
one. However, there are still several other channels that persist |
556 |
> |
through the lipid structure. |
557 |
> |
|
558 |
> |
\begin{figure} |
559 |
> |
\centering |
560 |
> |
\includegraphics[width=\linewidth]{iRodImage.eps} |
561 |
> |
\caption[Image sequence of the inverted rod aggregation]{Image sequence of the inverted rod aggregation. color scheme is the same as in Fig.~\ref{lipidFig:blImage}.} |
562 |
> |
\label{lipidFig:iRimage} |
563 |
> |
\end{figure} |
564 |
> |
|
565 |
> |
\section{\label{lipidSec:Conclusion}Conclusion} |
566 |
> |
|
567 |
> |
We have presented a phospholipid model capable of spontaneous |
568 |
> |
aggregation into a bilayer and an inverted rod structure. The time |
569 |
> |
scales of the macro-molecular aggregations are in excess of 24~ns. In |
570 |
> |
addition the model's bilayer properties have been explored over a |
571 |
> |
range of temperatures through prefabricated bilayers. No freezing |
572 |
> |
transition is seen in the temperature range of our current |
573 |
> |
simulations. However, structural information from the lowest |
574 |
> |
temperature may imply that a freezing event is on a much longer time |
575 |
> |
scale than that explored in this current research. Further studies of |
576 |
> |
this system could extend the time length of the simulations at the low |
577 |
> |
temperatures to observe whether lipid crystallization occurs within the |
578 |
> |
framework of this model. |