--- trunk/matt_papers/canidacy_paper/canidacy_paper.tex	2002/08/26 22:09:28	99
+++ trunk/matt_papers/canidacy_paper/canidacy_paper.tex	2002/09/11 03:51:11	107
@@ -1,6 +1,7 @@
 \documentclass[11pt]{article}
 
 \usepackage{graphicx}
+\usepackage{color}
 \usepackage{floatflt}
 \usepackage{amsmath}
 \usepackage{amssymb}
@@ -33,102 +34,148 @@ Notre Dame, Indiana 46556}
 
 \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 increaseing the number
-of molecules simulated without drastically increasing the
+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,
-charge distributions are replaced with dipoles, and unified atoms are
-used in place of water and phospholipid head groups.
+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 charge charge potential, with that of a relitively short
-range, $\frac{1}{r^{3}}$ for dipole - dipole potentials
-(Equations~\ref{eq:dipolePot} and \ref{eq:chargePot}). This allows us
-to use computaional simplifications algorithms such as Verlet neighbor
-lists,\cite{allen87:csl} which gives computaional scaling by $N$. This
-is in comparison to the Ewald sum\cite{leach01:mm} needed to compute
-the charge - charge interactions which scales at best by $N
-\ln N$.
+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$.
 
-\begin{equation}
-V^{\text{dp}}_{ij}(\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} \cdot \mathbf{r}_{ij}) %
-		(\boldsymbol{\mu} \cdot \mathbf{r}_{ij}) }{r^{5}_{ij}} \biggr]
-\label{eq:dipolePot}
-\end{equation}
-
-\begin{equation}
-V^{\text{ch}}_{ij}(\mathbf{r}_{ij}) = \frac{q_{i}q_{j}}%
-	{4\pi\epsilon_{0} r_{ij}}
-\label{eq:chargePot}
-\end{equation} 
-
-The second step taken to simplify the number of calculationsis to
+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 alkanes in the tails (Section~\ref{sec:lipidModel}).
+replace the alkyl units in the tails (Section~\ref{sec:lipidModel}).
 
-The time scale simplifications are taken so that the simulation can
-take long time steps. By incresing the time steps taken by the
-simulation, we are able to integrate the simulation trajectory with
-fewer calculations. However, care must be taken to conserve the energy
-of the simulation. This is a constraint placed upon the system by
-simulating in the microcanonical ensemble. 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 affect of the surface.
+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, typically the fastest periodical 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.
+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}}
-% The dipole magnitude is 2.35 D and the Lennard-Jones parameters are $\sigma = 3.051 \mbox{\AA}$ and $\epsilon = 0.152$ kcal/mol.}
-\label{fig:ssdModel}
-\end{floatingfigure} 
+%\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. Like the soft Stockmayer sphere, the SSD
+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.
 
-This SSD water's motion is then governed by the following potential:
+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},
+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 with $\sigma_{\text{w}}
-= 3.051 \mbox{ \AA}$ and $\epsilon_{\text{w}} = 0.152\text{
-kcal/mol}$. $V_{\text{dp}}$ is the dipole potential with
-$|\mu_{\text{w}}| = 2.35\text{ D}$.
+$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 iteraction, 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 potentail. Its form is as follows:
+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}) =
@@ -139,8 +186,7 @@ Where $v^\circ$ is responsible for scaling the strengt
 		\boldsymbol{\Omega}_{j})]
 \label{eq:spPot}
 \end{equation}
-Where $v^\circ$ is responsible for scaling the strength of the
-interaction.
+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})$
@@ -150,7 +196,7 @@ w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsy
 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:apPot2}
+\label{eq:spPot2}
 \end{equation}
 and
 \begin{equation}
@@ -161,13 +207,18 @@ The correction 
 \end{split}
 \label{eq:spCorrection}
 \end{equation}
-The correction 
+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$. 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 alligned with the dipole moment. 
+vanishes when $\theta_{ij}$ is $0^\circ$ or $180^\circ$.
 
-Finaly, the  sticky potentail is scaled by a cutoff function, $s(r_{ij})$ that scales smoothly between 0 and 1. It is represented by:
+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}
@@ -180,22 +231,30 @@ s(r_{ij}) = 
 \label{eq:spCutoff}
 \end{equation}
 
+Despite the apparent complexity of Equation \ref{eq:spPot}, the SSD
+model is still computationaly 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 extremly 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}
+%\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 unifed $\text{CH}_2$ and $\text{CH}_3$ atoms and are also modeled
+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}.
 
@@ -210,7 +269,7 @@ The non bonded interactions, $V_{\text{LJ}}$ and $V_{\
 \label{eq:lipidModelPot}
 \end{equation}
 
-The non bonded interactions, $V_{\text{LJ}}$ and $V_{\text{dp}}$, are
+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
@@ -220,6 +279,91 @@ V_{\text{bend}}(\theta_{ijk}) = k_{\theta}\frac{(\thet
 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$.
 
- \pagebreak \bibliographystyle{achemso}
+\section{Initial Simulation: 25 lipids in water}
+\label{sec:5x5}
+
+\subsection{Starting Configuration and Parameters}
+\label{sec:5x5Start}
+
+Our first simulation was an array of 25 single chained lipids in a sea
+of water (Figure \ref{fig:5x5Start}). The total number of water
+molecules was 1386, giving a final of water concentration of 70\%
+wt. The simulation box measured 34.5~$\mbox{\AA}$ x 39.4~$\mbox{\AA}$
+x 39.4~$\mbox{\AA}$ with periodic boundary conditions imposed. The
+system was 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 at 4.03~$\mbox{\AA}$ is the
+nearest neighbor separation of the heads of two lipids.
+
+
+
+
+\section{Second Simulation: 50 randomly oriented lipids in water}
+
+the second simulation
+
+\section{Future Directions}
+
+
+\pagebreak
+\bibliographystyle{achemso}
 \bibliography{canidacy_paper} \end{document}
