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 lenghths 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:Ensembles}.


\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 macrostructures: $\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:scd}$\text{S}_{\text{{\sc cd}}}$ order parameters}

The $\text{S}_{\text{{\sc cd}}}$ order parameter is often reported in
the experimental charecterizations 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<(3\cos\Theta_i\cos\Theta_j - \delta_{ij})\Bigr>
\label{lipidEq:scd2}
\end{equation}
Here, $\Theta_i$ is the angle the $i$th carbon atom frame axis 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 perpindicular to $z$ and
in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$;
$\mathbf{\hat{x}}\rightarrow$ vector perpindicular 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 perpindicular to the
bilayer axis. The {\sc cd} bond vector for carbons near the head group
are usually ordered perpindicular to the bilayer normal, with tails
farther away tending toward disorder. This makes the order paramter
negative for most carbons, and as such $|S_{\text{{\sc cd}}}|$ is more
commonly reported than $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}


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


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

\begin{table}
\caption[Structural properties of the bilayers]{Bilayer Structural properties as a function of temperature.}
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
Temperature (K) & $<L_{\perp}>$ ($\mbox{\AA}$) & %
        $<A_{\parallel}>$ ($\mbox{\AA}^2$) & $<P_2>_{\text{Lipid}}$ & %
        $<P_2>_{\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}
Revision:	1083
Committed:	Fri Mar 5 00:07:19 2004 UTC (20 years, 6 months ago) by mmeineke
Content type:	application/x-tex
File size:	19344 byte(s)
Log Message:	added several figures to the lipid paper, changed the ndthesis class around a little. some small fixes in oopse.