ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/matt_papers/canidacy_paper/canidacy_paper.tex
Revision: 507
Committed: Mon Apr 28 16:03:10 2003 UTC (21 years, 2 months ago) by mmeineke
Content type: application/x-tex
File size: 23895 byte(s)
Log Message:
added all the files

File Contents

# User Rev Content
1 mmeineke 95 \documentclass[11pt]{article}
2    
3     \usepackage{graphicx}
4 mmeineke 106 \usepackage{color}
5 mmeineke 98 \usepackage{floatflt}
6 mmeineke 95 \usepackage{amsmath}
7     \usepackage{amssymb}
8 mmeineke 111 \usepackage{subfigure}
9     \usepackage{palatino}
10 mmeineke 95 \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 mmeineke 98
27 mmeineke 95 \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 mmeineke 106 Simulations of phospholipid bilayers are, by necessity, quite
40     complex. The lipid molecules are large molecules containing many
41 mmeineke 110 atoms, and the head group of the lipid will typically contain charge
42 mmeineke 106 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 mmeineke 110 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 mmeineke 106 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 mmeineke 111 and the trajectories are integrated for times up to 10 ns.
54 mmeineke 106
55     These limitations make it difficult to study certain biologically
56     interesting phenomena that don't fit within the short time and length
57 mmeineke 111 scale requirements. One such phenomenon is the existence of the ripple
58 mmeineke 106 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 mmeineke 110 the numbers of atoms becomes reasonable.
67 mmeineke 106
68     Another system of interest would be drug molecule diffusion through
69 mmeineke 110 the membrane. Due to the fluid-like properties of a lipid membrane,
70 mmeineke 106 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 mmeineke 110 calculations needed at each time step or to increase the length of
77     each time step.
78 mmeineke 106
79    
80 mmeineke 95 \section{Methodology}
81    
82 mmeineke 96 \subsection{Length and Time Scale Simplifications}
83 mmeineke 95
84 mmeineke 106 The length scale simplifications are aimed at increasing the number of
85     molecules that can be simulated without drastically increasing the
86 mmeineke 110 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 mmeineke 111 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 mmeineke 95
93 mmeineke 111 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 mmeineke 110 ($\frac{1}{r^{3}}$ for dipole - dipole potentials). Combined with
97     Verlet neighbor lists,\cite{allen87:csl} this should result in an
98 mmeineke 111 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 mmeineke 95
103 mmeineke 106 The second step taken to simplify the number of calculations is to
104 mmeineke 95 incorporate unified models for groups of atoms. In the case of water,
105     we use the soft sticky dipole (SSD) model developed by
106 mmeineke 110 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 mmeineke 111 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 mmeineke 95
117 mmeineke 106 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 mmeineke 110 the simulation, we are able to integrate a given length of time using
120     fewer calculations. However, care must be taken that any
121 mmeineke 111 simplifications used still conserve the total energy of the
122 mmeineke 106 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 mmeineke 95
127 mmeineke 96 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 mmeineke 106 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 mmeineke 95
135 mmeineke 111 \subsection{The Soft Sticky Dipole Water Model}
136 mmeineke 95 \label{sec:ssdModel}
137    
138 mmeineke 110 \begin{figure}
139 mmeineke 111 \centering
140 mmeineke 110 \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 mmeineke 96
145 mmeineke 106 The water model used in our simulations is a modified soft
146     Stockmayer-sphere model.\cite{stevens95} Like the Stockmayer-sphere, the SSD
147 mmeineke 110 model consists of a Lennard-Jones interaction site and a
148 mmeineke 98 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 mmeineke 110 The SSD water potential for a pair of water molecules is then given by
153     the following equation:
154 mmeineke 96 \begin{equation}
155 mmeineke 106 V_{\text{SSD}} = V_{\text{LJ}}(r_{i\!j}) + V_{\text{dp}}(\mathbf{r}_{i\!j},
156 mmeineke 96 \boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
157     + V_{\text{sp}}(\mathbf{r}_{i\!j},\boldsymbol{\Omega}_{i},
158     \boldsymbol{\Omega}_{j})
159 mmeineke 98 \label{eq:ssdTotPot}
160 mmeineke 96 \end{equation}
161 mmeineke 110 where $\mathbf{r}_{ij}$ is the vector between molecules $i$ and $j$,
162 mmeineke 111 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 mmeineke 106 \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 mmeineke 111 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 mmeineke 106 \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 mmeineke 110 where $\boldsymbol{\mu}_i$ is the dipole vector of molecule $i$,
189 mmeineke 111 $\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 mmeineke 110 $\boldsymbol{\Omega}_i$. The SSD model specifies a dipole magnitude of
192     2.35~D for water.
193 mmeineke 96
194 mmeineke 110 The hydrogen bonding is modeled by the $V_{\text{sp}}$
195 mmeineke 106 term of the potential. Its form is as follows:
196 mmeineke 96 \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 mmeineke 98 \label{eq:spPot}
205 mmeineke 96 \end{equation}
206 mmeineke 106 Where $v^\circ$ scales strength of the interaction.
207 mmeineke 98 $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 mmeineke 96 \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 mmeineke 106 \label{eq:spPot2}
217 mmeineke 96 \end{equation}
218 mmeineke 98 and
219 mmeineke 96 \begin{equation}
220     \begin{split}
221 mmeineke 98 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 mmeineke 96 \end{split}
225 mmeineke 98 \label{eq:spCorrection}
226 mmeineke 96 \end{equation}
227 mmeineke 106 The angles $\theta_{ij}$ and $\phi_{ij}$ are defined by the spherical
228 mmeineke 110 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 mmeineke 106 The correction
231 mmeineke 98 $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 mmeineke 106 vanishes when $\theta_{ij}$ is $0^\circ$ or $180^\circ$.
235 mmeineke 96
236 mmeineke 111 Finally, the sticky potential is scaled by a switching function,
237 mmeineke 110 $s(r_{ij})$, that scales smoothly between 0 and 1. It is represented
238 mmeineke 101 by:
239 mmeineke 96 \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 mmeineke 98 \label{eq:spCutoff}
249 mmeineke 96 \end{equation}
250    
251 mmeineke 106 Despite the apparent complexity of Equation \ref{eq:spPot}, the SSD
252 mmeineke 108 model is still computationally inexpensive. This is due to Equation
253 mmeineke 107 \ref{eq:spCutoff}. With $r_{L}$ being 2.75~$\mbox{\AA}$ and $r_{U}$
254 mmeineke 106 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 mmeineke 108 over an extremely short range, and then only with other SSD
257 mmeineke 110 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 mmeineke 98
261 mmeineke 95 \subsection{The Phospholipid Model}
262     \label{sec:lipidModel}
263    
264 mmeineke 110 \begin{figure}
265 mmeineke 111 \centering
266 mmeineke 110 \includegraphics[angle=-90,width=80mm]{lipidModel.epsi}
267 mmeineke 111 \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 mmeineke 110 \label{fig:lipidModel}
269     \end{figure}
270 mmeineke 95
271 mmeineke 99 The lipid molecules in our simulations are unified atom models. Figure
272 mmeineke 110 \ref{fig:lipidModel} shows a schematic for one of our
273 mmeineke 111 lipids. The head group of the phospholipid is replaced by a single
274 mmeineke 99 Lennard-Jones sphere with a freely oriented dipole placed at it's
275 mmeineke 110 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 mmeineke 99
281     \begin{equation}
282 mmeineke 110 V_{\text{lipid}} =
283     \sum_{i}V_{i}^{\text{internal}}
284 mmeineke 111 + \sum_i \sum_{j>i} \sum_{\alpha_i}
285     \sum_{\beta_j}
286 mmeineke 110 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 mmeineke 99 \label{eq:lipidModelPot}
289     \end{equation}
290 mmeineke 110 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 mmeineke 507 + \sum_{\alpha_i} \sum_{\beta_i > (\alpha_i + 4)}V_{\text{LJ}}
296 mmeineke 111 (r_{\alpha_i \beta_i})
297 mmeineke 110 \label{eq:lipidModelPotInternal}
298     \end{equation}
299 mmeineke 99
300 mmeineke 106 The non-bonded interactions, $V_{\text{LJ}}$ and $V_{\text{dp}}$, are
301 mmeineke 99 the Lennard-Jones and dipole-dipole interactions respectively. For the
302 mmeineke 110 bonded potentials, only the bend and the torsional potentials are
303 mmeineke 99 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 mmeineke 110 V_{\text{bend}}(\theta_{\alpha\beta\gamma})
308     = k_{\theta}\frac{(\theta_{\alpha\beta\gamma} - \theta_0)^2}{2}
309 mmeineke 99 \label{eq:bendPot}
310     \end{equation}
311 mmeineke 100 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 mmeineke 110 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 mmeineke 100 \label{eq:torsPot}
319     \end{equation}
320 mmeineke 102 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 mmeineke 106 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 mmeineke 99
328 mmeineke 102 \section{Initial Simulation: 25 lipids in water}
329     \label{sec:5x5}
330 mmeineke 101
331 mmeineke 106 \subsection{Starting Configuration and Parameters}
332 mmeineke 102 \label{sec:5x5Start}
333 mmeineke 101
334 mmeineke 110 \begin{figure}
335 mmeineke 111 \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 mmeineke 110 \end{figure}
346    
347     Our first simulation is an array of 25 single chain lipids in a sea
348 mmeineke 102 of water (Figure \ref{fig:5x5Start}). The total number of water
349 mmeineke 108 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 mmeineke 107 x 39.4~$\mbox{\AA}$ with periodic boundary conditions imposed. The
352 mmeineke 108 system is simulated in the micro-canonical (NVE) ensemble with an
353 mmeineke 102 average temperature of 300~K.
354    
355     \subsection{Results}
356     \label{sec:5x5Results}
357    
358 mmeineke 111 \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 mmeineke 102 Figure \ref{fig:5x5Final} shows a snapshot of the system at
376 mmeineke 110 6.27~ns. Note that the system has spontaneously self assembled into a
377 mmeineke 107 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 mmeineke 111 size may relax or expand to keep the system at 1~atm.
385 mmeineke 102
386 mmeineke 108 The simulation was analyzed using the radial distribution function,
387     $g(r)$, which has the form:
388 mmeineke 106 \begin{equation}
389 mmeineke 107 g(r) = \frac{V}{N_{\text{pairs}}}\langle \sum_{i} \sum_{j > i}
390 mmeineke 106 \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 mmeineke 102
402 mmeineke 111 For the species containing dipoles, a second pair-wise distribution
403 mmeineke 106 function was used, $g_{\gamma}(r)$. It is of the form:
404     \begin{equation}
405 mmeineke 110 g_{\gamma}(r) = \langle \sum_i \sum_{j>i}
406     (\cos \gamma_{ij}) \delta(| \mathbf{r} - \mathbf{r}_{ij}|) \rangle
407 mmeineke 106 \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 mmeineke 102
416 mmeineke 106 Figure \ref{fig:5x5HHCorr} shows the two self correlation functions
417 mmeineke 108 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 mmeineke 110 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 mmeineke 108 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 mmeineke 102
433 mmeineke 108 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 mmeineke 102
445 mmeineke 108 \section{Second Simulation: 50 randomly oriented lipids in water}
446     \label{sec:r50}
447 mmeineke 102
448 mmeineke 108 \subsection{Starting Configuration and Parameters}
449     \label{sec:r50Start}
450 mmeineke 102
451 mmeineke 110 \begin{figure}
452 mmeineke 111 \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 mmeineke 110 \end{figure}
463    
464 mmeineke 108 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 mmeineke 110 \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 mmeineke 108 imposed. The simulation was run in the NVE ensemble with an average
470     temperature of 300~K.
471 mmeineke 102
472 mmeineke 108 \subsection{Results}
473     \label{sec:r50Results}
474 mmeineke 102
475 mmeineke 111 \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 mmeineke 108 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 mmeineke 111 aggregation and bilayer formation will be, in Marrink's\cite{Marrink01}
498     simulation, bilayer aggregation is predicted to occur around 10~ns.
499 mmeineke 108
500 mmeineke 111 Figures \ref{fig:r50HHCorr}, \ref{fig:r50HXCorr}, and \ref{fig:r50CCg} are
501 mmeineke 108 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 mmeineke 110 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 mmeineke 108
513 mmeineke 101 \section{Future Directions}
514    
515 mmeineke 108 Current simulations indicate that our model is a feasible one, however
516 mmeineke 110 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 mmeineke 111 process of parallelizing the code for an even greater speedup. This
521     will allow us to simulate large enough systems to examine phenomena
522 mmeineke 110 such as the ripple phase and drug molecule diffusion
523 mmeineke 101
524 mmeineke 110 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 mmeineke 108 characterizing how our model parameters affect the bilayer properties,
527 mmeineke 110 we will tailor our model to more closely match real biological
528     molecules. With this information, we will then incorporate
529 mmeineke 108 biologically relevant molecules into the system and observe their
530     transport properties across the membrane.
531    
532     \section{Acknowledgments}
533    
534 mmeineke 111 I would like to thank Dr.~J.~Daniel Gezelter for his guidance on this
535 mmeineke 108 project. I would also like to acknowledge the following for their help
536     and discussions during this project: Christopher Fennell, Charles
537 mmeineke 110 Vardeman, Teng Lin, Megan Sprague, Patrick Conforti, and Dan
538     Combest. Funding for this project came from the National Science
539     Foundation.
540 mmeineke 108
541     \pagebreak
542     \bibliographystyle{achemso}
543     \bibliography{canidacy_paper}
544     \end{document}