--- trunk/tengDissertation/Lipid.tex 2006/05/26 17:56:36 2778 +++ trunk/tengDissertation/Lipid.tex 2006/06/23 20:21:54 2881 @@ -11,33 +11,35 @@ ns. In the second solvation shell, some water molecule small number of water molecules are strongly held around the different parts of the headgroup and are oriented by them with residence times for the first hydration shell being around 0.5 - 1 -ns. In the second solvation shell, some water molecules are weakly -bound, but are still essential for determining the properties of the -system. Transport of various molecular species into living cells is -one of the major functions of membranes. A thorough understanding of -the underlying molecular mechanism for solute diffusion is crucial -to the further studies of other related biological processes. All -transport across cell membranes takes place by one of two -fundamental processes: Passive transport is driven by bulk or -inter-diffusion of the molecules being transported or by membrane -pores which facilitate crossing. Active transport depends upon the -expenditure of cellular energy in the form of ATP hydrolysis. As the -central processes of membrane assembly, translocation of -phospholipids across membrane bilayers requires the hydrophilic head -of the phospholipid to pass through the highly hydrophobic interior -of the membrane, and for the hydrophobic tails to be exposed to the -aqueous environment. A number of studies indicate that the flipping -of phospholipids occurs rapidly in the eukaryotic ER and the -bacterial cytoplasmic membrane via a bi-directional, facilitated -diffusion process requiring no metabolic energy input. Another -system of interest would be the distribution of sites occupied by -inhaled anesthetics in membrane. Although the physiological effects -of anesthetics have been extensively studied, the controversy over +ns\cite{Ho1992}. In the second solvation shell, some water molecules +are weakly bound, but are still essential for determining the +properties of the system. Transport of various molecular species +into living cells is one of the major functions of membranes. A +thorough understanding of the underlying molecular mechanism for +solute diffusion is crucial to the further studies of other related +biological processes. All transport across cell membranes takes +place by one of two fundamental processes: Passive transport is +driven by bulk or inter-diffusion of the molecules being transported +or by membrane pores which facilitate crossing. Active transport +depends upon the expenditure of cellular energy in the form of ATP +hydrolysis. As the central processes of membrane assembly, +translocation of phospholipids across membrane bilayers requires the +hydrophilic head of the phospholipid to pass through the highly +hydrophobic interior of the membrane, and for the hydrophobic tails +to be exposed to the aqueous environment\cite{Sasaki2004}. A number +of studies indicate that the flipping of phospholipids occurs +rapidly in the eukaryotic endoplasmic reticulum and the bacterial +cytoplasmic membrane via a bi-directional, facilitated diffusion +process requiring no metabolic energy input. Another system of +interest is the distribution of sites occupied by inhaled +anesthetics in membrane. Although the physiological effects of +anesthetics have been extensively studied, the controversy over their effects on lipid bilayers still continues. Recent deuterium -NMR measurements on halothane in POPC lipid bilayers suggest the -anesthetics are primarily located at the hydrocarbon chain region. -Infrared spectroscopy experiments suggest that halothane in DMPC -lipid bilayers lives near the membrane/water interface. +NMR measurements on halothane on POPC lipid bilayers suggest the +anesthetics are primarily located at the hydrocarbon chain +region\cite{Baber1995}. However, infrared spectroscopy experiments +suggest that halothane in DMPC lipid bilayers lives near the +membrane/water interface\cite{Lieb1982}. Molecular dynamics simulations have proven to be a powerful tool for studying the functions of biological systems, providing structural, @@ -46,20 +48,20 @@ proposed in this chapter to overcome these difficultie the range of current simulation technologies. Several schemes are proposed in this chapter to overcome these difficulties. -\section{\label{lipidSection:model}Model} +\section{\label{lipidSection:model}Model and Methodology} \subsection{\label{lipidSection:SSD}The Soft Sticky Dipole Water Model} In a typical bilayer simulation, the dominant portion of the computation time will be spent calculating water-water interactions. As an efficient solvent model, the Soft Sticky Dipole (SSD) water -model is used as the explicit solvent in this project. Unlike other -water models which have partial charges distributed throughout the -whole molecule, the SSD water model consists of a single site which -is a Lennard-Jones interaction site, as well as a point dipole. A -tetrahedral potential is added to correct for hydrogen bonding. The -following equation describes the interaction between two water -molecules: +model\cite{Chandra1999,Fennell2004} is used as the explicit solvent +in this project. Unlike other water models which have partial +charges distributed throughout the whole molecule, the SSD water +model consists of a single site which is a Lennard-Jones interaction +site, as well as a point dipole. A tetrahedral potential is added to +correct for hydrogen bonding. The following equation describes the +interaction between two water molecules: \[ V_{SSD} = V_{LJ} (r_{ij} ) + V_{dp} (r_{ij} ,\Omega _i ,\Omega _j ) + V_{sticky} (r_{ij} ,\Omega _i ,\Omega _j ) @@ -93,21 +95,35 @@ Although dipole-dipole and sticky interactions are mor w'(r_{ij} ,\Omega _i ,\Omega _j ) = (\cos \theta _{ij} - 0.6)^2 (\cos \theta _{ij} + 0.8)^2 - w_0 \\ \end{array} \] -Although dipole-dipole and sticky interactions are more +Although the dipole-dipole and sticky interactions are more mathematically complicated than Coulomb interactions, the number of pair interactions is reduced dramatically both because the model only contains a single-point as well as "short range" nature of the -higher order interaction. +more expensive interaction. \subsection{\label{lipidSection:coarseGrained}The Coarse-Grained Phospholipid Model} -Figure 1 shows a schematic for our coarse-grained phospholipid -model. The lipid head group is modeled by a linear rigid body which -consists of three Lennard-Jones spheres and a centrally located -point-dipole. The backbone atoms in the glycerol motif are modeled -by Lennard-Jones spheres with dipoles. Alkyl groups in hydrocarbon -chains are replaced with unified CH2 or CH3 atoms. +Fig.~\ref{lipidFigure:coarseGrained} shows a schematic for our +coarse-grained phospholipid model. The lipid head group is modeled +by a linear rigid body which consists of three Lennard-Jones spheres +and a centrally located point-dipole. The backbone atoms in the +glycerol motif are modeled by Lennard-Jones spheres with dipoles. +Alkyl groups in hydrocarbon chains are replaced with unified +$\text{{\sc CH}}_2$ or $\text{{\sc CH}}_3$ atoms. +\begin{figure} +\centering +\includegraphics[width=3in]{coarse_grained.eps} +\caption[A representation of coarse-grained phospholipid model]{A +representation of coarse-grained phospholipid model. The lipid +headgroup consists of $\text{{\sc PO}}_4$ group (yellow), +$\text{{\sc NC}}_4$ group (blue) and a united $\text{{\sc C}}$ atom +(gray) $ with a dipole, while the glycerol backbone includes dipolar +$\text{{\sc CE}}$ (red) and $\text{{\sc CK}}$ (pink) groups. Alkyl +groups in hydrocarbon chains are simply represented by gray united +atoms.} \label{lipidFigure:coarseGrained} +\end{figure} + Accurate and efficient computation of electrostatics is one of the most difficult tasks in molecular modeling. Traditionally, the electrostatic interaction between two molecular species is @@ -119,13 +135,12 @@ Ewald summation method mathematically transforms this \] where $N_A$ and $N_B$ are the number of point charges in the two molecular species. Originally developed to study ionic crystals, the -Ewald summation method mathematically transforms this -straightforward but conditionally convergent summation into two more -complicated but rapidly convergent sums. One summation is carried -out in reciprocal space while the other is carried out in real -space. An alternative approach is a multipole expansion, which is -based on electrostatic moments, such as charge (monopole), dipole, -quadruple etc. +Ewald sum method mathematically transforms this straightforward but +conditionally convergent summation into two more complicated but +rapidly convergent sums. One summation is carried out in reciprocal +space while the other is carried out in real space. An alternative +approach is the multipole expansion, which is based on electrostatic +moments, such as charge (monopole), dipole, quadrupole etc. Here we consider a linear molecule which consists of two point charges $\pm q$ located at $ ( \pm \frac{d}{2},0)$. The @@ -138,6 +153,15 @@ The basic assumption of the multipole expansion is $r \theta } }}} \right) \] +\begin{figure} +\centering +\includegraphics[width=3in]{charge_dipole.eps} +\caption[An illustration of split-dipole +approximation]{Electrostatic potential due to a linear molecule +comprising two point charges with opposite charges. } +\label{lipidFigure:chargeDipole} +\end{figure} + The basic assumption of the multipole expansion is $r \gg d$ , thus, $\frac{{d^2 }}{4}$ inside the square root of the denominator is neglected. This is a reasonable approximation in most cases. @@ -162,7 +186,7 @@ $i$ and molecule $j$, and $R_{ij{$ is given by, \] where $\mu _i$ and $\mu _j$ are the dipole moment of molecule $i$ and molecule $j$ respectively, $r_{ij}$ is vector between molecule -$i$ and molecule $j$, and $R_{ij{$ is given by, +$i$ and molecule $j$, and $R_{ij}$ is given by, \[ R_{ij} = \sqrt {r_{ij}^2 + \frac{{d_i^2 }}{4} + \frac{{d_j^2 }}{4}} @@ -171,14 +195,21 @@ another. and respectively. This approximation to the multipole expansion maintains the fast fall-off of the multipole potentials but lacks the normal divergences when two polar groups get close to one -another. - +another. The comparision between different electrostatic +approximation is shown in \ref{lipidFigure:splitDipole}. Due to the +divergence at the central region of the headgroup introduced by +dipole-dipole approximation, we discover that water molecules are +locked into the central region in the simulation. This artifact can +be corrected using split-dipole approximation or other accurate +methods. \begin{figure} \centering \includegraphics[width=\linewidth]{split_dipole.eps} -\caption[Comparison between electrostatic approximation]{Electron -density profile along the bilayer normal.} -\label{lipidFigure:splitDipole} +\caption[Comparison between electrostatic +approximation]{Electrostatic potential map for two pairs of charges +with different alignments: (a) illustration of different alignments; +(b) charge-charge interaction; (c) dipole-dipole approximation; (d) +split-dipole approximation.} \label{lipidFigure:splitDipole} \end{figure} %\section{\label{lipidSection:methods}Methods} @@ -187,13 +218,14 @@ To exclude the inter-headgroup interaction, atomistic \subsection{One Lipid in Sea of Water Molecules} -To exclude the inter-headgroup interaction, atomistic models of one -lipid (DMPC or DLPE) in sea of water molecules (TIP3P) were built -and studied using atomistic molecular dynamics. The simulation was -analyzed using a set of radial distribution functions, which give -the probability of finding a pair of molecular species a distance -apart, relative to the probability expected for a completely random -distribution function at the same density. +To tune our parameters without the inter-headgroup interactions, +atomistic models of one lipid (DMPC or DLPE) in sea of water +molecules (TIP3P) were built and studied using atomistic molecular +dynamics. The simulation was analyzed using a set of radial +distribution functions, which give the probability of finding a pair +of molecular species a distance apart, relative to the probability +expected for a completely random distribution function at the same +density. \begin{equation} g_{AB} (r) = \frac{1}{{\rho _B }}\frac{1}{{N_A }} < \sum\limits_{i @@ -205,23 +237,33 @@ From figure 4(a), we can identify the first solvation } \delta (\cos \theta _{ij} - \cos \theta ) > \end{equation} -From figure 4(a), we can identify the first solvation shell (3.5 -$\AA$) and the second solvation shell (5.0 $\AA$) from both plots. -However, the corresponding orientations are different. In DLPE, -water molecules prefer to sit around -NH3 group due to the hydrogen -bonding. In contrast, because of the hydrophobic effect of the --N(CH3)3 group, the preferred position of water molecules in DMPC is -around the -PO4 Group. When the water molecules are far from the -headgroup, the distribution of the two angles should be uniform. The -correlation close to center of the headgroup dipole (< 5 $\AA$) in -both plots tell us that in the closely-bound region, the dipoles of -the water molecules are preferentially anti-aligned with the dipole -of headgroup. +From Fig.~\ref{lipidFigure:PCFAtom}, we can identify the first +solvation shell (3.5 $\AA$) and the second solvation shell (5.0 +$\AA$) from both plots. However, the corresponding orientations are +different. In DLPE, water molecules prefer to sit around $\text{{\sc +NH}}_3$ group due to the hydrogen bonding. In contrast, because of +the hydrophobic effect of the $ {\rm{N(CH}}_{\rm{3}} +{\rm{)}}_{\rm{3}} $ group, the preferred position of water molecules +in DMPC is around the $\text{{\sc PO}}_4$ Group. When the water +molecules are far from the headgroup, the distribution of the two +angles should be uniform. The correlation close to center of the +headgroup dipole in both plots tells us that in the closely-bound +region, the dipoles of the water molecules are preferentially +anti-aligned with the dipole of headgroup. When the water molecules +are far from the headgroup, the distribution of the two angles +should be uniform. The correlation close to center of the headgroup +dipole in both plots tell us that in the closely-bound region, the +dipoles of the water molecules are preferentially aligned with the +dipole of headgroup. \begin{figure} \centering \includegraphics[width=\linewidth]{g_atom.eps} -\caption[The pair correlation functions for atomistic models]{} +\caption[The pair correlation functions for atomistic models]{The +pair correlation functions for atomistic models: (a)$g(r,\cos \theta +)$ for DMPC; (b) $g(r,\cos \theta )$ for DLPE; (c)$g(r,\cos \omega +)$ for DMPC; (d)$g(r,\cos \omega )$ for DLPE; (e)$g(\cos \theta,\cos +\omega)$ for DMPC; (f)$g(\cos \theta,\cos \omega)$ for DMLPE.} \label{lipidFigure:PCFAtom} \end{figure} @@ -240,15 +282,18 @@ atoms. \begin{figure} \centering \includegraphics[width=\linewidth]{g_coarse.eps} -\caption[The pair correlation functions for coarse-grained models]{} +\caption[The pair correlation functions for coarse-grained +models]{The pair correlation functions for coarse-grained models: +(a)$g(r,\cos \theta )$ for DMPC; (b) $g(r,\cos \theta )$ for DLPE.} \label{lipidFigure:PCFCoarse} \end{figure} \begin{figure} \centering \includegraphics[width=\linewidth]{EWD_coarse.eps} -\caption[Excess water density of coarse-grained phospholipids]{ } -\label{lipidFigure:EWDCoarse} +\caption[Excess water density of coarse-grained +phospholipids]{Excess water density of coarse-grained +phospholipids.} \label{lipidFigure:EWDCoarse} \end{figure} \begin{table} @@ -259,10 +304,9 @@ atoms. \hline % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ... Atom type & Mass & $\sigma$ & $\epsilon$ & charge & Dipole \\ - $\text{{\sc CH}}_2$ & 14.03 & 3.95 & 0.0914 & 0 & 0 \\ $\text{{\sc CH}}_3$ & 15.04 & 3.75 & 0.195 & 0 & 0 \\ - $\text{{\sc CE}}$ & 28.01 & 3.427& 0.294 & 0 & 1.693 + $\text{{\sc CE}}$ & 28.01 & 3.427& 0.294 & 0 & 1.693 \\ $\text{{\sc CK}}$ & 28.01 & 3.592& 0.311 & 0 & 2.478 \\ $\text{{\sc PO}}_4$ & 108.995& 3.9 & 1.88 & -1& 0 \\ $\text{{\sc HDP}}$ & 14.03 & 4.0 & 0.13 & 0 & 0 \\ @@ -276,16 +320,15 @@ molecules has been constructed from an atomistic coord \subsection{Bilayer Simulations Using Coarse-grained Model} A bilayer system consisting of 128 DMPC lipids and 3655 water -molecules has been constructed from an atomistic coordinate -file.[15] The MD simulation is performed at constant temperature, T -= 300K, and constant pressure, p = 1 atm, and consisted of an -equilibration period of 2 ns. During the equilibration period, the -system was initially simulated at constant volume for 1ns. Once the -system was equilibrated at constant volume, the cell dimensions of -the system was relaxed by performing under NPT conditions using -Nos¨¦-Hoover extended system isothermal-isobaric dynamics. After -equilibration, different properties were evaluated over a production -run of 5 ns. +molecules has been constructed from an atomistic coordinate file. +The MD simulation is performed at constant temperature, T = 300K, +and constant pressure, p = 1 atm, and consisted of an equilibration +period of 2 ns. During the equilibration period, the system was +initially simulated at constant volume for 1 ns. Once the system was +equilibrated at constant volume, the cell dimensions of the system +was relaxed by performing under NPT conditions using Nos\'{e}-Hoover +extended system isothermal-isobaric dynamics. After equilibration, +different properties were evaluated over a production run of 5 ns. \begin{figure} \centering @@ -296,13 +339,13 @@ molecules.} \label{lipidFigure:bilayer} \end{figure} -\subsubsection{Electron Density Profile (EDP)} +\subsubsection{\textbf{Electron Density Profile (EDP)}} Assuming a gaussian distribution of electrons on each atomic center with a variance estimated from the size of the van der Waals radius, the EDPs which are proportional to the density profiles measured along the bilayer normal obtained by x-ray scattering experiment, -can be expressed by +can be expressed by\cite{Tu1995} \begin{equation} \rho _{x - ray} (z)dz \propto \sum\limits_{i = 1}^N {\frac{{n_i }}{V}\frac{1}{{\sqrt {2\pi \sigma ^2 } }}e^{ - (z - z_i )^2 /2\sigma @@ -319,7 +362,7 @@ $\AA$. , is defined as the distance between two peaks in the electron density profile, calculated from our simulations to be 34.1 $\AA$. This value is close to the x-ray diffraction experimental value 34.4 -$\AA$. +$\AA$\cite{Petrache1998}. \begin{figure} \centering @@ -332,7 +375,7 @@ and total density due to DMPC in blue.} \label{lipidFigure:electronDensity} \end{figure} -\subsubsection{$\text{S}_{\text{{\sc cd}}}$ Order Parameter} +\subsubsection{\textbf{$\text{S}_{\text{{\sc cd}}}$ Order Parameter}} Measuring deuterium order parameters by NMR is a useful technique to study the orientation of hydrocarbon chains in phospholipids. The @@ -353,7 +396,7 @@ at each point of the chain \end{itemize} In coarse-grained model, although there are no explicit hydrogens, the order parameter can still be written in terms of carbon ordering -at each point of the chain +at each point of the chain\cite{Egberts1988} \begin{equation} S_{ij} = \frac{1}{2} < 3\cos \theta _i \cos \theta _j - \delta _{ij} >. @@ -361,7 +404,7 @@ shown are the experimental data of Tiburu. The fact th Fig.~\ref{lipidFigure:Scd} shows the order parameter profile calculated for our coarse-grained DMPC bilayer system at 300K. Also -shown are the experimental data of Tiburu. The fact that +shown are the experimental data of Tu\cite{Tu1995}. The fact that $\text{S}_{\text{{\sc cd}}}$ order parameters calculated from simulation are higher than the experimental ones is ascribed to the assumption of the locations of implicit hydrogen atoms which are @@ -373,6 +416,29 @@ of $|\text{S}_{\text{{\sc cd}}}|$ between coarse-grain \includegraphics[width=\linewidth]{scd.eps} \caption[$\text{S}_{\text{{\sc cd}}}$ order parameter]{A comparison of $|\text{S}_{\text{{\sc cd}}}|$ between coarse-grained model -(blue) and DMPC\cite{petrache00} (black) near 300~K.} +(blue) and DMPC\cite{Petrache2000} (black) near 300~K.} \label{lipidFigure:Scd} \end{figure} + +%\subsection{Bilayer Simulations Using Gay-Berne Ellipsoid Model} + +\section{\label{lipidSection:Conclusion}Conclusion} + +Atomistic simulations have been used in this study to determine the +preferred orientation and location of water molecules relative to +the location and orientation of the PC and PE lipid headgroups. +Based on the results from our all-atom simulations, we developed a +simple coarse-grained model which captures the essential features of +the headgroup solvation which is crucial to transport process in +membrane system. In addition, the model has been explored in a +bilayer system was shown to have reasonable electron density +profiles, $\text{S}_{\text{{\sc cd}}}$ order parameter and other +structural properties. The accuracy of this model is achieved by +matching atomistic result. It is also easy to represent different +phospholipids by changing a few parameters of the model. Another +important characteristic of this model distinguishing itself from +other models\cite{Goetz1998,Marrink2004} is the computational speed +gained by introducing a short range electrostatic approximation. +Further studies of this system using z-constraint method could +extend the time length of the simulations to study transport +phenomena in large-scale membrane systems.