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 ($L_{\beta}$) phase 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, direct Ewald sum or the slightly faster particle
mesh Ewald (PME) 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 technologies. One
such example is the observance of a ripple ($P_{\beta^{\prime}}$)
phase which appears between the $L_{\beta}$ and $L_{\alpha}$ phases of
certain phospholipid bilayers
(Fig.~\ref{lipidFig:phaseDiag}).\cite{katsaras00,sengupta00} These
ripples are known from x-ray diffraction data to have periodicities on
the order of 100-200~$\mbox{\AA}$.\cite{katsaras00} 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.

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

The time and length scale limitations are most striking in transport
phenomena.  Due to the fluid-like properties of lipid membranes, not
all small molecule diffusion across the membranes happens at pores.
Some molecules of interest may incorporate themselves directly into
the membrane.  Once there, they may exhibit appreciable waiting times
(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 Lipowsky\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 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. 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 simplification of atomistic detail.  The model uses
unified-atom 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 therefore
eliminates the need for the 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 simulated 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. This effect has been missing from previous
reduced models.

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 possible long range interactions between atomic groups in the
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}
Here $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 intra-molecular bend and torsion
interactions.  The bond between two beads on a chain is of fixed
length, and is maintained using 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 is a bend
potential which keeps the dipole roughly perpendicular to the
molecular body, 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}). This bend mimics the hinge between
the phosphatidyl part of the PC head group and the remainder of the
molecule.  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, long-range interactions for
pairs of atoms not directly involved in a bond, bend, or torsion are
calculated.

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

All simulations presented here use a two-chain 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 (approximately half the magnitude of
the dipole on the PC head group\cite{Cevc87}), 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|c|}
\hline
     & mass (amu) & $\sigma$($\mbox{\AA}$)  & $\epsilon$ (kcal/mol) %
        & $|\mathbf{\hat{\mu}}|$ (Debye) \\ \hline
{\sc head} & 72  & 4.0 & 0.185 & 10.6 \\ \hline
{\sc ch}\cite{Siepmann1998} & 13.02 & 4.0 & 0.0189 & 0.0 \\ \hline
$\text{{\sc ch}}_2$\cite{Siepmann1998} & 14.03 & 3.95 & 0.18 & 0.0 \\ \hline
$\text{{\sc ch}}_3$\cite{Siepmann1998} & 15.04 & 3.75 & 0.25 & 0.0 \\ \hline
{\sc ssd}\cite{liu96:new_model} & 18.03 & 3.051 & 0.152 & 0.0 \\ \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}
\caption{BEND PARAMETERS FOR THE TWO CHAIN PHOSPHOLIPIDS}
\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}
\begin{minipage}{\linewidth}
\begin{center}
\vspace{2mm}
All alkane parameters are based off of those from TraPPE.\cite{Siepmann1998}
\end{center}
\end{minipage}
\end{center}
\end{table}

\begin{table}
\caption{TORSION PARAMETERS FOR THE TWO CHAIN PHOSPHOLIPIDS}
\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}
\begin{minipage}{\linewidth}
\begin{center}
\vspace{2mm}
All alkane parameters are based off of those from TraPPE.\cite{Siepmann1998}
\end{center}
\end{minipage}
\end{center}
\end{table}


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

As mentioned previously, the water model used throughout these
simulations was the {\sc ssd/e} model of Fennell and
Gezelter,\cite{fennell04} earlier forms of this model can be found in
Ichiye \emph{et
al}.\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}$. integrator described in
Sec.~\ref{oopseSec:integrate}. With $\tau_B = 1.5$~ps and $\tau_T =
1.2$~ps, the box volume stabilizes after 50~ps, and fluctuates about
its equilibrium value by $\sim 0.6\%$, temperature fluctuations are
about $\sim 1.4\%$ of their set value, and pressure fluctuations are
the largest, varying as much as $\pm 250$~atm. However, such large
fluctuations in pressure are typical for liquid state simulations.\cite{leach01:mm}


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

Two main classes of starting configurations were used in this research:
random and ordered bilayers.  The ordered bilayer simulations were all
started from an equilibrated bilayer configuration 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 principal axis perpendicular to the plane.
The bottom leaf simply mirrored the top leaf, and the appropriate
number of water molecules 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 water molecules were then tested against
the lipid using a 5.0~$\mbox{\AA}$ overlap test with any atom in the
lipid.  This gave an estimate for the number of water molecules each
lipid would displace in a simulation box. A target number of water
molecules was then defined which included the number of water
molecules each lipid would displace, the number of water molecules
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 water molecules 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 placed.  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 water
to lipid ratio was achieved.  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 under the
NPTxyz integrator.

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

The first class of simulations were started from ordered bilayers. All
configurations consisted of 60 lipid molecules with 30 lipids on each
leaf, and were hydrated with 1620 {\sc ssd/e} 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 configuration 7~ns in to the 300~K
simulation. Their temperatures were modified with the thermostatting
algorithm in the NPTxyz 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/e} molecules simulated at 300~K. The second
($\text{R}_{\text{II}}$) was 90 lipids with 1350 {\sc ssd/e} molecules
simulated at 350~K. Both simulations were integrated for more than
20~ns to observe whether our model is capable of spontaneous
aggregation into known 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: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 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. 

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.

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.

\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. 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}}}$.

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:
\begin{itemize}
\item $\mathbf{\hat{z}}$ is the vector from $C_{n-1}$ to $C_{n+1}$.
\item $\mathbf{\hat{y}}$ is the vector that is perpendicular to $z$ and
in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$.
\item $\mathbf{\hat{x}}$ is the vector perpendicular to
$\mathbf{\hat{y}}$ and $\mathbf{\hat{z}}$.
\end{itemize}
This assumes that the hydrogen atoms are always in a plane
perpendicular to the $C_{n-1}-C_{n}-C_{n+1}$ plane.

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 model's ordering, and the
experimentally observed ordering 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. 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}}}$.

\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 bodies of the molecules making
up the bilayer. Each lipid molecule can be thought of as a cylindrical
rod with the head group at the top. If all of the rods are perfectly
aligned, the $P_2$ order parameter will be $1.0$. If the rods are
completely disordered, the $P_2$ order parameter will be 0. For a
collection of unit vectors pointing along the principal axes of the
rods, the $P_2$ order parameter can be solved via the following
method.\cite{zannoni94}

Define an ordering tensor $\overleftrightarrow{\mathsf{Q}}$, such that,
\begin{equation}
\overleftrightarrow{\mathsf{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 tensor to be written:
\begin{equation}
\overleftrightarrow{\mathsf{Q}} = \frac{1}{N}\sum_i^N \biggl[
        \mathbf{\hat{u}}_i \otimes \mathbf{\hat{u}}_i 
        - \frac{1}{3} \cdot \mathsf{1} \biggr].
\label{lipidEq:po2}
\end{equation}

After constructing the tensor, diagonalizing
$\overleftrightarrow{\mathsf{Q}}$ yields three eigenvalues and
eigenvectors. The eigenvector associated with the largest eigenvalue,
$\lambda_{\text{max}}$, is the director axis  for the system of unit
vectors. The director axis 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 for 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 (the long axis). 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 themselves
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. 

\begin{table}
\caption{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 Structural Data}

Also summarized in Table~\ref{lipidTab:blSummary}, are the bilayer
thickness ($\langle L_{\perp}\rangle$) and area per lipid ($\langle
A_{\parallel}\rangle$). 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 our
values are lower this is most likely due to the shorter chain length
of our model (8 versus 14 for DMPC).

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

The lateral diffusion constant, $D_L$, is the constant characterizing
the diffusive motion of the lipid molecules within the plane of the bilayer. It
is given by the following Einstein relation:\cite{allen87:csl}
\begin{equation}
D_L = \lim_{t\rightarrow\infty}\frac{1}{4t}\langle |\mathbf{r}(t) 
        - \mathbf{r}(0)|^2\rangle,
\end{equation}
where $\mathbf{r}(t)$ is the $xy$ position of the lipid at time $t$
(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.

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 two orders of magnitude larger than those
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
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
the {\sc ssd} water, $9.84\times 10^{-6}\,\text{cm}^2/\text{s}$, is
reasonably close to the bulk water diffusion constant ($2.2999\times
10^{-5}\,\text{cm}^2/\text{s}$).\cite{Holz00}

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

\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: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 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} 

\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 water molecules permeating the
lipid macro-structure. At 35.0~ns, the central water channel has
stabilized and several smaller water channels have been absorbed by
the main one. However, there is still an appreciable water
concentration throughout 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 simple unified-atom phospholipid model capable of
spontaneous aggregation into a bilayer and an inverted rod
structure. The time scales of the macro-molecular aggregations are
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
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 270~K
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 can occur within
the framework of this model.

Potential problems that may be obstacles in further research, is the
lack of detail in the head region. As the chains are almost directly
attached to the {\sc head} atom, there is no buffer between the
actions of the head group and the tails. Another disadvantage of the
model is the dipole approximation will alter results when details
concerning a charged solute's interactions with the bilayer. However,
it is important to keep in mind that the dipole approximation can be
kept an advantage by examining solutes that do not require point
charges, or at the least, require only dipole approximations
themselves. Other advantages of the model include the ability to alter
the size of the unified-atoms so that the size of the lipid can be
increased without adding to the number of interactions in the
system. However, what sets our model apart from other current
simplified models,\cite{goetz98,marrink04} is the information gained
by observing the ordering of the head groups dipole's in relation to
each other and the solvent without the need for point charges and the
Ewald sum.
Revision:	1115
Committed:	Thu Apr 15 19:19:06 2004 UTC (21 years ago) by mmeineke
Content type:	application/x-tex
File size:	34393 byte(s)
Log Message:	final copies sent to the printer