ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/matt_papers/canidacy_paper/canidacy_paper.tex
(Generate patch)

Comparing trunk/matt_papers/canidacy_paper/canidacy_paper.tex (file contents):
Revision 105 by mmeineke, Tue Aug 27 20:54:03 2002 UTC vs.
Revision 106 by mmeineke, Wed Sep 11 03:26:01 2002 UTC

# Line 1 | Line 1
1   \documentclass[11pt]{article}
2  
3   \usepackage{graphicx}
4 + \usepackage{color}
5   \usepackage{floatflt}
6   \usepackage{amsmath}
7   \usepackage{amssymb}
# Line 33 | Line 34 | Notre Dame, Indiana 46556}
34  
35   \section{Background and Research Goals}
36  
37 + Simulations of phospholipid bilayers are, by necessity, quite
38 + complex. The lipid molecules are large molecules containing many
39 + atoms,and the head group of the lipid will typically contain charge
40 + separated ions which set up a large dipole within the molecule. Adding
41 + to the complexity are the number of water molecules needed to properly
42 + solvate the lipid bilayer. Because of these factors, many current
43 + simulations are limited in both length and time scale due to to the
44 + sheer number of calculations performed at every time step and the
45 + lifetime of the researcher. A typical
46 + simulation\cite{saiz02,lindahl00,venable00,Marrink01} will have around
47 + 64 phospholipids forming a bilayer approximately 40~$\mbox{\AA}$ by
48 + 50~$\mbox{\AA}$ with roughly 25 waters for every lipid. This means
49 + there are on the order of 8,000 atoms needed to simulate these systems
50 + and the trajectories in turn are integrated for times up to 10 ns.
51 +
52 + These limitations make it difficult to study certain biologically
53 + interesting phenomena that don't fit within the short time and length
54 + scale requirements. One such phenomena is the existence of the ripple
55 + phase ($P_{\beta'}$) of the bilayer between the gel phase
56 + ($L_{\beta'}$) and the fluid phase ($L_{\alpha}$). The $P_{\beta'}$
57 + phase has been shown to have a ripple period of
58 + 100-200~$\mbox{\AA}$.\cite{katsaras00,sengupta00} A simulation of this
59 + length scale would require approximately 1,300 lipid molecules and
60 + roughly 25 waters for every lipid to fully solvate the bilayer. With
61 + the large number of atoms involved in a simulation of this magnitude,
62 + steps \emph{must} be taken to simplify the system to the point where
63 + these numbers are reasonable.
64 +
65 + Another system of interest would be drug molecule diffusion through
66 + the membrane. Due to the fluid like properties of a lipid membrane,
67 + not all diffusion takes place at membrane channels. It is of interest
68 + to study certain molecules that may incorporate themselves directly
69 + into the membrane. These molecules may then have an appreciable
70 + waiting time (on the order of nanoseconds) within the
71 + bilayer. Simulation of such a long time scale again requires
72 + simplification of the system in order to lower the number of
73 + calculations needed at each time step.
74 +
75 +
76   \section{Methodology}
77  
78   \subsection{Length and Time Scale Simplifications}
79  
80 < The length scale simplifications are aimed at increaseing the number
81 < of molecules simulated without drastically increasing the
80 > The length scale simplifications are aimed at increasing the number of
81 > molecules that can be simulated without drastically increasing the
82   computational cost of the system. This is done by a combination of
83   substituting less expensive interactions for expensive ones and
84   decreasing the number of interaction sites per molecule. Namely,
85 < charge distributions are replaced with dipoles, and unified atoms are
86 < used in place of water and phospholipid head groups.
85 > point charge distributions are replaced with dipoles, and unified atoms are
86 > used in place of water, phospholipid head groups, and alkyl groups.
87  
88   The replacement of charge distributions with dipoles allows us to
89 < replace an interaction that has a relatively long range, $\frac{1}{r}$
90 < for the charge charge potential, with that of a relitively short
91 < range, $\frac{1}{r^{3}}$ for dipole - dipole potentials
92 < (Equations~\ref{eq:dipolePot} and \ref{eq:chargePot}). This allows us
93 < to use computaional simplifications algorithms such as Verlet neighbor
94 < lists,\cite{allen87:csl} which gives computaional scaling by $N$. This
95 < is in comparison to the Ewald sum\cite{leach01:mm} needed to compute
96 < the charge - charge interactions which scales at best by $N
57 < \ln N$.
89 > replace an interaction that has a relatively long range,
90 > $(\frac{1}{r})$, for the coulomb potential, with that of a relatively
91 > short range, $(\frac{1}{r^{3}})$, for dipole - dipole
92 > potentials. Combined with a computational simplification algorithm
93 > such as a Verlet neighbor list,\cite{allen87:csl} this should give
94 > computational scaling by $N$. This is in comparison to the Ewald
95 > sum\cite{leach01:mm} needed to compute the coulomb interactions, which
96 > scales at best by $N \ln N$.
97  
98 < \begin{equation}
60 < V^{\text{dp}}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
61 <        \boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
62 <        \frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
63 <        -
64 <        \frac{3(\boldsymbol{\mu} \cdot \mathbf{r}_{ij}) %
65 <                (\boldsymbol{\mu} \cdot \mathbf{r}_{ij}) }{r^{5}_{ij}} \biggr]
66 < \label{eq:dipolePot}
67 < \end{equation}
68 <
69 < \begin{equation}
70 < V^{\text{ch}}_{ij}(\mathbf{r}_{ij}) = \frac{q_{i}q_{j}}%
71 <        {4\pi\epsilon_{0} r_{ij}}
72 < \label{eq:chargePot}
73 < \end{equation}
74 <
75 < The second step taken to simplify the number of calculationsis to
98 > The second step taken to simplify the number of calculations is to
99   incorporate unified models for groups of atoms. In the case of water,
100   we use the soft sticky dipole (SSD) model developed by
101   Ichiye\cite{Liu96} (Section~\ref{sec:ssdModel}). For the phospholipids, a
102   unified head atom with a dipole will replace the atoms in the head
103   group, while unified $\text{CH}_2$ and $\text{CH}_3$ atoms will
104 < replace the alkanes in the tails (Section~\ref{sec:lipidModel}).
104 > replace the alkyl units in the tails (Section~\ref{sec:lipidModel}).
105  
106 < The time scale simplifications are taken so that the simulation can
107 < take long time steps. By incresing the time steps taken by the
108 < simulation, we are able to integrate the simulation trajectory with
109 < fewer calculations. However, care must be taken to conserve the energy
110 < of the simulation. This is a constraint placed upon the system by
111 < simulating in the microcanonical ensemble. In practice, this means
112 < taking steps small enough to resolve all motion in the system without
113 < accidently moving an object too far along a repulsive energy surface
114 < before it feels the affect of the surface.
106 > The time scale simplifications are introduced so that we can take
107 > longer time steps. By increasing the size of the time steps taken by
108 > the simulation, we are able to integrate the simulation trajectory
109 > with fewer calculations. However, care must be taken that any
110 > simplifications used, still conserve the total energy of the
111 > simulation. In practice, this means taking steps small enough to
112 > resolve all motion in the system without accidently moving an object
113 > too far along a repulsive energy surface before it feels the effect of
114 > the surface.
115  
116   In our simulation we have chosen to constrain all bonds to be of fixed
117   length. This means the bonds are no longer allowed to vibrate about
118 < their equilibrium positions, typically the fastest periodical motion
119 < in a dynamics simulation. By taking this action, we are able to take
120 < time steps of 3 fs while still maintaining constant energy. This is in
121 < contrast to the 1 fs time step typically needed to conserve energy when
122 < bonds lengths are allowed to oscillate.
118 > their equilibrium positions. Bond vibrations are typically the fastest
119 > periodic motion in a dynamics simulation. By taking this action, we
120 > are able to take time steps of 3 fs while still maintaining constant
121 > energy. This is in contrast to the 1 fs time step typically needed to
122 > conserve energy when bonds lengths are allowed to oscillate.
123  
124   \subsection{The Soft Sticky Water Model}
125   \label{sec:ssdModel}
126  
127 < \begin{floatingfigure}{55mm}
128 < \includegraphics[width=45mm]{ssd.epsi}
129 < \caption{The SSD model with the oxygen and hydrogen atoms drawn in for reference. \vspace{5mm}}
130 < % The dipole magnitude is 2.35 D and the Lennard-Jones parameters are $\sigma = 3.051 \mbox{\AA}$ and $\epsilon = 0.152$ kcal/mol.}
131 < \label{fig:ssdModel}
109 < \end{floatingfigure}
127 > %\begin{floatingfigure}{55mm}
128 > %\includegraphics[width=45mm]{ssd.epsi}
129 > %\caption{The SSD model with the oxygen and hydrogen atoms drawn in for reference. \vspace{5mm}}
130 > %\label{fig:ssdModel}
131 > %\end{floatingfigure}
132  
133 < The water model used in our simulations is a modified soft Stockmayer
134 < sphere model. Like the soft Stockmayer sphere, the SSD
133 > The water model used in our simulations is a modified soft
134 > Stockmayer-sphere model.\cite{stevens95} Like the Stockmayer-sphere, the SSD
135   model\cite{Liu96} consists of a Lennard-Jones interaction site and a
136   dipole both located at the water's center of mass (Figure
137   \ref{fig:ssdModel}). However, the SSD model extends this by adding a
138   tetrahedral potential to correct for hydrogen bonding.
139  
140 < This SSD water's motion is then governed by the following potential:
140 > The SSD water potential is then given by the following equation:
141   \begin{equation}
142 < V_{\text{ssd}} = V_{\text{LJ}}(r_{i\!j}) + V_{\text{dp}}(\mathbf{r}_{i\!j},
142 > V_{\text{SSD}} = V_{\text{LJ}}(r_{i\!j}) + V_{\text{dp}}(\mathbf{r}_{i\!j},
143          \boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
144          + V_{\text{sp}}(\mathbf{r}_{i\!j},\boldsymbol{\Omega}_{i},
145          \boldsymbol{\Omega}_{j})
146   \label{eq:ssdTotPot}
147   \end{equation}
148 < $V_{\text{LJ}}$ is the Lennard-Jones potential with $\sigma_{\text{w}}
149 < = 3.051 \mbox{ \AA}$ and $\epsilon_{\text{w}} = 0.152\text{
150 < kcal/mol}$. $V_{\text{dp}}$ is the dipole potential with
151 < $|\mu_{\text{w}}| = 2.35\text{ D}$.
148 > $V_{\text{LJ}}$ is the Lennard-Jones potential:
149 > \begin{equation}
150 > V_{\text{LJ}} =
151 >        4\epsilon_{ij} \biggl[
152 >        \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
153 >        - \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
154 >        \biggr]
155 > \label{eq:lennardJonesPot}
156 > \end{equation}
157 > where $r_{ij}$ is the distance between two $ij$ pairs, $\sigma_{ij}$
158 > scales the length of the iteraction, and $\epsilon_{ij}$ scales the
159 > energy of the potential. For SSD, $\sigma_{\text{SSD}} = 3.051 \mbox{
160 > \AA}$ and $\epsilon_{\text{SSD}} = 0.152\text{ kcal/mol}$.
161 > $V_{\text{dp}}$ is the dipole potential:
162 > \begin{equation}
163 > V_{\text{dp}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
164 >        \boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
165 >        \frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
166 >        -
167 >        \frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
168 >                (\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
169 >                {r^{5}_{ij}} \biggr]
170 > \label{eq:dipolePot}
171 > \end{equation}
172 > where $\mathbf{r}_{ij}$ is the vector between $i$ and $j$,
173 > $\boldsymbol{\Omega}$ is the orientation of the species, and
174 > $\boldsymbol{\mu}$ is the dipole vector. The SSD model specifies a dipole
175 > magnitude of 2.35~D for water.
176  
177   The hydrogen bonding of the model is governed by the $V_{\text{sp}}$
178 < term of the potentail. Its form is as follows:
178 > term of the potential. Its form is as follows:
179   \begin{equation}
180   V_{\text{sp}}(\mathbf{r}_{i\!j},\boldsymbol{\Omega}_{i},
181          \boldsymbol{\Omega}_{j}) =
# Line 140 | Line 186 | Where $v^\circ$ is responsible for scaling the strengt
186                  \boldsymbol{\Omega}_{j})]
187   \label{eq:spPot}
188   \end{equation}
189 < Where $v^\circ$ is responsible for scaling the strength of the
144 < interaction.
189 > Where $v^\circ$ scales strength of the interaction.
190   $w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
191   and
192   $w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
# Line 151 | Line 196 | w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsy
196   w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) =
197          \sin\theta_{ij} \sin 2\theta_{ij} \cos 2\phi_{ij}
198          + \sin \theta_{ji} \sin 2\theta_{ji} \cos 2\phi_{ji}
199 < \label{eq:apPot2}
199 > \label{eq:spPot2}
200   \end{equation}
201   and
202   \begin{equation}
# Line 162 | Line 207 | The correction
207   \end{split}
208   \label{eq:spCorrection}
209   \end{equation}
210 < The correction
210 > The angles $\theta_{ij}$ and $\phi_{ij}$ are defined by the spherical
211 > polar coordinates of the position of sphere $j$ in the reference frame
212 > fixed on sphere $i$ with the z-axis aligned with the dipole moment.
213 > The correction
214   $w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
215   is needed because
216   $w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
217 < vanishes when $\theta_{ij}$ is $0^\circ$ or $180^\circ$. The angles
170 < $\theta_{ij}$ and $\phi_{ij}$ are defined by the spherical polar
171 < coordinates of the position of sphere $j$ in the reference frame fixed
172 < on sphere $i$ with the z-axis alligned with the dipole moment.
217 > vanishes when $\theta_{ij}$ is $0^\circ$ or $180^\circ$.
218  
219 < Finaly, the sticky potentail is scaled by a cutoff function,
219 > Finally, the sticky potential is scaled by a cutoff function,
220   $s(r_{ij})$ that scales smoothly between 0 and 1. It is represented
221   by:
222   \begin{equation}
# Line 186 | Line 231 | s(r_{ij}) =
231   \label{eq:spCutoff}
232   \end{equation}
233  
234 + Despite the apparent complexity of Equation \ref{eq:spPot}, the SSD
235 + model is still computationaly inexpensive. This is due to Equation
236 + \ref{eq:spCutoff}. With $r_{L}$ being 2.75~$\mbox\AA}$ and $r_{U}$
237 + being equal to either 3.35~$\mbox{\AA}$ for $s(r_{ij})$ or
238 + 4.0~$\mbox{\AA}$ for $s'(r_{ij})$, the sticky potential is only active
239 + over an extremly short range, and then only with other SSD
240 + molecules. Therefore, it's predominant interaction is through it's
241 + point dipole and Lennard-Jones sphere.
242  
243   \subsection{The Phospholipid Model}
244   \label{sec:lipidModel}
# Line 201 | Line 254 | are unifed $\text{CH}_2$ and $\text{CH}_3$ atoms and a
254   lipids. The Head group of the phospholipid is replaced by a single
255   Lennard-Jones sphere with a freely oriented dipole placed at it's
256   center. The magnitude of it's dipole moment is 20.6 D. The tail atoms
257 < are unifed $\text{CH}_2$ and $\text{CH}_3$ atoms and are also modeled
257 > are unified $\text{CH}_2$ and $\text{CH}_3$ atoms and are also modeled
258   as Lennard-Jones spheres. The total potential for the lipid is
259   represented by Equation \ref{eq:lipidModelPot}.
260  
# Line 216 | Line 269 | The non bonded interactions, $V_{\text{LJ}}$ and $V_{\
269   \label{eq:lipidModelPot}
270   \end{equation}
271  
272 < The non bonded interactions, $V_{\text{LJ}}$ and $V_{\text{dp}}$, are
272 > The non-bonded interactions, $V_{\text{LJ}}$ and $V_{\text{dp}}$, are
273   the Lennard-Jones and dipole-dipole interactions respectively. For the
274   non-bonded potentials, only the bend and the torsional potentials are
275   calculated. The bond potential is not calculated, and the bond lengths
# Line 229 | Line 282 | V_{\text{tors.}}(\phi_{ijkl})= c_1[1+\cos\phi]
282   Where $k_{\theta}$ sets the stiffness of the bend potential, and $\theta_0$
283   sets the equilibrium bend angle. The torsion potential was given by:
284   \begin{equation}
285 < V_{\text{tors.}}(\phi_{ijkl})= c_1[1+\cos\phi]
286 <        + c_2 [1 - \cos(2\phi)] + c_3[1 + \cos(3\phi)]
285 > V_{\text{tors.}}(\phi_{ijkl})= c_1[1+\cos\phi_{ijkl}]
286 >        + c_2 [1 - \cos(2\phi_{ijkl})] + c_3[1 + \cos(3\phi_{ijkl})]
287   \label{eq:torsPot}
288   \end{equation}
289   All parameters for bonded and non-bonded potentials in the tail atoms
290   were taken from TraPPE.\cite{Siepmann1998} The bonded interactions for
291   the head atom were also taken from TraPPE, however it's dipole moment
292 < and mass were based on the properties of ``DMPC?'''s head group. The
293 < Lennard-Jones parameter for the head group was chosen such that it was
294 < roughly twice the size of a $\text{CH}_3$ atom, and it's well depth
295 < was set to be aproximately equal to that of $\text{CH}_3$.
292 > and mass were based on the properties of the phosphatidylcholine head
293 > group. The Lennard-Jones parameter for the head group was chosen such
294 > that it was roughly twice the size of a $\text{CH}_3$ atom, and it's
295 > well depth was set to be approximately equal to that of $\text{CH}_3$.
296  
297   \section{Initial Simulation: 25 lipids in water}
298   \label{sec:5x5}
299  
300 < \subsection{Starting Configurtion and Parameters}
300 > \subsection{Starting Configuration and Parameters}
301   \label{sec:5x5Start}
302  
250 \begin{floatingfigure}{85mm}
251 \includegraphics[width=75mm]{5x5-initial.eps}
252 \caption{A snapshot of the initial configuration of the 25 lipid simulation.}
253 \label{fig:5x5Start}
254 \end{floatingfigure}
255
303   Our first simulation was an array of 25 single chained lipids in a sea
304   of water (Figure \ref{fig:5x5Start}). The total number of water
305   molecules was 1386, giving a final of water concentration of 70\%
306   wt. The simulation box measured 34.5~$\mbox{\AA}$ x 39.4~$\mbox{\AA}$
307 < x 39.4~$\mbox{\AA}$ with periodic boundry connditions invoked. The
307 > x 39.4~$\mbox{\AA}$ with periodic boundary conditions invoked. The
308   system was simulated in the micro-canonical (NVE) ensemble with an
309   average temperature of 300~K.
310  
264
265 yada
266
267 yada
268
269 yada
270 yya
271 d
272
273 uada
274
275 adsd
276
277 asfa
278
279 asfdads
280
311   \subsection{Results}
312   \label{sec:5x5Results}
313  
284
314   Figure \ref{fig:5x5Final} shows a snapshot of the system at
315 < 3.6~ns. Here the bylayer has formed with the lipid chains tilted to
316 < the bylayer normal.
315 > 3.6~ns. Here it is exciting to note that the system has spontaneously
316 > self assembled into a bilayer. Discussion of the length scales of the
317 > bilayer will follow in this section. However, it is interesting to
318 > note several qualitative properties of the system revealed by this
319 > snapshot. First, is how the tail chains are tilted to the bilayer
320 > normal. This is usually indicative of the gel ($L_{\beta'}$)
321 > phase. Likely, the box size is too small for the bilayer to relax to
322 > the fluid ($P_{\alpha}$ phase. Showing the need for a constant
323 > pressure simulation versus the constant volume of the current
324 > ensemble.
325  
326 < yada
326 > The trajectory of the simulation was analyzed using the pair-wise
327 > radial distribution function, $g(r)$, which has the form:
328 > \begin{equation}
329 > g(r) = \frac{V}{N_{\text{pairs}}}\langle \sum_{i} \sum_{j \neq i}
330 >        \delta(|\mathbf{r} - \mathbf{r}_{ij}|) \rangle
331 > \label{eq:gofr}
332 > \end{equation}
333 > Equation \ref{eq:gofr} gives us information about the spacing of two
334 > species as a function of radius. Essentially, if the observer were
335 > located at atom $i$ and were looking out in all directions, $g(r)$
336 > shows the relative density of atom $j$ at any given radius, $r$,
337 > normalized by the expected density of atom $j$ at $r$. In a
338 > homogeneously mixed fluid, $g(r)$ will approach 1 at large $r$, as a
339 > fluid contains no long range structure to contribute peaks in the
340 > number density.
341  
342 < yada
342 > For the species containing dipoles, a second pair wise distribution
343 > function was used, $g_{\gamma}(r)$. It is of the form:
344 > \begin{equation}
345 > g_{\gamma}(r) = foobar
346 > \label{eq:gammaofr}
347 > \end{equation}
348 > Where $\gamma_{ij}$ is the angle between the dipole of atom $j$ with
349 > respect to the dipole of atom $i$. This correlation will vary between
350 > $+1$ and $-1$ when the two dipoles are perfectly aligned and
351 > anti-aligned respectively. This then gives us information about how
352 > directional species are aligned with each other as a function of
353 > distance.
354  
355 < yada
355 > Figure \ref{fig:5x5HHCorr} shows the two self correlation functions
356 > for the Head groups of the lipids. The first peak at 4.03~$\mbox{\AA}$ is the
357 > nearest neighbor separation of the heads of two lipids.
358  
295 yada
359  
297 %\begin{floatingfigure}{85mm}
298 %\includegraphics[angle=-90,width=75mm]{5x5-3.6ns.epsi}
299 %\caption{A snapshot of the system at 3.6~ns.\vspace{5mm}}
300 %\label{fig:5x5Final}
301 %\end{floatingfigure}
360  
361  
304
305
362   \section{Second Simulation: 50 randomly oriented lipids in water}
363  
364   the second simulation
365  
310 \section{Preliminary Results}
311
312 \section{Discussion}
313
366   \section{Future Directions}
367  
368  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines