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.
%review of coarse-grained modeling
Several schemes are proposed in this chapter to overcome these
difficulties.

\section{\label{lipidSection:model}Model}

\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,Fennel2004} 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]{}
\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[Electrostatic potential due to a linear molecule comprising
two point charges]{Electrostatic potential due to a linear molecule
comprising two point 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.
%description of the comparsion
\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\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{$\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{petrache00} (black) near 300~K.}
\label{lipidFigure:Scd}
\end{figure}

%\subsection{Bilayer Simulations Using Gay-Berne Ellipsoid Model}
Revision:	2806
Committed:	Tue Jun 6 20:49:05 2006 UTC (19 years, 4 months ago) by tim
Content type:	application/x-tex
File size:	17937 byte(s)
Log Message:	change the size of figure
#	User	Rev	Content
1	tim	2685	\chapter{\label{chapt:lipid}LIPID MODELING}
2
3			\section{\label{lipidSection:introduction}Introduction}
4
5	tim	2731	Under biologically relevant conditions, phospholipids are solvated
6			in aqueous solutions at a roughly 25:1 ratio. Solvation can have a
7			tremendous impact on transport phenomena in biological membranes
8			since it can affect the dynamics of ions and molecules that are
9			transferred across membranes. Studies suggest that because of the
10			directional hydrogen bonding ability of the lipid headgroups, a
11			small number of water molecules are strongly held around the
12			different parts of the headgroup and are oriented by them with
13			residence times for the first hydration shell being around 0.5 - 1
14	tim	2786	ns\cite{Ho1992}. In the second solvation shell, some water molecules
15			are weakly bound, but are still essential for determining the
16			properties of the system. Transport of various molecular species
17			into living cells is one of the major functions of membranes. A
18			thorough understanding of the underlying molecular mechanism for
19			solute diffusion is crucial to the further studies of other related
20			biological processes. All transport across cell membranes takes
21			place by one of two fundamental processes: Passive transport is
22			driven by bulk or inter-diffusion of the molecules being transported
23			or by membrane pores which facilitate crossing. Active transport
24			depends upon the expenditure of cellular energy in the form of ATP
25			hydrolysis. As the central processes of membrane assembly,
26			translocation of phospholipids across membrane bilayers requires the
27			hydrophilic head of the phospholipid to pass through the highly
28			hydrophobic interior of the membrane, and for the hydrophobic tails
29			to be exposed to the aqueous environment\cite{Sasaki2004}. A number
30			of studies indicate that the flipping of phospholipids occurs
31			rapidly in the eukaryotic ER and the bacterial cytoplasmic membrane
32			via a bi-directional, facilitated diffusion process requiring no
33			metabolic energy input. Another system of interest would be the
34			distribution of sites occupied by inhaled anesthetics in membrane.
35			Although the physiological effects of anesthetics have been
36			extensively studied, the controversy over their effects on lipid
37			bilayers still continues. Recent deuterium NMR measurements on
38			halothane in POPC lipid bilayers suggest the anesthetics are
39			primarily located at the hydrocarbon chain region\cite{Baber1995}.
40	tim	2776	Infrared spectroscopy experiments suggest that halothane in DMPC
41	tim	2786	lipid bilayers lives near the membrane/water
42			interface\cite{Lieb1982}.
43	tim	2731
44	tim	2776	Molecular dynamics simulations have proven to be a powerful tool for
45			studying the functions of biological systems, providing structural,
46			thermodynamic and dynamical information. Unfortunately, much of
47			biological interest happens on time and length scales well beyond
48	tim	2786	the range of current simulation technologies.
49			%review of coarse-grained modeling
50			Several schemes are proposed in this chapter to overcome these
51			difficulties.
52	tim	2731
53	tim	2685	\section{\label{lipidSection:model}Model}
54
55	tim	2776	\subsection{\label{lipidSection:SSD}The Soft Sticky Dipole Water Model}
56	tim	2685
57	tim	2776	In a typical bilayer simulation, the dominant portion of the
58			computation time will be spent calculating water-water interactions.
59			As an efficient solvent model, the Soft Sticky Dipole (SSD) water
60	tim	2786	model\cite{Chandra1999,Fennel2004} is used as the explicit solvent
61			in this project. Unlike other water models which have partial
62			charges distributed throughout the whole molecule, the SSD water
63			model consists of a single site which is a Lennard-Jones interaction
64			site, as well as a point dipole. A tetrahedral potential is added to
65			correct for hydrogen bonding. The following equation describes the
66			interaction between two water molecules:
67	tim	2776	\[
68			V_{SSD} = V_{LJ} (r_{ij} ) + V_{dp} (r_{ij} ,\Omega _i ,\Omega _j )
69			+ V_{sticky} (r_{ij} ,\Omega _i ,\Omega _j )
70			\]
71			where $r_{ij}$ is the vector between molecule $i$ and molecule $j$,
72			$\Omega _i$ and $\Omega _j$ are the orientational degrees of freedom
73			for molecule $i$ and molecule $j$ respectively.
74			\[
75			V_{LJ} (r_{ij} ) = 4\varepsilon _{ij} \left[ {\left( {\frac{{\sigma
76			_{ij} }}{{r_{ij} }}} \right)^{12} - \left( {\frac{{\sigma _{ij}
77			}}{{r_{ij} }}} \right)^6 } \right]
78			\]
79			\[
80			V_{dp} (r_{ij} ,\Omega _i ,\Omega _j ) = \frac{1}{{4\pi \varepsilon
81			_0 }}\left[ {\frac{{\mu _i \cdot \mu _j }}{{r_{ij}^3 }} -
82			\frac{{3\left( {\mu _i \cdot r_{ij} } \right)\left( {\mu _i \cdot
83			r_{ij} } \right)}}{{r_{ij}^5 }}} \right]
84			\]
85			\[
86			V_{sticky} (r_{ij} ,\Omega _i ,\Omega _j ) = v_0 [s(r_{ij} )w(r_{ij}
87			,\Omega _i ,\Omega _j ) + s'(r_{ij} )w'(r_{ij} ,\Omega _i ,\Omega _j
88			)]
89			\]
90			where $v_0$ is a strength parameter, $s$ and $s'$ are cubic
91			switching functions, while $w$ and $w'$ are responsible for the
92			tetrahedral potential and the short-range correction to the dipolar
93			interaction respectively.
94			\[
95			\begin{array}{l}
96			w(r_{ij} ,\Omega _i ,\Omega _j ) = \sin \theta _{ij} \sin 2\theta _{ij} \cos 2\phi _{ij} \\
97			w'(r_{ij} ,\Omega _i ,\Omega _j ) = (\cos \theta _{ij} - 0.6)^2 (\cos \theta _{ij} + 0.8)^2 - w_0 \\
98			\end{array}
99			\]
100			Although dipole-dipole and sticky interactions are more
101			mathematically complicated than Coulomb interactions, the number of
102			pair interactions is reduced dramatically both because the model
103			only contains a single-point as well as "short range" nature of the
104			higher order interaction.
105
106			\subsection{\label{lipidSection:coarseGrained}The Coarse-Grained Phospholipid Model}
107
108	tim	2781	Fig.~\ref{lipidFigure:coarseGrained} shows a schematic for our
109			coarse-grained phospholipid model. The lipid head group is modeled
110			by a linear rigid body which consists of three Lennard-Jones spheres
111			and a centrally located point-dipole. The backbone atoms in the
112			glycerol motif are modeled by Lennard-Jones spheres with dipoles.
113			Alkyl groups in hydrocarbon chains are replaced with unified
114			$\text{{\sc CH}}_2$ or $\text{{\sc CH}}_3$ atoms.
115	tim	2776
116	tim	2781	\begin{figure}
117			\centering
118	tim	2806	\includegraphics[width=3in]{coarse_grained.eps}
119	tim	2781	\caption[A representation of coarse-grained phospholipid model]{}
120			\label{lipidFigure:coarseGrained}
121			\end{figure}
122
123	tim	2776	Accurate and efficient computation of electrostatics is one of the
124			most difficult tasks in molecular modeling. Traditionally, the
125			electrostatic interaction between two molecular species is
126			calculated as a sum of interactions between pairs of point charges,
127			using Coulomb's law:
128			\[
129			V = \sum\limits_{i = 1}^{N_A } {\sum\limits_{j = 1}^{N_B }
130			{\frac{{q_i q_j }}{{4\pi \varepsilon _0 r_{ij} }}} }
131			\]
132			where $N_A$ and $N_B$ are the number of point charges in the two
133			molecular species. Originally developed to study ionic crystals, the
134			Ewald summation method mathematically transforms this
135			straightforward but conditionally convergent summation into two more
136			complicated but rapidly convergent sums. One summation is carried
137			out in reciprocal space while the other is carried out in real
138			space. An alternative approach is a multipole expansion, which is
139			based on electrostatic moments, such as charge (monopole), dipole,
140			quadruple etc.
141
142			Here we consider a linear molecule which consists of two point
143			charges $\pm q$ located at $ ( \pm \frac{d}{2},0)$. The
144			electrostatic potential at point $P$ is given by:
145			\[
146			\frac{1}{{4\pi \varepsilon _0 }}\left( {\frac{{ - q}}{{r_ - }} +
147			\frac{{ + q}}{{r_ + }}} \right) = \frac{1}{{4\pi \varepsilon _0
148			}}\left( {\frac{{ - q}}{{\sqrt {r^2 + \frac{{d^2 }}{4} + rd\cos
149			\theta } }} + \frac{q}{{\sqrt {r^2 + \frac{{d^2 }}{4} - rd\cos
150			\theta } }}} \right)
151			\]
152
153	tim	2781	\begin{figure}
154			\centering
155	tim	2806	\includegraphics[width=3in]{charge_dipole.eps}
156	tim	2781	\caption[Electrostatic potential due to a linear molecule comprising
157			two point charges]{Electrostatic potential due to a linear molecule
158			comprising two point charges} \label{lipidFigure:chargeDipole}
159			\end{figure}
160
161	tim	2776	The basic assumption of the multipole expansion is $r \gg d$ , thus,
162			$\frac{{d^2 }}{4}$ inside the square root of the denominator is
163			neglected. This is a reasonable approximation in most cases.
164			Unfortunately, in our headgroup model, the distance of charge
165			separation $d$ (4.63 $\AA$ in PC headgroup) may be comparable to
166			$r$. Nevertheless, we could still assume $ \cos \theta \approx 0$
167			in the central region of the headgroup. Using Taylor expansion and
168			associating appropriate terms with electric moments will result in a
169			"split-dipole" approximation:
170			\[
171			V(r) = \frac{1}{{4\pi \varepsilon _0 }}\frac{{r\mu \cos \theta
172			}}{{R^3 }}
173			\]
174			where$R = \sqrt {r^2 + \frac{{d^2 }}{4}}$ Placing a dipole at point
175			$P$ and applying the same strategy, the interaction between two
176			split-dipoles is then given by:
177			\[
178			V_{sd} (r_{ij} ,\Omega _i ,\Omega _j ) = \frac{1}{{4\pi \varepsilon
179			_0 }}\left[ {\frac{{\mu _i \cdot \mu _j }}{{R_{ij}^3 }} -
180			\frac{{3\left( {\mu _i \cdot r_{ij} } \right)\left( {\mu _i \cdot
181			r_{ij} } \right)}}{{R_{ij}^5 }}} \right]
182			\]
183			where $\mu _i$ and $\mu _j$ are the dipole moment of molecule $i$
184			and molecule $j$ respectively, $r_{ij}$ is vector between molecule
185	tim	2800	$i$ and molecule $j$, and $R_{ij}$ is given by,
186	tim	2776	\[
187			R_{ij} = \sqrt {r_{ij}^2 + \frac{{d_i^2 }}{4} + \frac{{d_j^2
188			}}{4}}
189			\]
190			where $d_i$ and $d_j$ are the charge separation distance of dipole
191			and respectively. This approximation to the multipole expansion
192			maintains the fast fall-off of the multipole potentials but lacks
193			the normal divergences when two polar groups get close to one
194			another.
195	tim	2786	%description of the comparsion
196	tim	2776	\begin{figure}
197			\centering
198			\includegraphics[width=\linewidth]{split_dipole.eps}
199			\caption[Comparison between electrostatic approximation]{Electron
200			density profile along the bilayer normal.}
201			\label{lipidFigure:splitDipole}
202			\end{figure}
203
204			%\section{\label{lipidSection:methods}Methods}
205
206	tim	2685	\section{\label{lipidSection:resultDiscussion}Results and Discussion}
207	tim	2776
208			\subsection{One Lipid in Sea of Water Molecules}
209
210			To exclude the inter-headgroup interaction, atomistic models of one
211			lipid (DMPC or DLPE) in sea of water molecules (TIP3P) were built
212			and studied using atomistic molecular dynamics. The simulation was
213			analyzed using a set of radial distribution functions, which give
214			the probability of finding a pair of molecular species a distance
215			apart, relative to the probability expected for a completely random
216			distribution function at the same density.
217
218			\begin{equation}
219			g_{AB} (r) = \frac{1}{{\rho _B }}\frac{1}{{N_A }} < \sum\limits_{i
220			\in A} {\sum\limits_{j \in B} {\delta (r - r_{ij} )} } >
221			\end{equation}
222			\begin{equation}
223			g_{AB} (r,\cos \theta ) = \frac{1}{{\rho _B }}\frac{1}{{N_A }} <
224			\sum\limits_{i \in A} {\sum\limits_{j \in B} {\delta (r - r_{ij} )}
225			} \delta (\cos \theta _{ij} - \cos \theta ) >
226			\end{equation}
227
228			From figure 4(a), we can identify the first solvation shell (3.5
229			$\AA$) and the second solvation shell (5.0 $\AA$) from both plots.
230			However, the corresponding orientations are different. In DLPE,
231			water molecules prefer to sit around -NH3 group due to the hydrogen
232			bonding. In contrast, because of the hydrophobic effect of the
233			-N(CH3)3 group, the preferred position of water molecules in DMPC is
234			around the -PO4 Group. When the water molecules are far from the
235			headgroup, the distribution of the two angles should be uniform. The
236			correlation close to center of the headgroup dipole (< 5 $\AA$) in
237			both plots tell us that in the closely-bound region, the dipoles of
238			the water molecules are preferentially anti-aligned with the dipole
239			of headgroup.
240
241			\begin{figure}
242			\centering
243			\includegraphics[width=\linewidth]{g_atom.eps}
244			\caption[The pair correlation functions for atomistic models]{}
245			\label{lipidFigure:PCFAtom}
246			\end{figure}
247
248			The initial configurations of coarse-grained systems are constructed
249			from the previous atomistic ones. The parameters for the
250			coarse-grained model in Table~\ref{lipidTable:parameter} are
251			estimated and tuned using isothermal-isobaric molecular dynamics.
252			Pair distribution functions calculated from coarse-grained models
253			preserve the basic characteristics of the atomistic simulations. The
254			water density, measured in a head-group-fixed reference frame,
255			surrounding two phospholipid headgroups is shown in
256			Fig.~\ref{lipidFigure:PCFCoarse}. It is clear that the phosphate end
257			in DMPC and the amine end in DMPE are the two most heavily solvated
258			atoms.
259
260			\begin{figure}
261			\centering
262			\includegraphics[width=\linewidth]{g_coarse.eps}
263			\caption[The pair correlation functions for coarse-grained models]{}
264			\label{lipidFigure:PCFCoarse}
265			\end{figure}
266
267			\begin{figure}
268			\centering
269			\includegraphics[width=\linewidth]{EWD_coarse.eps}
270			\caption[Excess water density of coarse-grained phospholipids]{ }
271			\label{lipidFigure:EWDCoarse}
272			\end{figure}
273
274			\begin{table}
275			\caption{The Parameters For Coarse-grained Phospholipids}
276			\label{lipidTable:parameter}
277			\begin{center}
278			\begin{tabular}{\|l\|c\|c\|c\|c\|c\|}
279			\hline
280			% after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
281			Atom type & Mass & $\sigma$ & $\epsilon$ & charge & Dipole \\
282			$\text{{\sc CH}}_2$ & 14.03 & 3.95 & 0.0914 & 0 & 0 \\
283			$\text{{\sc CH}}_3$ & 15.04 & 3.75 & 0.195 & 0 & 0 \\
284	tim	2781	$\text{{\sc CE}}$ & 28.01 & 3.427& 0.294 & 0 & 1.693 \\
285	tim	2776	$\text{{\sc CK}}$ & 28.01 & 3.592& 0.311 & 0 & 2.478 \\
286			$\text{{\sc PO}}_4$ & 108.995& 3.9 & 1.88 & -1& 0 \\
287			$\text{{\sc HDP}}$ & 14.03 & 4.0 & 0.13 & 0 & 0 \\
288			$\text{{\sc NC}}_4$ & 73.137 & 4.9 & 0.88 & +1& 0 \\
289			$\text{{\sc NH}}_3$ & 17.03 & 3.25 & 0.17 & +1& 0\\
290			\hline
291			\end{tabular}
292			\end{center}
293			\end{table}
294
295			\subsection{Bilayer Simulations Using Coarse-grained Model}
296
297			A bilayer system consisting of 128 DMPC lipids and 3655 water
298			molecules has been constructed from an atomistic coordinate
299			file.[15] The MD simulation is performed at constant temperature, T
300			= 300K, and constant pressure, p = 1 atm, and consisted of an
301			equilibration period of 2 ns. During the equilibration period, the
302			system was initially simulated at constant volume for 1ns. Once the
303			system was equilibrated at constant volume, the cell dimensions of
304			the system was relaxed by performing under NPT conditions using
305			Nos��-Hoover extended system isothermal-isobaric dynamics. After
306			equilibration, different properties were evaluated over a production
307			run of 5 ns.
308
309			\begin{figure}
310			\centering
311			\includegraphics[width=\linewidth]{bilayer.eps}
312			\caption[Image of a coarse-grained bilayer system]{A coarse-grained
313			bilayer system consisting of 128 DMPC lipids and 3655 SSD water
314			molecules.}
315			\label{lipidFigure:bilayer}
316			\end{figure}
317
318			\subsubsection{Electron Density Profile (EDP)}
319
320			Assuming a gaussian distribution of electrons on each atomic center
321			with a variance estimated from the size of the van der Waals radius,
322			the EDPs which are proportional to the density profiles measured
323			along the bilayer normal obtained by x-ray scattering experiment,
324	tim	2786	can be expressed by\cite{Tu1995}
325	tim	2776	\begin{equation}
326			\rho _{x - ray} (z)dz \propto \sum\limits_{i = 1}^N {\frac{{n_i
327			}}{V}\frac{1}{{\sqrt {2\pi \sigma ^2 } }}e^{ - (z - z_i )^2 /2\sigma
328			^2 } dz},
329			\end{equation}
330			where $\sigma$ is the variance equal to the van der Waals radius,
331			$n_i$ is the atomic number of site $i$ and $V$ is the volume of the
332			slab between $z$ and $z+dz$ . The highest density of total EDP
333			appears at the position of lipid-water interface corresponding to
334			headgroup, glycerol, and carbonyl groups of the lipids and the
335			distribution of water locked near the head groups, while the lowest
336			electron density is in the hydrocarbon region. As a good
337			approximation to the thickness of the bilayer, the headgroup spacing
338			, is defined as the distance between two peaks in the electron
339			density profile, calculated from our simulations to be 34.1 $\AA$.
340			This value is close to the x-ray diffraction experimental value 34.4
341	tim	2786	$\AA$\cite{Petrache1998}.
342	tim	2776
343			\begin{figure}
344			\centering
345			\includegraphics[width=\linewidth]{electron_density.eps}
346			\caption[The density profile of the lipid bilayers]{Electron density
347			profile along the bilayer normal. The water density is shown in red,
348			the density due to the headgroups in green, the glycerol backbone in
349			brown, $\text{{\sc CH}}_2$ in yellow, $\text{{\sc CH}}_3$ in cyan,
350			and total density due to DMPC in blue.}
351			\label{lipidFigure:electronDensity}
352			\end{figure}
353
354			\subsubsection{$\text{S}_{\text{{\sc cd}}}$ Order Parameter}
355
356			Measuring deuterium order parameters by NMR is a useful technique to
357			study the orientation of hydrocarbon chains in phospholipids. The
358			order parameter tensor $S$ is defined by:
359			\begin{equation}
360			S_{ij} = \frac{1}{2} < 3\cos \theta _i \cos \theta _j - \delta
361			_{ij} >
362			\end{equation}
363			where $\theta$ is the angle between the $i$th molecular axis and
364			the bilayer normal ($z$ axis). The brackets denote an average over
365			time and molecules. The molecular axes are defined:
366			\begin{itemize}
367			\item $\mathbf{\hat{z}}$ is the vector from $C_{n-1}$ to $C_{n+1}$.
368			\item $\mathbf{\hat{y}}$ is the vector that is perpendicular to $z$ and
369			in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$.
370			\item $\mathbf{\hat{x}}$ is the vector perpendicular to
371			$\mathbf{\hat{y}}$ and $\mathbf{\hat{z}}$.
372			\end{itemize}
373			In coarse-grained model, although there are no explicit hydrogens,
374			the order parameter can still be written in terms of carbon ordering
375	tim	2786	at each point of the chain\cite{Egberts1988}
376	tim	2776	\begin{equation}
377			S_{ij} = \frac{1}{2} < 3\cos \theta _i \cos \theta _j - \delta
378			_{ij} >.
379			\end{equation}
380
381			Fig.~\ref{lipidFigure:Scd} shows the order parameter profile
382			calculated for our coarse-grained DMPC bilayer system at 300K. Also
383	tim	2786	shown are the experimental data of Tu\cite{Tu1995}. The fact that
384	tim	2776	$\text{S}_{\text{{\sc cd}}}$ order parameters calculated from
385			simulation are higher than the experimental ones is ascribed to the
386			assumption of the locations of implicit hydrogen atoms which are
387			fixed in coarse-grained models at positions relative to the CC
388			vector.
389
390			\begin{figure}
391			\centering
392			\includegraphics[width=\linewidth]{scd.eps}
393			\caption[$\text{S}_{\text{{\sc cd}}}$ order parameter]{A comparison
394			of $\|\text{S}_{\text{{\sc cd}}}\|$ between coarse-grained model
395			(blue) and DMPC\cite{petrache00} (black) near 300~K.}
396			\label{lipidFigure:Scd}
397			\end{figure}
398	tim	2781
399			%\subsection{Bilayer Simulations Using Gay-Berne Ellipsoid Model}