trunk/tengDissertation/Lipid.tex

\chapter{\label{chapt:lipid}LIPID MODELING}

\section{\label{lipidSection:introduction}Introduction}

Under biologically relevant conditions, phospholipids are solvated
in aqueous solutions at a roughly 25:1 ratio. Solvation can have a
tremendous impact on transport phenomena in biological membranes
since it can affect the dynamics of ions and molecules that are
transferred across membranes. Studies suggest that because of the
directional hydrogen bonding ability of the lipid headgroups, a
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\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 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,
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
introduced in this chapter to overcome these difficulties.

\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\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:
\begin{equation}
V_{SSD}  = V_{LJ} (r_{ij} ) + V_{dp} (r_{ij} ,\Omega _i ,\Omega _j )
+ V_{sticky} (r_{ij} ,\Omega _i ,\Omega _j )
\label{lipidSection:ssdEquation}
\end{equation}
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. The potential terms
in Eq.~\ref{lipidSection:ssdEquation} are given by :
\begin{eqnarray}
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{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[
\hat{u}_{i} \cdot \hat{u}_{j} - 3(\hat{u}_i \cdot \hat{\mathbf{r}}_{ij}) %
(\hat{u}_j \cdot \hat{\mathbf{r}}_{ij}) \biggr],\\
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 )]
\end{eqnarray}
where $v_0$ is a strength parameter, $s$ and $s'$ are cubic
switching functions and $w$ and $w'$  are responsible for the
tetrahedral potential and the short-range correction to the dipolar
interaction respectively:
\begin{eqnarray}
 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{eqnarray}
Here $\theta _{ij}$ and $\phi _{ij}$ are the spherical polar angles
representing relative orientations of molecule $j$ in the body-fixed
frame of molecule $i$. 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 and
because of the "short range" nature of the more expensive
interaction.

\subsection{\label{lipidSection:coarseGrained}The Coarse-Grained Phospholipid Model}

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 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
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 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
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)
\]
\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.
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. 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 have discovered 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]{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:resultDiscussion}Results and Discussion}

\subsection{One Lipid in Sea of Water Molecules}

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{eqnarray}
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} )} }  >, \\
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{eqnarray}
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]{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}

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]{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]{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}{lccccc}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  Atom type & Mass & $\sigma$ & $\epsilon$ & charge & Dipole \\
  \hline
  $\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.
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
\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{\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\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 _i ^2 } }}e^{ - (z - z_i )^2
/2\sigma _i ^2 } dz},
\end{equation}
where $\sigma _i$ 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\cite{Petrache1998}.

\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{\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
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 _i$ is the angle between the  $i$th molecular axis and
the bilayer normal ($z$ axis). The brackets denote an average over
time and lipid 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 our 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\cite{Egberts1988}
\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 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
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{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.
Revision:	2909
Committed:	Thu Jun 29 23:00:35 2006 UTC (18 years ago) by tim
Content type:	application/x-tex
File size:	21171 byte(s)
Log Message:	more corrections