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}

Figure 1 shows a schematic for our coarse-grained phospholipid
model. The lipid head group is modeled by a linear rigid body which
consists of three Lennard-Jones spheres and a centrally located
point-dipole. The backbone atoms in the glycerol motif are modeled
by Lennard-Jones spheres with dipoles. Alkyl groups in hydrocarbon
chains are replaced with unified CH2 or CH3 atoms.

Accurate and efficient computation of electrostatics is one of the
most difficult tasks in molecular modeling. Traditionally, the
electrostatic interaction between two molecular species is
calculated as a sum of interactions between pairs of point charges,
using Coulomb's law:
\[
V = \sum\limits_{i = 1}^{N_A } {\sum\limits_{j = 1}^{N_B }
{\frac{{q_i q_j }}{{4\pi \varepsilon _0 r_{ij} }}} }
\]
where $N_A$ and $N_B$ are the number of point charges in the two
molecular species. Originally developed to study ionic crystals, the
Ewald summation method mathematically transforms this
straightforward but conditionally convergent summation into two more
complicated but rapidly convergent sums. One summation is carried
out in reciprocal space while the other is carried out in real
space. An alternative approach is a multipole expansion, which is
based on electrostatic moments, such as charge (monopole), dipole,
quadruple etc.

Here we consider a linear molecule which consists of two point
charges $\pm q$ located at $ ( \pm \frac{d}{2},0)$. The
electrostatic potential at point $P$ is given by:
\[
\frac{1}{{4\pi \varepsilon _0 }}\left( {\frac{{ - q}}{{r_ -  }} +
\frac{{ + q}}{{r_ +  }}} \right) = \frac{1}{{4\pi \varepsilon _0
}}\left( {\frac{{ - q}}{{\sqrt {r^2  + \frac{{d^2 }}{4} + rd\cos
\theta } }} + \frac{q}{{\sqrt {r^2  + \frac{{d^2 }}{4} - rd\cos
\theta } }}} \right)
\]

The basic assumption of the multipole expansion is $r \gg d$ , thus,
$\frac{{d^2 }}{4}$ inside the square root of the denominator is
neglected. This is a reasonable approximation in most cases.
Unfortunately, in our headgroup model, the distance of charge
separation $d$ (4.63 $\AA$ in PC headgroup) may be comparable to
$r$. Nevertheless, we could still assume  $ \cos \theta  \approx 0$
in the central region of the headgroup. Using Taylor expansion and
associating appropriate terms with electric moments will result in a
"split-dipole" approximation:
\[
V(r) = \frac{1}{{4\pi \varepsilon _0 }}\frac{{r\mu \cos \theta
}}{{R^3 }}
\]
where$R = \sqrt {r^2  + \frac{{d^2 }}{4}}$ Placing a dipole at point
$P$ and applying the same strategy, the interaction between two
split-dipoles is then given by:
\[
V_{sd} (r_{ij} ,\Omega _i ,\Omega _j ) = \frac{1}{{4\pi \varepsilon
_0 }}\left[ {\frac{{\mu _i  \cdot \mu _j }}{{R_{ij}^3 }} -
\frac{{3\left( {\mu _i  \cdot r_{ij} } \right)\left( {\mu _i  \cdot
r_{ij} } \right)}}{{R_{ij}^5 }}} \right]
\]
where $\mu _i$  and  $\mu _j$ are the dipole moment of molecule $i$
and molecule $j$ respectively, $r_{ij}$ is vector between molecule
$i$ and molecule $j$, and $R_{ij{$ is given by,
\[
R_{ij}  = \sqrt {r_{ij}^2  + \frac{{d_i^2 }}{4} + \frac{{d_j^2
}}{4}}
\]
where $d_i$ and $d_j$ are the charge separation distance of dipole
and respectively. This approximation to the multipole expansion
maintains the fast fall-off of the multipole potentials but lacks
the normal divergences when two polar groups get close to one
another.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{split_dipole.eps}
\caption[Comparison between electrostatic approximation]{Electron
density profile along the bilayer normal.}
\label{lipidFigure:splitDipole}
\end{figure}

%\section{\label{lipidSection:methods}Methods}

\section{\label{lipidSection:resultDiscussion}Results and Discussion}

\subsection{One Lipid in Sea of Water Molecules}

To exclude the inter-headgroup interaction, atomistic models of one
lipid (DMPC or DLPE) in sea of water molecules (TIP3P) were built
and studied using atomistic molecular dynamics. The simulation was
analyzed using a set of radial distribution functions, which give
the probability of finding a pair of molecular species a distance
apart, relative to the probability expected for a completely random
distribution function at the same density.

\begin{equation}
g_{AB} (r) = \frac{1}{{\rho _B }}\frac{1}{{N_A }} < \sum\limits_{i
\in A} {\sum\limits_{j \in B} {\delta (r - r_{ij} )} }  >
\end{equation}
\begin{equation}
g_{AB} (r,\cos \theta ) = \frac{1}{{\rho _B }}\frac{1}{{N_A }} <
\sum\limits_{i \in A} {\sum\limits_{j \in B} {\delta (r - r_{ij} )}
} \delta (\cos \theta _{ij}  - \cos \theta ) >
\end{equation}

From figure 4(a), we can identify the first solvation shell (3.5
$\AA$) and the second solvation shell (5.0 $\AA$) from both plots.
However, the corresponding orientations are different. In DLPE,
water molecules prefer to sit around -NH3 group due to the hydrogen
bonding. In contrast, because of the hydrophobic effect of the
-N(CH3)3 group, the preferred position of water molecules in DMPC is
around the -PO4 Group. When the water molecules are far from the
headgroup, the distribution of the two angles should be uniform. The
correlation close to center of the headgroup dipole (< 5 $\AA$) in
both plots tell us that in the closely-bound region, the dipoles of
the water molecules are preferentially anti-aligned with the dipole
of headgroup.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{g_atom.eps}
\caption[The pair correlation functions for atomistic models]{}
\label{lipidFigure:PCFAtom}
\end{figure}

The initial configurations of coarse-grained systems are constructed
from the previous atomistic ones. The parameters for the
coarse-grained model in Table~\ref{lipidTable:parameter} are
estimated and tuned using isothermal-isobaric molecular dynamics.
Pair distribution functions calculated from coarse-grained models
preserve the basic characteristics of the atomistic simulations. The
water density, measured in a head-group-fixed reference frame,
surrounding two phospholipid headgroups is shown in
Fig.~\ref{lipidFigure:PCFCoarse}. It is clear that the phosphate end
in DMPC and the amine end in DMPE are the two most heavily solvated
atoms.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{g_coarse.eps}
\caption[The pair correlation functions for coarse-grained models]{}
\label{lipidFigure:PCFCoarse}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{EWD_coarse.eps}
\caption[Excess water density of coarse-grained phospholipids]{ }
\label{lipidFigure:EWDCoarse}
\end{figure}

\begin{table}
\caption{The Parameters For Coarse-grained Phospholipids}
\label{lipidTable:parameter}
\begin{center}
\begin{tabular}{|l|c|c|c|c|c|}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  Atom type & Mass & $\sigma$ & $\epsilon$ & charge & Dipole \\

  $\text{{\sc CH}}_2$ & 14.03  & 3.95 & 0.0914 & 0 & 0 \\
  $\text{{\sc CH}}_3$ & 15.04  & 3.75 & 0.195  & 0 & 0 \\
  $\text{{\sc CE}}$   & 28.01  & 3.427& 0.294  & 0 & 1.693
  $\text{{\sc CK}}$   & 28.01  & 3.592& 0.311  & 0 & 2.478 \\
  $\text{{\sc PO}}_4$ & 108.995& 3.9  & 1.88   & -1&  0 \\
  $\text{{\sc HDP}}$  & 14.03  & 4.0  & 0.13   & 0 &  0 \\
  $\text{{\sc NC}}_4$ & 73.137 & 4.9  & 0.88   & +1&  0 \\
  $\text{{\sc NH}}_3$ & 17.03  & 3.25 & 0.17   & +1&  0\\
  \hline
\end{tabular}
\end{center}
\end{table}

\subsection{Bilayer Simulations Using Coarse-grained Model}

A bilayer system consisting of 128 DMPC lipids and 3655 water
molecules has been constructed from an atomistic coordinate
file.[15] The MD simulation is performed at constant temperature, T
= 300K, and constant pressure, p = 1 atm, and consisted of an
equilibration period of 2 ns. During the equilibration period, the
system was initially simulated at constant volume for 1ns. Once the
system was equilibrated at constant volume, the cell dimensions of
the system was relaxed by performing under NPT conditions using
Nos��-Hoover extended system isothermal-isobaric dynamics. After
equilibration, different properties were evaluated over a production
run of 5 ns.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{bilayer.eps}
\caption[Image of a coarse-grained bilayer system]{A coarse-grained
bilayer system consisting of 128 DMPC lipids and 3655 SSD water
molecules.}
\label{lipidFigure:bilayer}
\end{figure}

\subsubsection{Electron Density Profile (EDP)}

Assuming a gaussian distribution of electrons on each atomic center
with a variance estimated from the size of the van der Waals radius,
the EDPs which are proportional to the density profiles measured
along the bilayer normal obtained by x-ray scattering experiment,
can be expressed by
\begin{equation}
\rho _{x - ray} (z)dz \propto \sum\limits_{i = 1}^N {\frac{{n_i
}}{V}\frac{1}{{\sqrt {2\pi \sigma ^2 } }}e^{ - (z - z_i )^2 /2\sigma
^2 } dz},
\end{equation}
where $\sigma$ is the variance equal to the van der Waals radius,
$n_i$ is the atomic number of site $i$ and $V$ is the volume of the
slab between $z$ and $z+dz$ . The highest density of total EDP
appears at the position of lipid-water interface corresponding to
headgroup, glycerol, and carbonyl groups of the lipids and the
distribution of water locked near the head groups, while the lowest
electron density is in the hydrocarbon region. As a good
approximation to the thickness of the bilayer, the headgroup spacing
, is defined as the distance between two peaks in the electron
density profile, calculated from our simulations to be 34.1 $\AA$.
This value is close to the x-ray diffraction experimental value 34.4
$\AA$

\begin{figure}
\centering
\includegraphics[width=\linewidth]{electron_density.eps}
\caption[The density profile of the lipid bilayers]{Electron density
profile along the bilayer normal. The water density is shown in red,
the density due to the headgroups in green, the glycerol backbone in
brown, $\text{{\sc CH}}_2$ in yellow, $\text{{\sc CH}}_3$ in cyan,
and total density due to DMPC in blue.}
\label{lipidFigure:electronDensity}
\end{figure}

\subsubsection{$\text{S}_{\text{{\sc cd}}}$ Order Parameter}

Measuring deuterium order parameters by NMR is a useful technique to
study the orientation of hydrocarbon chains in phospholipids. The
order parameter tensor $S$ is defined by:
\begin{equation}
S_{ij}  = \frac{1}{2} < 3\cos \theta _i \cos \theta _j  - \delta
_{ij}  >
\end{equation}
where $\theta$ is the angle between the  $i$th molecular axis and
the bilayer normal ($z$ axis). The brackets denote an average over
time and molecules. The molecular axes are defined:
\begin{itemize}
\item $\mathbf{\hat{z}}$ is the vector from $C_{n-1}$ to $C_{n+1}$.
\item $\mathbf{\hat{y}}$ is the vector that is perpendicular to $z$ and
in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$.
\item $\mathbf{\hat{x}}$ is the vector perpendicular to
$\mathbf{\hat{y}}$ and $\mathbf{\hat{z}}$.
\end{itemize}
In coarse-grained model, although there are no explicit hydrogens,
the order parameter can still be written in terms of carbon ordering
at each point of the chain
\begin{equation}
S_{ij}  = \frac{1}{2} < 3\cos \theta _i \cos \theta _j  - \delta
_{ij}  >.
\end{equation}

Fig.~\ref{lipidFigure:Scd} shows the order parameter profile
calculated for our coarse-grained DMPC bilayer system at 300K. Also
shown are the experimental data of Tiburu. The fact that
$\text{S}_{\text{{\sc cd}}}$ order parameters calculated from
simulation are higher than the experimental ones is ascribed to the
assumption of the locations of implicit hydrogen atoms which are
fixed in coarse-grained models at positions relative to the CC
vector.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{scd.eps}
\caption[$\text{S}_{\text{{\sc cd}}}$ order parameter]{A comparison
of $|\text{S}_{\text{{\sc cd}}}|$ between coarse-grained model
(blue) and DMPC\cite{petrache00} (black) near 300~K.}
\label{lipidFigure:Scd}
\end{figure}
Revision:	2776
Committed:	Thu May 25 21:32:14 2006 UTC (18 years, 1 month ago) by tim
Content type:	application/x-tex
File size:	17118 byte(s)
Log Message:	Adding Lipid 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			Figure 1 shows a schematic for our coarse-grained phospholipid
105			model. The lipid head group is modeled by a linear rigid body which
106			consists of three Lennard-Jones spheres and a centrally located
107			point-dipole. The backbone atoms in the glycerol motif are modeled
108			by Lennard-Jones spheres with dipoles. Alkyl groups in hydrocarbon
109			chains are replaced with unified CH2 or CH3 atoms.
110
111			Accurate and efficient computation of electrostatics is one of the
112			most difficult tasks in molecular modeling. Traditionally, the
113			electrostatic interaction between two molecular species is
114			calculated as a sum of interactions between pairs of point charges,
115			using Coulomb's law:
116			\[
117			V = \sum\limits_{i = 1}^{N_A } {\sum\limits_{j = 1}^{N_B }
118			{\frac{{q_i q_j }}{{4\pi \varepsilon _0 r_{ij} }}} }
119			\]
120			where $N_A$ and $N_B$ are the number of point charges in the two
121			molecular species. Originally developed to study ionic crystals, the
122			Ewald summation method mathematically transforms this
123			straightforward but conditionally convergent summation into two more
124			complicated but rapidly convergent sums. One summation is carried
125			out in reciprocal space while the other is carried out in real
126			space. An alternative approach is a multipole expansion, which is
127			based on electrostatic moments, such as charge (monopole), dipole,
128			quadruple etc.
129
130			Here we consider a linear molecule which consists of two point
131			charges $\pm q$ located at $ ( \pm \frac{d}{2},0)$. The
132			electrostatic potential at point $P$ is given by:
133			\[
134			\frac{1}{{4\pi \varepsilon _0 }}\left( {\frac{{ - q}}{{r_ - }} +
135			\frac{{ + q}}{{r_ + }}} \right) = \frac{1}{{4\pi \varepsilon _0
136			}}\left( {\frac{{ - q}}{{\sqrt {r^2 + \frac{{d^2 }}{4} + rd\cos
137			\theta } }} + \frac{q}{{\sqrt {r^2 + \frac{{d^2 }}{4} - rd\cos
138			\theta } }}} \right)
139			\]
140
141			The basic assumption of the multipole expansion is $r \gg d$ , thus,
142			$\frac{{d^2 }}{4}$ inside the square root of the denominator is
143			neglected. This is a reasonable approximation in most cases.
144			Unfortunately, in our headgroup model, the distance of charge
145			separation $d$ (4.63 $\AA$ in PC headgroup) may be comparable to
146			$r$. Nevertheless, we could still assume $ \cos \theta \approx 0$
147			in the central region of the headgroup. Using Taylor expansion and
148			associating appropriate terms with electric moments will result in a
149			"split-dipole" approximation:
150			\[
151			V(r) = \frac{1}{{4\pi \varepsilon _0 }}\frac{{r\mu \cos \theta
152			}}{{R^3 }}
153			\]
154			where$R = \sqrt {r^2 + \frac{{d^2 }}{4}}$ Placing a dipole at point
155			$P$ and applying the same strategy, the interaction between two
156			split-dipoles is then given by:
157			\[
158			V_{sd} (r_{ij} ,\Omega _i ,\Omega _j ) = \frac{1}{{4\pi \varepsilon
159			_0 }}\left[ {\frac{{\mu _i \cdot \mu _j }}{{R_{ij}^3 }} -
160			\frac{{3\left( {\mu _i \cdot r_{ij} } \right)\left( {\mu _i \cdot
161			r_{ij} } \right)}}{{R_{ij}^5 }}} \right]
162			\]
163			where $\mu _i$ and $\mu _j$ are the dipole moment of molecule $i$
164			and molecule $j$ respectively, $r_{ij}$ is vector between molecule
165			$i$ and molecule $j$, and $R_{ij{$ is given by,
166			\[
167			R_{ij} = \sqrt {r_{ij}^2 + \frac{{d_i^2 }}{4} + \frac{{d_j^2
168			}}{4}}
169			\]
170			where $d_i$ and $d_j$ are the charge separation distance of dipole
171			and respectively. This approximation to the multipole expansion
172			maintains the fast fall-off of the multipole potentials but lacks
173			the normal divergences when two polar groups get close to one
174			another.
175
176			\begin{figure}
177			\centering
178			\includegraphics[width=\linewidth]{split_dipole.eps}
179			\caption[Comparison between electrostatic approximation]{Electron
180			density profile along the bilayer normal.}
181			\label{lipidFigure:splitDipole}
182			\end{figure}
183
184			%\section{\label{lipidSection:methods}Methods}
185
186	tim	2685	\section{\label{lipidSection:resultDiscussion}Results and Discussion}
187	tim	2776
188			\subsection{One Lipid in Sea of Water Molecules}
189
190			To exclude the inter-headgroup interaction, atomistic models of one
191			lipid (DMPC or DLPE) in sea of water molecules (TIP3P) were built
192			and studied using atomistic molecular dynamics. The simulation was
193			analyzed using a set of radial distribution functions, which give
194			the probability of finding a pair of molecular species a distance
195			apart, relative to the probability expected for a completely random
196			distribution function at the same density.
197
198			\begin{equation}
199			g_{AB} (r) = \frac{1}{{\rho _B }}\frac{1}{{N_A }} < \sum\limits_{i
200			\in A} {\sum\limits_{j \in B} {\delta (r - r_{ij} )} } >
201			\end{equation}
202			\begin{equation}
203			g_{AB} (r,\cos \theta ) = \frac{1}{{\rho _B }}\frac{1}{{N_A }} <
204			\sum\limits_{i \in A} {\sum\limits_{j \in B} {\delta (r - r_{ij} )}
205			} \delta (\cos \theta _{ij} - \cos \theta ) >
206			\end{equation}
207
208			From figure 4(a), we can identify the first solvation shell (3.5
209			$\AA$) and the second solvation shell (5.0 $\AA$) from both plots.
210			However, the corresponding orientations are different. In DLPE,
211			water molecules prefer to sit around -NH3 group due to the hydrogen
212			bonding. In contrast, because of the hydrophobic effect of the
213			-N(CH3)3 group, the preferred position of water molecules in DMPC is
214			around the -PO4 Group. When the water molecules are far from the
215			headgroup, the distribution of the two angles should be uniform. The
216			correlation close to center of the headgroup dipole (< 5 $\AA$) in
217			both plots tell us that in the closely-bound region, the dipoles of
218			the water molecules are preferentially anti-aligned with the dipole
219			of headgroup.
220
221			\begin{figure}
222			\centering
223			\includegraphics[width=\linewidth]{g_atom.eps}
224			\caption[The pair correlation functions for atomistic models]{}
225			\label{lipidFigure:PCFAtom}
226			\end{figure}
227
228			The initial configurations of coarse-grained systems are constructed
229			from the previous atomistic ones. The parameters for the
230			coarse-grained model in Table~\ref{lipidTable:parameter} are
231			estimated and tuned using isothermal-isobaric molecular dynamics.
232			Pair distribution functions calculated from coarse-grained models
233			preserve the basic characteristics of the atomistic simulations. The
234			water density, measured in a head-group-fixed reference frame,
235			surrounding two phospholipid headgroups is shown in
236			Fig.~\ref{lipidFigure:PCFCoarse}. It is clear that the phosphate end
237			in DMPC and the amine end in DMPE are the two most heavily solvated
238			atoms.
239
240			\begin{figure}
241			\centering
242			\includegraphics[width=\linewidth]{g_coarse.eps}
243			\caption[The pair correlation functions for coarse-grained models]{}
244			\label{lipidFigure:PCFCoarse}
245			\end{figure}
246
247			\begin{figure}
248			\centering
249			\includegraphics[width=\linewidth]{EWD_coarse.eps}
250			\caption[Excess water density of coarse-grained phospholipids]{ }
251			\label{lipidFigure:EWDCoarse}
252			\end{figure}
253
254			\begin{table}
255			\caption{The Parameters For Coarse-grained Phospholipids}
256			\label{lipidTable:parameter}
257			\begin{center}
258			\begin{tabular}{\|l\|c\|c\|c\|c\|c\|}
259			\hline
260			% after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
261			Atom type & Mass & $\sigma$ & $\epsilon$ & charge & Dipole \\
262
263			$\text{{\sc CH}}_2$ & 14.03 & 3.95 & 0.0914 & 0 & 0 \\
264			$\text{{\sc CH}}_3$ & 15.04 & 3.75 & 0.195 & 0 & 0 \\
265			$\text{{\sc CE}}$ & 28.01 & 3.427& 0.294 & 0 & 1.693
266			$\text{{\sc CK}}$ & 28.01 & 3.592& 0.311 & 0 & 2.478 \\
267			$\text{{\sc PO}}_4$ & 108.995& 3.9 & 1.88 & -1& 0 \\
268			$\text{{\sc HDP}}$ & 14.03 & 4.0 & 0.13 & 0 & 0 \\
269			$\text{{\sc NC}}_4$ & 73.137 & 4.9 & 0.88 & +1& 0 \\
270			$\text{{\sc NH}}_3$ & 17.03 & 3.25 & 0.17 & +1& 0\\
271			\hline
272			\end{tabular}
273			\end{center}
274			\end{table}
275
276			\subsection{Bilayer Simulations Using Coarse-grained Model}
277
278			A bilayer system consisting of 128 DMPC lipids and 3655 water
279			molecules has been constructed from an atomistic coordinate
280			file.[15] The MD simulation is performed at constant temperature, T
281			= 300K, and constant pressure, p = 1 atm, and consisted of an
282			equilibration period of 2 ns. During the equilibration period, the
283			system was initially simulated at constant volume for 1ns. Once the
284			system was equilibrated at constant volume, the cell dimensions of
285			the system was relaxed by performing under NPT conditions using
286			Nos��-Hoover extended system isothermal-isobaric dynamics. After
287			equilibration, different properties were evaluated over a production
288			run of 5 ns.
289
290			\begin{figure}
291			\centering
292			\includegraphics[width=\linewidth]{bilayer.eps}
293			\caption[Image of a coarse-grained bilayer system]{A coarse-grained
294			bilayer system consisting of 128 DMPC lipids and 3655 SSD water
295			molecules.}
296			\label{lipidFigure:bilayer}
297			\end{figure}
298
299			\subsubsection{Electron Density Profile (EDP)}
300
301			Assuming a gaussian distribution of electrons on each atomic center
302			with a variance estimated from the size of the van der Waals radius,
303			the EDPs which are proportional to the density profiles measured
304			along the bilayer normal obtained by x-ray scattering experiment,
305			can be expressed by
306			\begin{equation}
307			\rho _{x - ray} (z)dz \propto \sum\limits_{i = 1}^N {\frac{{n_i
308			}}{V}\frac{1}{{\sqrt {2\pi \sigma ^2 } }}e^{ - (z - z_i )^2 /2\sigma
309			^2 } dz},
310			\end{equation}
311			where $\sigma$ is the variance equal to the van der Waals radius,
312			$n_i$ is the atomic number of site $i$ and $V$ is the volume of the
313			slab between $z$ and $z+dz$ . The highest density of total EDP
314			appears at the position of lipid-water interface corresponding to
315			headgroup, glycerol, and carbonyl groups of the lipids and the
316			distribution of water locked near the head groups, while the lowest
317			electron density is in the hydrocarbon region. As a good
318			approximation to the thickness of the bilayer, the headgroup spacing
319			, is defined as the distance between two peaks in the electron
320			density profile, calculated from our simulations to be 34.1 $\AA$.
321			This value is close to the x-ray diffraction experimental value 34.4
322			$\AA$
323
324			\begin{figure}
325			\centering
326			\includegraphics[width=\linewidth]{electron_density.eps}
327			\caption[The density profile of the lipid bilayers]{Electron density
328			profile along the bilayer normal. The water density is shown in red,
329			the density due to the headgroups in green, the glycerol backbone in
330			brown, $\text{{\sc CH}}_2$ in yellow, $\text{{\sc CH}}_3$ in cyan,
331			and total density due to DMPC in blue.}
332			\label{lipidFigure:electronDensity}
333			\end{figure}
334
335			\subsubsection{$\text{S}_{\text{{\sc cd}}}$ Order Parameter}
336
337			Measuring deuterium order parameters by NMR is a useful technique to
338			study the orientation of hydrocarbon chains in phospholipids. The
339			order parameter tensor $S$ is defined by:
340			\begin{equation}
341			S_{ij} = \frac{1}{2} < 3\cos \theta _i \cos \theta _j - \delta
342			_{ij} >
343			\end{equation}
344			where $\theta$ is the angle between the $i$th molecular axis and
345			the bilayer normal ($z$ axis). The brackets denote an average over
346			time and molecules. The molecular axes are defined:
347			\begin{itemize}
348			\item $\mathbf{\hat{z}}$ is the vector from $C_{n-1}$ to $C_{n+1}$.
349			\item $\mathbf{\hat{y}}$ is the vector that is perpendicular to $z$ and
350			in the plane through $C_{n-1}$, $C_{n}$, and $C_{n+1}$.
351			\item $\mathbf{\hat{x}}$ is the vector perpendicular to
352			$\mathbf{\hat{y}}$ and $\mathbf{\hat{z}}$.
353			\end{itemize}
354			In coarse-grained model, although there are no explicit hydrogens,
355			the order parameter can still be written in terms of carbon ordering
356			at each point of the chain
357			\begin{equation}
358			S_{ij} = \frac{1}{2} < 3\cos \theta _i \cos \theta _j - \delta
359			_{ij} >.
360			\end{equation}
361
362			Fig.~\ref{lipidFigure:Scd} shows the order parameter profile
363			calculated for our coarse-grained DMPC bilayer system at 300K. Also
364			shown are the experimental data of Tiburu. The fact that
365			$\text{S}_{\text{{\sc cd}}}$ order parameters calculated from
366			simulation are higher than the experimental ones is ascribed to the
367			assumption of the locations of implicit hydrogen atoms which are
368			fixed in coarse-grained models at positions relative to the CC
369			vector.
370
371			\begin{figure}
372			\centering
373			\includegraphics[width=\linewidth]{scd.eps}
374			\caption[$\text{S}_{\text{{\sc cd}}}$ order parameter]{A comparison
375			of $\|\text{S}_{\text{{\sc cd}}}\|$ between coarse-grained model
376			(blue) and DMPC\cite{petrache00} (black) near 300~K.}
377			\label{lipidFigure:Scd}
378			\end{figure}