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 (21 years, 2 months ago) by mmeineke
Content type:	application/x-tex
File size:	11584 byte(s)
Log Message:	typed up the lipid section
#	Content
1
2
3	\chapter{\label{chapt:lipid}Phospholipid Simulations}
4
5	\section{\label{lipidSec:Intro}Introduction}
6
7	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	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
74	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	\section{\label{lipidSec:Methods}Methods}
95
96
97
98	\subsection{\label{lipidSec:lipidModel}The Lipid Model}
99
100	\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	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
172	\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.