--- trunk/matt_papers/canidacy_paper/canidacy_paper.tex 2002/09/11 03:26:01 106 +++ trunk/matt_papers/canidacy_paper/canidacy_paper.tex 2002/09/16 22:10:09 110 @@ -36,18 +36,19 @@ atoms,and the head group of the lipid will typically c Simulations of phospholipid bilayers are, by necessity, quite complex. The lipid molecules are large molecules containing many -atoms,and the head group of the lipid will typically contain charge +atoms, and the head group of the lipid will typically contain charge separated ions which set up a large dipole within the molecule. Adding to the complexity are the number of water molecules needed to properly -solvate the lipid bilayer. Because of these factors, many current -simulations are limited in both length and time scale due to to the -sheer number of calculations performed at every time step and the -lifetime of the researcher. A typical +solvate the lipid bilayer, typically 25 water molecules for every +lipid molecule. Because of these factors, many current simulations are +limited in both length and time scale due to to the sheer number of +calculations performed at every time step and the lifetime of the +researcher. A typical simulation\cite{saiz02,lindahl00,venable00,Marrink01} will have around 64 phospholipids forming a bilayer approximately 40~$\mbox{\AA}$ by 50~$\mbox{\AA}$ with roughly 25 waters for every lipid. This means there are on the order of 8,000 atoms needed to simulate these systems -and the trajectories in turn are integrated for times up to 10 ns. +and the trajectories are integrated for times up to 10 ns. These limitations make it difficult to study certain biologically interesting phenomena that don't fit within the short time and length @@ -60,17 +61,18 @@ these numbers are reasonable. roughly 25 waters for every lipid to fully solvate the bilayer. With the large number of atoms involved in a simulation of this magnitude, steps \emph{must} be taken to simplify the system to the point where -these numbers are reasonable. +the numbers of atoms becomes reasonable. Another system of interest would be drug molecule diffusion through -the membrane. Due to the fluid like properties of a lipid membrane, +the membrane. Due to the fluid-like properties of a lipid membrane, not all diffusion takes place at membrane channels. It is of interest to study certain molecules that may incorporate themselves directly into the membrane. These molecules may then have an appreciable waiting time (on the order of nanoseconds) within the bilayer. Simulation of such a long time scale again requires simplification of the system in order to lower the number of -calculations needed at each time step. +calculations needed at each time step or to increase the length of +each time step. \section{Methodology} @@ -79,34 +81,36 @@ computational cost of the system. This is done by a co The length scale simplifications are aimed at increasing the number of molecules that can be simulated without drastically increasing the -computational cost of the system. This is done by a combination of -substituting less expensive interactions for expensive ones and -decreasing the number of interaction sites per molecule. Namely, -point charge distributions are replaced with dipoles, and unified atoms are -used in place of water, phospholipid head groups, and alkyl groups. +computational cost of the simulation. This is done through a +combination of substituting less expensive interactions for expensive +ones and decreasing the number of interaction sites per +molecule. Namely, point charge distributions are replaced with +dipoles, and unified atoms are used in place of water, phospholipid +head groups, and alkyl groups. The replacement of charge distributions with dipoles allows us to -replace an interaction that has a relatively long range, -$(\frac{1}{r})$, for the coulomb potential, with that of a relatively -short range, $(\frac{1}{r^{3}})$, for dipole - dipole -potentials. Combined with a computational simplification algorithm -such as a Verlet neighbor list,\cite{allen87:csl} this should give -computational scaling by $N$. This is in comparison to the Ewald -sum\cite{leach01:mm} needed to compute the coulomb interactions, which -scales at best by $N \ln N$. +replace an interaction that has a relatively long range ($\frac{1}{r}$ +for the coulomb potential) with that of a relatively short range +($\frac{1}{r^{3}}$ for dipole - dipole potentials). Combined with +Verlet neighbor lists,\cite{allen87:csl} this should result in an +algorithm wich scales linearly with increasing system size. This is in +comparison to the Ewald sum\cite{leach01:mm} needed to compute +periodic replicas of the coulomb interactions, which scales at best by +$N \ln N$. The second step taken to simplify the number of calculations is to incorporate unified models for groups of atoms. In the case of water, we use the soft sticky dipole (SSD) model developed by -Ichiye\cite{Liu96} (Section~\ref{sec:ssdModel}). For the phospholipids, a -unified head atom with a dipole will replace the atoms in the head -group, while unified $\text{CH}_2$ and $\text{CH}_3$ atoms will -replace the alkyl units in the tails (Section~\ref{sec:lipidModel}). +Ichiye\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md} +(Section~\ref{sec:ssdModel}). For the phospholipids, a unified head +atom with a dipole will replace the atoms in the head group, while +unified $\text{CH}_2$ and $\text{CH}_3$ atoms will replace the alkyl +units in the tails (Section~\ref{sec:lipidModel}). The time scale simplifications are introduced so that we can take longer time steps. By increasing the size of the time steps taken by -the simulation, we are able to integrate the simulation trajectory -with fewer calculations. However, care must be taken that any +the simulation, we are able to integrate a given length of time using +fewer calculations. However, care must be taken that any simplifications used, still conserve the total energy of the simulation. In practice, this means taking steps small enough to resolve all motion in the system without accidently moving an object @@ -124,20 +128,23 @@ conserve energy when bonds lengths are allowed to osci \subsection{The Soft Sticky Water Model} \label{sec:ssdModel} -%\begin{floatingfigure}{55mm} -%\includegraphics[width=45mm]{ssd.epsi} -%\caption{The SSD model with the oxygen and hydrogen atoms drawn in for reference. \vspace{5mm}} -%\label{fig:ssdModel} -%\end{floatingfigure} +\begin{figure} +\begin{center} +\includegraphics[width=50mm]{ssd.epsi} +\caption{The SSD model with the oxygen and hydrogen atoms drawn in for reference.} +\end{center} +\label{fig:ssdModel} +\end{figure} The water model used in our simulations is a modified soft Stockmayer-sphere model.\cite{stevens95} Like the Stockmayer-sphere, the SSD -model\cite{Liu96} consists of a Lennard-Jones interaction site and a +model consists of a Lennard-Jones interaction site and a dipole both located at the water's center of mass (Figure \ref{fig:ssdModel}). However, the SSD model extends this by adding a tetrahedral potential to correct for hydrogen bonding. -The SSD water potential is then given by the following equation: +The SSD water potential for a pair of water molecules is then given by +the following equation: \begin{equation} V_{\text{SSD}} = V_{\text{LJ}}(r_{i\!j}) + V_{\text{dp}}(\mathbf{r}_{i\!j}, \boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) @@ -145,7 +152,9 @@ $V_{\text{LJ}}$ is the Lennard-Jones potential: \boldsymbol{\Omega}_{j}) \label{eq:ssdTotPot} \end{equation} -$V_{\text{LJ}}$ is the Lennard-Jones potential: +where $\mathbf{r}_{ij}$ is the vector between molecules $i$ and $j$, +and $\boldsymbol{\Omega}$ is the orientation of molecule $i$ or $j$ +respectively. $V_{\text{LJ}}$ is the Lennard-Jones potential: \begin{equation} V_{\text{LJ}} = 4\epsilon_{ij} \biggl[ @@ -154,8 +163,8 @@ where $r_{ij}$ is the distance between two $ij$ pairs, \biggr] \label{eq:lennardJonesPot} \end{equation} -where $r_{ij}$ is the distance between two $ij$ pairs, $\sigma_{ij}$ -scales the length of the iteraction, and $\epsilon_{ij}$ scales the +here $\sigma_{ij}$ +scales the length of the interaction, and $\epsilon_{ij}$ scales the energy of the potential. For SSD, $\sigma_{\text{SSD}} = 3.051 \mbox{ \AA}$ and $\epsilon_{\text{SSD}} = 0.152\text{ kcal/mol}$. $V_{\text{dp}}$ is the dipole potential: @@ -169,12 +178,12 @@ where $\mathbf{r}_{ij}$ is the vector between $i$ and {r^{5}_{ij}} \biggr] \label{eq:dipolePot} \end{equation} -where $\mathbf{r}_{ij}$ is the vector between $i$ and $j$, -$\boldsymbol{\Omega}$ is the orientation of the species, and -$\boldsymbol{\mu}$ is the dipole vector. The SSD model specifies a dipole -magnitude of 2.35~D for water. +where $\boldsymbol{\mu}_i$ is the dipole vector of molecule $i$, +$\boldsymbol{\mu}_i$ takes its orientation from +$\boldsymbol{\Omega}_i$. The SSD model specifies a dipole magnitude of +2.35~D for water. -The hydrogen bonding of the model is governed by the $V_{\text{sp}}$ +The hydrogen bonding is modeled by the $V_{\text{sp}}$ term of the potential. Its form is as follows: \begin{equation} V_{\text{sp}}(\mathbf{r}_{i\!j},\boldsymbol{\Omega}_{i}, @@ -208,8 +217,8 @@ polar coordinates of the position of sphere $j$ in the \label{eq:spCorrection} \end{equation} The angles $\theta_{ij}$ and $\phi_{ij}$ are defined by the spherical -polar coordinates of the position of sphere $j$ in the reference frame -fixed on sphere $i$ with the z-axis aligned with the dipole moment. +coordinates of the position of molecule $j$ in the reference frame +fixed on molecule $i$ with the z-axis aligned with the dipole moment. The correction $w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$ is needed because @@ -217,7 +226,7 @@ $s(r_{ij})$ that scales smoothly between 0 and 1. It i vanishes when $\theta_{ij}$ is $0^\circ$ or $180^\circ$. Finally, the sticky potential is scaled by a cutoff function, -$s(r_{ij})$ that scales smoothly between 0 and 1. It is represented +$s(r_{ij})$, that scales smoothly between 0 and 1. It is represented by: \begin{equation} s(r_{ij}) = @@ -232,58 +241,72 @@ model is still computationaly inexpensive. This is due \end{equation} Despite the apparent complexity of Equation \ref{eq:spPot}, the SSD -model is still computationaly inexpensive. This is due to Equation -\ref{eq:spCutoff}. With $r_{L}$ being 2.75~$\mbox\AA}$ and $r_{U}$ +model is still computationally inexpensive. This is due to Equation +\ref{eq:spCutoff}. With $r_{L}$ being 2.75~$\mbox{\AA}$ and $r_{U}$ being equal to either 3.35~$\mbox{\AA}$ for $s(r_{ij})$ or 4.0~$\mbox{\AA}$ for $s'(r_{ij})$, the sticky potential is only active -over an extremly short range, and then only with other SSD -molecules. Therefore, it's predominant interaction is through it's -point dipole and Lennard-Jones sphere. +over an extremely short range, and then only with other SSD +molecules. Therefore, it's predominant interaction is through the +point dipole and the Lennard-Jones sphere. Of these, only the dipole +interaction can be considered ``long-range''. \subsection{The Phospholipid Model} \label{sec:lipidModel} -\begin{floatingfigure}{90mm} +\begin{figure} +\begin{center} \includegraphics[angle=-90,width=80mm]{lipidModel.epsi} \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}} +\end{center} \label{fig:lipidModel} -\end{floatingfigure} +\end{figure} The lipid molecules in our simulations are unified atom models. Figure -\ref{fig:lipidModel} shows a template drawing for one of our +\ref{fig:lipidModel} shows a schematic for one of our lipids. The Head group of the phospholipid is replaced by a single Lennard-Jones sphere with a freely oriented dipole placed at it's -center. The magnitude of it's dipole moment is 20.6 D. The tail atoms -are unified $\text{CH}_2$ and $\text{CH}_3$ atoms and are also modeled -as Lennard-Jones spheres. The total potential for the lipid is -represented by Equation \ref{eq:lipidModelPot}. +center. The magnitude of the dipole moment is 20.6 D, chosen to match +that of DPPC\cite{Cevc87}. The tail atoms are unified $\text{CH}_2$ +and $\text{CH}_3$ atoms and are also modeled as Lennard-Jones +spheres. The total potential for the lipid is represented by Equation +\ref{eq:lipidModelPot}. \begin{equation} -V_{\mbox{lipid}} = \overbrace{% - V_{\text{bend}}(\theta_{ijk}) +% - V_{\text{tors.}}(\phi_{ijkl})}^{bonded} - + \overbrace{% - V_{\text{LJ}}(r_{i\!j}) + - V_{\text{dp}}(r_{i\!j},\Omega_{i},\Omega_{j})% - }^{non-bonded} +V_{\text{lipid}} = + \sum_{i}V_{i}^{\text{internal}} + + \sum_i \sum_{j>i} \sum_{\text{$\alpha$ in $i$}} + \sum_{\text{$\beta$ in $j$}} + V_{\text{LJ}}(r_{\alpha_{i}\beta_{j}}) + +\sum_i\sum_{j>i}V_{\text{dp}}(r_{1_i,1_j},\Omega_{1_i},\Omega_{1_j}) \label{eq:lipidModelPot} \end{equation} +where, +\begin{equation} +V_{i}^{\text{internal}} = + \sum_{\text{bends}}V_{\text{bend}}(\theta_{\alpha\beta\gamma}) + + \sum_{\text{torsions}}V_{\text{tors.}}(\phi_{\alpha\beta\gamma\zeta}) + + \sum_{\alpha} \sum_{\beta>\alpha}V_{\text{LJ}}(r_{\alpha \beta}) +\label{eq:lipidModelPotInternal} +\end{equation} The non-bonded interactions, $V_{\text{LJ}}$ and $V_{\text{dp}}$, are the Lennard-Jones and dipole-dipole interactions respectively. For the -non-bonded potentials, only the bend and the torsional potentials are +bonded potentials, only the bend and the torsional potentials are calculated. The bond potential is not calculated, and the bond lengths are constrained via RATTLE.\cite{leach01:mm} The bend potential is of the form: \begin{equation} -V_{\text{bend}}(\theta_{ijk}) = k_{\theta}\frac{(\theta_{ijk} - \theta_0)^2}{2} +V_{\text{bend}}(\theta_{\alpha\beta\gamma}) + = k_{\theta}\frac{(\theta_{\alpha\beta\gamma} - \theta_0)^2}{2} \label{eq:bendPot} \end{equation} Where $k_{\theta}$ sets the stiffness of the bend potential, and $\theta_0$ sets the equilibrium bend angle. The torsion potential was given by: \begin{equation} -V_{\text{tors.}}(\phi_{ijkl})= c_1[1+\cos\phi_{ijkl}] - + c_2 [1 - \cos(2\phi_{ijkl})] + c_3[1 + \cos(3\phi_{ijkl})] +V_{\text{tors.}}(\phi_{\alpha\beta\gamma\zeta}) + = c_1 [1+\cos\phi_{\alpha\beta\gamma\zeta}] + + c_2 [1 - \cos(2\phi_{\alpha\beta\gamma\zeta})] + + c_3 [1 + \cos(3\phi_{\alpha\beta\gamma\zeta})] \label{eq:torsPot} \end{equation} All parameters for bonded and non-bonded potentials in the tail atoms @@ -300,33 +323,48 @@ Our first simulation was an array of 25 single chained \subsection{Starting Configuration and Parameters} \label{sec:5x5Start} -Our first simulation was an array of 25 single chained lipids in a sea +\begin{figure} +\begin{center} +\includegraphics[width=70mm]{5x5-initial.eps} +\caption{The starting configuration of the 25 lipid system. A box is drawn around the periodic image.} +\end{center} +\label{fig:5x5Start} +\end{figure} + +\begin{figure} +\begin{center} +\includegraphics[width=70mm]{5x5-6.27ns.eps} +\caption{The 25 lipid system at 6.27~ns} +\end{center} +\label{fig:5x5Final} +\end{figure} + +Our first simulation is an array of 25 single chain lipids in a sea of water (Figure \ref{fig:5x5Start}). The total number of water -molecules was 1386, giving a final of water concentration of 70\% -wt. The simulation box measured 34.5~$\mbox{\AA}$ x 39.4~$\mbox{\AA}$ -x 39.4~$\mbox{\AA}$ with periodic boundary conditions invoked. The -system was simulated in the micro-canonical (NVE) ensemble with an +molecules is 1386, giving a final of water concentration of 70\% +wt. The simulation box measures 34.5~$\mbox{\AA}$ x 39.4~$\mbox{\AA}$ +x 39.4~$\mbox{\AA}$ with periodic boundary conditions imposed. The +system is simulated in the micro-canonical (NVE) ensemble with an average temperature of 300~K. \subsection{Results} \label{sec:5x5Results} Figure \ref{fig:5x5Final} shows a snapshot of the system at -3.6~ns. Here it is exciting to note that the system has spontaneously -self assembled into a bilayer. Discussion of the length scales of the -bilayer will follow in this section. However, it is interesting to -note several qualitative properties of the system revealed by this -snapshot. First, is how the tail chains are tilted to the bilayer -normal. This is usually indicative of the gel ($L_{\beta'}$) -phase. Likely, the box size is too small for the bilayer to relax to -the fluid ($P_{\alpha}$ phase. Showing the need for a constant -pressure simulation versus the constant volume of the current -ensemble. +6.27~ns. Note that the system has spontaneously self assembled into a +bilayer. Discussion of the length scales of the bilayer will follow in +this section. However, it is interesting to note a key qualitative +property of the system revealed by this snapshot, the tail chains are +tilted to the bilayer normal. This is usually indicative of the gel +($L_{\beta'}$) phase. In this system, the box size is probably too +small for the bilayer to relax to the fluid ($P_{\alpha}$) phase. This +demonstrates a need for an isobaric-isothermal ensemble where the box +size may relax or expand to keep the system at a 1~atm. -The trajectory of the simulation was analyzed using the pair-wise -radial distribution function, $g(r)$, which has the form: +The simulation was analyzed using the radial distribution function, +$g(r)$, which has the form: \begin{equation} -g(r) = \frac{V}{N_{\text{pairs}}}\langle \sum_{i} \sum_{j \neq i} +g(r) = \frac{V}{N_{\text{pairs}}}\langle \sum_{i} \sum_{j > i} \delta(|\mathbf{r} - \mathbf{r}_{ij}|) \rangle \label{eq:gofr} \end{equation} @@ -342,7 +380,8 @@ g_{\gamma}(r) = foobar For the species containing dipoles, a second pair wise distribution function was used, $g_{\gamma}(r)$. It is of the form: \begin{equation} -g_{\gamma}(r) = foobar +g_{\gamma}(r) = \langle \sum_i \sum_{j>i} + (\cos \gamma_{ij}) \delta(| \mathbf{r} - \mathbf{r}_{ij}|) \rangle \label{eq:gammaofr} \end{equation} Where $\gamma_{ij}$ is the angle between the dipole of atom $j$ with @@ -353,19 +392,115 @@ for the Head groups of the lipids. The first peak at 4 distance. Figure \ref{fig:5x5HHCorr} shows the two self correlation functions -for the Head groups of the lipids. The first peak at 4.03~$\mbox{\AA}$ is the -nearest neighbor separation of the heads of two lipids. +for the Head groups of the lipids. The first peak in the $g(r)$ at +4.03~$\mbox{\AA}$ is the nearest neighbor separation of the heads of +two lipids. This corresponds to a maximum in the $g_{\gamma}(r)$ which +means that the two neighbors on the same leaf have their dipoles +aligned. The broad peak at 6.5~$\mbox{\AA}$ is the inter-surface +spacing. Here, there is a corresponding anti-alignment in the angular +correlation. This means that although the dipoles are aligned on the +same monolayer, the dipoles will orient themselves to be anti-aligned +on a opposite facing monolayer. With this information, the two peaks +at 7.0~$\mbox{\AA}$ and 7.4~$\mbox{\AA}$ are head groups on the same +monolayer, and they are the second nearest neighbors to the head +group. The peak is likely a split peak because of the small statistics +of this system. Finally, the peak at 8.0~$\mbox{\AA}$ is likely the +second nearest neighbor on the opposite monolayer because of the +anti-alignment evident in the angular correlation. +Figure \ref{fig:5x5CCg} shows the radial distribution function for the +$\text{CH}_2$ unified atoms. The spacing of the atoms along the tail +chains accounts for the regularly spaced sharp peaks, but the broad +underlying peak with its maximum at 4.6~$\mbox{\AA}$ is the +distribution of chain-chain distances between two lipids. The final +Figure, Figure \ref{fig:5x5HXCorr}, includes the correlation functions +between the Head group and the SSD atoms. The peak in $g(r)$ at +3.6~$\mbox{\AA}$ is the distance of closest approach between the two, +and $g_{\gamma}(r)$ shows that the SSD atoms will align their dipoles +with the head groups at close distance. However, as one increases the +distance, the SSD atoms are no longer aligned. +\section{Second Simulation: 50 randomly oriented lipids in water} +\label{sec:r50} +\subsection{Starting Configuration and Parameters} +\label{sec:r50Start} -\section{Second Simulation: 50 randomly oriented lipids in water} +\begin{figure} +\begin{center} +\includegraphics[width=70mm]{r50-initial.eps} +\caption{The starting configuration of the 50 lipid system.} +\end{center} +\label{fig:r50Start} +\end{figure} -the second simulation +\begin{figure} +\begin{center} +\includegraphics[width=70mm]{r50-2.21ns.eps} +\caption{The 50 lipid system at 2.21~ns} +\end{center} +\label{fig:r50Final} +\end{figure} +The second simulation consists of 50 single chained lipid molecules +embedded in a sea of 1384 SSD waters (54\% wt.). The lipids in this +simulation were started with random orientation and location (Figure +\ref{fig:r50Start} ) The simulation box measured 26.6~$\mbox{\AA}$ x +26.6~$\mbox{\AA}$ x 108.4~$\mbox{\AA}$ with periodic boundary conditions +imposed. The simulation was run in the NVE ensemble with an average +temperature of 300~K. + +\subsection{Results} +\label{sec:r50Results} + +Figure \ref{fig:r50Final} is a snapshot of the system at 2.0~ns. Here +we see that the system has already aggregated into several micelles +and two are even starting to merge. It will be interesting to watch as +this simulation continues what the total time scale for the micelle +aggregation and bilayer formation will be. + +Figures \ref{fig:r50HHCorr}, \ref{fig:r50CCg}, and \ref{fig:r50} are +the same correlation functions for the random 50 simulation as for the +previous simulation of 25 lipids. What is most interesting to note, is +the high degree of similarity between the correlation functions +between systems. Even though the 25 lipid simulation formed a bilayer +and the random 50 simulation is still in the micelle stage, both have +an inter-surface spacing of 6.5~$\mbox{\AA}$ with the same +characteristic anti-alignment between surfaces. Not as surprising, is +the consistency of the closest packing statistics between +systems. Namely, a head-head closest approach distance of +4~$\mbox{\AA}$, and similar findings for the chain-chain and +head-water distributions as in the 25 lipid system. + \section{Future Directions} +Current simulations indicate that our model is a feasible one, however +improvements will need to be made to allow the system to be simulated +in the isobaric-isothermal ensemble. This will relax the system to an +equilibrium configuration at room temperature and pressure allowing us +to compare our model to experimental results. Also, we are in the +process of parallelizeing the code for an even greater speedup. This +will allow us to simulate the size systems needed to examine phenomena +such as the ripple phase and drug molecule diffusion -\pagebreak -\bibliographystyle{achemso} -\bibliography{canidacy_paper} \end{document} +Once the work has been completed on the simulation engine, we will +then use it to explore the phase diagram for our model. By +characterizing how our model parameters affect the bilayer properties, +we will tailor our model to more closely match real biological +molecules. With this information, we will then incorporate +biologically relevant molecules into the system and observe their +transport properties across the membrane. + +\section{Acknowledgments} + +I would like to thank Dr. J.Daniel Gezelter for his guidance on this +project. I would also like to acknowledge the following for their help +and discussions during this project: Christopher Fennell, Charles +Vardeman, Teng Lin, Megan Sprague, Patrick Conforti, and Dan +Combest. Funding for this project came from the National Science +Foundation. + +\pagebreak +\bibliographystyle{achemso} +\bibliography{canidacy_paper} +\end{document}