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 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\cite{Baber1995}.
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
proposed 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:
\[
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}

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}}$ (read) 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 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)
\]

\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 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=3in]{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: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 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 tell 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}{|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.
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��-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 ^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$\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$ 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\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 are 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 result from all-atom simulations, we developed a simple
coarse-grained model capturing 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 which
is shown to have reasonable electron density profile,
$\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
phosphorlipids 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
gaining by introducing 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:	2847
Committed:	Fri Jun 9 19:25:07 2006 UTC (18 years ago) by tim
Content type:	application/x-tex
File size:	20860 byte(s)
Log Message:	minor fix