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. 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}

\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}

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=\linewidth]{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=\linewidth]{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.

\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}

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