ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/matt_papers/canidacy_paper/canidacy_paper.tex
Revision: 111
Committed: Wed Sep 18 15:53:30 2002 UTC (21 years, 9 months ago) by mmeineke
Content type: application/x-tex
File size: 23889 byte(s)
Log Message:
the final copy went out this morning. YAYYYYYYYYY!!!!!!!!!!!!!!!

File Contents

# Content
1 \documentclass[11pt]{article}
2
3 \usepackage{graphicx}
4 \usepackage{color}
5 \usepackage{floatflt}
6 \usepackage{amsmath}
7 \usepackage{amssymb}
8 \usepackage{subfigure}
9 \usepackage{palatino}
10 \usepackage[ref]{overcite}
11
12
13
14 \pagestyle{plain}
15 \pagenumbering{arabic}
16 \oddsidemargin 0.0cm \evensidemargin 0.0cm
17 \topmargin -21pt \headsep 10pt
18 \textheight 9.0in \textwidth 6.5in
19 \brokenpenalty=10000
20 \renewcommand{\baselinestretch}{1.2}
21 \renewcommand\citemid{\ } % no comma in optional reference note
22
23
24 \begin{document}
25
26
27 \title{A Mesoscale Model for Phospholipid Simulations}
28
29 \author{Matthew A. Meineke\\
30 Department of Chemistry and Biochemistry\\
31 University of Notre Dame\\
32 Notre Dame, Indiana 46556}
33
34 \date{\today}
35 \maketitle
36
37 \section{Background and Research Goals}
38
39 Simulations of phospholipid bilayers are, by necessity, quite
40 complex. The lipid molecules are large molecules containing many
41 atoms, and the head group of the lipid will typically contain charge
42 separated ions which set up a large dipole within the molecule. Adding
43 to the complexity are the number of water molecules needed to properly
44 solvate the lipid bilayer, typically 25 water molecules for every
45 lipid molecule. Because of these factors, many current simulations are
46 limited in both length and time scale due to to the sheer number of
47 calculations performed at every time step and the lifetime of the
48 researcher. A typical
49 simulation\cite{saiz02,lindahl00,venable00,Marrink01} will have around
50 64 phospholipids forming a bilayer approximately 40~$\mbox{\AA}$ by
51 50~$\mbox{\AA}$ with roughly 25 waters for every lipid. This means
52 there are on the order of 8,000 atoms needed to simulate these systems
53 and the trajectories are integrated for times up to 10 ns.
54
55 These limitations make it difficult to study certain biologically
56 interesting phenomena that don't fit within the short time and length
57 scale requirements. One such phenomenon is the existence of the ripple
58 phase ($P_{\beta'}$) of the bilayer between the gel phase
59 ($L_{\beta'}$) and the fluid phase ($L_{\alpha}$). The $P_{\beta'}$
60 phase has been shown to have a ripple period of
61 100-200~$\mbox{\AA}$.\cite{katsaras00,sengupta00} A simulation of this
62 length scale would require approximately 1,300 lipid molecules and
63 roughly 25 waters for every lipid to fully solvate the bilayer. With
64 the large number of atoms involved in a simulation of this magnitude,
65 steps \emph{must} be taken to simplify the system to the point where
66 the numbers of atoms becomes reasonable.
67
68 Another system of interest would be drug molecule diffusion through
69 the membrane. Due to the fluid-like properties of a lipid membrane,
70 not all diffusion takes place at membrane channels. It is of interest
71 to study certain molecules that may incorporate themselves directly
72 into the membrane. These molecules may then have an appreciable
73 waiting time (on the order of nanoseconds) within the
74 bilayer. Simulation of such a long time scale again requires
75 simplification of the system in order to lower the number of
76 calculations needed at each time step or to increase the length of
77 each time step.
78
79
80 \section{Methodology}
81
82 \subsection{Length and Time Scale Simplifications}
83
84 The length scale simplifications are aimed at increasing the number of
85 molecules that can be simulated without drastically increasing the
86 computational cost of the simulation. This is done through a
87 combination of substituting less expensive interactions for expensive
88 ones and decreasing the number of interaction sites per
89 molecule. Namely, groups of point charges are replaced with single
90 point-dipoles, and unified atoms are used in place of water,
91 phospholipid head groups, and alkyl groups.
92
93 The replacement of charges with dipoles allows us to replace an
94 interaction that has a relatively long range ($\frac{1}{r}$ for the
95 coulomb potential) with that of a relatively short range
96 ($\frac{1}{r^{3}}$ for dipole - dipole potentials). Combined with
97 Verlet neighbor lists,\cite{allen87:csl} this should result in an
98 algorithm which scales linearly with increasing system size. This is
99 in comparison to the Ewald sum\cite{leach01:mm} needed to compute
100 periodic replicas of the coulomb interactions, which scales at
101 best\cite{darden93:pme} by $N \ln N$.
102
103 The second step taken to simplify the number of calculations is to
104 incorporate unified models for groups of atoms. In the case of water,
105 we use the soft sticky dipole (SSD) model developed by
106 Ichiye\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md}
107 (Section~\ref{sec:ssdModel}). For the phospholipids, a unified head
108 atom with a dipole will replace the atoms in the head group, while
109 unified $\text{CH}_2$ and $\text{CH}_3$ atoms will replace the alkyl
110 units in the tails (Section~\ref{sec:lipidModel}). This model is
111 similar in practice to that of Lipowsky and Goetz\cite{goetz98}, where
112 the whole system is reduced to attractive and repulsive Lennard-Jones
113 spheres. However, our model gives a greater level of detail to each
114 unified molecule, namely through dipole, bend, and torsion
115 interactions.
116
117 The time scale simplifications are introduced so that we can take
118 longer time steps. By increasing the size of the time steps taken by
119 the simulation, we are able to integrate a given length of time using
120 fewer calculations. However, care must be taken that any
121 simplifications used still conserve the total energy of the
122 simulation. In practice, this means taking steps small enough to
123 resolve all motion in the system without accidently moving an object
124 too far along a repulsive energy surface before it feels the effect of
125 the surface.
126
127 In our simulation we have chosen to constrain all bonds to be of fixed
128 length. This means the bonds are no longer allowed to vibrate about
129 their equilibrium positions. Bond vibrations are typically the fastest
130 periodic motion in a dynamics simulation. By taking this action, we
131 are able to take time steps of 3 fs while still maintaining constant
132 energy. This is in contrast to the 1 fs time step typically needed to
133 conserve energy when bonds lengths are allowed to oscillate.
134
135 \subsection{The Soft Sticky Dipole Water Model}
136 \label{sec:ssdModel}
137
138 \begin{figure}
139 \centering
140 \includegraphics[width=50mm]{ssd.epsi}
141 \caption{The SSD model with the oxygen and hydrogen atoms drawn in for reference.}
142 \label{fig:ssdModel}
143 \end{figure}
144
145 The water model used in our simulations is a modified soft
146 Stockmayer-sphere model.\cite{stevens95} Like the Stockmayer-sphere, the SSD
147 model consists of a Lennard-Jones interaction site and a
148 dipole both located at the water's center of mass (Figure
149 \ref{fig:ssdModel}). However, the SSD model extends this by adding a
150 tetrahedral potential to correct for hydrogen bonding.
151
152 The SSD water potential for a pair of water molecules is then given by
153 the following equation:
154 \begin{equation}
155 V_{\text{SSD}} = V_{\text{LJ}}(r_{i\!j}) + V_{\text{dp}}(\mathbf{r}_{i\!j},
156 \boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
157 + V_{\text{sp}}(\mathbf{r}_{i\!j},\boldsymbol{\Omega}_{i},
158 \boldsymbol{\Omega}_{j})
159 \label{eq:ssdTotPot}
160 \end{equation}
161 where $\mathbf{r}_{ij}$ is the vector between molecules $i$ and $j$,
162 and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ are the
163 Euler angles of molecule $i$ or $j$ respectively. $V_{\text{LJ}}$ is
164 the Lennard-Jones potential:
165 \begin{equation}
166 V_{\text{LJ}} =
167 4\epsilon_{ij} \biggl[
168 \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
169 - \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
170 \biggr]
171 \label{eq:lennardJonesPot}
172 \end{equation}
173 where $\sigma_{ij}$ scales the length of the interaction, and
174 $\epsilon_{ij}$ scales the energy of the potential. For SSD,
175 $\sigma_{\text{SSD}} = 3.051 \mbox{
176 \AA}$ and $\epsilon_{\text{SSD}} = 0.152\text{ kcal/mol}$.
177 $V_{\text{dp}}$ is the dipole potential:
178 \begin{equation}
179 V_{\text{dp}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
180 \boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
181 \frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
182 -
183 \frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
184 (\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
185 {r^{5}_{ij}} \biggr]
186 \label{eq:dipolePot}
187 \end{equation}
188 where $\boldsymbol{\mu}_i$ is the dipole vector of molecule $i$,
189 $\boldsymbol{\mu}_i$ points along the z-axis in the body fixed
190 frame. This frame is then oriented in space by the Euler angles,
191 $\boldsymbol{\Omega}_i$. The SSD model specifies a dipole magnitude of
192 2.35~D for water.
193
194 The hydrogen bonding is modeled by the $V_{\text{sp}}$
195 term of the potential. Its form is as follows:
196 \begin{equation}
197 V_{\text{sp}}(\mathbf{r}_{i\!j},\boldsymbol{\Omega}_{i},
198 \boldsymbol{\Omega}_{j}) =
199 v^{\circ}[s(r_{ij})w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
200 \boldsymbol{\Omega}_{j})
201 +
202 s'(r_{ij})w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
203 \boldsymbol{\Omega}_{j})]
204 \label{eq:spPot}
205 \end{equation}
206 Where $v^\circ$ scales strength of the interaction.
207 $w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
208 and
209 $w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
210 are responsible for the tetrahedral potential and a correction to the
211 tetrahedral potential respectively. They are,
212 \begin{equation}
213 w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) =
214 \sin\theta_{ij} \sin 2\theta_{ij} \cos 2\phi_{ij}
215 + \sin \theta_{ji} \sin 2\theta_{ji} \cos 2\phi_{ji}
216 \label{eq:spPot2}
217 \end{equation}
218 and
219 \begin{equation}
220 \begin{split}
221 w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) =
222 &(\cos\theta_{ij}-0.6)^2(\cos\theta_{ij} + 0.8)^2 \\
223 &+ (\cos\theta_{ji}-0.6)^2(\cos\theta_{ji} + 0.8)^2 - 2w^{\circ}
224 \end{split}
225 \label{eq:spCorrection}
226 \end{equation}
227 The angles $\theta_{ij}$ and $\phi_{ij}$ are defined by the spherical
228 coordinates of the position of molecule $j$ in the reference frame
229 fixed on molecule $i$ with the z-axis aligned with the dipole moment.
230 The correction
231 $w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
232 is needed because
233 $w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
234 vanishes when $\theta_{ij}$ is $0^\circ$ or $180^\circ$.
235
236 Finally, the sticky potential is scaled by a switching function,
237 $s(r_{ij})$, that scales smoothly between 0 and 1. It is represented
238 by:
239 \begin{equation}
240 s(r_{ij}) =
241 \begin{cases}
242 1& \text{if $r_{ij} < r_{L}$}, \\
243 \frac{(r_{U} - r_{ij})^2 (r_{U} + 2r_{ij}
244 - 3r_{L})}{(r_{U}-r_{L})^3}&
245 \text{if $r_{L} \leq r_{ij} \leq r_{U}$},\\
246 0& \text{if $r_{ij} \geq r_{U}$}.
247 \end{cases}
248 \label{eq:spCutoff}
249 \end{equation}
250
251 Despite the apparent complexity of Equation \ref{eq:spPot}, the SSD
252 model is still computationally inexpensive. This is due to Equation
253 \ref{eq:spCutoff}. With $r_{L}$ being 2.75~$\mbox{\AA}$ and $r_{U}$
254 being equal to either 3.35~$\mbox{\AA}$ for $s(r_{ij})$ or
255 4.0~$\mbox{\AA}$ for $s'(r_{ij})$, the sticky potential is only active
256 over an extremely short range, and then only with other SSD
257 molecules. Therefore, it's predominant interaction is through the
258 point dipole and the Lennard-Jones sphere. Of these, only the dipole
259 interaction can be considered ``long-range''.
260
261 \subsection{The Phospholipid Model}
262 \label{sec:lipidModel}
263
264 \begin{figure}
265 \centering
266 \includegraphics[angle=-90,width=80mm]{lipidModel.epsi}
267 \caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length.}
268 \label{fig:lipidModel}
269 \end{figure}
270
271 The lipid molecules in our simulations are unified atom models. Figure
272 \ref{fig:lipidModel} shows a schematic for one of our
273 lipids. The head group of the phospholipid is replaced by a single
274 Lennard-Jones sphere with a freely oriented dipole placed at it's
275 center. The magnitude of the dipole moment is 20.6 D, chosen to match
276 that of DPPC\cite{Cevc87}. The tail atoms are unified $\text{CH}_2$
277 and $\text{CH}_3$ atoms and are also modeled as Lennard-Jones
278 spheres. The total potential for the lipid is represented by Equation
279 \ref{eq:lipidModelPot}.
280
281 \begin{equation}
282 V_{\text{lipid}} =
283 \sum_{i}V_{i}^{\text{internal}}
284 + \sum_i \sum_{j>i} \sum_{\alpha_i}
285 \sum_{\beta_j}
286 V_{\text{LJ}}(r_{\alpha_{i}\beta_{j}})
287 +\sum_i\sum_{j>i}V_{\text{dp}}(r_{1_i,1_j},\Omega_{1_i},\Omega_{1_j})
288 \label{eq:lipidModelPot}
289 \end{equation}
290 where,
291 \begin{equation}
292 V_{i}^{\text{internal}} =
293 \sum_{\text{bends}}V_{\text{bend}}(\theta_{\alpha\beta\gamma})
294 + \sum_{\text{torsions}}V_{\text{tors.}}(\phi_{\alpha\beta\gamma\zeta})
295 + \sum_{\alpha_i} \sum_{\beta_i > \alpha_i}V_{\text{LJ}}
296 (r_{\alpha_i \beta_i})
297 \label{eq:lipidModelPotInternal}
298 \end{equation}
299
300 The non-bonded interactions, $V_{\text{LJ}}$ and $V_{\text{dp}}$, are
301 the Lennard-Jones and dipole-dipole interactions respectively. For the
302 bonded potentials, only the bend and the torsional potentials are
303 calculated. The bond potential is not calculated, and the bond lengths
304 are constrained via RATTLE.\cite{leach01:mm} The bend potential is of
305 the form:
306 \begin{equation}
307 V_{\text{bend}}(\theta_{\alpha\beta\gamma})
308 = k_{\theta}\frac{(\theta_{\alpha\beta\gamma} - \theta_0)^2}{2}
309 \label{eq:bendPot}
310 \end{equation}
311 Where $k_{\theta}$ sets the stiffness of the bend potential, and $\theta_0$
312 sets the equilibrium bend angle. The torsion potential was given by:
313 \begin{equation}
314 V_{\text{tors.}}(\phi_{\alpha\beta\gamma\zeta})
315 = c_1 [1+\cos\phi_{\alpha\beta\gamma\zeta}]
316 + c_2 [1 - \cos(2\phi_{\alpha\beta\gamma\zeta})]
317 + c_3 [1 + \cos(3\phi_{\alpha\beta\gamma\zeta})]
318 \label{eq:torsPot}
319 \end{equation}
320 All parameters for bonded and non-bonded potentials in the tail atoms
321 were taken from TraPPE.\cite{Siepmann1998} The bonded interactions for
322 the head atom were also taken from TraPPE, however it's dipole moment
323 and mass were based on the properties of the phosphatidylcholine head
324 group. The Lennard-Jones parameter for the head group was chosen such
325 that it was roughly twice the size of a $\text{CH}_3$ atom, and it's
326 well depth was set to be approximately equal to that of $\text{CH}_3$.
327
328 \section{Initial Simulation: 25 lipids in water}
329 \label{sec:5x5}
330
331 \subsection{Starting Configuration and Parameters}
332 \label{sec:5x5Start}
333
334 \begin{figure}
335 \centering
336 \mbox{
337 \subfigure[The starting configuration of the 25 lipid system.]{%
338 \label{fig:5x5Start}%
339 \includegraphics[width=70mm]{5x5-initial.eps}}\quad
340 \subfigure[The 25 lipid system at 6.27~ns.]{%
341 \label{fig:5x5Final}%
342 \includegraphics[width=70mm]{5x5-6.27ns.eps}}
343 }
344 \caption{Snapshots of the 25 lipid system. Carbon tail atoms are drawn in gray, the phospholipid head groups are colored blue, and the waters are scaled down for visibility. A box has been drawn around the periodic image.}
345 \end{figure}
346
347 Our first simulation is an array of 25 single chain lipids in a sea
348 of water (Figure \ref{fig:5x5Start}). The total number of water
349 molecules is 1386, giving a final of water concentration of 70\%
350 wt. The simulation box measures 34.5~$\mbox{\AA}$ x 39.4~$\mbox{\AA}$
351 x 39.4~$\mbox{\AA}$ with periodic boundary conditions imposed. The
352 system is simulated in the micro-canonical (NVE) ensemble with an
353 average temperature of 300~K.
354
355 \subsection{Results}
356 \label{sec:5x5Results}
357
358 \begin{figure}
359 \centering
360 \subfigure[The self correlation of the phospholipid head groups. $g(r)$ is on the top, the bottom chart is the $g_\gamma(r)$.]{%
361 \label{fig:5x5HHCorr}%
362 \includegraphics[angle=-90,width=80mm]{all5x5-HEAD-HEAD.epsi}%
363 }
364 \subfigure[The $g(r)$ for the $\text{CH}_2$ molecules in the chain tails]{%
365 \label{fig:5x5CCg}%
366 \includegraphics[angle=-90,width=80mm]{all5x5-CH2-CH2.epsi}}
367 \subfigure%
368 [The pair correlations between the head groups and the water]{%
369 \label{fig:5x5HXCorr}%
370 \includegraphics[angle=-90,width=80mm]{all5x5-HEAD-X.epsi}}
371 \caption{The pair correlation functions for the 25 lipid system}
372 \end{figure}
373
374
375 Figure \ref{fig:5x5Final} shows a snapshot of the system at
376 6.27~ns. Note that the system has spontaneously self assembled into a
377 bilayer. Discussion of the length scales of the bilayer will follow in
378 this section. However, it is interesting to note a key qualitative
379 property of the system revealed by this snapshot, the tail chains are
380 tilted to the bilayer normal. This is usually indicative of the gel
381 ($L_{\beta'}$) phase. In this system, the box size is probably too
382 small for the bilayer to relax to the fluid ($P_{\alpha}$) phase. This
383 demonstrates a need for an isobaric-isothermal ensemble where the box
384 size may relax or expand to keep the system at 1~atm.
385
386 The simulation was analyzed using the radial distribution function,
387 $g(r)$, which has the form:
388 \begin{equation}
389 g(r) = \frac{V}{N_{\text{pairs}}}\langle \sum_{i} \sum_{j > i}
390 \delta(|\mathbf{r} - \mathbf{r}_{ij}|) \rangle
391 \label{eq:gofr}
392 \end{equation}
393 Equation \ref{eq:gofr} gives us information about the spacing of two
394 species as a function of radius. Essentially, if the observer were
395 located at atom $i$ and were looking out in all directions, $g(r)$
396 shows the relative density of atom $j$ at any given radius, $r$,
397 normalized by the expected density of atom $j$ at $r$. In a
398 homogeneously mixed fluid, $g(r)$ will approach 1 at large $r$, as a
399 fluid contains no long range structure to contribute peaks in the
400 number density.
401
402 For the species containing dipoles, a second pair-wise distribution
403 function was used, $g_{\gamma}(r)$. It is of the form:
404 \begin{equation}
405 g_{\gamma}(r) = \langle \sum_i \sum_{j>i}
406 (\cos \gamma_{ij}) \delta(| \mathbf{r} - \mathbf{r}_{ij}|) \rangle
407 \label{eq:gammaofr}
408 \end{equation}
409 Where $\gamma_{ij}$ is the angle between the dipole of atom $j$ with
410 respect to the dipole of atom $i$. This correlation will vary between
411 $+1$ and $-1$ when the two dipoles are perfectly aligned and
412 anti-aligned respectively. This then gives us information about how
413 directional species are aligned with each other as a function of
414 distance.
415
416 Figure \ref{fig:5x5HHCorr} shows the two self correlation functions
417 for the Head groups of the lipids. The first peak in the $g(r)$ at
418 4.03~$\mbox{\AA}$ is the nearest neighbor separation of the heads of
419 two lipids. This corresponds to a maximum in the $g_{\gamma}(r)$ which
420 means that the two neighbors on the same leaf have their dipoles
421 aligned. The broad peak at 6.5~$\mbox{\AA}$ is the inter-surface
422 spacing. Here, there is a corresponding anti-alignment in the angular
423 correlation. This means that although the dipoles are aligned on the
424 same monolayer, the dipoles will orient themselves to be anti-aligned
425 on a opposite facing monolayer. With this information, the two peaks
426 at 7.0~$\mbox{\AA}$ and 7.4~$\mbox{\AA}$ are head groups on the same
427 monolayer, and they are the second nearest neighbors to the head
428 group. The peak is likely a split peak because of the small statistics
429 of this system. Finally, the peak at 8.0~$\mbox{\AA}$ is likely the
430 second nearest neighbor on the opposite monolayer because of the
431 anti-alignment evident in the angular correlation.
432
433 Figure \ref{fig:5x5CCg} shows the radial distribution function for the
434 $\text{CH}_2$ unified atoms. The spacing of the atoms along the tail
435 chains accounts for the regularly spaced sharp peaks, but the broad
436 underlying peak with its maximum at 4.6~$\mbox{\AA}$ is the
437 distribution of chain-chain distances between two lipids. The final
438 Figure, Figure \ref{fig:5x5HXCorr}, includes the correlation functions
439 between the Head group and the SSD atoms. The peak in $g(r)$ at
440 3.6~$\mbox{\AA}$ is the distance of closest approach between the two,
441 and $g_{\gamma}(r)$ shows that the SSD atoms will align their dipoles
442 with the head groups at close distance. However, as one increases the
443 distance, the SSD atoms are no longer aligned.
444
445 \section{Second Simulation: 50 randomly oriented lipids in water}
446 \label{sec:r50}
447
448 \subsection{Starting Configuration and Parameters}
449 \label{sec:r50Start}
450
451 \begin{figure}
452 \centering
453 \mbox{
454 \subfigure[The starting configuration of the 50 lipid system.]{%
455 \label{fig:r50Start}%
456 \includegraphics[width=70mm]{r50-initial.eps}}\quad
457 \subfigure[The 50 lipid system at 2.21~ns]{%
458 \label{fig:r50Final}%
459 \includegraphics[width=70mm]{r50-2.21ns.eps}}
460 }
461 \caption{Snapshots of the 50 lipid system}
462 \end{figure}
463
464 The second simulation consists of 50 single chained lipid molecules
465 embedded in a sea of 1384 SSD waters (54\% wt.). The lipids in this
466 simulation were started with random orientation and location (Figure
467 \ref{fig:r50Start} ) The simulation box measured 26.6~$\mbox{\AA}$ x
468 26.6~$\mbox{\AA}$ x 108.4~$\mbox{\AA}$ with periodic boundary conditions
469 imposed. The simulation was run in the NVE ensemble with an average
470 temperature of 300~K.
471
472 \subsection{Results}
473 \label{sec:r50Results}
474
475 \begin{figure}
476 \centering
477 \subfigure[The self correlation of the phospholipid head groups.]{%
478 \label{fig:r50HHCorr}%
479 \includegraphics[angle=-90,width=80mm]{r50-HEAD-HEAD.epsi}%
480 }
481 \subfigure%
482 [The pair correlations between the head groups and the water]{%
483 \label{fig:r50HXCorr}%
484 \includegraphics[angle=-90,width=80mm]{r50-HEAD-X.epsi}%
485 }
486 \subfigure[The $g(r)$ for the $\text{CH}_2$ molecules in the chain tails]{%
487 \label{fig:r50CCg}%
488 \includegraphics[angle=-90,width=80mm]{r50-CH2-CH2.epsi}}
489
490 \caption{The pair correlation functions for the 50 lipid system}
491 \end{figure}
492
493 Figure \ref{fig:r50Final} is a snapshot of the system at 2.21~ns. Here
494 we see that the system has already aggregated into several micelles
495 and two are even starting to merge. It will be interesting to watch as
496 this simulation continues what the total time scale for the micelle
497 aggregation and bilayer formation will be, in Marrink's\cite{Marrink01}
498 simulation, bilayer aggregation is predicted to occur around 10~ns.
499
500 Figures \ref{fig:r50HHCorr}, \ref{fig:r50HXCorr}, and \ref{fig:r50CCg} are
501 the same correlation functions for the random 50 simulation as for the
502 previous simulation of 25 lipids. What is most interesting to note, is
503 the high degree of similarity between the correlation functions
504 between systems. Even though the 25 lipid simulation formed a bilayer
505 and the random 50 simulation is still in the micelle stage, both have
506 an inter-surface spacing of 6.5~$\mbox{\AA}$ with the same
507 characteristic anti-alignment between surfaces. Not as surprising, is
508 the consistency of the closest packing statistics between
509 systems. Namely, a head-head closest approach distance of
510 4~$\mbox{\AA}$, and similar findings for the chain-chain and
511 head-water distributions as in the 25 lipid system.
512
513 \section{Future Directions}
514
515 Current simulations indicate that our model is a feasible one, however
516 improvements will need to be made to allow the system to be simulated
517 in the isobaric-isothermal ensemble. This will relax the system to an
518 equilibrium configuration at room temperature and pressure allowing us
519 to compare our model to experimental results. Also, we are in the
520 process of parallelizing the code for an even greater speedup. This
521 will allow us to simulate large enough systems to examine phenomena
522 such as the ripple phase and drug molecule diffusion
523
524 Once the work has been completed on the simulation engine, we will
525 then use it to explore the phase diagram for our model. By
526 characterizing how our model parameters affect the bilayer properties,
527 we will tailor our model to more closely match real biological
528 molecules. With this information, we will then incorporate
529 biologically relevant molecules into the system and observe their
530 transport properties across the membrane.
531
532 \section{Acknowledgments}
533
534 I would like to thank Dr.~J.~Daniel Gezelter for his guidance on this
535 project. I would also like to acknowledge the following for their help
536 and discussions during this project: Christopher Fennell, Charles
537 Vardeman, Teng Lin, Megan Sprague, Patrick Conforti, and Dan
538 Combest. Funding for this project came from the National Science
539 Foundation.
540
541 \pagebreak
542 \bibliographystyle{achemso}
543 \bibliography{canidacy_paper}
544 \end{document}