trunk/mattDisertation/lipid.tex



\chapter{\label{chapt:lipid}Phospholipid Simulations}

\section{\label{lipidSec:Intro}Introduction}

In the past 10 years, increasing computer speeds have allowed for the
atomistic simulation of phospholipid bilayers for increasingly
relevant lengths of time.  These simulations have ranged from
simulation of the gel phase ($L_{\beta}$) of
dipalmitoylphosphatidylcholine (DPPC),\cite{lindahl00} to the
spontaneous aggregation of DPPC molecules into fluid phase
($L_{\alpha}$) bilayers.\cite{marrink01} With the exception of a few
ambitious
simulations,\cite{marrink01:undulation,marrink:2002,lindahl00} most
investigations are limited to a range of 64 to 256
phospholipids.\cite{lindahl00,sum:2003,venable00,gomez:2003,smondyrev:1999,marrink01}
The expense of the force calculations involved when performing these
simulations limits the system size. To properly hydrate a bilayer, one
typically needs 25 water molecules for every lipid, bringing the total
number of atoms simulated to roughly 8,000 for a system of 64 DPPC
molecules. Added to the difficulty is the electrostatic nature of the
phospholipid head groups and water, requiring either the
computationally expensive Ewald sum or the faster, particle mesh Ewald
sum.\cite{nina:2002,norberg:2000,patra:2003} These factors all limit
the system size and time scales of bilayer simulations.

Unfortunately, much of biological interest happens on time and length
scales well beyond the range of current simulation technology. One
such example is the observance of a ripple phase
($P_{\beta^{\prime}}$) between the $L_{\beta}$ and $L_{\alpha}$ phases
of certain phospholipid bilayers.\cite{katsaras00,sengupta00} These
ripples are shown to have periodicity on the order of
100-200~$\mbox{\AA}$. A simulation on this length scale would have
approximately 1,300 lipid molecules with an additional 25 water
molecules per lipid to fully solvate the bilayer. A simulation of this
size is impractical with current atomistic models.

The time and length scale limitations are most striking in transport
phenomena.  Due to the fluid-like properties of a lipid membrane, not
all diffusion across the membrane happens at pores.  Some molecules of
interest may incorporate themselves directly into the membrane.  Once
here, they may possess an appreciable waiting time (on the order of
10's to 100's of nanoseconds) within the bilayer. Such long simulation
times are nearly impossible to obtain when integrating the system with
atomistic detail.

To address these issues, several schemes have been proposed.  One
approach by Goetz and Liposky\cite{goetz98} is to model the entire
system as Lennard-Jones spheres. Phospholipids are represented by
chains of beads with the top most beads identified as the head
atoms. Polar and non-polar interactions are mimicked through
attractive and soft-repulsive potentials respectively.  A similar
model proposed by Marrinck \emph{et. al}.\cite{marrink04}~uses a
similar technique for modeling polar and non-polar interactions with
Lennard-Jones spheres. However, they also include charges on the head
group spheres to mimic the electrostatic interactions of the
bilayer. While the solvent spheres are kept charge-neutral and
interact with the bilayer solely through an attractive Lennard-Jones
potential.

The model used in this investigation adds more information to the
interactions than the previous two models, while still balancing the
need for simplifications over atomistic detail.  The model uses
Lennard-Jones spheres for the head and tail groups of the
phospholipids, allowing for the ability to scale the parameters to
reflect various sized chain configurations while keeping the number of
interactions small.  What sets this model apart, however, is the use
of dipoles to represent the electrostatic nature of the
phospholipids. The dipole electrostatic interaction is shorter range
than Coulombic ($\frac{1}{r^3}$ versus $\frac{1}{r}$), and eliminates
the need for a costly Ewald sum.

Another key feature of this model, is the use of a dipolar water model
to represent the solvent. The soft sticky dipole ({\sc ssd}) water
\cite{liu96:new_model} relies on the dipole for long range electrostatic
effects, but also contains a short range correction for hydrogen
bonding. In this way the systems in this research mimic the entropic
contribution to the hydrophobic effect due to hydrogen-bond network
deformation around a non-polar entity, \emph{i.e.}~the phospholipid
molecules.

The following is an outline of this chapter.
Sec.~\ref{lipidSec:Methods} is an introduction to the lipid model used
in these simulations, as well as clarification about the water model
and integration techniques. The various simulations explored in this
research are outlined in
Sec.~\ref{lipidSec:ExpSetup}. Sec.~\ref{lipidSec:resultsDis} gives a
summary and interpretation of the results.  Finally, the conclusions
of this chapter are presented in Sec.~\ref{lipidSec:Conclusion}.

\section{\label{lipidSec:Methods}Methods}

\subsection{\label{lipidSec:lipidModel}The Lipid Model}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{twoChainFig.eps}
\caption[The two chained lipid model]{Schematic diagram of the double chain phospholipid model. The head group (in red) has a point dipole, $\boldsymbol{\mu}$, located at its center of mass. The two chains are eight methylene groups in length.}
\label{lipidFig:lipidModel}
\end{figure}

The phospholipid model used in these simulations is based on the
design illustrated in Fig.~\ref{lipidFig:lipidModel}.  The head group
of the phospholipid is replaced by a single Lennard-Jones sphere of
diameter $\sigma_{\text{head}}$, with $\epsilon_{\text{head}}$ scaling
the well depth of its van der Walls interaction.  This sphere also
contains a single dipole of magnitude $|\boldsymbol{\mu}|$, where
$|\boldsymbol{\mu}|$ can be varied to mimic the charge separation of a
given phospholipid head group.  The atoms of the tail region are
modeled by beads representing multiple methyl groups.  They are free
of partial charges or dipoles, and contain only Lennard-Jones
interaction sites at their centers of mass.  As with the head groups,
their potentials can be scaled by $\sigma_{\text{tail}}$ and
$\epsilon_{\text{tail}}$.

The long range interactions between lipids are given by the following:
\begin{equation}
V_{\text{LJ}}(r_{ij}) = 
        4\epsilon_{ij} \biggl[
        \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
        - \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
        \biggr]
\label{lipidEq:LJpot}
\end{equation}
and
\begin{equation}
V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
        \boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[
        \boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j}
        -
        3(\boldsymbol{\hat{u}}_i \cdot \mathbf{\hat{r}}_{ij}) %
                (\boldsymbol{\hat{u}}_j \cdot \mathbf{\hat{r}}_{ij} \biggr]
\label{lipidEq:dipolePot}
\end{equation}
Where $V_{\text{LJ}}$ is the Lennard-Jones potential and
$V_{\text{dipole}}$ is the dipole-dipole potential.  As previously
stated, $\sigma_{ij}$ and $\epsilon_{ij}$ are the Lennard-Jones
parameters which scale the length and depth of the interaction
respectively, and $r_{ij}$ is the distance between beads $i$ and $j$.
In $V_{\text{dipole}}$, $\mathbf{r}_{ij}$ is the vector starting at
bead $i$ and pointing towards bead $j$.  Vectors $\mathbf{\Omega}_i$
and $\mathbf{\Omega}_j$ are the orientational degrees of freedom for
beads $i$ and $j$.  $|\mu_i|$ is the magnitude of the dipole moment of
$i$, and $\boldsymbol{\hat{u}}_i$ is the standard unit orientation
vector rotated with Euler angles: $\boldsymbol{\Omega}_i$.

The model also allows for the bonded interactions bends, and torsions.
The bond between two beads on a chain is of fixed length, and is
maintained according to the {\sc rattle} algorithm.\cite{andersen83}
The bends are subject to a harmonic potential:
\begin{equation}
V_{\text{bend}}(\theta) = k_{\theta}( \theta - \theta_0 )^2
\label{lipidEq:bendPot}
\end{equation}
where $k_{\theta}$ scales the strength of the harmonic well, and
$\theta$ is the angle between bond vectors
(Fig.~\ref{lipidFig:lipidModel}). In addition, we have placed a
``ghost'' bend on the phospholipid head. The ghost bend adds a
potential to keep the dipole pointed along the bilayer surface, where
$\theta$ is now the angle the dipole makes with respect to the {\sc
head}-$\text{{\sc ch}}_2$ bond vector
(Fig.~\ref{lipidFig:ghostBend}). The torsion potential is given by:
\begin{equation}
V_{\text{torsion}}(\phi) =  
        k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
\label{lipidEq:torsionPot}
\end{equation}
Here, the parameters $k_0$, $k_1$, $k_2$, and $k_3$ fit the cosine
power series to the desired torsion potential surface, and $\phi$ is
the angle the two end atoms have rotated about the middle bond
(Fig.:\ref{lipidFig:lipidModel}).  Long range interactions such as the
Lennard-Jones potential are excluded for atom pairs involved in the
same bond, bend, or torsion.  However, internal interactions not
directly involved in a bonded pair are calculated.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{ghostBendFig.eps}
\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.}
\label{lipidFig:ghostBend}
\end{figure}

All simulations presented here use a two chained lipid as pictured in
Fig.~\ref{lipidFig:lipidModel}.  The chains are both eight beads long,
and their mass and Lennard Jones parameters are summarized in
Table~\ref{lipidTable:tcLJParams}. The magnitude of the dipole moment
for the head bead is 10.6~Debye, and the bend and torsion parameters
are summarized in Table~\ref{lipidTable:tcBendParams} and
\ref{lipidTable:tcTorsionParams}.

\begin{table}
\caption{The Lennard Jones Parameters for the two chain phospholipids.}
\label{lipidTable:tcLJParams}
\begin{center}
\begin{tabular}{|l|c|c|c|}
\hline
     & mass (amu) & $\sigma$($\mbox{\AA}$)  & $\epsilon$ (kcal/mol) \\ \hline
{\sc head} & 72  & 4.0 & 0.185 \\ \hline
{\sc ch}\cite{Siepmann1998} & 13.02 & 4.0 & 0.0189 \\ \hline
$\text{{\sc ch}}_2$\cite{Siepmann1998} & 14.03 & 3.95 & 0.18 \\ \hline
$\text{{\sc ch}}_3$\cite{Siepmann1998} & 15.04 & 3.75 & 0.25 \\ \hline
{\sc ssd}\cite{liu96:new_model} & 18.03 & 3.051 & 0.152 \\ \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}
\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}}
\label{lipidTable:tcBendParams}
\begin{center}
\begin{tabular}{|l|c|c|}
\hline
   & $k_{\theta}$ ( kcal/($\text{mol deg}^2$) ) & $\theta_0$ ( deg ) \\ \hline
{\sc ghost}-{\sc head}-$\text{{\sc ch}}_2$ & 0.00177 & 129.78 \\ \hline
$x$-{\sc ch}-$y$ & 58.84 & 112.0 \\ \hline
$x$-$\text{{\sc ch}}_2$-$y$ & 58.84 & 114.0 \\ \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}
\caption[Torsion Parameters for the two chain phospholipids]{Torsion Parameters for the two chain phospholipids. Alkane parameters based on TraPPE.\cite{Siepmann1998}}
\label{lipidTable:tcTorsionParams}
\begin{center}
\begin{tabular}{|l|c|c|c|c|}
\hline
All are in kcal/mol $\rightarrow$ & $k_3$ & $k_2$ & $k_1$ & $k_0$ \\ \hline
$x$-{\sc ch}-$y$-$z$ & 3.3254 & -0.4215 & -1.686 & 1.1661 \\ \hline
$x$-$\text{{\sc ch}}_2$-$\text{{\sc ch}}_2$-$y$ & 5.9602 & -0.568 & -3.802 & 2.1586 \\ \hline
\end{tabular}
\end{center}
\end{table}


\section{\label{lipidSec:furtherMethod}Further Methodology}

As mentioned previously, the water model used throughout these
simulations was the {\sc ssd} model of
Ichiye.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md} A
discussion of the model can be found in Sec.~\ref{oopseSec:SSD}. As
for the integration of the equations of motion, all simulations were
performed in an orthorhombic periodic box with a thermostat on
velocities, and an independent barostat on each Cartesian axis $x$,
$y$, and $z$.  This is the $\text{NPT}_{xyz}$. ensemble described in
Sec.~\ref{oopseSec:integrate}.


\subsection{\label{lipidSec:ExpSetup}Experimental Setup}

Two main starting configuration classes were used in this research:
random and ordered bilayers.  The ordered bilayer starting
configurations were all started from an equilibrated bilayer at
300~K. The original configuration for the first 300~K run was
assembled by placing the phospholipids centers of mass on a planar
hexagonal lattice.  The lipids were oriented with their long axis
perpendicular to the plane.  The second leaf simply mirrored the first
leaf, and the appropriate number of waters were then added above and
below the bilayer.

The random configurations took more work to generate.  To begin, a
test lipid was placed in a simulation box already containing water at
the intended density.  The waters were then tested for overlap with
the lipid using a 5.0~$\mbox{\AA}$ buffer distance.  This gave an
estimate for the number of waters each lipid would displace in a
simulation box. A target number of waters was then defined which
included the number of waters each lipid would displace, the number of
waters desired to solvate each lipid, and a factor to pad the
initial box with a few extra water molecules. 

Next, a cubic simulation box was created that contained at least the
target number of waters in an FCC lattice (the lattice was for ease of
placement).  What followed was a RSA simulation similar to those of
Chapt.~\ref{chapt:RSA}. The lipids were sequentially given a random
position and orientation within the box.  If a lipid's position caused
atomic overlap with any previously placed lipid, its position and
orientation were rejected, and a new random placement site was
attempted. The RSA simulation proceeded until all phospholipids had
been adsorbed.  After placement of all lipid molecules, water
molecules with locations that overlapped with the atomic coordinates
of the lipids were removed.

Finally, water molecules were removed at random until the desired
number of waters per lipid was reached.  The typical low final density
for these initial configurations was not a problem, as the box shrinks
to an appropriate size within the first 50~ps of a simulation in the
$\text{NPT}_{xyz}$ ensemble.

\subsection{\label{lipidSec:Configs}Configurations}

The first class of simulations were started from ordered
bilayers. They were all configurations consisting of 60 lipid
molecules with 30 lipids on each leaf, and were hydrated with 1620
{\sc ssd} molecules. The original configuration was assembled
according to Sec.~\ref{lipidSec:ExpSetup} and simulated for a length
of 10~ns at 300~K. The other temperature runs were started from a
frame 7~ns into the 300~K simulation. Their temperatures were reset
with the thermostating algorithm in the $\text{NPT}_{xyz}$
integrator. All of the temperature variants were also run for 10~ns,
with only the last 5~ns being used for accumulation of statistics.

The second class of simulations were two configurations started from
randomly dispersed lipids in a ``gas'' of water. The first
($\text{R}_{\text{I}}$) was a simulation containing 72 lipids with
1800 {\sc ssd} molecules simulated at 300~K. The second
($\text{R}_{\text{II}}$) was 90 lipids with 1350 {\sc ssd} molecules
simulated at 350~K. Both simulations were integrated for more than
20~ns, and illustrate the spontaneous aggregation of the lipid model
into phospholipid macro-structures: $\text{R}_{\text{I}}$ into a
bilayer, and $\text{R}_{\text{II}}$ into a inverted rod.

\section{\label{lipidSec:resultsDis}Results and Discussion}

\subsection{\label{lipidSec:diffusion}Lateral Diffusion Constants}

The lateral diffusion constant, $D_L$, is the constant characterizing
the diffusive motion of the lipid within the plane of the bilayer. It
is given by the following Einstein relation valid at long
times:\cite{allen87:csl}
\begin{equation}
2tD_L = \frac{1}{2}\langle |\mathbf{r}(t) - \mathbf{r}(0)|^2\rangle
\end{equation}
Where $\mathbf{r}(t)$ is the position of the lipid at time $t$, and is
constrained to lie within a plane. For the bilayer simulations the
plane of constrained motion was that perpendicular to the bilayer
normal, namely the $xy$-plane.

Fig.~\ref{lipidFig:diffusionFig} shows the lateral diffusion constants
as a function of temperature. There is a definite increase in the
lateral diffusion with higher temperatures, which is exactly what one
would expect with greater fluidity of the chains. However, the
diffusion constants are all two orders of magnitude smaller than those
typical of DPPC.\cite{Cevc87} This is counter-intuitive as the DPPC
molecule is sterically larger and heavier than our model. This could
be an indication that our model's chains are too interwoven and hinder
the motion of the lipid, or that a simplification in the model's
forces has led to a slowing of diffusive behavior within the
bilayer. In contrast, the diffusion constant of the {\sc ssd} water,
$9.84\times 10^{-6}\,\text{cm}^2/\text{s}$, compares favorably with
that of bulk water ($2.2999\times
10^{-5}\,\text{cm}^2/\text{s}$\cite{Holz00}).

\begin{figure}
\centering
\includegraphics[width=\linewidth]{diffusionFig.eps}
\caption[The lateral diffusion constants versus temperature]{The lateral diffusion constants for the bilayers as a function of temperature.}
\label{lipidFig:diffusionFig}
\end{figure}

\subsection{\label{lipidSec:densProf}Density Profile}

Fig.~\ref{lipidFig:densityProfile} illustrates the densities of the
atoms in the bilayer systems normalized by the bulk density as a
function of distance from the center of the box. The profile is taken
along the bilayer normal, in this case the $z$ axis. The profile at
270~K shows several structural features that are largely smoothed out
by 300~K. The left peak for the {\sc head} atoms is split at 270~K,
implying that some freezing of the structure might already be occurring
at this temperature. From the dynamics, the tails at this temperature
are very much fluid, but the profile could indicate that a phase
transition may simply be beyond the length scale of the current
simulation. In all profiles, the water penetrates almost
5~$\mbox{\AA}$ into the bilayer, completely solvating the {\sc head}
atoms. The $\text{{\sc ch}}_3$ atoms although mainly centered at the
middle of the bilayer, show appreciable penetration into the head
group region. This indicates that the chains have enough mobility to
bend back upward to allow the ends to explore areas around the {\sc
head} atoms. It is unlikely that this is penetration from a lipid of
the opposite face, as the lipids are only 12~$\mbox{\AA}$ in length,
and the typical leaf spacing as measured from the {\sc head-head}
spacing in the profile is 17.5~$\mbox{\AA}$.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{densityProfile.eps}
\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.}
\label{lipidFig:densityProfile}
\end{figure}


\subsection{\label{lipidSec:scd}$\text{S}_{\text{{\sc cd}}}$ Order Parameters}

The $\text{S}_{\text{{\sc cd}}}$ order parameter is often reported in
the experimental characterizations of phospholipids. It is obtained
through deuterium NMR, and measures the ordering of the carbon
deuterium bond in relation to the bilayer normal at various points
along the chains. In our model, there are no explicit hydrogens, but
the order parameter can be written in terms of the carbon ordering at
each point in the chain:\cite{egberts88}
\begin{equation}
S_{\text{{\sc cd}}}  = \frac{2}{3}S_{xx} + \frac{1}{3}S_{yy}
\label{lipidEq:scd1}
\end{equation}
Where $S_{ij}$ is given by:
\begin{equation}
S_{ij} = \frac{1}{2}\Bigl\langle(3\cos\Theta_i\cos\Theta_j 
        - \delta_{ij})\Bigr\rangle
\label{lipidEq:scd2}
\end{equation}
Here, $\Theta_i$ is the angle the $i$th axis in the reference frame of
the carbon atom makes with the bilayer normal. The brackets denote an
average over time and molecules. The carbon atom axes are defined:
$\mathbf{\hat{z}}\rightarrow$ vector from $C_{n-1}$ to $C_{n+1}$;
$\mathbf{\hat{y}}\rightarrow$ vector that is perpendicular to $z$ and
in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$;
$\mathbf{\hat{x}}\rightarrow$ vector perpendicular to
$\mathbf{\hat{y}}$ and $\mathbf{\hat{z}}$.

The order parameter has a range of $[1,-\frac{1}{2}]$. A value of 1
implies full order aligned to the bilayer axis, 0 implies full
disorder, and $-\frac{1}{2}$ implies full order perpendicular to the
bilayer axis. The {\sc cd} bond vector for carbons near the head group
are usually ordered perpendicular to the bilayer normal, with tails
farther away tending toward disorder. This makes the order parameter
negative for most carbons, and as such $|S_{\text{{\sc cd}}}|$ is more
commonly reported than $S_{\text{{\sc cd}}}$.

Fig.~\ref{lipidFig:scdFig} shows the $S_{\text{{\sc cd}}}$ order
parameters for the bilayer system at 300~K. There is no appreciable
difference in the plots for the various temperatures, however, there
is a larger difference between our models ordering, and that of
DMPC. As our values are closer to $-\frac{1}{2}$, this implies more
ordering perpendicular to the normal than in a real system. This is
due to the model having only one carbon group separating the chains
from the top of the lipid. In DMPC, with the flexibility inherent in a
multiple atom head group, as well as a glycerol linkage between the
head group and the acyl chains, there is more loss of ordering by the
point when the chains start.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{scdFig.eps}
\caption[$\text{S}_{\text{{\sc cd}}}$ order parameter for our model]{A comparison of $|\text{S}_{\text{{\sc cd}}}|$ between our model (blue) and DMPC\cite{petrache00} (black) near 300~K.}
\label{lipidFig:scdFig}
\end{figure}

\subsection{\label{lipidSec:p2Order}$P_2$ Order Parameter}

The $P_2$ order parameter allows us to measure the amount of
directional ordering that exists in the bilayer. Each lipid molecule
can be thought of as a cylindrical tube with the head group at the
top. If all of the cylinders are perfectly aligned, the $P_2$ order
parameter will be $1.0$. If the cylinders are completely dispersed,
the $P_2$ order parameter will be 0. For a collection of unit vectors,
the $P_2$ order parameter can be solved via the following
method.\cite{zannoni94}

Define an ordering matrix $\mathbf{Q}$, such that,
\begin{equation}
\mathbf{Q} = \frac{1}{N}\sum_i^N %
        \begin{pmatrix} %
        u_{ix}u_{ix}-\frac{1}{3} & u_{ix}u_{iy} & u_{ix}u_{iz} \\
        u_{iy}u_{ix} & u_{iy}u_{iy}-\frac{1}{3} & u_{iy}u_{iz} \\
        u_{iz}u_{ix} & u_{iz}u_{iy} & u_{iz}u_{iz}-\frac{1}{3} %
        \end{pmatrix}
\label{lipidEq:po1}
\end{equation}
Where the $u_{i\alpha}$ is the $\alpha$ element of the unit vector
$\mathbf{\hat{u}}_i$, and the sum over $i$ averages over the whole
collection of unit vectors. This allows the matrix element
$Q_{\alpha\beta}$ to be written:
\begin{equation}
Q_{\alpha\beta} = \langle u_{\alpha}u_{\beta} - 
        \frac{1}{3}\delta_{\alpha\beta} \rangle
\label{lipidEq:po2}
\end{equation}

Having constructed the matrix, diagonalizing $\mathbf{Q}$ yields three
eigenvalues and eigenvectors. The eigenvector associated with the
largest eigenvalue, $\lambda_{\text{max}}$, is the director for the
system of unit vectors. The director is the average direction all of
the unit vectors are pointing. The $P_2$ order parameter is then
simply
\begin{equation}
\langle P_2 \rangle = \frac{3}{2}\lambda_{\text{max}}
\label{lipidEq:po3}
\end{equation}

Table~\ref{lipidTab:blSummary} summarizes the $P_2$ values for the
bilayers, as well as the dipole orientations. The unit vector for the
lipid molecules was defined by finding the moment of inertia for each
lipid, then setting $\mathbf{\hat{u}}$ to point along the axis of
minimum inertia. For the {\sc head} atoms, the unit vector simply
pointed in the same direction as the dipole moment. For the lipid
molecules, the ordering was consistent across all temperatures, with
the director pointed along the $z$ axis of the box. More
interestingly, is the high degree of ordering the dipoles impose on
the {\sc head} atoms. The directors for the dipoles consistently
pointed along the plane of the bilayer, with the directors
anti-aligned on the top and bottom leaf.

\begin{table}
\caption[Structural properties of the bilayers]{Bilayer Structural properties as a function of temperature.}
\label{lipidTab:blSummary}
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
Temperature (K) & $\langle L_{\perp}\rangle$ ($\mbox{\AA}$) & %
        $\langle A_{\parallel}\rangle$ ($\mbox{\AA}^2$) & %
        $\langle P_2\rangle_{\text{Lipid}}$ & %
        $\langle P_2\rangle_{\text{{\sc head}}}$ \\ \hline
270 & 18.1 & 58.1 & 0.253 & 0.494 \\ \hline
275 & 17.2 & 56.7 & 0.295 & 0.514 \\ \hline
277 & 16.9 & 58.0 & 0.301 & 0.541 \\ \hline
280 & 17.4 & 58.0 & 0.274 & 0.488 \\ \hline
285 & 16.9 & 57.6 & 0.270 & 0.616 \\ \hline
290 & 17.0 & 57.6 & 0.263 & 0.534 \\ \hline
293 & 17.5 & 58.0 & 0.227 & 0.643 \\ \hline
300 & 16.9 & 57.6 & 0.315 & 0.536 \\ \hline
\end{tabular}
\end{center}
\end{table}

\subsection{\label{lipidSec:miscData}Further Bilayer Data}

Also summarized in Table~\ref{lipidTab:blSummary}, are the bilayer
thickness and area per lipid. The bilayer thickness was measured from
the peak to peak {\sc head} atom distance in the density profiles. The
area per lipid data compares favorably with values typically seen for
DMPC (60.0~$\mbox{\AA}^2$ at 303~K)\cite{petrache00}. Although are
values are lower this is most likely due to the shorter chain length
of our model (8 versus 14 for DMPC).

\subsection{\label{lipidSec:randBilayer}Bilayer Aggregation}

A very important accomplishment for our model is its ability to
spontaneously form bilayers from a randomly dispersed starting
configuration. Fig.~\ref{lipidFig:blImage} shows an image sequence for
the bilayer aggregation. After 3.0~ns, the basic form of the bilayer
can already be seen. By 7.0~ns, the bilayer has a lipid bridge
stretched across the simulation box to itself that will turn out to be
very long lived ($\sim$20~ns), as well as a water pore, that will
persist for the length of the current simulation. At 24~ns, the lipid
bridge is dispersed, and the bilayer is still integrating the lipid
molecules from the bridge into itself, and has still been unable to
disperse the water pore.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{bLayerImage.eps}
\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.}
\label{lipidFig:blImage}
\end{figure}

\subsection{\label{lipidSec:randIrod}Inverted Rod Aggregation}

Fig.~\ref{lipidFig:iRimage} shows a second aggregation sequence
simulated in this research. Here the fraction of water had been
significantly decreased to observe how the model would respond. After
1.5~ns, The main body of water in the system has already collected
into a central water channel. By 10.0~ns, the channel has widened
slightly, but there are still many sub channels permeating the lipid
macro-structure. At 23.0~ns, the central water channel has stabilized
and several smaller water channels have been absorbed to main
one. However, there are still several other channels that persist
through the lipid structure.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{iRodImage.eps}
\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}.}
\label{lipidFig:iRimage}
\end{figure}

\section{\label{lipidSec:Conclusion}Conclusion}

We have presented a phospholipid model capable of spontaneous
aggregation into a bilayer and an inverted rod structure. The time
scales of the macro-molecular aggregations are in excess of 24~ns. In
addition the model's bilayer properties have been explored over a
range of temperatures through prefabricated bilayers. No freezing
transition is seen in the temperature range of our current
simulations. However, structural information from the lowest
temperature may imply that a freezing event is on a much longer time
scale than that explored in this current research. Further studies of
this system could extend the time length of the simulations at the low
temperatures to observe whether lipid crystallization occurs within the
framework of this model.
Revision:	1087
Committed:	Fri Mar 5 22:16:34 2004 UTC (20 years, 4 months ago) by mmeineke
Content type:	application/x-tex
File size:	28993 byte(s)
Log Message:	semi-final draft