--- trunk/tengDissertation/Lipid.tex	2006/05/25 20:31:10	2775
+++ trunk/tengDissertation/Lipid.tex	2006/05/25 21:32:14	2776
@@ -11,38 +11,368 @@ ns.[14] In the second solvation shell, some water mole
 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.[14] 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.[18] 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 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.[16] Infrared spectroscopy experiments
-suggest that halothane in DMPC lipid bilayers lives near the
-membrane/water interface. [17]
+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
+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.
 
+Molecular dynamics simulations have proven to be a powerful tool for
+studying the functions of biological systems, providing structural,
+thermodynamic and dynamical information. Unfortunately, much of
+biological interest happens on time and length scales well beyond
+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:methods}Methods}
+\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:
+\[
+V_{SSD}  = V_{LJ} (r_{ij} ) + V_{dp} (r_{ij} ,\Omega _i ,\Omega _j )
++ V_{sticky} (r_{ij} ,\Omega _i ,\Omega _j )
+\]
+where $r_{ij}$ is the vector between molecule $i$ and molecule $j$,
+$\Omega _i$ and $\Omega _j$ are the orientational degrees of freedom
+for molecule $i$ and molecule $j$ respectively.
+\[
+V_{LJ} (r_{ij} ) = 4\varepsilon _{ij} \left[ {\left( {\frac{{\sigma
+_{ij} }}{{r_{ij} }}} \right)^{12}  - \left( {\frac{{\sigma _{ij}
+}}{{r_{ij} }}} \right)^6 } \right]
+\]
+\[
+V_{dp} (r_{ij} ,\Omega _i ,\Omega _j ) = \frac{1}{{4\pi \varepsilon
+_0 }}\left[ {\frac{{\mu _i  \cdot \mu _j }}{{r_{ij}^3 }} -
+\frac{{3\left( {\mu _i  \cdot r_{ij} } \right)\left( {\mu _i  \cdot
+r_{ij} } \right)}}{{r_{ij}^5 }}} \right]
+\]
+\[
+V_{sticky} (r_{ij} ,\Omega _i ,\Omega _j ) = v_0 [s(r_{ij} )w(r_{ij}
+,\Omega _i ,\Omega _j ) + s'(r_{ij} )w'(r_{ij} ,\Omega _i ,\Omega _j
+)]
+\]
+where $v_0$ is a strength parameter, $s$ and $s'$ are cubic
+switching functions, while $w$   and $w'$  are responsible for the
+tetrahedral potential and the short-range correction to the dipolar
+interaction respectively.
+\[
+\begin{array}{l}
+ w(r_{ij} ,\Omega _i ,\Omega _j ) = \sin \theta _{ij} \sin 2\theta _{ij} \cos 2\phi _{ij}  \\
+ 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
+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.
+
+\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.
+
+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
+calculated as a sum of interactions between pairs of point charges,
+using Coulomb's law:
+\[
+V = \sum\limits_{i = 1}^{N_A } {\sum\limits_{j = 1}^{N_B }
+{\frac{{q_i q_j }}{{4\pi \varepsilon _0 r_{ij} }}} }
+\]
+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.
+
+Here we consider a linear molecule which consists of two point
+charges $\pm q$ located at $ ( \pm \frac{d}{2},0)$. The
+electrostatic potential at point $P$ is given by:
+\[
+\frac{1}{{4\pi \varepsilon _0 }}\left( {\frac{{ - q}}{{r_ -  }} +
+\frac{{ + q}}{{r_ +  }}} \right) = \frac{1}{{4\pi \varepsilon _0
+}}\left( {\frac{{ - q}}{{\sqrt {r^2  + \frac{{d^2 }}{4} + rd\cos
+\theta } }} + \frac{q}{{\sqrt {r^2  + \frac{{d^2 }}{4} - rd\cos
+\theta } }}} \right)
+\]
+
+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.
+Unfortunately, in our headgroup model, the distance of charge
+separation $d$ (4.63 $\AA$ in PC headgroup) may be comparable to
+$r$. Nevertheless, we could still assume  $ \cos \theta  \approx 0$
+in the central region of the headgroup. Using Taylor expansion and
+associating appropriate terms with electric moments will result in a
+"split-dipole" approximation:
+\[
+V(r) = \frac{1}{{4\pi \varepsilon _0 }}\frac{{r\mu \cos \theta
+}}{{R^3 }}
+\]
+where$R = \sqrt {r^2  + \frac{{d^2 }}{4}}$ Placing a dipole at point
+$P$ and applying the same strategy, the interaction between two
+split-dipoles is then given by:
+\[
+V_{sd} (r_{ij} ,\Omega _i ,\Omega _j ) = \frac{1}{{4\pi \varepsilon
+_0 }}\left[ {\frac{{\mu _i  \cdot \mu _j }}{{R_{ij}^3 }} -
+\frac{{3\left( {\mu _i  \cdot r_{ij} } \right)\left( {\mu _i  \cdot
+r_{ij} } \right)}}{{R_{ij}^5 }}} \right]
+\]
+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,
+\[
+R_{ij}  = \sqrt {r_{ij}^2  + \frac{{d_i^2 }}{4} + \frac{{d_j^2
+}}{4}}
+\]
+where $d_i$ and $d_j$ are the charge separation distance of dipole
+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.
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{split_dipole.eps}
+\caption[Comparison between electrostatic approximation]{Electron
+density profile along the bilayer normal.}
+\label{lipidFigure:splitDipole}
+\end{figure}
+
+%\section{\label{lipidSection:methods}Methods}
+
 \section{\label{lipidSection:resultDiscussion}Results and Discussion}
+
+\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.
+
+\begin{equation}
+g_{AB} (r) = \frac{1}{{\rho _B }}\frac{1}{{N_A }} < \sum\limits_{i
+\in A} {\sum\limits_{j \in B} {\delta (r - r_{ij} )} }  >
+\end{equation}
+\begin{equation}
+g_{AB} (r,\cos \theta ) = \frac{1}{{\rho _B }}\frac{1}{{N_A }} <
+\sum\limits_{i \in A} {\sum\limits_{j \in B} {\delta (r - r_{ij} )}
+} \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.
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{g_atom.eps}
+\caption[The pair correlation functions for atomistic models]{}
+\label{lipidFigure:PCFAtom}
+\end{figure}
+
+The initial configurations of coarse-grained systems are constructed
+from the previous atomistic ones. The parameters for the
+coarse-grained model in Table~\ref{lipidTable:parameter} are
+estimated and tuned using isothermal-isobaric molecular dynamics.
+Pair distribution functions calculated from coarse-grained models
+preserve the basic characteristics of the atomistic simulations. The
+water density, measured in a head-group-fixed reference frame,
+surrounding two phospholipid headgroups is shown in
+Fig.~\ref{lipidFigure:PCFCoarse}. It is clear that the phosphate end
+in DMPC and the amine end in DMPE are the two most heavily solvated
+atoms.
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{g_coarse.eps}
+\caption[The pair correlation functions for coarse-grained models]{}
+\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}
+\end{figure}
+
+\begin{table}
+\caption{The Parameters For Coarse-grained Phospholipids}
+\label{lipidTable:parameter}
+\begin{center}
+\begin{tabular}{|l|c|c|c|c|c|}
+  \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 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 \\
+  $\text{{\sc NC}}_4$ & 73.137 & 4.9  & 0.88   & +1&  0 \\
+  $\text{{\sc NH}}_3$ & 17.03  & 3.25 & 0.17   & +1&  0\\
+  \hline
+\end{tabular}
+\end{center}
+\end{table}
+
+\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.
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{bilayer.eps}
+\caption[Image of a coarse-grained bilayer system]{A coarse-grained
+bilayer system consisting of 128 DMPC lipids and 3655 SSD water
+molecules.}
+\label{lipidFigure:bilayer}
+\end{figure}
+
+\subsubsection{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
+\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
+^2 } dz},
+\end{equation}
+where $\sigma$ is the variance equal to the van der Waals radius,
+$n_i$ is the atomic number of site $i$ and $V$ is the volume of the
+slab between $z$ and $z+dz$ . The highest density of total EDP
+appears at the position of lipid-water interface corresponding to
+headgroup, glycerol, and carbonyl groups of the lipids and the
+distribution of water locked near the head groups, while the lowest
+electron density is in the hydrocarbon region. As a good
+approximation to the thickness of the bilayer, the headgroup spacing
+, 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$
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{electron_density.eps}
+\caption[The density profile of the lipid bilayers]{Electron density
+profile along the bilayer normal. The water density is shown in red,
+the density due to the headgroups in green, the glycerol backbone in
+brown, $\text{{\sc CH}}_2$ in yellow, $\text{{\sc CH}}_3$ in cyan,
+and total density due to DMPC in blue.}
+\label{lipidFigure:electronDensity}
+\end{figure}
+
+\subsubsection{$\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
+order parameter tensor $S$ is defined by:
+\begin{equation}
+S_{ij}  = \frac{1}{2} < 3\cos \theta _i \cos \theta _j  - \delta
+_{ij}  >
+\end{equation}
+where $\theta$ is the angle between the  $i$th molecular axis and
+the bilayer normal ($z$ axis). The brackets denote an average over
+time and molecules. The molecular axes are defined:
+\begin{itemize}
+\item $\mathbf{\hat{z}}$ is the vector from $C_{n-1}$ to $C_{n+1}$.
+\item $\mathbf{\hat{y}}$ is the vector that is perpendicular to $z$ and
+in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$.
+\item $\mathbf{\hat{x}}$ is the vector perpendicular to
+$\mathbf{\hat{y}}$ and $\mathbf{\hat{z}}$.
+\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
+\begin{equation}
+S_{ij}  = \frac{1}{2} < 3\cos \theta _i \cos \theta _j  - \delta
+_{ij}  >.
+\end{equation}
+
+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
+$\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
+fixed in coarse-grained models at positions relative to the CC
+vector.
+
+\begin{figure}
+\centering
+\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.}
+\label{lipidFigure:Scd}
+\end{figure}