| 5 |
|
\usepackage{floatflt} |
| 6 |
|
\usepackage{amsmath} |
| 7 |
|
\usepackage{amssymb} |
| 8 |
+ |
\usepackage{subfigure} |
| 9 |
+ |
\usepackage{palatino} |
| 10 |
|
\usepackage[ref]{overcite} |
| 11 |
|
|
| 12 |
|
|
| 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. |
| 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 phenomena is the existence of the ripple |
| 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 |
| 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, point charge distributions are replaced with |
| 90 |
< |
dipoles, and unified atoms are used in place of water, phospholipid |
| 91 |
< |
head groups, and alkyl groups. |
| 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 charge distributions with dipoles allows us to |
| 94 |
< |
replace an interaction that has a relatively long range ($\frac{1}{r}$ |
| 95 |
< |
for the coulomb potential) with that of a relatively short range |
| 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 wich scales linearly with increasing system size. This is in |
| 99 |
< |
comparison to the Ewald sum\cite{leach01:mm} needed to compute |
| 100 |
< |
periodic replicas of the coulomb interactions, which scales at best by |
| 101 |
< |
$N \ln N$. |
| 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, |
| 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}). |
| 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 |
| 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 |
| 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 Water Model} |
| 135 |
> |
\subsection{The Soft Sticky Dipole Water Model} |
| 136 |
|
\label{sec:ssdModel} |
| 137 |
|
|
| 138 |
|
\begin{figure} |
| 139 |
< |
\begin{center} |
| 139 |
> |
\centering |
| 140 |
|
\includegraphics[width=50mm]{ssd.epsi} |
| 141 |
|
\caption{The SSD model with the oxygen and hydrogen atoms drawn in for reference.} |
| 135 |
– |
\end{center} |
| 142 |
|
\label{fig:ssdModel} |
| 143 |
|
\end{figure} |
| 144 |
|
|
| 159 |
|
\label{eq:ssdTotPot} |
| 160 |
|
\end{equation} |
| 161 |
|
where $\mathbf{r}_{ij}$ is the vector between molecules $i$ and $j$, |
| 162 |
< |
and $\boldsymbol{\Omega}$ is the orientation of molecule $i$ or $j$ |
| 163 |
< |
respectively. $V_{\text{LJ}}$ is the Lennard-Jones potential: |
| 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[ |
| 170 |
|
\biggr] |
| 171 |
|
\label{eq:lennardJonesPot} |
| 172 |
|
\end{equation} |
| 173 |
< |
here $\sigma_{ij}$ |
| 174 |
< |
scales the length of the interaction, and $\epsilon_{ij}$ scales the |
| 175 |
< |
energy of the potential. For SSD, $\sigma_{\text{SSD}} = 3.051 \mbox{ |
| 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} |
| 186 |
|
\label{eq:dipolePot} |
| 187 |
|
\end{equation} |
| 188 |
|
where $\boldsymbol{\mu}_i$ is the dipole vector of molecule $i$, |
| 189 |
< |
$\boldsymbol{\mu}_i$ takes its orientation from |
| 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 |
|
|
| 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 cutoff function, |
| 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} |
| 262 |
|
\label{sec:lipidModel} |
| 263 |
|
|
| 264 |
|
\begin{figure} |
| 265 |
< |
\begin{center} |
| 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. \vspace{5mm}} |
| 260 |
< |
\end{center} |
| 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 |
| 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$ |
| 281 |
|
\begin{equation} |
| 282 |
|
V_{\text{lipid}} = |
| 283 |
|
\sum_{i}V_{i}^{\text{internal}} |
| 284 |
< |
+ \sum_i \sum_{j>i} \sum_{\text{$\alpha$ in $i$}} |
| 285 |
< |
\sum_{\text{$\beta$ in $j$}} |
| 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} |
| 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} \sum_{\beta>\alpha}V_{\text{LJ}}(r_{\alpha \beta}) |
| 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 |
|
|
| 332 |
|
\label{sec:5x5Start} |
| 333 |
|
|
| 334 |
|
\begin{figure} |
| 335 |
< |
\begin{center} |
| 336 |
< |
\includegraphics[width=70mm]{5x5-initial.eps} |
| 337 |
< |
\caption{The starting configuration of the 25 lipid system. A box is drawn around the periodic image.} |
| 338 |
< |
\end{center} |
| 339 |
< |
\label{fig:5x5Start} |
| 340 |
< |
\end{figure} |
| 341 |
< |
|
| 342 |
< |
\begin{figure} |
| 343 |
< |
\begin{center} |
| 344 |
< |
\includegraphics[width=70mm]{5x5-6.27ns.eps} |
| 337 |
< |
\caption{The 25 lipid system at 6.27~ns} |
| 338 |
< |
\end{center} |
| 339 |
< |
\label{fig:5x5Final} |
| 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 |
| 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 |
| 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 a 1~atm. |
| 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: |
| 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 |
| 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} |
| 449 |
|
\label{sec:r50Start} |
| 450 |
|
|
| 451 |
|
\begin{figure} |
| 452 |
< |
\begin{center} |
| 453 |
< |
\includegraphics[width=70mm]{r50-initial.eps} |
| 454 |
< |
\caption{The starting configuration of the 50 lipid system.} |
| 455 |
< |
\end{center} |
| 456 |
< |
\label{fig:r50Start} |
| 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 |
|
|
| 437 |
– |
\begin{figure} |
| 438 |
– |
\begin{center} |
| 439 |
– |
\includegraphics[width=70mm]{r50-2.21ns.eps} |
| 440 |
– |
\caption{The 50 lipid system at 2.21~ns} |
| 441 |
– |
\end{center} |
| 442 |
– |
\label{fig:r50Final} |
| 443 |
– |
\end{figure} |
| 444 |
– |
|
| 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 |
| 472 |
|
\subsection{Results} |
| 473 |
|
\label{sec:r50Results} |
| 474 |
|
|
| 475 |
< |
Figure \ref{fig:r50Final} is a snapshot of the system at 2.0~ns. Here |
| 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. |
| 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:r50CCg}, and \ref{fig:r50} are |
| 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 |
| 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 parallelizeing the code for an even greater speedup. This |
| 521 |
< |
will allow us to simulate the size systems needed to examine phenomena |
| 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 |
| 531 |
|
|
| 532 |
|
\section{Acknowledgments} |
| 533 |
|
|
| 534 |
< |
I would like to thank Dr. J.Daniel Gezelter for his guidance on this |
| 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 |
