| 1 |
|
|
| 2 |
|
|
| 3 |
< |
\chapter{\label{chapt:lipid}Phospholipid Simulations} |
| 3 |
> |
\chapter{\label{chapt:lipid}PHOSPHOLIPID SIMULATIONS} |
| 4 |
|
|
| 5 |
|
\section{\label{lipidSec:Intro}Introduction} |
| 6 |
|
|
| 7 |
|
In the past 10 years, increasing computer speeds have allowed for the |
| 8 |
|
atomistic simulation of phospholipid bilayers for increasingly |
| 9 |
< |
relevant lenghths of time. These simulations have ranged from |
| 10 |
< |
simulation of the gel phase ($L_{\beta}$) of |
| 9 |
> |
relevant lengths of time. These simulations have ranged from |
| 10 |
> |
simulation of the gel ($L_{\beta}$) phase 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 |
| 13 |
> |
($L_{\alpha}$) bilayers.\cite{Marrink01} With the exception of a few |
| 14 |
|
ambitious |
| 15 |
< |
simulations,\cite{marrink01:undulation,marrink:2002,lindahl00} most |
| 15 |
> |
simulations, \cite{marrink01:undulation,marrink:2002,lindahl00} most |
| 16 |
|
investigations are limited to a range of 64 to 256 |
| 17 |
< |
phospholipids.\cite{lindahl00,sum:2003,venable00,gomez:2003,smondyrev:1999,marrink01} |
| 17 |
> |
phospholipids.\cite{lindahl00,sum:2003,venable00,gomez:2003,smondyrev:1999,Marrink01} |
| 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 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. |
| 24 |
> |
computationally expensive, direct Ewald sum or the slightly faster particle |
| 25 |
> |
mesh Ewald (PME) sum.\cite{nina:2002,norberg:2000,patra:2003} These factors |
| 26 |
> |
all limit the system size and time scales of bilayer simulations. |
| 27 |
|
|
| 28 |
|
Unfortunately, much of biological interest happens on time and length |
| 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. |
| 29 |
> |
scales well beyond the range of current simulation technologies. One |
| 30 |
> |
such example is the observance of a ripple ($P_{\beta^{\prime}}$) |
| 31 |
> |
phase which appears between the $L_{\beta}$ and $L_{\alpha}$ phases of |
| 32 |
> |
certain phospholipid bilayers |
| 33 |
> |
(Fig.~\ref{lipidFig:phaseDiag}).\cite{katsaras00,sengupta00} These |
| 34 |
> |
ripples are known from x-ray diffraction data to have periodicities on |
| 35 |
> |
the order of 100-200~$\mbox{\AA}$.\cite{katsaras00} A simulation on |
| 36 |
> |
this length scale would have approximately 1,300 lipid molecules with |
| 37 |
> |
an additional 25 water molecules per lipid to fully solvate the |
| 38 |
> |
bilayer. A simulation of this size is impractical with current |
| 39 |
> |
atomistic models. |
| 40 |
|
|
| 41 |
+ |
\begin{figure} |
| 42 |
+ |
\centering |
| 43 |
+ |
\includegraphics[width=\linewidth]{ripple.eps} |
| 44 |
+ |
\caption[Diagram of the bilayer gel to fluid phase transition]{Diagram showing the $P_{\beta^{\prime}}$ phase as a transition between the $L_{\beta}$ and $L_{\alpha}$ phases} |
| 45 |
+ |
\label{lipidFig:phaseDiag} |
| 46 |
+ |
\end{figure} |
| 47 |
+ |
|
| 48 |
|
The time and length scale limitations are most striking in transport |
| 49 |
< |
phenomena. Due to the fluid-like properties of a lipid membrane, not |
| 50 |
< |
all diffusion across the membrane happens at pores. Some molecules of |
| 51 |
< |
interest may incorporate themselves directly into the membrane. Once |
| 52 |
< |
here, they may possess an appreciable waiting time (on the order of |
| 53 |
< |
10's to 100's of nanoseconds) within the bilayer. Such long simulation |
| 54 |
< |
times are nearly impossible to obtain when integrating the system with |
| 55 |
< |
atomistic detail. |
| 49 |
> |
phenomena. Due to the fluid-like properties of lipid membranes, not |
| 50 |
> |
all small molecule diffusion across the membranes happens at pores. |
| 51 |
> |
Some molecules of interest may incorporate themselves directly into |
| 52 |
> |
the membrane. Once there, they may exhibit appreciable waiting times |
| 53 |
> |
(on the order of 10's to 100's of nanoseconds) within the |
| 54 |
> |
bilayer. Such long simulation times are nearly impossible to obtain |
| 55 |
> |
when integrating the system with atomistic detail. |
| 56 |
|
|
| 57 |
|
To address these issues, several schemes have been proposed. One |
| 58 |
< |
approach by Goetz and Liposky\cite{goetz98} is to model the entire |
| 58 |
> |
approach by Goetz and Lipowsky\cite{goetz98} is to model the entire |
| 59 |
|
system as Lennard-Jones spheres. Phospholipids are represented by |
| 60 |
|
chains of beads with the top most beads identified as the head |
| 61 |
|
atoms. Polar and non-polar interactions are mimicked through |
| 62 |
< |
attractive and soft-repulsive potentials respectively. A similar |
| 63 |
< |
model proposed by Marrinck \emph{et. al}.\cite{marrink04}~uses a |
| 64 |
< |
similar technique for modeling polar and non-polar interactions with |
| 62 |
> |
attractive and soft-repulsive potentials respectively. A model |
| 63 |
> |
proposed by Marrinck \emph{et. al}.\cite{marrink04}~uses a similar |
| 64 |
> |
technique for modeling polar and non-polar interactions with |
| 65 |
|
Lennard-Jones spheres. However, they also include charges on the head |
| 66 |
|
group spheres to mimic the electrostatic interactions of the |
| 67 |
< |
bilayer. While the solvent spheres are kept charge-neutral and |
| 67 |
> |
bilayer. The solvent spheres are kept charge-neutral and |
| 68 |
|
interact with the bilayer solely through an attractive Lennard-Jones |
| 69 |
|
potential. |
| 70 |
|
|
| 71 |
|
The model used in this investigation adds more information to the |
| 72 |
|
interactions than the previous two models, while still balancing the |
| 73 |
< |
need for simplifications over atomistic detail. The model uses |
| 74 |
< |
Lennard-Jones spheres for the head and tail groups of the |
| 73 |
> |
need for simplification of atomistic detail. The model uses |
| 74 |
> |
unified-atom Lennard-Jones spheres for the head and tail groups of the |
| 75 |
|
phospholipids, allowing for the ability to scale the parameters to |
| 76 |
|
reflect various sized chain configurations while keeping the number of |
| 77 |
|
interactions small. What sets this model apart, however, is the use |
| 78 |
|
of dipoles to represent the electrostatic nature of the |
| 79 |
|
phospholipids. The dipole electrostatic interaction is shorter range |
| 80 |
< |
than Coulombic ($\frac{1}{r^3}$ versus $\frac{1}{r}$), and eliminates |
| 81 |
< |
the need for a costly Ewald sum. |
| 80 |
> |
than Coulombic ($\frac{1}{r^3}$ versus $\frac{1}{r}$), and therefore |
| 81 |
> |
eliminates the need for the costly Ewald sum. |
| 82 |
|
|
| 83 |
< |
Another key feature of this model, is the use of a dipolar water model |
| 83 |
> |
Another key feature of this model is the use of a dipolar water model |
| 84 |
|
to represent the solvent. The soft sticky dipole ({\sc ssd}) water |
| 85 |
|
\cite{liu96:new_model} relies on the dipole for long range electrostatic |
| 86 |
|
effects, but also contains a short range correction for hydrogen |
| 87 |
< |
bonding. In this way the systems in this research mimic the entropic |
| 88 |
< |
contribution to the hydrophobic effect due to hydrogen-bond network |
| 89 |
< |
deformation around a non-polar entity, \emph{i.e.}~the phospholipid |
| 90 |
< |
molecules. |
| 87 |
> |
bonding. In this way the simulated systems in this research mimic the |
| 88 |
> |
entropic contribution to the hydrophobic effect due to hydrogen-bond |
| 89 |
> |
network deformation around a non-polar entity, \emph{i.e.}~the |
| 90 |
> |
phospholipid molecules. This effect has been missing from previous |
| 91 |
> |
reduced models. |
| 92 |
|
|
| 93 |
|
The following is an outline of this chapter. |
| 94 |
|
Sec.~\ref{lipidSec:Methods} is an introduction to the lipid model used |
| 124 |
|
their potentials can be scaled by $\sigma_{\text{tail}}$ and |
| 125 |
|
$\epsilon_{\text{tail}}$. |
| 126 |
|
|
| 127 |
< |
The long range interactions between lipids are given by the following: |
| 127 |
> |
The possible long range interactions between atomic groups in the |
| 128 |
> |
lipids are given by the following: |
| 129 |
|
\begin{equation} |
| 130 |
|
V_{\text{LJ}}(r_{ij}) = |
| 131 |
|
4\epsilon_{ij} \biggl[ |
| 132 |
|
\biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12} |
| 133 |
|
- \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6} |
| 134 |
< |
\biggr] |
| 134 |
> |
\biggr], |
| 135 |
|
\label{lipidEq:LJpot} |
| 136 |
|
\end{equation} |
| 137 |
|
and |
| 141 |
|
\boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j} |
| 142 |
|
- |
| 143 |
|
3(\boldsymbol{\hat{u}}_i \cdot \mathbf{\hat{r}}_{ij}) % |
| 144 |
< |
(\boldsymbol{\hat{u}}_j \cdot \mathbf{\hat{r}}_{ij} \biggr] |
| 144 |
> |
(\boldsymbol{\hat{u}}_j \cdot \mathbf{\hat{r}}_{ij} )\biggr]. |
| 145 |
|
\label{lipidEq:dipolePot} |
| 146 |
|
\end{equation} |
| 147 |
< |
Where $V_{\text{LJ}}$ is the Lennard-Jones potential and |
| 147 |
> |
Here $V_{\text{LJ}}$ is the Lennard-Jones potential and |
| 148 |
|
$V_{\text{dipole}}$ is the dipole-dipole potential. As previously |
| 149 |
|
stated, $\sigma_{ij}$ and $\epsilon_{ij}$ are the Lennard-Jones |
| 150 |
|
parameters which scale the length and depth of the interaction |
| 154 |
|
and $\mathbf{\Omega}_j$ are the orientational degrees of freedom for |
| 155 |
|
beads $i$ and $j$. $|\mu_i|$ is the magnitude of the dipole moment of |
| 156 |
|
$i$, and $\boldsymbol{\hat{u}}_i$ is the standard unit orientation |
| 157 |
< |
vector rotated with euler angles: $\boldsymbol{\Omega}_i$. |
| 157 |
> |
vector rotated with Euler angles: $\boldsymbol{\Omega}_i$. |
| 158 |
|
|
| 159 |
< |
The model also allows for the bonded interactions bends, and torsions. |
| 160 |
< |
The bond between two beads on a chain is of fixed length, and is |
| 161 |
< |
maintained according to the {\sc rattle} algorithm.\cite{andersen83} |
| 162 |
< |
The bends are subject to a harmonic potential: |
| 159 |
> |
The model also allows for the intra-molecular bend and torsion |
| 160 |
> |
interactions. The bond between two beads on a chain is of fixed |
| 161 |
> |
length, and is maintained using the {\sc rattle} |
| 162 |
> |
algorithm.\cite{andersen83} The bends are subject to a harmonic |
| 163 |
> |
potential: |
| 164 |
|
\begin{equation} |
| 165 |
< |
V_{\text{bend}}(\theta) = k_{\theta}( \theta - \theta_0 )^2 |
| 165 |
> |
V_{\text{bend}}(\theta) = k_{\theta}( \theta - \theta_0 )^2, |
| 166 |
|
\label{lipidEq:bendPot} |
| 167 |
|
\end{equation} |
| 168 |
|
where $k_{\theta}$ scales the strength of the harmonic well, and |
| 169 |
|
$\theta$ is the angle between bond vectors |
| 170 |
|
(Fig.~\ref{lipidFig:lipidModel}). In addition, we have placed a |
| 171 |
< |
``ghost'' bend on the phospholipid head. The ghost bend adds a |
| 172 |
< |
potential to keep the dipole pointed along the bilayer surface, where |
| 173 |
< |
$\theta$ is now the angle the dipole makes with respect to the {\sc |
| 174 |
< |
head}-$\text{{\sc ch}}_2$ bond vector |
| 175 |
< |
(Fig.~\ref{lipidFig:ghostBend}). The torsion potential is given by: |
| 171 |
> |
``ghost'' bend on the phospholipid head. The ghost bend is a bend |
| 172 |
> |
potential which keeps the dipole roughly perpendicular to the |
| 173 |
> |
molecular body, where $\theta$ is now the angle the dipole makes with |
| 174 |
> |
respect to the {\sc head}-$\text{{\sc ch}}_2$ bond vector |
| 175 |
> |
(Fig.~\ref{lipidFig:ghostBend}). This bend mimics the hinge between |
| 176 |
> |
the phosphatidyl part of the PC head group and the remainder of the |
| 177 |
> |
molecule. The torsion potential is given by: |
| 178 |
|
\begin{equation} |
| 179 |
|
V_{\text{torsion}}(\phi) = |
| 180 |
< |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0 |
| 180 |
> |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0. |
| 181 |
|
\label{lipidEq:torsionPot} |
| 182 |
|
\end{equation} |
| 183 |
|
Here, the parameters $k_0$, $k_1$, $k_2$, and $k_3$ fit the cosine |
| 185 |
|
the angle the two end atoms have rotated about the middle bond |
| 186 |
|
(Fig.:\ref{lipidFig:lipidModel}). Long range interactions such as the |
| 187 |
|
Lennard-Jones potential are excluded for atom pairs involved in the |
| 188 |
< |
same bond, bend, or torsion. However, internal interactions not |
| 189 |
< |
directly involved in a bonded pair are calculated. |
| 188 |
> |
same bond, bend, or torsion. However, long-range interactions for |
| 189 |
> |
pairs of atoms not directly involved in a bond, bend, or torsion are |
| 190 |
> |
calculated. |
| 191 |
|
|
| 192 |
|
\begin{figure} |
| 193 |
|
\centering |
| 194 |
< |
\includegraphics[width=\linewidth]{ghostBendFig.eps} |
| 195 |
< |
\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.} |
| 194 |
> |
\includegraphics[width=0.5\linewidth]{ghostBendFig.eps} |
| 195 |
> |
\caption[Depiction of the ``ghost'' bend]{The ``ghost'' bend is a bend potential added to restrain 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.} |
| 196 |
|
\label{lipidFig:ghostBend} |
| 197 |
|
\end{figure} |
| 198 |
|
|
| 199 |
< |
All simulations presented here use a two chained lipid as pictured in |
| 199 |
> |
All simulations presented here use a two-chain lipid as pictured in |
| 200 |
|
Fig.~\ref{lipidFig:lipidModel}. The chains are both eight beads long, |
| 201 |
|
and their mass and Lennard Jones parameters are summarized in |
| 202 |
|
Table~\ref{lipidTable:tcLJParams}. The magnitude of the dipole moment |
| 203 |
< |
for the head bead is 10.6~Debye, and the bend and torsion parameters |
| 204 |
< |
are summarized in Table~\ref{lipidTable:tcBendParams} and |
| 203 |
> |
for the head bead is 10.6~Debye (approximately half the magnitude of |
| 204 |
> |
the dipole on the PC head group\cite{Cevc87}), and the bend and |
| 205 |
> |
torsion parameters are summarized in |
| 206 |
> |
Table~\ref{lipidTable:tcBendParams} and |
| 207 |
|
\ref{lipidTable:tcTorsionParams}. |
| 208 |
|
|
| 209 |
|
\begin{table} |
| 210 |
< |
\caption{The Lennard Jones Parameters for the two chain phospholipids.} |
| 210 |
> |
\caption{THE LENNARD JONES PARAMETERS FOR THE TWO CHAIN PHOSPHOLIPIDS} |
| 211 |
|
\label{lipidTable:tcLJParams} |
| 212 |
|
\begin{center} |
| 213 |
< |
\begin{tabular}{|l|c|c|c|} |
| 213 |
> |
\begin{tabular}{|l|c|c|c|c|} |
| 214 |
|
\hline |
| 215 |
< |
& mass (amu) & $\sigma$($\mbox{\AA}$) & $\epsilon$ (kcal/mol) \\ \hline |
| 216 |
< |
{\sc head} & 72 & 4.0 & 0.185 \\ \hline |
| 217 |
< |
{\sc ch}\cite{Siepmann1998} & 13.02 & 4.0 & 0.0189 \\ \hline |
| 218 |
< |
$\text{{\sc ch}}_2$\cite{Siepmann1998} & 14.03 & 3.95 & 0.18 \\ \hline |
| 219 |
< |
$\text{{\sc ch}}_3$\cite{Siepmann1998} & 15.04 & 3.75 & 0.25 \\ \hline |
| 220 |
< |
{\sc ssd}\cite{liu96:new_model} & 18.03 & 3.051 & 0.152 \\ \hline |
| 215 |
> |
& mass (amu) & $\sigma$($\mbox{\AA}$) & $\epsilon$ (kcal/mol) % |
| 216 |
> |
& $|\mathbf{\hat{\mu}}|$ (Debye) \\ \hline |
| 217 |
> |
{\sc head} & 72 & 4.0 & 0.185 & 10.6 \\ \hline |
| 218 |
> |
{\sc ch}\cite{Siepmann1998} & 13.02 & 4.0 & 0.0189 & 0.0 \\ \hline |
| 219 |
> |
$\text{{\sc ch}}_2$\cite{Siepmann1998} & 14.03 & 3.95 & 0.18 & 0.0 \\ \hline |
| 220 |
> |
$\text{{\sc ch}}_3$\cite{Siepmann1998} & 15.04 & 3.75 & 0.25 & 0.0 \\ \hline |
| 221 |
> |
{\sc ssd}\cite{liu96:new_model} & 18.03 & 3.051 & 0.152 & 0.0 \\ \hline |
| 222 |
|
\end{tabular} |
| 223 |
|
\end{center} |
| 224 |
|
\end{table} |
| 225 |
|
|
| 226 |
|
\begin{table} |
| 227 |
< |
\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}} |
| 227 |
> |
\caption{BEND PARAMETERS FOR THE TWO CHAIN PHOSPHOLIPIDS} |
| 228 |
|
\label{lipidTable:tcBendParams} |
| 229 |
|
\begin{center} |
| 230 |
|
\begin{tabular}{|l|c|c|} |
| 234 |
|
$x$-{\sc ch}-$y$ & 58.84 & 112.0 \\ \hline |
| 235 |
|
$x$-$\text{{\sc ch}}_2$-$y$ & 58.84 & 114.0 \\ \hline |
| 236 |
|
\end{tabular} |
| 237 |
+ |
\begin{minipage}{\linewidth} |
| 238 |
+ |
\begin{center} |
| 239 |
+ |
\vspace{2mm} |
| 240 |
+ |
All alkane parameters are based off of those from TraPPE.\cite{Siepmann1998} |
| 241 |
|
\end{center} |
| 242 |
+ |
\end{minipage} |
| 243 |
+ |
\end{center} |
| 244 |
|
\end{table} |
| 245 |
|
|
| 246 |
|
\begin{table} |
| 247 |
< |
\caption[Torsion Parameters for the two chain phospholipids]{Torsion Parameters for the two chain phospholipids. Alkane parameters based on TraPPE.\cite{Siepmann1998}} |
| 247 |
> |
\caption{TORSION PARAMETERS FOR THE TWO CHAIN PHOSPHOLIPIDS} |
| 248 |
|
\label{lipidTable:tcTorsionParams} |
| 249 |
|
\begin{center} |
| 250 |
|
\begin{tabular}{|l|c|c|c|c|} |
| 253 |
|
$x$-{\sc ch}-$y$-$z$ & 3.3254 & -0.4215 & -1.686 & 1.1661 \\ \hline |
| 254 |
|
$x$-$\text{{\sc ch}}_2$-$\text{{\sc ch}}_2$-$y$ & 5.9602 & -0.568 & -3.802 & 2.1586 \\ \hline |
| 255 |
|
\end{tabular} |
| 256 |
+ |
\begin{minipage}{\linewidth} |
| 257 |
+ |
\begin{center} |
| 258 |
+ |
\vspace{2mm} |
| 259 |
+ |
All alkane parameters are based off of those from TraPPE.\cite{Siepmann1998} |
| 260 |
|
\end{center} |
| 261 |
+ |
\end{minipage} |
| 262 |
+ |
\end{center} |
| 263 |
|
\end{table} |
| 264 |
|
|
| 265 |
|
|
| 266 |
|
\section{\label{lipidSec:furtherMethod}Further Methodology} |
| 267 |
|
|
| 268 |
|
As mentioned previously, the water model used throughout these |
| 269 |
< |
simulations was the {\sc ssd} model of |
| 270 |
< |
Ichiye.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md} A |
| 269 |
> |
simulations was the {\sc ssd/e} model of Fennell and |
| 270 |
> |
Gezelter,\cite{fennell04} earlier forms of this model can be found in |
| 271 |
> |
Ichiye \emph{et |
| 272 |
> |
al}.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md} A |
| 273 |
|
discussion of the model can be found in Sec.~\ref{oopseSec:SSD}. As |
| 274 |
|
for the integration of the equations of motion, all simulations were |
| 275 |
|
performed in an orthorhombic periodic box with a thermostat on |
| 276 |
|
velocities, and an independent barostat on each Cartesian axis $x$, |
| 277 |
< |
$y$, and $z$. This is the $\text{NPT}_{xyz}$. ensemble described in |
| 278 |
< |
Sec.~\ref{oopseSec:Ensembles}. |
| 277 |
> |
$y$, and $z$. This is the $\text{NPT}_{xyz}$. integrator described in |
| 278 |
> |
Sec.~\ref{oopseSec:integrate}. With $\tau_B = 1.5$~ps and $\tau_T = |
| 279 |
> |
1.2$~ps, the box volume stabilizes after 50~ps, and fluctuates about |
| 280 |
> |
its equilibrium value by $\sim 0.6\%$, temperature fluctuations are |
| 281 |
> |
about $\sim 1.4\%$ of their set value, and pressure fluctuations are |
| 282 |
> |
the largest, varying as much as $\pm 250$~atm. However, such large |
| 283 |
> |
fluctuations in pressure are typical for liquid state simulations.\cite{leach01:mm} |
| 284 |
|
|
| 285 |
|
|
| 286 |
|
\subsection{\label{lipidSec:ExpSetup}Experimental Setup} |
| 287 |
|
|
| 288 |
< |
Two main starting configuration classes were used in this research: |
| 289 |
< |
random and ordered bilayers. The ordered bilayer starting |
| 290 |
< |
configurations were all started from an equilibrated bilayer at |
| 291 |
< |
300~K. The original configuration for the first 300~K run was |
| 292 |
< |
assembled by placing the phospholipids centers of mass on a planar |
| 293 |
< |
hexagonal lattice. The lipids were oriented with their long axis |
| 294 |
< |
perpendicular to the plane. The second leaf simply mirrored the first |
| 295 |
< |
leaf, and the appropriate number of waters were then added above and |
| 259 |
< |
below the bilayer. |
| 288 |
> |
Two main classes of starting configurations were used in this research: |
| 289 |
> |
random and ordered bilayers. The ordered bilayer simulations were all |
| 290 |
> |
started from an equilibrated bilayer configuration at 300~K. The original |
| 291 |
> |
configuration for the first 300~K run was assembled by placing the |
| 292 |
> |
phospholipids centers of mass on a planar hexagonal lattice. The |
| 293 |
> |
lipids were oriented with their principal axis perpendicular to the plane. |
| 294 |
> |
The bottom leaf simply mirrored the top leaf, and the appropriate |
| 295 |
> |
number of water molecules were then added above and below the bilayer. |
| 296 |
|
|
| 297 |
|
The random configurations took more work to generate. To begin, a |
| 298 |
|
test lipid was placed in a simulation box already containing water at |
| 299 |
< |
the intended density. The waters were then tested for overlap with |
| 300 |
< |
the lipid using a 5.0~$\mbox{\AA}$ buffer distance. This gave an |
| 301 |
< |
estimate for the number of waters each lipid would displace in a |
| 302 |
< |
simulation box. A target number of waters was then defined which |
| 303 |
< |
included the number of waters each lipid would displace, the number of |
| 304 |
< |
waters desired to solvate each lipid, and a factor to pad the |
| 305 |
< |
initial box with a few extra water molecules. |
| 299 |
> |
the intended density. The water molecules were then tested against |
| 300 |
> |
the lipid using a 5.0~$\mbox{\AA}$ overlap test with any atom in the |
| 301 |
> |
lipid. This gave an estimate for the number of water molecules each |
| 302 |
> |
lipid would displace in a simulation box. A target number of water |
| 303 |
> |
molecules was then defined which included the number of water |
| 304 |
> |
molecules each lipid would displace, the number of water molecules |
| 305 |
> |
desired to solvate each lipid, and a factor to pad the initial box |
| 306 |
> |
with a few extra water molecules. |
| 307 |
|
|
| 308 |
|
Next, a cubic simulation box was created that contained at least the |
| 309 |
< |
target number of waters in an FCC lattice (the lattice was for ease of |
| 309 |
> |
target number of water molecules in an FCC lattice (the lattice was for ease of |
| 310 |
|
placement). What followed was a RSA simulation similar to those of |
| 311 |
|
Chapt.~\ref{chapt:RSA}. The lipids were sequentially given a random |
| 312 |
|
position and orientation within the box. If a lipid's position caused |
| 313 |
|
atomic overlap with any previously placed lipid, its position and |
| 314 |
|
orientation were rejected, and a new random placement site was |
| 315 |
|
attempted. The RSA simulation proceeded until all phospholipids had |
| 316 |
< |
been adsorbed. After placement of all lipid molecules, water |
| 316 |
> |
been placed. After placement of all lipid molecules, water |
| 317 |
|
molecules with locations that overlapped with the atomic coordinates |
| 318 |
|
of the lipids were removed. |
| 319 |
|
|
| 320 |
< |
Finally, water molecules were removed at random until the desired |
| 321 |
< |
number of waters per lipid was reached. The typical low final density |
| 322 |
< |
for these initial configurations was not a problem, as the box shrinks |
| 323 |
< |
to an appropriate size within the first 50~ps of a simulation in the |
| 324 |
< |
$\text{NPT}_{xyz}$ ensemble. |
| 320 |
> |
Finally, water molecules were removed at random until the desired water |
| 321 |
> |
to lipid ratio was achieved. The typical low final density for these |
| 322 |
> |
initial configurations was not a problem, as the box shrinks to an |
| 323 |
> |
appropriate size within the first 50~ps of a simulation under the |
| 324 |
> |
NPTxyz integrator. |
| 325 |
|
|
| 326 |
|
\subsection{\label{lipidSec:Configs}Configurations} |
| 327 |
|
|
| 328 |
< |
The first class of simulations were started from ordered |
| 329 |
< |
bilayers. They were all configurations consisting of 60 lipid |
| 330 |
< |
molecules with 30 lipids on each leaf, and were hydrated with 1620 |
| 331 |
< |
{\sc ssd} molecules. The original configuration was assembled |
| 332 |
< |
according to Sec.~\ref{lipidSec:ExpSetup} and simulated for a length |
| 333 |
< |
of 10~ns at 300~K. The other temperature runs were started from a |
| 334 |
< |
frame 7~ns into the 300~K simulation. Their temperatures were reset |
| 335 |
< |
with the thermostating algorithm in the $\text{NPT}_{xyz}$ |
| 336 |
< |
integrator. All of the temperature variants were also run for 10~ns, |
| 337 |
< |
with only the last 5~ns being used for accumulation of statistics. |
| 328 |
> |
The first class of simulations were started from ordered bilayers. All |
| 329 |
> |
configurations consisted of 60 lipid molecules with 30 lipids on each |
| 330 |
> |
leaf, and were hydrated with 1620 {\sc ssd/e} molecules. The original |
| 331 |
> |
configuration was assembled according to Sec.~\ref{lipidSec:ExpSetup} |
| 332 |
> |
and simulated for a length of 10~ns at 300~K. The other temperature |
| 333 |
> |
runs were started from a configuration 7~ns in to the 300~K |
| 334 |
> |
simulation. Their temperatures were modified with the thermostatting |
| 335 |
> |
algorithm in the NPTxyz integrator. All of the temperature variants |
| 336 |
> |
were also run for 10~ns, with only the last 5~ns being used for |
| 337 |
> |
accumulation of statistics. |
| 338 |
|
|
| 339 |
|
The second class of simulations were two configurations started from |
| 340 |
|
randomly dispersed lipids in a ``gas'' of water. The first |
| 341 |
|
($\text{R}_{\text{I}}$) was a simulation containing 72 lipids with |
| 342 |
< |
1800 {\sc ssd} molecules simulated at 300~K. The second |
| 343 |
< |
($\text{R}_{\text{II}}$) was 90 lipids with 1350 {\sc ssd} molecules |
| 342 |
> |
1800 {\sc ssd/e} molecules simulated at 300~K. The second |
| 343 |
> |
($\text{R}_{\text{II}}$) was 90 lipids with 1350 {\sc ssd/e} molecules |
| 344 |
|
simulated at 350~K. Both simulations were integrated for more than |
| 345 |
< |
20~ns, and illustrate the spontaneous aggregation of the lipid model |
| 346 |
< |
into phospholipid macrostructures: $\text{R}_{\text{I}}$ into a |
| 347 |
< |
bilayer, and $\text{R}_{\text{II}}$ into a inverted rod. |
| 345 |
> |
20~ns to observe whether our model is capable of spontaneous |
| 346 |
> |
aggregation into known phospholipid macro-structures: |
| 347 |
> |
$\text{R}_{\text{I}}$ into a bilayer, and $\text{R}_{\text{II}}$ into |
| 348 |
> |
a inverted rod. |
| 349 |
|
|
| 350 |
|
\section{\label{lipidSec:resultsDis}Results and Discussion} |
| 351 |
|
|
| 352 |
< |
\subsection{\label{lipidSec:scd}$\text{S}_{\text{{\sc cd}}}$ order parameters} |
| 352 |
> |
\subsection{\label{lipidSec:densProf}Density Profile} |
| 353 |
> |
|
| 354 |
> |
Fig.~\ref{lipidFig:densityProfile} illustrates the densities of the |
| 355 |
> |
atoms in the bilayer systems normalized by the bulk density as a |
| 356 |
> |
function of distance from the center of the box. The profile is taken |
| 357 |
> |
along the bilayer normal (in this case the $z$ axis). The first interesting point to note is the penetration of water into the membrane. Water penetrates about 5~$\mbox{\AA}$ into the bilayer, completely solvating the head groups. This is common in atomistic and some coarse grain simulations of phospholipid bilayers.\cite{Marrink01,marrink04,klein01} It is an indication that the water molecules are very attracted to the polar head region, yet there is still enough of a hydrophobic effect to exclude water from the inside of the bilayer. |
| 358 |
|
|
| 359 |
+ |
Another interesting point is the fluidity of the chains. Although the ends of the tails, the $\text{{\sc ch}}_3$ atoms, are mostly concentrated at the centers of the bilayers, they have a significant density around the head regions. This indicates that there is much freedom of movement in the chains of our model. Typical atomistic simulations of DPPC show the terminal groups concentrated at the center of the bilayer.\cite{marrink03:vesicles} This is most likely an indication that our chain lengths are too small, and given longer chains, the tail groups would stay more deeply buried in the bilayer. |
| 360 |
+ |
|
| 361 |
+ |
The last point to consider is the splitting in the density peak of the {\sc head} atom at 270~K. This implies that there is some freezing of structure at this temperature. By 280~K, this feature is smoothed out, demonstrating a more fluid phase in the bilayer. Within the time scale of the simulation, the gel phase has not formed at 270~K, so this splitting in the peak is likely a glassy transition in the head groups, and could possibly indicate that we are simulating in a super cooled region of our phospholipid model. |
| 362 |
+ |
|
| 363 |
+ |
\begin{figure} |
| 364 |
+ |
\centering |
| 365 |
+ |
\includegraphics[width=\linewidth]{densityProfile.eps} |
| 366 |
+ |
\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.} |
| 367 |
+ |
\label{lipidFig:densityProfile} |
| 368 |
+ |
\end{figure} |
| 369 |
+ |
|
| 370 |
+ |
|
| 371 |
+ |
\subsection{\label{lipidSec:scd}$\text{S}_{\text{{\sc cd}}}$ Order Parameters} |
| 372 |
+ |
|
| 373 |
|
The $\text{S}_{\text{{\sc cd}}}$ order parameter is often reported in |
| 374 |
< |
the experimental charecterizations of phospholipids. It is obtained |
| 374 |
> |
the experimental characterizations of phospholipids. It is obtained |
| 375 |
|
through deuterium NMR, and measures the ordering of the carbon |
| 376 |
|
deuterium bond in relation to the bilayer normal at various points |
| 377 |
< |
along the chains. In our model, there are no explicit hydrogens, but |
| 377 |
> |
along the chains. The order parameter has a range of $[1,-\frac{1}{2}]$. A value of 1 |
| 378 |
> |
implies full order aligned to the bilayer axis, 0 implies full |
| 379 |
> |
disorder, and $-\frac{1}{2}$ implies full order perpendicular to the |
| 380 |
> |
bilayer axis. The {\sc cd} bond vector for carbons near the head group |
| 381 |
> |
are usually ordered perpendicular to the bilayer normal, with tails |
| 382 |
> |
farther away tending toward disorder. This makes the order parameter |
| 383 |
> |
negative for most carbons, and as such $|S_{\text{{\sc cd}}}|$ is more |
| 384 |
> |
commonly reported than $S_{\text{{\sc cd}}}$. |
| 385 |
> |
|
| 386 |
> |
In our model, there are no explicit hydrogens, but |
| 387 |
|
the order parameter can be written in terms of the carbon ordering at |
| 388 |
|
each point in the chain:\cite{egberts88} |
| 389 |
|
\begin{equation} |
| 390 |
< |
S_{\text{{\sc cd}}} = \frac{2}{3}S_{xx} + \frac{1}{3}S_{yy} |
| 390 |
> |
S_{\text{{\sc cd}}} = \frac{2}{3}S_{xx} + \frac{1}{3}S_{yy}, |
| 391 |
|
\label{lipidEq:scd1} |
| 392 |
|
\end{equation} |
| 393 |
< |
Where $S_{ij}$ is given by: |
| 393 |
> |
where $S_{ij}$ is given by: |
| 394 |
|
\begin{equation} |
| 395 |
< |
S_{ij} = \frac{1}{2}\Bigl<(3\cos\Theta_i\cos\Theta_j - \delta_{ij})\Bigr> |
| 395 |
> |
S_{ij} = \frac{1}{2}\Bigl\langle(3\cos\Theta_i\cos\Theta_j |
| 396 |
> |
- \delta_{ij})\Bigr\rangle. |
| 397 |
|
\label{lipidEq:scd2} |
| 398 |
|
\end{equation} |
| 399 |
< |
Here, $\Theta_i$ is the angle the $i$th carbon atom frame axis makes |
| 400 |
< |
with the bilayer normal. The brackets denote an average over time and |
| 401 |
< |
molecules. The carbon atom axes are defined: |
| 402 |
< |
$\mathbf{\hat{z}}\rightarrow$ vector from $C_{n-1}$ to $C_{n+1}$; |
| 403 |
< |
$\mathbf{\hat{y}}\rightarrow$ vector that is perpindicular to $z$ and |
| 404 |
< |
in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$; |
| 405 |
< |
$\mathbf{\hat{x}}\rightarrow$ vector perpindicular to |
| 399 |
> |
Here, $\Theta_i$ is the angle the $i$th axis in the reference frame of |
| 400 |
> |
the carbon atom makes with the bilayer normal. The brackets denote an |
| 401 |
> |
average over time and molecules. The carbon atom axes are defined: |
| 402 |
> |
\begin{itemize} |
| 403 |
> |
\item $\mathbf{\hat{z}}$ is the vector from $C_{n-1}$ to $C_{n+1}$. |
| 404 |
> |
\item $\mathbf{\hat{y}}$ is the vector that is perpendicular to $z$ and |
| 405 |
> |
in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$. |
| 406 |
> |
\item $\mathbf{\hat{x}}$ is the vector perpendicular to |
| 407 |
|
$\mathbf{\hat{y}}$ and $\mathbf{\hat{z}}$. |
| 408 |
+ |
\end{itemize} |
| 409 |
+ |
This assumes that the hydrogen atoms are always in a plane |
| 410 |
+ |
perpendicular to the $C_{n-1}-C_{n}-C_{n+1}$ plane. |
| 411 |
|
|
| 412 |
< |
The order parameter has a range of $[1,-\frac{1}{2}]$. A value of 1 |
| 413 |
< |
implies full order aligned to the bilayer axis, 0 implies full |
| 414 |
< |
disorder, and $-\frac{1}{2}$ implies full order perpindicular to the |
| 415 |
< |
bilayer axis. The {\sc cd} bond vector for carbons near the head group |
| 416 |
< |
are usually ordered perpindicular to the bilayer normal, with tails |
| 417 |
< |
farther away tending toward disorder. This makes the order paramter |
| 418 |
< |
negative for most carbons, and as such $|S_{\text{{\sc cd}}}|$ is more |
| 419 |
< |
commonly reported than $S_{\text{{\sc cd}}}$. |
| 412 |
> |
Fig.~\ref{lipidFig:scdFig} shows the $S_{\text{{\sc cd}}}$ order |
| 413 |
> |
parameters for the bilayer system at 300~K. There is no appreciable |
| 414 |
> |
difference in the plots for the various temperatures, however, there |
| 415 |
> |
is a larger difference between our model's ordering, and the |
| 416 |
> |
experimentally observed ordering of DMPC. As our values are closer to |
| 417 |
> |
$-\frac{1}{2}$, this implies more ordering perpendicular to the normal |
| 418 |
> |
than in a real system. This is due to the model having only one carbon |
| 419 |
> |
group separating the chains from the top of the lipid. In DMPC, with |
| 420 |
> |
the flexibility inherent in a multiple atom head group, as well as a |
| 421 |
> |
glycerol linkage between the head group and the acyl chains, there is |
| 422 |
> |
more loss of ordering by the point when the chains start. Also, there is more ordering in the model due to the our assumptions about the locations of the hydrogen atoms. Our method assumes a rigid location for each hydrogen atom based on the carbon positions. This does not allow for any small fluctuations in their positions that would be inherent in an atomistic simulation or in experiments. These small fluctuations would serve to lower the ordering measured in the $S_{\text{{\sc cd}}}$. |
| 423 |
|
|
| 350 |
– |
|
| 351 |
– |
|
| 352 |
– |
|
| 424 |
|
\begin{figure} |
| 425 |
|
\centering |
| 426 |
|
\includegraphics[width=\linewidth]{scdFig.eps} |
| 428 |
|
\label{lipidFig:scdFig} |
| 429 |
|
\end{figure} |
| 430 |
|
|
| 431 |
+ |
\subsection{\label{lipidSec:p2Order}$P_2$ Order Parameter} |
| 432 |
|
|
| 433 |
< |
\begin{figure} |
| 434 |
< |
\centering |
| 435 |
< |
\includegraphics[width=\linewidth]{densityProfile.eps} |
| 436 |
< |
\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.} |
| 437 |
< |
\label{lipidFig:densityProfile} |
| 438 |
< |
\end{figure} |
| 433 |
> |
The $P_2$ order parameter allows us to measure the amount of |
| 434 |
> |
directional ordering that exists in the bodies of the molecules making |
| 435 |
> |
up the bilayer. Each lipid molecule can be thought of as a cylindrical |
| 436 |
> |
rod with the head group at the top. If all of the rods are perfectly |
| 437 |
> |
aligned, the $P_2$ order parameter will be $1.0$. If the rods are |
| 438 |
> |
completely disordered, the $P_2$ order parameter will be 0. For a |
| 439 |
> |
collection of unit vectors pointing along the principal axes of the |
| 440 |
> |
rods, the $P_2$ order parameter can be solved via the following |
| 441 |
> |
method.\cite{zannoni94} |
| 442 |
|
|
| 443 |
+ |
Define an ordering tensor $\overleftrightarrow{\mathsf{Q}}$, such that, |
| 444 |
+ |
\begin{equation} |
| 445 |
+ |
\overleftrightarrow{\mathsf{Q}} = \frac{1}{N}\sum_i^N % |
| 446 |
+ |
\begin{pmatrix} % |
| 447 |
+ |
u_{ix}u_{ix}-\frac{1}{3} & u_{ix}u_{iy} & u_{ix}u_{iz} \\ |
| 448 |
+ |
u_{iy}u_{ix} & u_{iy}u_{iy}-\frac{1}{3} & u_{iy}u_{iz} \\ |
| 449 |
+ |
u_{iz}u_{ix} & u_{iz}u_{iy} & u_{iz}u_{iz}-\frac{1}{3} % |
| 450 |
+ |
\end{pmatrix}, |
| 451 |
+ |
\label{lipidEq:po1} |
| 452 |
+ |
\end{equation} |
| 453 |
+ |
where the $u_{i\alpha}$ is the $\alpha$ element of the unit vector |
| 454 |
+ |
$\mathbf{\hat{u}}_i$, and the sum over $i$ averages over the whole |
| 455 |
+ |
collection of unit vectors. This allows the tensor to be written: |
| 456 |
+ |
\begin{equation} |
| 457 |
+ |
\overleftrightarrow{\mathsf{Q}} = \frac{1}{N}\sum_i^N \biggl[ |
| 458 |
+ |
\mathbf{\hat{u}}_i \otimes \mathbf{\hat{u}}_i |
| 459 |
+ |
- \frac{1}{3} \cdot \mathsf{1} \biggr]. |
| 460 |
+ |
\label{lipidEq:po2} |
| 461 |
+ |
\end{equation} |
| 462 |
|
|
| 463 |
+ |
After constructing the tensor, diagonalizing |
| 464 |
+ |
$\overleftrightarrow{\mathsf{Q}}$ yields three eigenvalues and |
| 465 |
+ |
eigenvectors. The eigenvector associated with the largest eigenvalue, |
| 466 |
+ |
$\lambda_{\text{max}}$, is the director axis for the system of unit |
| 467 |
+ |
vectors. The director axis is the average direction all of the unit vectors |
| 468 |
+ |
are pointing. The $P_2$ order parameter is then simply |
| 469 |
+ |
\begin{equation} |
| 470 |
+ |
\langle P_2 \rangle = \frac{3}{2}\lambda_{\text{max}}. |
| 471 |
+ |
\label{lipidEq:po3} |
| 472 |
+ |
\end{equation} |
| 473 |
|
|
| 474 |
< |
\begin{figure} |
| 475 |
< |
\centering |
| 476 |
< |
\includegraphics[width=\linewidth]{diffusionFig.eps} |
| 477 |
< |
\caption[The lateral difusion constants versus temperature]{The lateral diffusion constants for the bilayers as a function of temperature.} |
| 478 |
< |
\label{lipidFig:diffusionFig} |
| 479 |
< |
\end{figure} |
| 474 |
> |
Table~\ref{lipidTab:blSummary} summarizes the $P_2$ values for the |
| 475 |
> |
bilayers, as well as for the dipole orientations. The unit vector for the |
| 476 |
> |
lipid molecules was defined by finding the moment of inertia for each |
| 477 |
> |
lipid, then setting $\mathbf{\hat{u}}$ to point along the axis of |
| 478 |
> |
minimum inertia (the long axis). For the {\sc head} atoms, the unit vector simply |
| 479 |
> |
pointed in the same direction as the dipole moment. For the lipid |
| 480 |
> |
molecules, the ordering was consistent across all temperatures, with |
| 481 |
> |
the director pointed along the $z$ axis of the box. More |
| 482 |
> |
interestingly, is the high degree of ordering the dipoles impose on |
| 483 |
> |
the {\sc head} atoms. The directors for the dipoles themselves |
| 484 |
> |
consistently pointed along the plane of the bilayer, with head groups lining up in rows of alternating alignment. The ordering implies that the dipole interaction is a little too strong, or that perhaps the dipoles are allowed to approach each other a bit too closely. A possible change in future models would alter the size or shape of the head group to discourage too rigid ordering of the dipoles. |
| 485 |
|
|
| 486 |
|
\begin{table} |
| 487 |
< |
\caption[Structural properties of the bilayers]{Bilayer Structural properties as a function of temperature.} |
| 487 |
> |
\caption{BILAYER STRUCTURAL PROPERTIES AS A FUNCTION OF TEMPERATURE} |
| 488 |
> |
\label{lipidTab:blSummary} |
| 489 |
|
\begin{center} |
| 490 |
|
\begin{tabular}{|c|c|c|c|c|} |
| 491 |
|
\hline |
| 492 |
< |
Temperature (K) & $<L_{\perp}>$ ($\mbox{\AA}$) & % |
| 493 |
< |
$<A_{\parallel}>$ ($\mbox{\AA}^2$) & $<P_2>_{\text{Lipid}}$ & % |
| 494 |
< |
$<P_2>_{\text{{\sc head}}}$ \\ \hline |
| 492 |
> |
Temperature (K) & $\langle L_{\perp}\rangle$ ($\mbox{\AA}$) & % |
| 493 |
> |
$\langle A_{\parallel}\rangle$ ($\mbox{\AA}^2$) & % |
| 494 |
> |
$\langle P_2\rangle_{\text{Lipid}}$ & % |
| 495 |
> |
$\langle P_2\rangle_{\text{{\sc head}}}$ \\ \hline |
| 496 |
|
270 & 18.1 & 58.1 & 0.253 & 0.494 \\ \hline |
| 497 |
|
275 & 17.2 & 56.7 & 0.295 & 0.514 \\ \hline |
| 498 |
|
277 & 16.9 & 58.0 & 0.301 & 0.541 \\ \hline |
| 504 |
|
\end{tabular} |
| 505 |
|
\end{center} |
| 506 |
|
\end{table} |
| 507 |
+ |
|
| 508 |
+ |
\subsection{\label{lipidSec:miscData}Further Structural Data} |
| 509 |
+ |
|
| 510 |
+ |
Also summarized in Table~\ref{lipidTab:blSummary}, are the bilayer |
| 511 |
+ |
thickness ($\langle L_{\perp}\rangle$) and area per lipid ($\langle |
| 512 |
+ |
A_{\parallel}\rangle$). The bilayer thickness was measured from the |
| 513 |
+ |
peak to peak {\sc head} atom distance in the density profiles. The |
| 514 |
+ |
area per lipid data compares favorably with values typically seen for |
| 515 |
+ |
DMPC (60.0~$\mbox{\AA}^2$ at 303~K)\cite{petrache00}. Although our |
| 516 |
+ |
values are lower this is most likely due to the shorter chain length |
| 517 |
+ |
of our model (8 versus 14 for DMPC). |
| 518 |
+ |
|
| 519 |
+ |
\subsection{\label{lipidSec:diffusion}Lateral Diffusion Constants} |
| 520 |
+ |
|
| 521 |
+ |
The lateral diffusion constant, $D_L$, is the constant characterizing |
| 522 |
+ |
the diffusive motion of the lipid molecules within the plane of the bilayer. It |
| 523 |
+ |
is given by the following Einstein relation:\cite{allen87:csl} |
| 524 |
+ |
\begin{equation} |
| 525 |
+ |
D_L = \lim_{t\rightarrow\infty}\frac{1}{4t}\langle |\mathbf{r}(t) |
| 526 |
+ |
- \mathbf{r}(0)|^2\rangle, |
| 527 |
+ |
\end{equation} |
| 528 |
+ |
where $\mathbf{r}(t)$ is the $xy$ position of the lipid at time $t$ |
| 529 |
+ |
(assuming the $z$-axis is parallel to the bilayer normal). Calculating the $D_L$ involves first plotting the mean square displacement, $\langle |\mathbf{r}(t) - \mathbf{r}(0)|^2\rangle$, finding the slope at long times, and dividing the slope by 4 to give the diffusion constant (Fig.~\ref{lipidFig:msdFig}). When finding the slope only the 1~ns to 3~ns times are considered. Points at the longer times are not included due to the lack of good statistics at long time intervals. |
| 530 |
+ |
|
| 531 |
+ |
Fig.~\ref{lipidFig:diffusionFig} shows the lateral diffusion constants |
| 532 |
+ |
as a function of temperature. There is a definite increase in the |
| 533 |
+ |
lateral diffusion with higher temperatures, which is exactly what one |
| 534 |
+ |
would expect with greater fluidity of the chains. However, the |
| 535 |
+ |
diffusion constants are two orders of magnitude larger than those |
| 536 |
+ |
typical of DPPC ($\sim10^{-9}\text{cm}^2/\text{s}$ over this temperature range).\cite{Cevc87} This is what one would expect as the DPPC |
| 537 |
+ |
molecule is sterically larger and heavier than our model, indicating that further modifications to the model should increase the lengths of the tail chains, and perhaps explore larger, more massive head groups. In contrast, the diffusion constant of |
| 538 |
+ |
the {\sc ssd} water, $9.84\times 10^{-6}\,\text{cm}^2/\text{s}$, is |
| 539 |
+ |
reasonably close to the bulk water diffusion constant ($2.2999\times |
| 540 |
+ |
10^{-5}\,\text{cm}^2/\text{s}$).\cite{Holz00} |
| 541 |
+ |
|
| 542 |
+ |
\begin{figure} |
| 543 |
+ |
\centering |
| 544 |
+ |
\includegraphics[width=\linewidth]{msdFig.eps} |
| 545 |
+ |
\caption[Lateral mean square displacement for the phospholipid at 300~K]{This is a representative lateral mean square displacement for the center of mass motion of the phospholipid model. This particular example is from the 300~K simulation. The box is drawn about the region used in the calculation of the diffusion constant.} |
| 546 |
+ |
\label{lipidFig:msdFig} |
| 547 |
+ |
\end{figure} |
| 548 |
+ |
|
| 549 |
+ |
\begin{figure} |
| 550 |
+ |
\centering |
| 551 |
+ |
\includegraphics[width=\linewidth]{diffusionFig.eps} |
| 552 |
+ |
\caption[The lateral diffusion constants versus temperature]{The lateral diffusion constants for the bilayers as a function of temperature.} |
| 553 |
+ |
\label{lipidFig:diffusionFig} |
| 554 |
+ |
\end{figure} |
| 555 |
+ |
|
| 556 |
+ |
\subsection{\label{lipidSec:randBilayer}Bilayer Aggregation} |
| 557 |
+ |
|
| 558 |
+ |
A very important accomplishment for our model is its ability to |
| 559 |
+ |
spontaneously form bilayers from a randomly dispersed starting |
| 560 |
+ |
configuration. Fig.~\ref{lipidFig:blImage} shows an image sequence for |
| 561 |
+ |
the bilayer aggregation. After 1.0~ns, bulk aggregation has occured. By 5.0~ns, the basic bilayer aggregation can be seen, however there is a vertical lipid bridge connecting the periodic image of the bilayer to itself. At 15.0~ns, the lipid bridge has finally broken up, and the lipid molecules are starting to re-incorporate themselves into the bilayer. A water pore is still present through the membrane. In the last frame at 42.0~ns, the water pore is still present, although does show some signs of rejoining the bulk water section. These behaviors are typical for coarse grain model simulations, which can have lipid bridge lifetimes of up to 20~ns, and water pores typically lasting 3 to 25~ns.\cite{marrink04} |
| 562 |
+ |
|
| 563 |
+ |
\begin{figure} |
| 564 |
+ |
\centering |
| 565 |
+ |
\includegraphics[width=\linewidth]{bLayerImage.eps} |
| 566 |
+ |
\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.} |
| 567 |
+ |
\label{lipidFig:blImage} |
| 568 |
+ |
\end{figure} |
| 569 |
+ |
|
| 570 |
+ |
\subsection{\label{lipidSec:randIrod}Inverted Rod Aggregation} |
| 571 |
+ |
|
| 572 |
+ |
Fig.~\ref{lipidFig:iRimage} shows a second aggregation sequence |
| 573 |
+ |
simulated in this research. Here the fraction of water had been |
| 574 |
+ |
significantly decreased to observe how the model would respond. After |
| 575 |
+ |
1.5~ns, The main body of water in the system has already collected |
| 576 |
+ |
into a central water channel. By 10.0~ns, the channel has widened |
| 577 |
+ |
slightly, but there are still many water molecules permeating the |
| 578 |
+ |
lipid macro-structure. At 35.0~ns, the central water channel has |
| 579 |
+ |
stabilized and several smaller water channels have been absorbed by |
| 580 |
+ |
the main one. However, there is still an appreciable water |
| 581 |
+ |
concentration throughout the lipid structure. |
| 582 |
+ |
|
| 583 |
+ |
\begin{figure} |
| 584 |
+ |
\centering |
| 585 |
+ |
\includegraphics[width=\linewidth]{iRodImage.eps} |
| 586 |
+ |
\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}.} |
| 587 |
+ |
\label{lipidFig:iRimage} |
| 588 |
+ |
\end{figure} |
| 589 |
+ |
|
| 590 |
+ |
\section{\label{lipidSec:Conclusion}Conclusion} |
| 591 |
+ |
|
| 592 |
+ |
We have presented a simple unified-atom phospholipid model capable of |
| 593 |
+ |
spontaneous aggregation into a bilayer and an inverted rod |
| 594 |
+ |
structure. The time scales of the macro-molecular aggregations are |
| 595 |
+ |
approximately 24~ns, with water permeation of the structures persisting for times longer than the scope of both aggregations. In addition the model's properties have been |
| 596 |
+ |
explored over a range of temperatures through prefabricated |
| 597 |
+ |
bilayers. No freezing transition is seen in the temperature range of |
| 598 |
+ |
our current simulations. However, structural information from 270~K |
| 599 |
+ |
may imply that a freezing event is on a much longer time scale than |
| 600 |
+ |
that explored in this current research. Further studies of this system |
| 601 |
+ |
could extend the time length of the simulations at the low |
| 602 |
+ |
temperatures to observe whether lipid crystallization can occur within |
| 603 |
+ |
the framework of this model. |
| 604 |
+ |
|
| 605 |
+ |
Potential problems that may be obstacles in further research, is the |
| 606 |
+ |
lack of detail in the head region. As the chains are almost directly |
| 607 |
+ |
attached to the {\sc head} atom, there is no buffer between the |
| 608 |
+ |
actions of the head group and the tails. Another disadvantage of the |
| 609 |
+ |
model is the dipole approximation will alter results when details |
| 610 |
+ |
concerning a charged solute's interactions with the bilayer. However, |
| 611 |
+ |
it is important to keep in mind that the dipole approximation can be |
| 612 |
+ |
kept an advantage by examining solutes that do not require point |
| 613 |
+ |
charges, or at the least, require only dipole approximations |
| 614 |
+ |
themselves. Other advantages of the model include the ability to alter |
| 615 |
+ |
the size of the unified-atoms so that the size of the lipid can be |
| 616 |
+ |
increased without adding to the number of interactions in the |
| 617 |
+ |
system. However, what sets our model apart from other current |
| 618 |
+ |
simplified models,\cite{goetz98,marrink04} is the information gained |
| 619 |
+ |
by observing the ordering of the head groups dipole's in relation to |
| 620 |
+ |
each other and the solvent without the need for point charges and the |
| 621 |
+ |
Ewald sum. |
