trunk/mattDisertation/lipid.tex



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

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

In the past 10 years, computer speeds have allowed for the atomistic
simulation of phospholipid bilayers.  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 64 to 256
phospholipids.\cite{lindahl00,sum:2003,venable00,gomez:2003,smondyrev:1999,marrink01}
This is due to the expense of the computer calculations involved when
performing these simulations.  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 the computationally
expensive Ewald sum or its slightly faster derivative particle mesh
Ewald sum.\cite{nina:2002,norberg:2000,patra:2003} These factors all
limit the potential size and time lengths of bilayer simulations.

Unfortunately, much of biological interest happens on time and length
scales infeasible with current simulation. 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.

Another class of simulations to consider, are those dealing with the
diffusion of molecules through a bilayer.  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
difficulty to obtain when integrating the system with atomistic
detail.

Addressing 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}$), eliminating 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.

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 simulation setups
explored in this research are outlined in
Sec.~\ref{lipidSec:ExpSetup}. Sec.~\ref{lipidSec:Results} and
Sec.~\ref{lipidSec:Discussion} give a summary of the results and
interpretation of those results respectively.  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 unified atom beads.  They are free of partial charges or
dipoles, containing 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}
        -
        \frac{3(\boldsymbol{\hat{u}}_i \cdot \mathbf{r}_{ij}) %
                (\boldsymbol{\hat{u}}_j \cdot \mathbf{r}_{ij}) }
                {r^{2}_{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 of $\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. 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.

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 fudge factor to pad the
initialization. 

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 adsorbed lipid, its position and
orientation were rejected, and a new random adsorption site was
attempted. The RSA simulation proceeded until all phospholipids had
been adsorbed.  After adsorption, all water molecules with locations
that overlapped with the atomic coordinates of the lipids were
removed. 

Finally, water molecules were removed one by one 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
would shrink to an appropriate size within the first 50~ps of a
simulation in the $\text{NPT}_{xyz}$ ensemble.

\subsection{\label{lipidSec:Configs}The simulation configurations}

Table ~\ref{lipidTable:simNames} summarizes the names and important
details of the simulations.  The B set of simulations were all started
in an ordered bilayer and observed over a period of 10~ns. Simulation
RL was integrated for approximately 20~ns starting from a random
configuration as an example of spontaneous bilayer aggregation.
Lastly, simulation RH was also started from a random configuration,
but with a lesser water content and higher temperature to show the
spontaneous aggregation of an inverted hexagonal lamellar phase.
Revision:	1061
Committed:	Fri Feb 20 17:52:57 2004 UTC (21 years, 9 months ago) by mmeineke
Content type:	application/x-tex
File size:	14554 byte(s)
Log Message:	cleaned up the references and the tables in the lipid section. also spell checked the thing.
#	User	Rev	Content
1	mmeineke	971
2
3			\chapter{\label{chapt:lipid}Phospholipid Simulations}
4
5			\section{\label{lipidSec:Intro}Introduction}
6
7	mmeineke	1001	In the past 10 years, computer speeds have allowed for the atomistic
8			simulation of phospholipid bilayers. These simulations have ranged
9			from simulation of the gel phase ($L_{\beta}$) of
10	mmeineke	1061	dipalmitoylphosphatidylcholine (DPPC),\cite{lindahl00} to the
11	mmeineke	1001	spontaneous aggregation of DPPC molecules into fluid phase
12	mmeineke	1061	($L_{\alpha}$) bilayers.\cite{marrink01} With the exception of a few
13			ambitious
14			simulations,\cite{marrink01:undulation,marrink:2002,lindahl00} most
15	mmeineke	1001	investigations are limited to 64 to 256
16	mmeineke	1061	phospholipids.\cite{lindahl00,sum:2003,venable00,gomez:2003,smondyrev:1999,marrink01}
17	mmeineke	1001	This is due to the expense of the computer calculations involved when
18			performing these simulations. To properly hydrate a bilayer, one
19			typically needs 25 water molecules for every lipid, bringing the total
20			number of atoms simulated to roughly 8,000 for a system of 64 DPPC
21	mmeineke	1061	molecules. Added to the difficulty is the electrostatic nature of the
22	mmeineke	1001	phospholipid head groups and water, requiring the computationally
23			expensive Ewald sum or its slightly faster derivative particle mesh
24	mmeineke	1061	Ewald sum.\cite{nina:2002,norberg:2000,patra:2003} These factors all
25			limit the potential size and time lengths of bilayer simulations.
26	mmeineke	1001
27			Unfortunately, much of biological interest happens on time and length
28	mmeineke	1061	scales infeasible with current simulation. One such example is the
29			observance of a ripple phase ($P_{\beta^{\prime}}$) between the
30			$L_{\beta}$ and $L_{\alpha}$ phases of certain phospholipid
31			bilayers.\cite{katsaras00,sengupta00} These ripples are shown to
32	mmeineke	1001	have periodicity on the order of 100-200~$\mbox{\AA}$. A simulation on
33			this length scale would have approximately 1,300 lipid molecules with
34			an additional 25 water molecules per lipid to fully solvate the
35			bilayer. A simulation of this size is impractical with current
36			atomistic models.
37
38			Another class of simulations to consider, are those dealing with the
39			diffusion of molecules through a bilayer. Due to the fluid-like
40			properties of a lipid membrane, not all diffusion across the membrane
41			happens at pores. Some molecules of interest may incorporate
42			themselves directly into the membrane. Once here, they may possess an
43			appreciable waiting time (on the order of 10's to 100's of
44	mmeineke	1061	nanoseconds) within the bilayer. Such long simulation times are
45	mmeineke	1001	difficulty to obtain when integrating the system with atomistic
46			detail.
47
48			Addressing these issues, several schemes have been proposed. One
49	mmeineke	1061	approach by Goetz and Liposky\cite{goetz98} is to model the entire
50	mmeineke	1001	system as Lennard-Jones spheres. Phospholipids are represented by
51			chains of beads with the top most beads identified as the head
52			atoms. Polar and non-polar interactions are mimicked through
53			attractive and soft-repulsive potentials respectively. A similar
54	mmeineke	1061	model proposed by Marrinck \emph{et. al}.\cite{marrink04}~uses a
55	mmeineke	1001	similar technique for modeling polar and non-polar interactions with
56			Lennard-Jones spheres. However, they also include charges on the head
57			group spheres to mimic the electrostatic interactions of the
58			bilayer. While the solvent spheres are kept charge-neutral and
59			interact with the bilayer solely through an attractive Lennard-Jones
60			potential.
61
62			The model used in this investigation adds more information to the
63	mmeineke	1026	interactions than the previous two models, while still balancing the
64			need for simplifications over atomistic detail. The model uses
65			Lennard-Jones spheres for the head and tail groups of the
66	mmeineke	1061	phospholipids, allowing for the ability to scale the parameters to
67	mmeineke	1026	reflect various sized chain configurations while keeping the number of
68			interactions small. What sets this model apart, however, is the use
69	mmeineke	1061	of dipoles to represent the electrostatic nature of the
70	mmeineke	1026	phospholipids. The dipole electrostatic interaction is shorter range
71	mmeineke	1061	than Coulombic ($\frac{1}{r^3}$ versus $\frac{1}{r}$), eliminating the
72	mmeineke	1026	need for a costly Ewald sum.
73	mmeineke	1001
74	mmeineke	1026	Another key feature of this model, is the use of a dipolar water model
75	mmeineke	1061	to represent the solvent. The soft sticky dipole ({\sc ssd}) water
76			\cite{liu96:new_model} relies on the dipole for long range electrostatic
77			effects, but also contains a short range correction for hydrogen
78			bonding. In this way the systems in this research mimic the entropic
79			contribution to the hydrophobic effect due to hydrogen-bond network
80			deformation around a non-polar entity, \emph{i.e.}~the phospholipid.
81	mmeineke	1026
82			The following is an outline of this chapter.
83	mmeineke	1061	Sec.~\ref{lipidSec:Methods} is an introduction to the lipid model
84	mmeineke	1026	used in these simulations. As well as clarification about the water
85			model and integration techniques. The various simulation setups
86			explored in this research are outlined in
87			Sec.~\ref{lipidSec:ExpSetup}. Sec.~\ref{lipidSec:Results} and
88			Sec.~\ref{lipidSec:Discussion} give a summary of the results and
89			interpretation of those results respectively. Finally, the
90			conclusions of this chapter are presented in
91			Sec.~\ref{lipidSec:Conclusion}.
92
93	mmeineke	971	\section{\label{lipidSec:Methods}Methods}
94
95	mmeineke	1026	\subsection{\label{lipidSec:lipidModel}The Lipid Model}
96
97	mmeineke	971	\begin{figure}
98	mmeineke	1061	\centering
99			\includegraphics[width=\linewidth]{twoChainFig.eps}
100			\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.}
101	mmeineke	971	\label{lipidFig:lipidModel}
102			\end{figure}
103
104			The phospholipid model used in these simulations is based on the
105			design illustrated in Fig.~\ref{lipidFig:lipidModel}. The head group
106			of the phospholipid is replaced by a single Lennard-Jones sphere of
107	mmeineke	1061	diameter $\sigma_{\text{head}}$, with $\epsilon_{\text{head}}$ scaling
108			the well depth of its van der Walls interaction. This sphere also
109			contains a single dipole of magnitude $\|\boldsymbol{\mu}\|$, where
110			$\|\boldsymbol{\mu}\|$ can be varied to mimic the charge separation of a
111	mmeineke	971	given phospholipid head group. The atoms of the tail region are
112			modeled by unified atom beads. They are free of partial charges or
113			dipoles, containing only Lennard-Jones interaction sites at their
114			centers of mass. As with the head groups, their potentials can be
115	mmeineke	1061	scaled by $\sigma_{\text{tail}}$ and $\epsilon_{\text{tail}}$.
116	mmeineke	971
117			The long range interactions between lipids are given by the following:
118			\begin{equation}
119	mmeineke	1061	V_{\text{LJ}}(r_{ij}) =
120			4\epsilon_{ij} \biggl[
121			\biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
122			- \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
123			\biggr]
124	mmeineke	971	\label{lipidEq:LJpot}
125			\end{equation}
126			and
127			\begin{equation}
128	mmeineke	1061	V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
129			\boldsymbol{\Omega}_{j}) = \frac{\|\mu_i\|\|\mu_j\|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[
130			\boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j}
131			-
132			\frac{3(\boldsymbol{\hat{u}}_i \cdot \mathbf{r}_{ij}) %
133			(\boldsymbol{\hat{u}}_j \cdot \mathbf{r}_{ij}) }
134			{r^{2}_{ij}} \biggr]
135	mmeineke	971	\label{lipidEq:dipolePot}
136			\end{equation}
137			Where $V_{\text{LJ}}$ is the Lennard-Jones potential and
138			$V_{\text{dipole}}$ is the dipole-dipole potential. As previously
139			stated, $\sigma_{ij}$ and $\epsilon_{ij}$ are the Lennard-Jones
140			parameters which scale the length and depth of the interaction
141			respectively, and $r_{ij}$ is the distance between beads $i$ and $j$.
142			In $V_{\text{dipole}}$, $\mathbf{r}_{ij}$ is the vector starting at
143			bead$i$ and pointing towards bead $j$. Vectors $\mathbf{\Omega}_i$
144			and $\mathbf{\Omega}_j$ are the orientational degrees of freedom for
145			beads $i$ and $j$. $\|\mu_i\|$ is the magnitude of the dipole moment of
146			$i$, and $\boldsymbol{\hat{u}}_i$ is the standard unit orientation
147			vector of $\boldsymbol{\Omega}_i$.
148
149	mmeineke	1061	The model also allows for the bonded interactions bends, and torsions.
150			The bond between two beads on a chain is of fixed length, and is
151			maintained according to the {\sc rattle} algorithm.\cite{andersen83}
152	mmeineke	971	The bends are subject to a harmonic potential:
153			\begin{equation}
154	mmeineke	1061	V_{\text{bend}}(\theta) = k_{\theta}( \theta - \theta_0 )^2
155	mmeineke	971	\label{lipidEq:bendPot}
156			\end{equation}
157	mmeineke	1061	where $k_{\theta}$ scales the strength of the harmonic well, and
158			$\theta$ is the angle between bond vectors
159			(Fig.~\ref{lipidFig:lipidModel}). In addition, we have placed a
160			``ghost'' bend on the phospholipid head. The ghost bend adds a
161			potential to keep the dipole pointed along the bilayer surface, where
162			$theta$ is now the angle the dipole makes with respect to the {\sc
163			head}-$\text{{\sc ch}}_2$ bond vector. The torsion potential is given
164			by:
165	mmeineke	971	\begin{equation}
166	mmeineke	1061	V_{\text{torsion}}(\phi) =
167			k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
168	mmeineke	971	\label{lipidEq:torsionPot}
169			\end{equation}
170			Here, the parameters $k_0$, $k_1$, $k_2$, and $k_3$ fit the cosine
171			power series to the desired torsion potential surface, and $\phi$ is
172	mmeineke	1061	the angle the two end atoms have rotated about the middle bond
173			(Fig.:\ref{lipidFig:lipidModel}). Long range interactions such as the
174			Lennard-Jones potential are excluded for atom pairs involved in the
175			same bond, bend, or torsion. However, internal interactions not
176	mmeineke	971	directly involved in a bonded pair are calculated.
177
178	mmeineke	1026	All simulations presented here use a two chained lipid as pictured in
179	mmeineke	1061	Fig.~\ref{lipidFig:lipidModel}. The chains are both eight beads long,
180	mmeineke	1026	and their mass and Lennard Jones parameters are summarized in
181			Table~\ref{lipidTable:tcLJParams}. The magnitude of the dipole moment
182			for the head bead is 10.6~Debye, and the bend and torsion parameters
183	mmeineke	1061	are summarized in Table~\ref{lipidTable:tcBendParams} and
184			\ref{lipidTable:tcTorsionParams}.
185	mmeineke	971
186	mmeineke	1061	\begin{table}
187			\caption{The Lennard Jones Parameters for the two chain phospholipids.}
188			\label{lipidTable:tcLJParams}
189			\begin{center}
190			\begin{tabular}{\|l\|c\|c\|c\|}
191			\hline
192			& mass (amu) & $\sigma$($\mbox{\AA}$) & $\epsilon$ (kcal/mol) \\ \hline
193			{\sc head} & 72 & 4.0 & 0.185 \\ \hline
194			{\sc ch}\cite{Siepmann1998} & 13.02 & 4.0 & 0.0189 \\ \hline
195			$\text{{\sc ch}}_2$\cite{Siepmann1998} & 14.03 & 3.95 & 0.18 \\ \hline
196			$\text{{\sc ch}}_3$\cite{Siepmann1998} & 15.04 & 3.75 & 0.25 \\ \hline
197			{\sc ssd}\cite{liu96:new_model} & 18.03 & 3.051 & 0.152 \\ \hline
198			\end{tabular}
199			\end{center}
200			\end{table}
201	mmeineke	1026
202	mmeineke	1061	\begin{table}
203			\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}}
204			\label{lipidTable:tcBendParams}
205			\begin{center}
206			\begin{tabular}{\|l\|c\|c\|}
207			\hline
208			& $k_{\theta}$ ( kcal/($\text{mol deg}^2$) ) & $\theta_0$ ( deg ) \\ \hline
209			{\sc ghost}-{\sc head}-$\text{{\sc ch}}_2$ & 0.00177 & 129.78 \\ \hline
210			$x$-{\sc ch}-$y$ & 58.84 & 112.0 \\ \hline
211			$x$-$\text{{\sc ch}}_2$-$y$ & 58.84 & 114.0 \\ \hline
212			\end{tabular}
213			\end{center}
214			\end{table}
215
216			\begin{table}
217			\caption[Torsion Parameters for the two chain phospholipids]{Torsion Parameters for the two chain phospholipids. Alkane parameters based on TraPPE.\cite{Siepmann1998}}
218			\label{lipidTable:tcTorsionParams}
219			\begin{center}
220			\begin{tabular}{\|l\|c\|c\|c\|c\|}
221			\hline
222			All are in kcal/mol $\rightarrow$ & $k_3$ & $k_2$ & $k_1$ & $k_0$ \\ \hline
223			$x$-{\sc ch}-$y$-$z$ & 3.3254 & -0.4215 & -1.686 & 1.1661 \\ \hline
224			$x$-$\text{{\sc ch}}_2$-$\text{{\sc ch}}_2$-$y$ & 5.9602 & -0.568 & -3.802 & 2.1586 \\ \hline
225			\end{tabular}
226			\end{center}
227			\end{table}
228
229
230			\section{\label{lipidSec:furtherMethod}Further Methodology}
231
232	mmeineke	1026	As mentioned previously, the water model used throughout these
233	mmeineke	1061	simulations was the {\sc ssd} model of
234			Ichiye.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md} A
235			discussion of the model can be found in Sec.~\ref{oopseSec:SSD}. As
236			for the integration of the equations of motion, all simulations were
237			performed in an orthorhombic periodic box with a thermostat on
238			velocities, and an independent barostat on each Cartesian axis $x$,
239			$y$, and $z$. This is the $\text{NPT}_{xyz}$. ensemble described in
240			Sec.~\ref{oopseSec:Ensembles}.
241	mmeineke	1026
242
243			\subsection{\label{lipidSec:ExpSetup}Experimental Setup}
244
245			Two main starting configuration classes were used in this research:
246			random and ordered bilayers. The ordered bilayer starting
247			configurations were all started from an equilibrated bilayer at
248			300~K. The original configuration for the first 300~K run was
249			assembled by placing the phospholipids centers of mass on a planar
250			hexagonal lattice. The lipids were oriented with their long axis
251			perpendicular to the plane. The second leaf simply mirrored the first
252			leaf, and the appropriate number of waters were then added above and
253			below the bilayer.
254
255			The random configurations took more work to generate. To begin, a
256			test lipid was placed in a simulation box already containing water at
257			the intended density. The waters were then tested for overlap with
258			the lipid using a 5.0~$\mbox{\AA}$ buffer distance. This gave an
259			estimate for the number of waters each lipid would displace in a
260			simulation box. A target number of waters was then defined which
261			included the number of waters each lipid would displace, the number of
262			waters desired to solvate each lipid, and a fudge factor to pad the
263			initialization.
264
265			Next, a cubic simulation box was created that contained at least the
266			target number of waters in an FCC lattice (the lattice was for ease of
267			placement). What followed was a RSA simulation similar to those of
268			Chapt.~\ref{chapt:RSA}. The lipids were sequentially given a random
269			position and orientation within the box. If a lipid's position caused
270			atomic overlap with any previously adsorbed lipid, its position and
271			orientation were rejected, and a new random adsorption site was
272			attempted. The RSA simulation proceeded until all phospholipids had
273			been adsorbed. After adsorption, all water molecules with locations
274			that overlapped with the atomic coordinates of the lipids were
275			removed.
276
277			Finally, water molecules were removed one by one at random until the
278			desired number of waters per lipid was reached. The typical low final
279			density for these initial configurations was not a problem, as the box
280			would shrink to an appropriate size within the first 50~ps of a
281			simulation in the $\text{NPT}_{xyz}$ ensemble.
282
283			\subsection{\label{lipidSec:Configs}The simulation configurations}
284
285			Table ~\ref{lipidTable:simNames} summarizes the names and important
286			details of the simulations. The B set of simulations were all started
287	mmeineke	1061	in an ordered bilayer and observed over a period of 10~ns. Simulation
288	mmeineke	1026	RL was integrated for approximately 20~ns starting from a random
289			configuration as an example of spontaneous bilayer aggregation.
290			Lastly, simulation RH was also started from a random configuration,
291			but with a lesser water content and higher temperature to show the
292	mmeineke	1061	spontaneous aggregation of an inverted hexagonal lamellar phase.