| 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 |
| 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} |
| 18 |
|
The expense of the force calculations involved when performing these |
| 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[ |
| 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 |
| 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 |
| 166 |
|
\label{lipidEq:bendPot} |
| 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 |
| 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[Lennard-Jones parameters for the two chain phospholipids]{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]{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]{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. |
| 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 profile at |
| 358 |
> |
270~K shows several structural features that are largely smoothed out |
| 359 |
> |
at 300~K. The left peak for the {\sc head} atoms is split at 270~K, |
| 360 |
> |
implying that some freezing of the structure into a gel phase might |
| 361 |
> |
already be occurring at this temperature. However, movies of the |
| 362 |
> |
trajectories at this temperature show that the tails are very fluid, |
| 363 |
> |
and have not gelled. But this profile could indicate that a phase |
| 364 |
> |
transition may simply be beyond the time length of the current |
| 365 |
> |
simulation, and that given more time the system may tend towards a gel |
| 366 |
> |
phase. In all profiles, the water penetrates almost 5~$\mbox{\AA}$ |
| 367 |
> |
into the bilayer, completely solvating the {\sc head} atoms. The |
| 368 |
> |
$\text{{\sc ch}}_3$ atoms, although mainly centered at the middle of |
| 369 |
> |
the bilayer, show appreciable penetration into the head group |
| 370 |
> |
region. This indicates that the chains have enough flexibility to bend |
| 371 |
> |
back upward to allow the ends to explore areas around the {\sc head} |
| 372 |
> |
atoms. It is unlikely that this is penetration from a lipid of the |
| 373 |
> |
opposite face, as the lipids are only 12~$\mbox{\AA}$ in length, and |
| 374 |
> |
the typical leaf spacing as measured from the {\sc head-head} spacing |
| 375 |
> |
in the profile is 17.5~$\mbox{\AA}$. |
| 376 |
|
|
| 377 |
+ |
\begin{figure} |
| 378 |
+ |
\centering |
| 379 |
+ |
\includegraphics[width=\linewidth]{densityProfile.eps} |
| 380 |
+ |
\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.} |
| 381 |
+ |
\label{lipidFig:densityProfile} |
| 382 |
+ |
\end{figure} |
| 383 |
+ |
|
| 384 |
+ |
|
| 385 |
+ |
\subsection{\label{lipidSec:scd}$\text{S}_{\text{{\sc cd}}}$ Order Parameters} |
| 386 |
+ |
|
| 387 |
|
The $\text{S}_{\text{{\sc cd}}}$ order parameter is often reported in |
| 388 |
< |
the experimental charecterizations of phospholipids. It is obtained |
| 388 |
> |
the experimental characterizations of phospholipids. It is obtained |
| 389 |
|
through deuterium NMR, and measures the ordering of the carbon |
| 390 |
|
deuterium bond in relation to the bilayer normal at various points |
| 391 |
|
along the chains. In our model, there are no explicit hydrogens, but |
| 397 |
|
\end{equation} |
| 398 |
|
Where $S_{ij}$ is given by: |
| 399 |
|
\begin{equation} |
| 400 |
< |
S_{ij} = \frac{1}{2}\Bigl<(3\cos\Theta_i\cos\Theta_j - \delta_{ij})\Bigr> |
| 400 |
> |
S_{ij} = \frac{1}{2}\Bigl\langle(3\cos\Theta_i\cos\Theta_j |
| 401 |
> |
- \delta_{ij})\Bigr\rangle |
| 402 |
|
\label{lipidEq:scd2} |
| 403 |
|
\end{equation} |
| 404 |
< |
Here, $\Theta_i$ is the angle the $i$th carbon atom frame axis makes |
| 405 |
< |
with the bilayer normal. The brackets denote an average over time and |
| 406 |
< |
molecules. The carbon atom axes are defined: |
| 407 |
< |
$\mathbf{\hat{z}}\rightarrow$ vector from $C_{n-1}$ to $C_{n+1}$; |
| 408 |
< |
$\mathbf{\hat{y}}\rightarrow$ vector that is perpindicular to $z$ and |
| 409 |
< |
in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$; |
| 410 |
< |
$\mathbf{\hat{x}}\rightarrow$ vector perpindicular to |
| 404 |
> |
Here, $\Theta_i$ is the angle the $i$th axis in the reference frame of |
| 405 |
> |
the carbon atom makes with the bilayer normal. The brackets denote an |
| 406 |
> |
average over time and molecules. The carbon atom axes are defined: |
| 407 |
> |
\begin{itemize} |
| 408 |
> |
\item $\mathbf{\hat{z}}\rightarrow$ vector from $C_{n-1}$ to $C_{n+1}$ |
| 409 |
> |
\item $\mathbf{\hat{y}}\rightarrow$ vector that is perpendicular to $z$ and |
| 410 |
> |
in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$ |
| 411 |
> |
\item $\mathbf{\hat{x}}\rightarrow$ vector perpendicular to |
| 412 |
|
$\mathbf{\hat{y}}$ and $\mathbf{\hat{z}}$. |
| 413 |
+ |
\end{itemize} |
| 414 |
+ |
This assumes that the hydrogen atoms are always in a plane |
| 415 |
+ |
perpendicular to the $C_{n-1}-C_{n}-C_{n+1}$ plane. |
| 416 |
|
|
| 417 |
|
The order parameter has a range of $[1,-\frac{1}{2}]$. A value of 1 |
| 418 |
|
implies full order aligned to the bilayer axis, 0 implies full |
| 419 |
< |
disorder, and $-\frac{1}{2}$ implies full order perpindicular to the |
| 419 |
> |
disorder, and $-\frac{1}{2}$ implies full order perpendicular to the |
| 420 |
|
bilayer axis. The {\sc cd} bond vector for carbons near the head group |
| 421 |
< |
are usually ordered perpindicular to the bilayer normal, with tails |
| 422 |
< |
farther away tending toward disorder. This makes the order paramter |
| 421 |
> |
are usually ordered perpendicular to the bilayer normal, with tails |
| 422 |
> |
farther away tending toward disorder. This makes the order parameter |
| 423 |
|
negative for most carbons, and as such $|S_{\text{{\sc cd}}}|$ is more |
| 424 |
|
commonly reported than $S_{\text{{\sc cd}}}$. |
| 425 |
|
|
| 426 |
+ |
Fig.~\ref{lipidFig:scdFig} shows the $S_{\text{{\sc cd}}}$ order |
| 427 |
+ |
parameters for the bilayer system at 300~K. There is no appreciable |
| 428 |
+ |
difference in the plots for the various temperatures, however, there |
| 429 |
+ |
is a larger difference between our model's ordering, and the |
| 430 |
+ |
experimentally observed ordering of DMPC. As our values are closer to |
| 431 |
+ |
$-\frac{1}{2}$, this implies more ordering perpendicular to the normal |
| 432 |
+ |
than in a real system. This is due to the model having only one carbon |
| 433 |
+ |
group separating the chains from the top of the lipid. In DMPC, with |
| 434 |
+ |
the flexibility inherent in a multiple atom head group, as well as a |
| 435 |
+ |
glycerol linkage between the head group and the acyl chains, there is |
| 436 |
+ |
more loss of ordering by the point when the chains start. |
| 437 |
|
|
| 351 |
– |
|
| 352 |
– |
|
| 438 |
|
\begin{figure} |
| 439 |
|
\centering |
| 440 |
|
\includegraphics[width=\linewidth]{scdFig.eps} |
| 442 |
|
\label{lipidFig:scdFig} |
| 443 |
|
\end{figure} |
| 444 |
|
|
| 445 |
+ |
\subsection{\label{lipidSec:p2Order}$P_2$ Order Parameter} |
| 446 |
|
|
| 447 |
< |
\begin{figure} |
| 448 |
< |
\centering |
| 449 |
< |
\includegraphics[width=\linewidth]{densityProfile.eps} |
| 450 |
< |
\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.} |
| 451 |
< |
\label{lipidFig:densityProfile} |
| 452 |
< |
\end{figure} |
| 447 |
> |
The $P_2$ order parameter allows us to measure the amount of |
| 448 |
> |
directional ordering that exists in the bodies of the molecules making |
| 449 |
> |
up the bilayer. Each lipid molecule can be thought of as a cylindrical |
| 450 |
> |
rod with the head group at the top. If all of the rods are perfectly |
| 451 |
> |
aligned, the $P_2$ order parameter will be $1.0$. If the rods are |
| 452 |
> |
completely disordered, the $P_2$ order parameter will be 0. For a |
| 453 |
> |
collection of unit vectors pointing along the principal axes of the |
| 454 |
> |
rods, the $P_2$ order parameter can be solved via the following |
| 455 |
> |
method.\cite{zannoni94} |
| 456 |
|
|
| 457 |
+ |
Define an ordering tensor $\overleftrightarrow{\mathsf{Q}}$, such that, |
| 458 |
+ |
\begin{equation} |
| 459 |
+ |
\overleftrightarrow{\mathsf{Q}} = \frac{1}{N}\sum_i^N % |
| 460 |
+ |
\begin{pmatrix} % |
| 461 |
+ |
u_{ix}u_{ix}-\frac{1}{3} & u_{ix}u_{iy} & u_{ix}u_{iz} \\ |
| 462 |
+ |
u_{iy}u_{ix} & u_{iy}u_{iy}-\frac{1}{3} & u_{iy}u_{iz} \\ |
| 463 |
+ |
u_{iz}u_{ix} & u_{iz}u_{iy} & u_{iz}u_{iz}-\frac{1}{3} % |
| 464 |
+ |
\end{pmatrix} |
| 465 |
+ |
\label{lipidEq:po1} |
| 466 |
+ |
\end{equation} |
| 467 |
+ |
Where the $u_{i\alpha}$ is the $\alpha$ element of the unit vector |
| 468 |
+ |
$\mathbf{\hat{u}}_i$, and the sum over $i$ averages over the whole |
| 469 |
+ |
collection of unit vectors. This allows the tensor to be written: |
| 470 |
+ |
\begin{equation} |
| 471 |
+ |
\overleftrightarrow{\mathsf{Q}} = \frac{1}{N}\sum_i^N \biggl[ |
| 472 |
+ |
\mathbf{\hat{u}}_i \otimes \mathbf{\hat{u}}_i |
| 473 |
+ |
- \frac{1}{3} \cdot \mathsf{1} \biggr] |
| 474 |
+ |
\label{lipidEq:po2} |
| 475 |
+ |
\end{equation} |
| 476 |
|
|
| 477 |
+ |
After constructing the tensor, diagonalizing |
| 478 |
+ |
$\overleftrightarrow{\mathsf{Q}}$ yields three eigenvalues and |
| 479 |
+ |
eigenvectors. The eigenvector associated with the largest eigenvalue, |
| 480 |
+ |
$\lambda_{\text{max}}$, is the director axis for the system of unit |
| 481 |
+ |
vectors. The director axis is the average direction all of the unit vectors |
| 482 |
+ |
are pointing. The $P_2$ order parameter is then simply |
| 483 |
+ |
\begin{equation} |
| 484 |
+ |
\langle P_2 \rangle = \frac{3}{2}\lambda_{\text{max}} |
| 485 |
+ |
\label{lipidEq:po3} |
| 486 |
+ |
\end{equation} |
| 487 |
|
|
| 488 |
< |
\begin{figure} |
| 489 |
< |
\centering |
| 490 |
< |
\includegraphics[width=\linewidth]{diffusionFig.eps} |
| 491 |
< |
\caption[The lateral difusion constants versus temperature]{The lateral diffusion constants for the bilayers as a function of temperature.} |
| 492 |
< |
\label{lipidFig:diffusionFig} |
| 493 |
< |
\end{figure} |
| 488 |
> |
Table~\ref{lipidTab:blSummary} summarizes the $P_2$ values for the |
| 489 |
> |
bilayers, as well as the dipole orientations. The unit vector for the |
| 490 |
> |
lipid molecules was defined by finding the moment of inertia for each |
| 491 |
> |
lipid, then setting $\mathbf{\hat{u}}$ to point along the axis of |
| 492 |
> |
minimum inertia. For the {\sc head} atoms, the unit vector simply |
| 493 |
> |
pointed in the same direction as the dipole moment. For the lipid |
| 494 |
> |
molecules, the ordering was consistent across all temperatures, with |
| 495 |
> |
the director pointed along the $z$ axis of the box. More |
| 496 |
> |
interestingly, is the high degree of ordering the dipoles impose on |
| 497 |
> |
the {\sc head} atoms. The directors for the dipoles themselves |
| 498 |
> |
consistently pointed along the plane of the bilayer, with the |
| 499 |
> |
directors anti-aligned on the top and bottom leaf. |
| 500 |
|
|
| 501 |
|
\begin{table} |
| 502 |
< |
\caption[Structural properties of the bilayers]{Bilayer Structural properties as a function of temperature.} |
| 502 |
> |
\caption[Structural properties of the bilayers]{BILAYER STRUCTURAL PROPERTIES AS A FUNCTION OF TEMPERATURE} |
| 503 |
> |
\label{lipidTab:blSummary} |
| 504 |
|
\begin{center} |
| 505 |
|
\begin{tabular}{|c|c|c|c|c|} |
| 506 |
|
\hline |
| 507 |
< |
Temperature (K) & $<L_{\perp}>$ ($\mbox{\AA}$) & % |
| 508 |
< |
$<A_{\parallel}>$ ($\mbox{\AA}^2$) & $<P_2>_{\text{Lipid}}$ & % |
| 509 |
< |
$<P_2>_{\text{{\sc head}}}$ \\ \hline |
| 507 |
> |
Temperature (K) & $\langle L_{\perp}\rangle$ ($\mbox{\AA}$) & % |
| 508 |
> |
$\langle A_{\parallel}\rangle$ ($\mbox{\AA}^2$) & % |
| 509 |
> |
$\langle P_2\rangle_{\text{Lipid}}$ & % |
| 510 |
> |
$\langle P_2\rangle_{\text{{\sc head}}}$ \\ \hline |
| 511 |
|
270 & 18.1 & 58.1 & 0.253 & 0.494 \\ \hline |
| 512 |
|
275 & 17.2 & 56.7 & 0.295 & 0.514 \\ \hline |
| 513 |
|
277 & 16.9 & 58.0 & 0.301 & 0.541 \\ \hline |
| 519 |
|
\end{tabular} |
| 520 |
|
\end{center} |
| 521 |
|
\end{table} |
| 522 |
+ |
|
| 523 |
+ |
\subsection{\label{lipidSec:miscData}Further Structural Data} |
| 524 |
+ |
|
| 525 |
+ |
Also summarized in Table~\ref{lipidTab:blSummary}, are the bilayer |
| 526 |
+ |
thickness ($\langle L_{\perp}\rangle$) and area per lipid ($\langle |
| 527 |
+ |
A_{\parallel}\rangle$). The bilayer thickness was measured from the |
| 528 |
+ |
peak to peak {\sc head} atom distance in the density profiles. The |
| 529 |
+ |
area per lipid data compares favorably with values typically seen for |
| 530 |
+ |
DMPC (60.0~$\mbox{\AA}^2$ at 303~K)\cite{petrache00}. Although our |
| 531 |
+ |
values are lower this is most likely due to the shorter chain length |
| 532 |
+ |
of our model (8 versus 14 for DMPC). |
| 533 |
+ |
|
| 534 |
+ |
\subsection{\label{lipidSec:diffusion}Lateral Diffusion Constants} |
| 535 |
+ |
|
| 536 |
+ |
The lateral diffusion constant, $D_L$, is the constant characterizing |
| 537 |
+ |
the diffusive motion of the lipid molecules within the plane of the bilayer. It |
| 538 |
+ |
is given by the following Einstein relation:\cite{allen87:csl} |
| 539 |
+ |
\begin{equation} |
| 540 |
+ |
D_L = \lim_{t\rightarrow\infty}\frac{1}{4t}\langle |\mathbf{r}(t) |
| 541 |
+ |
- \mathbf{r}(0)|^2\rangle |
| 542 |
+ |
\end{equation} |
| 543 |
+ |
Where $\mathbf{r}(t)$ is the $xy$ position of the lipid at time $t$ |
| 544 |
+ |
(assuming the $z$-axis is parallel to the bilayer normal). |
| 545 |
+ |
|
| 546 |
+ |
Fig.~\ref{lipidFig:diffusionFig} shows the lateral diffusion constants |
| 547 |
+ |
as a function of temperature. There is a definite increase in the |
| 548 |
+ |
lateral diffusion with higher temperatures, which is exactly what one |
| 549 |
+ |
would expect with greater fluidity of the chains. However, the |
| 550 |
+ |
diffusion constants are two orders of magnitude smaller than those |
| 551 |
+ |
typical of DPPC.\cite{Cevc87} This is counter-intuitive as the DPPC |
| 552 |
+ |
molecule is sterically larger and heavier than our model. This could |
| 553 |
+ |
be an indication that our model's chains are too interwoven and hinder |
| 554 |
+ |
the motion of the lipid or that the dipolar head groups are too |
| 555 |
+ |
tightly bound to each other. In contrast, the diffusion constant of |
| 556 |
+ |
the {\sc ssd} water, $9.84\times 10^{-6}\,\text{cm}^2/\text{s}$, is |
| 557 |
+ |
reasonably close to the bulk water diffusion constant ($2.2999\times |
| 558 |
+ |
10^{-5}\,\text{cm}^2/\text{s}$).\cite{Holz00} |
| 559 |
+ |
|
| 560 |
+ |
\begin{figure} |
| 561 |
+ |
\centering |
| 562 |
+ |
\includegraphics[width=\linewidth]{diffusionFig.eps} |
| 563 |
+ |
\caption[The lateral diffusion constants versus temperature]{The lateral diffusion constants for the bilayers as a function of temperature.} |
| 564 |
+ |
\label{lipidFig:diffusionFig} |
| 565 |
+ |
\end{figure} |
| 566 |
+ |
|
| 567 |
+ |
\subsection{\label{lipidSec:randBilayer}Bilayer Aggregation} |
| 568 |
+ |
|
| 569 |
+ |
A very important accomplishment for our model is its ability to |
| 570 |
+ |
spontaneously form bilayers from a randomly dispersed starting |
| 571 |
+ |
configuration. Fig.~\ref{lipidFig:blImage} shows an image sequence for |
| 572 |
+ |
the bilayer aggregation. After 3.0~ns, the basic form of the bilayer |
| 573 |
+ |
can already be seen. By 7.0~ns, the bilayer has a lipid bridge |
| 574 |
+ |
stretched across the simulation box to itself that will turn out to be |
| 575 |
+ |
very long lived ($\sim$20~ns), as well as a water pore, that will |
| 576 |
+ |
persist for the length of the current simulation. At 24~ns, the lipid |
| 577 |
+ |
bridge has broken, and the bilayer is still integrating the lipid |
| 578 |
+ |
molecules from the bridge into itself. However, the water pore is |
| 579 |
+ |
still present at 24~ns. |
| 580 |
+ |
|
| 581 |
+ |
\begin{figure} |
| 582 |
+ |
\centering |
| 583 |
+ |
\includegraphics[width=\linewidth]{bLayerImage.eps} |
| 584 |
+ |
\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.} |
| 585 |
+ |
\label{lipidFig:blImage} |
| 586 |
+ |
\end{figure} |
| 587 |
+ |
|
| 588 |
+ |
\subsection{\label{lipidSec:randIrod}Inverted Rod Aggregation} |
| 589 |
+ |
|
| 590 |
+ |
Fig.~\ref{lipidFig:iRimage} shows a second aggregation sequence |
| 591 |
+ |
simulated in this research. Here the fraction of water had been |
| 592 |
+ |
significantly decreased to observe how the model would respond. After |
| 593 |
+ |
1.5~ns, The main body of water in the system has already collected |
| 594 |
+ |
into a central water channel. By 10.0~ns, the channel has widened |
| 595 |
+ |
slightly, but there are still many water molecules permeating the |
| 596 |
+ |
lipid macro-structure. At 23.0~ns, the central water channel has |
| 597 |
+ |
stabilized and several smaller water channels have been absorbed by |
| 598 |
+ |
the main one. However, there is still an appreciable water |
| 599 |
+ |
concentration throughout the lipid structure. |
| 600 |
+ |
|
| 601 |
+ |
\begin{figure} |
| 602 |
+ |
\centering |
| 603 |
+ |
\includegraphics[width=\linewidth]{iRodImage.eps} |
| 604 |
+ |
\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}.} |
| 605 |
+ |
\label{lipidFig:iRimage} |
| 606 |
+ |
\end{figure} |
| 607 |
+ |
|
| 608 |
+ |
\section{\label{lipidSec:Conclusion}Conclusion} |
| 609 |
+ |
|
| 610 |
+ |
We have presented a simple unified-atom phospholipid model capable of |
| 611 |
+ |
spontaneous aggregation into a bilayer and an inverted rod |
| 612 |
+ |
structure. The time scales of the macro-molecular aggregations are |
| 613 |
+ |
approximately 24~ns. In addition the model's properties have been |
| 614 |
+ |
explored over a range of temperatures through prefabricated |
| 615 |
+ |
bilayers. No freezing transition is seen in the temperature range of |
| 616 |
+ |
our current simulations. However, structural information from 270~K |
| 617 |
+ |
may imply that a freezing event is on a much longer time scale than |
| 618 |
+ |
that explored in this current research. Further studies of this system |
| 619 |
+ |
could extend the time length of the simulations at the low |
| 620 |
+ |
temperatures to observe whether lipid crystallization can occur within |
| 621 |
+ |
the framework of this model. |
| 622 |
+ |
|
| 623 |
+ |
Potential problems that may be obstacles in further research, is the |
| 624 |
+ |
lack of detail in the head region. As the chains are almost directly |
| 625 |
+ |
attached to the {\sc head} atom, there is no buffer between the |
| 626 |
+ |
actions of the head group and the tails. Another disadvantage of the |
| 627 |
+ |
model is the dipole approximation will alter results when details |
| 628 |
+ |
concerning a charged solute's interactions with the bilayer. However, |
| 629 |
+ |
it is important to keep in mind that the dipole approximation can be |
| 630 |
+ |
kept an advantage by examining solutes that do not require point |
| 631 |
+ |
charges, or at the least, require only dipole approximations |
| 632 |
+ |
themselves. Other advantages of the model include the ability to alter |
| 633 |
+ |
the size of the unified-atoms so that the size of the lipid can be |
| 634 |
+ |
increased without adding to the number of interactions in the |
| 635 |
+ |
system. However, what sets our model apart from other current |
| 636 |
+ |
simplified models,\cite{goetz98,marrink04} is the information gained |
| 637 |
+ |
by observing the ordering of the head groups dipole's in relation to |
| 638 |
+ |
each other and the solvent without the need for point charges and the |
| 639 |
+ |
Ewald sum. |
