matt_papers/canidacy_paper/canidacy_paper.tex

\documentclass[11pt]{article}

\usepackage{graphicx}
\usepackage{color}
\usepackage{floatflt}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage[ref]{overcite}


\pagestyle{plain}
\pagenumbering{arabic}
\oddsidemargin 0.0cm \evensidemargin 0.0cm
\topmargin -21pt \headsep 10pt
\textheight 9.0in \textwidth 6.5in
\brokenpenalty=10000
\renewcommand{\baselinestretch}{1.2}
\renewcommand\citemid{\ } % no comma in optional reference note


\begin{document}


\title{A Mesoscale Model for Phospholipid Simulations}

\author{Matthew A. Meineke\\
Department of Chemistry and Biochemistry\\
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}
\maketitle

\section{Background and Research Goals}

Simulations of phospholipid bilayers are, by necessity, quite
complex. The lipid molecules are large molecules containing many
atoms,and the head group of the lipid will typically contain charge
separated ions which set up a large dipole within the molecule. Adding
to the complexity are the number of water molecules needed to properly
solvate the lipid bilayer. Because of these factors, many current
simulations are limited in both length and time scale due to to the
sheer number of calculations performed at every time step and the
lifetime of the researcher. A typical
simulation\cite{saiz02,lindahl00,venable00,Marrink01} will have around
64 phospholipids forming a bilayer approximately 40~$\mbox{\AA}$ by
50~$\mbox{\AA}$ with roughly 25 waters for every lipid. This means
there are on the order of 8,000 atoms needed to simulate these systems
and the trajectories in turn are integrated for times up to 10 ns.

These limitations make it difficult to study certain biologically
interesting phenomena that don't fit within the short time and length
scale requirements. One such phenomena is the existence of the ripple
phase ($P_{\beta'}$) of the bilayer between the gel phase
($L_{\beta'}$) and the fluid phase ($L_{\alpha}$). The $P_{\beta'}$
phase has been shown to have a ripple period of
100-200~$\mbox{\AA}$.\cite{katsaras00,sengupta00} A simulation of this
length scale would require approximately 1,300 lipid molecules and
roughly 25 waters for every lipid to fully solvate the bilayer. With
the large number of atoms involved in a simulation of this magnitude,
steps \emph{must} be taken to simplify the system to the point where
these numbers are reasonable.

Another system of interest would be drug molecule diffusion through
the membrane. Due to the fluid like properties of a lipid membrane,
not all diffusion takes place at membrane channels. It is of interest
to study certain molecules that may incorporate themselves directly
into the membrane. These molecules may then have an appreciable
waiting time (on the order of nanoseconds) within the
bilayer. Simulation of such a long time scale again requires
simplification of the system in order to lower the number of
calculations needed at each time step.


\section{Methodology}

\subsection{Length and Time Scale Simplifications}

The length scale simplifications are aimed at increasing the number of
molecules that can be simulated without drastically increasing the
computational cost of the system. This is done by a combination of
substituting less expensive interactions for expensive ones and
decreasing the number of interaction sites per molecule. Namely,
point charge distributions are replaced with dipoles, and unified atoms are
used in place of water, phospholipid head groups, and alkyl groups.

The replacement of charge distributions with dipoles allows us to
replace an interaction that has a relatively long range,
$(\frac{1}{r})$, for the coulomb potential, with that of a relatively
short range, $(\frac{1}{r^{3}})$, for dipole - dipole
potentials. Combined with a computational simplification algorithm
such as a Verlet neighbor list,\cite{allen87:csl} this should give
computational scaling by $N$. This is in comparison to the Ewald
sum\cite{leach01:mm} needed to compute the coulomb interactions, which
scales at best by $N \ln N$.

The second step taken to simplify the number of calculations is to
incorporate unified models for groups of atoms. In the case of water,
we use the soft sticky dipole (SSD) model developed by
Ichiye\cite{Liu96} (Section~\ref{sec:ssdModel}). For the phospholipids, a
unified head atom with a dipole will replace the atoms in the head
group, while unified $\text{CH}_2$ and $\text{CH}_3$ atoms will
replace the alkyl units in the tails (Section~\ref{sec:lipidModel}).

The time scale simplifications are introduced so that we can take
longer time steps. By increasing the size of the time steps taken by
the simulation, we are able to integrate the simulation trajectory
with fewer calculations. However, care must be taken that any
simplifications used, still conserve the total energy of the
simulation. In practice, this means taking steps small enough to
resolve all motion in the system without accidently moving an object
too far along a repulsive energy surface before it feels the effect of
the surface.

In our simulation we have chosen to constrain all bonds to be of fixed
length. This means the bonds are no longer allowed to vibrate about
their equilibrium positions. Bond vibrations are typically the fastest
periodic motion in a dynamics simulation. By taking this action, we
are able to take time steps of 3 fs while still maintaining constant
energy. This is in contrast to the 1 fs time step typically needed to
conserve energy when bonds lengths are allowed to oscillate.

\subsection{The Soft Sticky Water Model}
\label{sec:ssdModel}

%\begin{floatingfigure}{55mm}
%\includegraphics[width=45mm]{ssd.epsi}
%\caption{The SSD model with the oxygen and hydrogen atoms drawn in for reference. \vspace{5mm}}
%\label{fig:ssdModel}
%\end{floatingfigure} 

The water model used in our simulations is a modified soft
Stockmayer-sphere model.\cite{stevens95} Like the Stockmayer-sphere, the SSD
model\cite{Liu96} consists of a Lennard-Jones interaction site and a
dipole both located at the water's center of mass (Figure
\ref{fig:ssdModel}). However, the SSD model extends this by adding a
tetrahedral potential to correct for hydrogen bonding.

The SSD water potential is then given by the following equation:
\begin{equation}
V_{\text{SSD}} = V_{\text{LJ}}(r_{i\!j}) + V_{\text{dp}}(\mathbf{r}_{i\!j},
        \boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
        + V_{\text{sp}}(\mathbf{r}_{i\!j},\boldsymbol{\Omega}_{i},
        \boldsymbol{\Omega}_{j})
\label{eq:ssdTotPot}
\end{equation}
$V_{\text{LJ}}$ is the Lennard-Jones potential:
\begin{equation}
V_{\text{LJ}} = 
        4\epsilon_{ij} \biggl[
        \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
        - \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
        \biggr]
\label{eq:lennardJonesPot}
\end{equation}
where $r_{ij}$ is the distance between two $ij$ pairs, $\sigma_{ij}$
scales the length of the interaction, and $\epsilon_{ij}$ scales the
energy of the potential. For SSD, $\sigma_{\text{SSD}} = 3.051 \mbox{
\AA}$ and $\epsilon_{\text{SSD}} = 0.152\text{ kcal/mol}$.
$V_{\text{dp}}$ is the dipole potential:
\begin{equation}
V_{\text{dp}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
        \boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
        \frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
        -
        \frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
                (\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
                {r^{5}_{ij}} \biggr]
\label{eq:dipolePot}
\end{equation}
where $\mathbf{r}_{ij}$ is the vector between $i$ and $j$,
$\boldsymbol{\Omega}$ is the orientation of the species, and
$\boldsymbol{\mu}$ is the dipole vector. The SSD model specifies a dipole
magnitude of 2.35~D for water.

The hydrogen bonding of the model is governed by the $V_{\text{sp}}$
term of the potential. Its form is as follows:
\begin{equation}
V_{\text{sp}}(\mathbf{r}_{i\!j},\boldsymbol{\Omega}_{i},
        \boldsymbol{\Omega}_{j}) =
        v^{\circ}[s(r_{ij})w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
                \boldsymbol{\Omega}_{j})
        +
        s'(r_{ij})w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
                \boldsymbol{\Omega}_{j})]
\label{eq:spPot}
\end{equation}
Where $v^\circ$ scales strength of the interaction.
$w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
and
$w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
are responsible for the tetrahedral potential and a correction to the
tetrahedral potential respectively. They are,
\begin{equation}
w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) =
        \sin\theta_{ij} \sin 2\theta_{ij} \cos 2\phi_{ij}
        + \sin \theta_{ji} \sin 2\theta_{ji} \cos 2\phi_{ji}
\label{eq:spPot2}
\end{equation}
and
\begin{equation}
\begin{split}
w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) =
        &(\cos\theta_{ij}-0.6)^2(\cos\theta_{ij} + 0.8)^2 \\
        &+ (\cos\theta_{ji}-0.6)^2(\cos\theta_{ji} + 0.8)^2 - 2w^{\circ}
\end{split}
\label{eq:spCorrection}
\end{equation}
The angles $\theta_{ij}$ and $\phi_{ij}$ are defined by the spherical
polar coordinates of the position of sphere $j$ in the reference frame
fixed on sphere $i$ with the z-axis aligned with the dipole moment.
The correction
$w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
is needed because
$w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
vanishes when $\theta_{ij}$ is $0^\circ$ or $180^\circ$.

Finally, the sticky potential is scaled by a cutoff function,
$s(r_{ij})$ that scales smoothly between 0 and 1. It is represented
by:
\begin{equation}
s(r_{ij}) = 
        \begin{cases}
        1&      \text{if $r_{ij} < r_{L}$}, \\
        \frac{(r_{U} - r_{ij})^2 (r_{U} + 2r_{ij} 
                - 3r_{L})}{(r_{U}-r_{L})^3}&
                \text{if $r_{L} \leq r_{ij} \leq r_{U}$},\\
        0&      \text{if $r_{ij} \geq r_{U}$}.
        \end{cases}
\label{eq:spCutoff}
\end{equation}

Despite the apparent complexity of Equation \ref{eq:spPot}, the SSD
model is still computationally inexpensive. This is due to Equation
\ref{eq:spCutoff}. With $r_{L}$ being 2.75~$\mbox{\AA}$ and $r_{U}$
being equal to either 3.35~$\mbox{\AA}$ for $s(r_{ij})$ or
4.0~$\mbox{\AA}$ for $s'(r_{ij})$, the sticky potential is only active
over an extremely short range, and then only with other SSD
molecules. Therefore, it's predominant interaction is through it's
point dipole and Lennard-Jones sphere.

\subsection{The Phospholipid Model}
\label{sec:lipidModel}

%\begin{floatingfigure}{90mm}
%\includegraphics[angle=-90,width=80mm]{lipidModel.epsi}
%\caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length. \vspace{5mm}}
%\label{fig:lipidModel}
%\end{floatingfigure}

The lipid molecules in our simulations are unified atom models. Figure
\ref{fig:lipidModel} shows a template drawing for one of our
lipids. The Head group of the phospholipid is replaced by a single
Lennard-Jones sphere with a freely oriented dipole placed at it's
center. The magnitude of it's dipole moment is 20.6 D. The tail atoms
are unified $\text{CH}_2$ and $\text{CH}_3$ atoms and are also modeled
as Lennard-Jones spheres. The total potential for the lipid is
represented by Equation \ref{eq:lipidModelPot}.

\begin{equation}
V_{\mbox{lipid}} = \overbrace{%
        V_{\text{bend}}(\theta_{ijk}) +% 
        V_{\text{tors.}}(\phi_{ijkl})}^{bonded}
        + \overbrace{%
        V_{\text{LJ}}(r_{i\!j}) + 
        V_{\text{dp}}(r_{i\!j},\Omega_{i},\Omega_{j})%
        }^{non-bonded}
\label{eq:lipidModelPot}
\end{equation}

The non-bonded interactions, $V_{\text{LJ}}$ and $V_{\text{dp}}$, are
the Lennard-Jones and dipole-dipole interactions respectively. For the
non-bonded potentials, only the bend and the torsional potentials are
calculated. The bond potential is not calculated, and the bond lengths
are constrained via RATTLE.\cite{leach01:mm} The bend potential is of
the form: 
\begin{equation}
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}\frac{(\theta_{ijk} - \theta_0)^2}{2}
\label{eq:bendPot}
\end{equation}
Where $k_{\theta}$ sets the stiffness of the bend potential, and $\theta_0$
sets the equilibrium bend angle. The torsion potential was given by:
\begin{equation}
V_{\text{tors.}}(\phi_{ijkl})= c_1[1+\cos\phi_{ijkl}]
        + c_2 [1 - \cos(2\phi_{ijkl})] + c_3[1 + \cos(3\phi_{ijkl})]
\label{eq:torsPot}
\end{equation}
All parameters for bonded and non-bonded potentials in the tail atoms
were taken from TraPPE.\cite{Siepmann1998} The bonded interactions for
the head atom were also taken from TraPPE, however it's dipole moment
and mass were based on the properties of the phosphatidylcholine head
group. The Lennard-Jones parameter for the head group was chosen such
that it was roughly twice the size of a $\text{CH}_3$ atom, and it's
well depth was set to be approximately equal to that of $\text{CH}_3$.

\section{Initial Simulation: 25 lipids in water}
\label{sec:5x5}

\subsection{Starting Configuration and Parameters}
\label{sec:5x5Start}

Our first simulation is an array of 25 single chained lipids in a sea
of water (Figure \ref{fig:5x5Start}). The total number of water
molecules is 1386, giving a final of water concentration of 70\%
wt. The simulation box measures 34.5~$\mbox{\AA}$ x 39.4~$\mbox{\AA}$
x 39.4~$\mbox{\AA}$ with periodic boundary conditions imposed. The
system is simulated in the micro-canonical (NVE) ensemble with an
average temperature of 300~K.

\subsection{Results}
\label{sec:5x5Results}

Figure \ref{fig:5x5Final} shows a snapshot of the system at
3.6~ns. Note that the system has spontaneously self assembled into a
bilayer. Discussion of the length scales of the bilayer will follow in
this section. However, it is interesting to note a key qualitative
property of the system revealed by this snapshot, the tail chains are
tilted to the bilayer normal. This is usually indicative of the gel
($L_{\beta'}$) phase. In this system, the box size is probably too
small for the bilayer to relax to the fluid ($P_{\alpha}$) phase. This
demonstrates a need for an isobaric-isothermal ensemble where the box
size may relax or expand to keep the system at a 1~atm.

The simulation was analyzed using the radial distribution function,
$g(r)$, which has the form:
\begin{equation}
g(r) = \frac{V}{N_{\text{pairs}}}\langle \sum_{i} \sum_{j > i} 
        \delta(|\mathbf{r} - \mathbf{r}_{ij}|) \rangle
\label{eq:gofr}
\end{equation}
Equation \ref{eq:gofr} gives us information about the spacing of two
species as a function of radius. Essentially, if the observer were
located at atom $i$ and were looking out in all directions, $g(r)$
shows the relative density of atom $j$ at any given radius, $r$,
normalized by the expected density of atom $j$ at $r$. In a
homogeneously mixed fluid, $g(r)$ will approach 1 at large $r$, as a
fluid contains no long range structure to contribute peaks in the
number density.

For the species containing dipoles, a second pair wise distribution
function was used, $g_{\gamma}(r)$. It is of the form:
\begin{equation}
g_{\gamma}(r) = foobar
\label{eq:gammaofr}
\end{equation}
Where $\gamma_{ij}$ is the angle between the dipole of atom $j$ with
respect to the dipole of atom $i$. This correlation will vary between
$+1$ and $-1$ when the two dipoles are perfectly aligned and
anti-aligned respectively. This then gives us information about how
directional species are aligned with each other as a function of
distance.

Figure \ref{fig:5x5HHCorr} shows the two self correlation functions
for the Head groups of the lipids. The first peak in the $g(r)$ at
4.03~$\mbox{\AA}$ is the nearest neighbor separation of the heads of
two lipids. This corresponds to a maximum in the $g_{\gamma}(r)$ which
means that the two neighbors on the same monolayer have their dipoles
aligned. The broad peak at 6.5~$\mbox{\AA}$ is the inter-bilayer
spacing. Here, there is a corresponding anti-alignment in the angular
correlation. This means that although the dipoles are aligned on the
same monolayer, the dipoles will orient themselves to be anti-aligned
on a opposite facing monolayer. With this information, the two peaks
at 7.0~$\mbox{\AA}$ and 7.4~$\mbox{\AA}$ are head groups on the same
monolayer, and they are the second nearest neighbors to the head
group. The peak is likely a split peak because of the small statistics
of this system. Finally, the peak at 8.0~$\mbox{\AA}$ is likely the
second nearest neighbor on the opposite monolayer because of the
anti-alignment evident in the angular correlation.

Figure \ref{fig:5x5CCg} shows the radial distribution function for the
$\text{CH}_2$ unified atoms. The spacing of the atoms along the tail
chains accounts for the regularly spaced sharp peaks, but the broad
underlying peak with its maximum at 4.6~$\mbox{\AA}$ is the
distribution of chain-chain distances between two lipids. The final
Figure, Figure \ref{fig:5x5HXCorr}, includes the correlation functions
between the Head group and the SSD atoms. The peak in $g(r)$ at
3.6~$\mbox{\AA}$ is the distance of closest approach between the two,
and $g_{\gamma}(r)$ shows that the SSD atoms will align their dipoles
with the head groups at close distance. However, as one increases the
distance, the SSD atoms are no longer aligned.

\section{Second Simulation: 50 randomly oriented lipids in water}
\label{sec:r50}

\subsection{Starting Configuration and Parameters}
\label{sec:r50Start}

The second simulation consists of 50 single chained lipid molecules
embedded in a sea of 1384 SSD waters (54\% wt.). The lipids in this
simulation were started with random orientation and location (Figure
\ref{fig:r50Start} ) The simulation box measured 34.5~$\mbox{\AA}$ x
39.4~$\mbox{\AA}$ x 39.4~$\mbox{\AA}$ with periodic boundary conditions
imposed. The simulation was run in the NVE ensemble with an average
temperature of 300~K.

\subsection{Results}
\label{sec:r50Results}

Figure \ref{fig:r50Final} is a snapshot of the system at 2.0~ns. Here
we see that the system has already aggregated into several micelles
and two are even starting to merge. It will be interesting to watch as
this simulation continues what the total time scale for the micelle
aggregation and bilayer formation will be.

Figures \ref{fig:r50HHCorr}, \ref{fig:r50CCg}, and \ref{fig:r50} are
the same correlation functions for the random 50 simulation as for the
previous simulation of 25 lipids. What is most interesting to note, is
the high degree of similarity between the correlation functions for
each system. Even though the 25 lipid simulation formed a bilayer and
the random 50 simulation is still in the micelle stage, both have a
inter surface spacing of 6.5~$\mbox{\AA}$ with the same characteristic
anti-alignment between surfaces. Not as surprising, is the consistency
of the closest packing statistics between systems. Namely, a head-head
closest approach distance of 4~$\mbox{\AA}$, and similar findings for
the chain-chain and head-water distributions as in the 25 lipid
system.

\section{Future Directions}

Current simulations indicate that our model is a feasible one, however
improvements will need to be made to allow the system to simulate an
isobaric-isothermal ensemble. This will allow the system to relax to
an equilibrium configuration at room temperature and pressure allowing
us to compare our model to experimental results. Also, we plan to
parallelize the code for an even greater speedup. This will allow us
to simulate the size systems needed to examine phenomena such as the
ripple phase and drug molecule diffusion

Once the work has completed on the simulation engine, we would then
like to use it to explore phase diagram for our model. By
characterizing how our model parameters affect the bilayer properties,
we hope to tailor our model to more closely match real biological
molecules. With this information, we then hope to incorporate
biologically relevant molecules into the system and observe their
transport properties across the membrane.

\section{Acknowledgments}

I would like to thank Dr. J.Daniel Gezelter for his guidance on this
project. I would also like to acknowledge the following for their help
and discussions during this project: Christopher Fennell, Charles
Vardeman, Teng Lin, Megan Sprague, Patrick Conforti, and Dan Combest.

\pagebreak 
\bibliographystyle{achemso} 
\bibliography{canidacy_paper}
\end{document}
Revision:	108
Committed:	Wed Sep 11 15:58:41 2002 UTC (22 years ago) by mmeineke
Content type:	application/x-tex
File size:	20290 byte(s)
Log Message:	text is finished YAY!!!