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{Lindahl:2000} to the
spontaneous aggregation of DPPC molecules into fluid phase
($L_{\alpha}$ bilayers. \cite{Marrinck:2001} With the exception of a
few ambitious
simulations,\cite{Marrinch:2001b,Marrinck:2002,Lindahl:2000} most
investigations are limited to 64 to 256
phospholipids.\cite{Lindal:2000,Sum:2003,Venable:2000,Gomez:2003,Smondyrev:1999,Marrinck:2001a}
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 difficluty 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 lenghts of bilayer simulations.

Unfortunately, much of biological interest happens on time and length
scales unfeasible with current simulation. One such example is the
observance of a ripple phase ($P_{\beta'}$) between the $L_{\beta}$
and $L_{\alpha}$ phases of certain phospholipid
bilayers.\cite{Katsaras:2000,Sengupta:2000} 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{Goetz:1998} 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{Marrinck:2004}~ 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
phopholipids, 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 electrosttaic 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 ({\scssd})
water \cite{Liu:1996a} relies on the dipole for long range
electrostatic effects, butalso 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{lipoidSec: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}

\caption{Schematic diagram of the single chain phospholipid model}

\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 $fix$, with $fix$ scaling the well depth of its van der Walls
interaction.  This sphere also contains a single dipole of magnitude
$fix$, where $fix$ 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 $fix$ and $fix$.

The long range interactions between lipids are given by the following:
\begin{equation}
EQ Here
\label{lipidEq:LJpot}
\end{equation}
and
\begin{equation}
EQ Here
\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 of bonds, bends, and
torsions.  The bonds between two beads on a chain are of fixed length,
and are maintained according to the {\sc rattle} algorithm. \cite{fix}
The bends are subject to a harmonic potential:
\begin{equation}
eq here
\label{lipidEq:bendPot}
\end{equation}
where $fix$ scales the strength of the harmonic well, and $fix$ is the
angle between bond vectors $fix$ and $fix$.  The torsion potential is
given by:
\begin{equation}
eq here
\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 between bondvectors $fix$ and $fix$ along the vector $fix$
(see Fig.:\ref{lipidFig:lipidModel}).  Long range interactions such as
the Lennard-Jones potential are excluded for bead 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:twochain}.  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:teBTParams}.

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

As mentioned previously, the water model used throughout these
simulations was the {\scssd} model of
Ichiye.\cite{liu:1996a,Liu:1996b,Chandra:1999} 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. Simulution
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:	1026
Committed:	Thu Feb 5 00:37:36 2004 UTC (20 years, 7 months ago) by mmeineke
Content type:	application/x-tex
File size:	11584 byte(s)
Log Message:	typed up the lipid section
#	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			dipalmitoylphosphatidylcholine (DPPC), \cite{Lindahl:2000} to the
11			spontaneous aggregation of DPPC molecules into fluid phase
12			($L_{\alpha}$ bilayers. \cite{Marrinck:2001} With the exception of a
13			few ambitious
14			simulations,\cite{Marrinch:2001b,Marrinck:2002,Lindahl:2000} most
15			investigations are limited to 64 to 256
16			phospholipids.\cite{Lindal:2000,Sum:2003,Venable:2000,Gomez:2003,Smondyrev:1999,Marrinck:2001a}
17			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			molecules. Added to the difficluty is the electrostatic nature of the
22			phospholipid head groups and water, requiring the computationally
23			expensive Ewald sum or its slightly faster derivative particle mesh
24			Ewald sum.\cite{Nina:2002,Norberg:2000,Patra:2003} These factors all
25			limit the potential size and time lenghts of bilayer simulations.
26
27			Unfortunately, much of biological interest happens on time and length
28			scales unfeasible with current simulation. One such example is the
29			observance of a ripple phase ($P_{\beta'}$) between the $L_{\beta}$
30			and $L_{\alpha}$ phases of certain phospholipid
31			bilayers.\cite{Katsaras:2000,Sengupta:2000} These ripples are shown to
32			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			nanoseconds) within the bilayer. Such long simulation times are
45			difficulty to obtain when integrating the system with atomistic
46			detail.
47
48			Addressing these issues, several schemes have been proposed. One
49			approach by Goetz and Liposky\cite{Goetz:1998} is to model the entire
50			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			model proposed by Marrinck \emph{et. al.}\cite{Marrinck:2004}~ uses a
55			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			phopholipids, allowing for the ability to scale the parameters to
67			reflect various sized chain configurations while keeping the number of
68			interactions small. What sets this model apart, however, is the use
69			of dipoles to represent the electrosttaic nature of the
70			phospholipids. The dipole electrostatic interaction is shorter range
71			than coulombic ($\frac{1}{r^3}$ versus $\frac{1}{r}$), eliminating the
72			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			to represent the solvent. The soft sticky dipole ({\scssd})
76			water \cite{Liu:1996a} relies on the dipole for long range
77			electrostatic effects, butalso contains a short range correction for
78			hydrogen bonding. In this way the systems in this research mimic the
79			entropic contribution to the hydrophobic effect due to hydrogen-bond
80			network deformation around a non-polar entity, \emph{i.e.}~ the
81			phospholipid.
82
83			The following is an outline of this chapter.
84			Sec.~\ref{lipoidSec:Methods} is an introduction to the lipid model
85			used in these simulations. As well as clarification about the water
86			model and integration techniques. The various simulation setups
87			explored in this research are outlined in
88			Sec.~\ref{lipidSec:ExpSetup}. Sec.~\ref{lipidSec:Results} and
89			Sec.~\ref{lipidSec:Discussion} give a summary of the results and
90			interpretation of those results respectively. Finally, the
91			conclusions of this chapter are presented in
92			Sec.~\ref{lipidSec:Conclusion}.
93
94	mmeineke	971	\section{\label{lipidSec:Methods}Methods}
95
96
97	mmeineke	1026
98			\subsection{\label{lipidSec:lipidModel}The Lipid Model}
99
100	mmeineke	971	\begin{figure}
101
102			\caption{Schematic diagram of the single chain phospholipid model}
103
104			\label{lipidFig:lipidModel}
105
106			\end{figure}
107
108			The phospholipid model used in these simulations is based on the
109			design illustrated in Fig.~\ref{lipidFig:lipidModel}. The head group
110			of the phospholipid is replaced by a single Lennard-Jones sphere of
111			diameter $fix$, with $fix$ scaling the well depth of its van der Walls
112			interaction. This sphere also contains a single dipole of magnitude
113			$fix$, where $fix$ can be varied to mimic the charge separation of a
114			given phospholipid head group. The atoms of the tail region are
115			modeled by unified atom beads. They are free of partial charges or
116			dipoles, containing only Lennard-Jones interaction sites at their
117			centers of mass. As with the head groups, their potentials can be
118			scaled by $fix$ and $fix$.
119
120			The long range interactions between lipids are given by the following:
121			\begin{equation}
122			EQ Here
123			\label{lipidEq:LJpot}
124			\end{equation}
125			and
126			\begin{equation}
127			EQ Here
128			\label{lipidEq:dipolePot}
129			\end{equation}
130			Where $V_{\text{LJ}}$ is the Lennard-Jones potential and
131			$V_{\text{dipole}}$ is the dipole-dipole potential. As previously
132			stated, $\sigma_{ij}$ and $\epsilon_{ij}$ are the Lennard-Jones
133			parameters which scale the length and depth of the interaction
134			respectively, and $r_{ij}$ is the distance between beads $i$ and $j$.
135			In $V_{\text{dipole}}$, $\mathbf{r}_{ij}$ is the vector starting at
136			bead$i$ and pointing towards bead $j$. Vectors $\mathbf{\Omega}_i$
137			and $\mathbf{\Omega}_j$ are the orientational degrees of freedom for
138			beads $i$ and $j$. $\|\mu_i\|$ is the magnitude of the dipole moment of
139			$i$, and $\boldsymbol{\hat{u}}_i$ is the standard unit orientation
140			vector of $\boldsymbol{\Omega}_i$.
141
142			The model also allows for the bonded interactions of bonds, bends, and
143			torsions. The bonds between two beads on a chain are of fixed length,
144			and are maintained according to the {\sc rattle} algorithm. \cite{fix}
145			The bends are subject to a harmonic potential:
146			\begin{equation}
147			eq here
148			\label{lipidEq:bendPot}
149			\end{equation}
150			where $fix$ scales the strength of the harmonic well, and $fix$ is the
151			angle between bond vectors $fix$ and $fix$. The torsion potential is
152			given by:
153			\begin{equation}
154			eq here
155			\label{lipidEq:torsionPot}
156			\end{equation}
157			Here, the parameters $k_0$, $k_1$, $k_2$, and $k_3$ fit the cosine
158			power series to the desired torsion potential surface, and $\phi$ is
159			the angle between bondvectors $fix$ and $fix$ along the vector $fix$
160			(see Fig.:\ref{lipidFig:lipidModel}). Long range interactions such as
161			the Lennard-Jones potential are excluded for bead pairs involved in
162			the same bond, bend, or torsion. However, internal interactions not
163			directly involved in a bonded pair are calculated.
164
165	mmeineke	1026	All simulations presented here use a two chained lipid as pictured in
166			Fig.~\ref{lipidFig:twochain}. The chains are both eight beads long,
167			and their mass and Lennard Jones parameters are summarized in
168			Table~\ref{lipidTable:tcLJParams}. The magnitude of the dipole moment
169			for the head bead is 10.6~Debye, and the bend and torsion parameters
170			are summarized in Table~\ref{lipidTable:teBTParams}.
171	mmeineke	971
172	mmeineke	1026	\section{label{lipidSec:furtherMethod}Further Methodology}
173
174			As mentioned previously, the water model used throughout these
175			simulations was the {\scssd} model of
176			Ichiye.\cite{liu:1996a,Liu:1996b,Chandra:1999} A discussion of the
177			model can be found in Sec.~\ref{oopseSec:SSD}. As for the integration
178			of the equations of motion, all simulations were performed in an
179			orthorhombic periodic box with a thermostat on velocities, and an
180			independent barostat on each cartesian axis $x$, $y$, and $z$. This
181			is the $\text{NPT}_{xyz}$. ensemble described in Sec.~\ref{oopseSec:Ensembles}.
182
183
184			\subsection{\label{lipidSec:ExpSetup}Experimental Setup}
185
186			Two main starting configuration classes were used in this research:
187			random and ordered bilayers. The ordered bilayer starting
188			configurations were all started from an equilibrated bilayer at
189			300~K. The original configuration for the first 300~K run was
190			assembled by placing the phospholipids centers of mass on a planar
191			hexagonal lattice. The lipids were oriented with their long axis
192			perpendicular to the plane. The second leaf simply mirrored the first
193			leaf, and the appropriate number of waters were then added above and
194			below the bilayer.
195
196			The random configurations took more work to generate. To begin, a
197			test lipid was placed in a simulation box already containing water at
198			the intended density. The waters were then tested for overlap with
199			the lipid using a 5.0~$\mbox{\AA}$ buffer distance. This gave an
200			estimate for the number of waters each lipid would displace in a
201			simulation box. A target number of waters was then defined which
202			included the number of waters each lipid would displace, the number of
203			waters desired to solvate each lipid, and a fudge factor to pad the
204			initialization.
205
206			Next, a cubic simulation box was created that contained at least the
207			target number of waters in an FCC lattice (the lattice was for ease of
208			placement). What followed was a RSA simulation similar to those of
209			Chapt.~\ref{chapt:RSA}. The lipids were sequentially given a random
210			position and orientation within the box. If a lipid's position caused
211			atomic overlap with any previously adsorbed lipid, its position and
212			orientation were rejected, and a new random adsorption site was
213			attempted. The RSA simulation proceeded until all phospholipids had
214			been adsorbed. After adsorption, all water molecules with locations
215			that overlapped with the atomic coordinates of the lipids were
216			removed.
217
218			Finally, water molecules were removed one by one at random until the
219			desired number of waters per lipid was reached. The typical low final
220			density for these initial configurations was not a problem, as the box
221			would shrink to an appropriate size within the first 50~ps of a
222			simulation in the $\text{NPT}_{xyz}$ ensemble.
223
224			\subsection{\label{lipidSec:Configs}The simulation configurations}
225
226			Table ~\ref{lipidTable:simNames} summarizes the names and important
227			details of the simulations. The B set of simulations were all started
228			in an ordered bilayer and observed over a period of 10~ns. Simulution
229			RL was integrated for approximately 20~ns starting from a random
230			configuration as an example of spontaneous bilayer aggregation.
231			Lastly, simulation RH was also started from a random configuration,
232			but with a lesser water content and higher temperature to show the
233			spontaneous aggregation of an inverted hexagonal lamellar phase.