trunk/mdRipple/mdripple.tex

%\documentclass[aps,pre,twocolumn,amssymb,showpacs,floatfix]{revtex4}
\documentclass[aps,pre,preprint,amssymb,showpacs]{revtex4}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}

\begin{document}
\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}

%\bibliographystyle{aps}

\title{Dipolar Ordering in Molecular-scale Models of the Ripple Phase
in Lipid Membranes} 
\author{Xiuquan Sun and J. Daniel Gezelter}
\email[E-mail:]{gezelter@nd.edu}
\affiliation{Department of Chemistry and Biochemistry,\\
University of Notre Dame, \\ 
Notre Dame, Indiana 46556}

\date{\today} 

\begin{abstract}
The ripple phase in phosphatidylcholine (PC) bilayers has never been
completely explained. 
\end{abstract}

\pacs{}
\maketitle

\section{Introduction}
\label{sec:Int}
Fully hydrated lipids will aggregate spontaneously to form bilayers
which exhibit a variety of phases depending on their temperatures and
compositions. Among these phases, a periodic rippled phase
($P_{\beta'}$) appears as an intermediate phase between the gel
($L_\beta$) and fluid ($L_{\alpha}$) phases for relatively pure
phosphatidylcholine (PC) bilayers.  The ripple phase has attracted
substantial experimental interest over the past 30 years. Most
structural information of the ripple phase has been obtained by the
X-ray diffraction~\cite{Sun96,Katsaras00} and freeze-fracture electron
microscopy (FFEM).~\cite{Copeland80,Meyer96} Recently, Kaasgaard {\it
et al.} used atomic force microscopy (AFM) to observe ripple phase
morphology in bilayers supported on mica.~\cite{Kaasgaard03} The
experimental results provide strong support for a 2-dimensional
hexagonal packing lattice of the lipid molecules within the ripple
phase.  This is a notable change from the observed lipid packing
within the gel phase.~\cite{Cevc87}

A number of theoretical models have been presented to explain the
formation of the ripple phase. Marder {\it et al.} used a
curvature-dependent Landau-de Gennes free-energy functional to predict
a rippled phase.~\cite{Marder84} This model and other related continuum
models predict higher fluidity in convex regions and that concave
portions of the membrane correspond to more solid-like regions.
Carlson and Sethna used a packing-competition model (in which head
groups and chains have competing packing energetics) to predict the
formation of a ripple-like phase.  Their model predicted that the
high-curvature portions have lower-chain packing and correspond to
more fluid-like regions.  Goldstein and Leibler used a mean-field
approach with a planar model for {\em inter-lamellar} interactions to
predict rippling in multilamellar phases.~\cite{Goldstein88} McCullough
and Scott proposed that the {\em anisotropy of the nearest-neighbor 
interactions} coupled to hydrophobic constraining forces which
restrict height differences between nearest neighbors is the origin of
the ripple phase.~\cite{McCullough90} Lubensky and MacKintosh
introduced a Landau theory for tilt order and curvature of a single
membrane and concluded that {\em coupling of molecular tilt to membrane
curvature} is responsible for the production of
ripples.~\cite{Lubensky93} Misbah, Duplat and Houchmandzadeh proposed
that {\em inter-layer dipolar interactions} can lead to ripple
instabilities.~\cite{Misbah98} Heimburg presented a {\em coexistence
model} for ripple formation in which he postulates that fluid-phase
line defects cause sharp curvature between relatively flat gel-phase
regions.~\cite{Heimburg00} Kubica has suggested that a lattice model of
polar head groups could be valuable in trying to understand bilayer
phase formation.~\cite{Kubica02} Bannerjee used Monte Carlo simulations
of lamellar stacks of hexagonal lattices to show that large headgroups
and molecular tilt with respect to the membrane normal vector can
cause bulk rippling.~\cite{Bannerjee02}

In contrast, few large-scale molecular modelling studies have been
done due to the large size of the resulting structures and the time
required for the phases of interest to develop.  With all-atom (and
even unified-atom) simulations, only one period of the ripple can be
observed and only for timescales in the range of 10-100 ns.  One of
the most interesting molecular simulations was carried out by De Vries
{\it et al.}~\cite{deVries05}. According to their simulation results,
the ripple consists of two domains, one resembling the gel bilayer,
while in the other, the two leaves of the bilayer are fully
interdigitated.  The mechanism for the formation of the ripple phase
suggested by their work is a packing competition between the head
groups and the tails of the lipid molecules.\cite{Carlson87} Recently,
the ripple phase has also been studied by Lenz and Schmid using Monte
Carlo simulations.\cite{Lenz07} Their structures are similar to the De
Vries {\it et al.} structures except that the connection between the
two leaves of the bilayer is a narrow interdigitated line instead of
the fully interdigitated domain.  The symmetric ripple phase was also
observed by Lenz {\it et al.}, and their work supports other claims
that the mismatch between the size of the head group and tail of the
lipid molecules is the driving force for the formation of the ripple
phase. Ayton and Voth have found significant undulations in
zero-surface-tension states of membranes simulated via dissipative
particle dynamics, but their results are consistent with purely
thermal undulations.~\cite{Ayton02}

Although the organization of the tails of lipid molecules are
addressed by these molecular simulations and the packing competition
between headgroups and tails is strongly implicated as the primary
driving force for ripple formation, questions about the ordering of
the head groups in ripple phase has not been settled.

In a recent paper, we presented a simple ``web of dipoles'' spin
lattice model which provides some physical insight into relationship
between dipolar ordering and membrane buckling.\cite{Sun2007} We found
that dipolar elastic membranes can spontaneously buckle, forming
ripple-like topologies.  The driving force for the buckling in dipolar
elastic membranes the antiferroelectric ordering of the dipoles, and
this was evident in the ordering of the dipole director axis
perpendicular to the wave vector of the surface ripples.  A similiar
phenomenon has also been observed by Tsonchev {\it et al.} in their
work on the spontaneous formation of dipolar peptide chains into
curved nano-structures.\cite{Tsonchev04,Tsonchev04II}

In this paper, we construct a somewhat more realistic molecular-scale
lipid model than our previous ``web of dipoles'' and use molecular
dynamics simulations to elucidate the role of the head group dipoles
in the formation and morphology of the ripple phase.  We describe our
model and computational methodology in section \ref{sec:method}.
Details on the simulations are presented in section
\ref{sec:experiment}, with results following in section
\ref{sec:results}.  A final discussion of the role of dipolar heads in
the ripple formation can be found in section
\ref{sec:discussion}.

\section{Computational Model}
\label{sec:method}

\begin{figure}[htb]
\centering
\includegraphics[width=4in]{lipidModels} 
\caption{Three different representations of DPPC lipid molecules,
including the chemical structure, an atomistic model, and the
head-body ellipsoidal coarse-grained model used in this
work.\label{fig:lipidModels}}
\end{figure}

Our simple molecular-scale lipid model for studying the ripple phase
is based on two facts: one is that the most essential feature of lipid
molecules is their amphiphilic structure with polar head groups and
non-polar tails. Another fact is that the majority of lipid molecules
in the ripple phase are relatively rigid (i.e. gel-like) which makes
some fraction of the details of the chain dynamics negligible.  Figure
\ref{fig:lipidModels} shows the molecular strucure of a DPPC
molecule, as well as atomistic and molecular-scale representations of
a DPPC molecule.  The hydrophilic character of the head group is
largely due to the separation of charge between the nitrogen and
phosphate groups.  The zwitterionic nature of the PC headgroups leads
to abnormally large dipole moments (as high as 20.6 D), and this
strongly polar head group interacts strongly with the solvating water
layers immediately surrounding the membrane.  The hydrophobic tail
consists of fatty acid chains.  In our molecular scale model, lipid
molecules have been reduced to these essential features; the fatty
acid chains are represented by an ellipsoid with a dipolar ball
perched on one end to represent the effects of the charge-separated
head group.  In real PC lipids, the direction of the dipole is
nearly perpendicular to the tail, so we have fixed the direction of
the point dipole rigidly in this orientation.  

The ellipsoidal portions of the model interact via the Gay-Berne
potential which has seen widespread use in the liquid crystal
community.  Ayton and Voth have also used Gay-Berne ellipsoids for
modelling large length-scale properties of lipid
bilayers.\cite{Ayton01} In its original form, the Gay-Berne potential
was a single site model for the interactions of rigid ellipsoidal
molecules.\cite{Gay81} It can be thought of as a modification of the
Gaussian overlap model originally described by Berne and
Pechukas.\cite{Berne72} The potential is constructed in the familiar
form of the Lennard-Jones function using orientation-dependent
$\sigma$ and $\epsilon$ parameters,
\begin{eqnarray*}
V_{ij}({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat
r}_{ij}}) = 4\epsilon ({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
{\mathbf{\hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i},
{\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^{12}
-\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
{\mathbf{\hat r}_{ij}})+\sigma_0}\right)^6\right]
\label{eq:gb}
\end{eqnarray*}

The range $(\sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
\hat{r}}_{ij}))$, and strength $(\epsilon({\bf \hat{u}}_{i},{\bf
\hat{u}}_{j},{\bf \hat{r}}_{ij}))$ parameters 
are dependent on the relative orientations of the two molecules (${\bf
\hat{u}}_{i},{\bf \hat{u}}_{j}$) as well as the direction of the
intermolecular separation (${\bf \hat{r}}_{ij}$).  $\sigma$ and
$\sigma_0$ are also governed by shape mixing and anisotropy variables,
\begin {equation}
\begin{array}{rcl}
\sigma_0 & = & \sqrt{d_i^2 + d_j^2} \\ \\
\chi & = & \left[ \frac{\left( l_i^2 - d_i^2 \right) \left(l_j^2 -
d_j^2 \right)}{\left( l_j^2 + d_i^2 \right) \left(l_i^2 +
d_j^2 \right)}\right]^{1/2} \\ \\
\alpha^2 & = & \left[ \frac{\left( l_i^2 - d_i^2 \right) \left(l_j^2 +
d_i^2 \right)}{\left( l_j^2 - d_j^2 \right) \left(l_i^2 +
d_j^2 \right)}\right]^{1/2},
\end{array}
\end{equation}
where $l$ and $d$ describe the length and width of each uniaxial
ellipsoid.  These shape anisotropy parameters can then be used to
calculate the range function,
\begin {equation}
\sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf \hat{r}}_{ij}) = \sigma_{0}
 \left[ 1- \left\{ \frac{ \chi \alpha^2 ({\bf \hat{u}}_i \cdot {\bf
\hat{r}}_{ij} ) + \chi \alpha^{-2} ({\bf \hat{u}}_j \cdot {\bf
\hat{r}}_{ij} ) - 2 \chi^2 ({\bf \hat{u}}_i \cdot {\bf
\hat{r}}_{ij} )({\bf \hat{u}}_j \cdot {\bf
\hat{r}}_{ij} ) ({\bf \hat{u}}_i \cdot {\bf \hat{u}}_j)}{1 - \chi^2
\left({\bf \hat{u}}_i \cdot {\bf \hat{u}}_j\right)^2} \right\}
\right]^{-1/2}
\end{equation}

Gay-Berne ellipsoids also have an energy scaling parameter,
$\epsilon^s$, which describes the well depth for two identical
ellipsoids in a {\it side-by-side} configuration.  Additionaly, a well
depth aspect ratio, $\epsilon^r = \epsilon^e / \epsilon^s$, describes
the ratio between the well depths in the {\it end-to-end} and
side-by-side configurations.  As in the range parameter, a set of
mixing and anisotropy variables can be used to describe the well
depths for dissimilar particles,
\begin {eqnarray*}
\epsilon_0 & = & \sqrt{\epsilon^s_i  * \epsilon^s_j} \\ \\
\epsilon^r & = & \sqrt{\epsilon^r_i  * \epsilon^r_j} \\ \\
\chi' & = & \frac{1 - (\epsilon^r)^{1/\mu}}{1 + (\epsilon^r)^{1/\mu}}
\\ \\
\alpha'^2 & = & \left[1 + (\epsilon^r)^{1/\mu}\right]^{-1}
\end{eqnarray*}
The form of the strength function is somewhat complicated,
\begin {eqnarray*}
\epsilon({\bf \hat{u}}_{i}, {\bf \hat{u}}_{j},{\bf \hat{r}}_{ij}) & = &
\epsilon_{0}  \epsilon_{1}^{\nu}({\bf \hat{u}}_{i}.{\bf \hat{u}}_{j})
 \epsilon_{2}^{\mu}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
\hat{r}}_{ij}) \\ \\
\epsilon_{1}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j}) & = &
\left[1-\chi^{2}({\bf \hat{u}}_{i}.{\bf
\hat{u}}_{j})^{2}\right]^{-1/2} \\ \\
\epsilon_{2}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf \hat{r}}_{ij}) &
= &
 1 - \left\{ \frac{ \chi' \alpha'^2 ({\bf \hat{u}}_i \cdot {\bf
\hat{r}}_{ij} ) + \chi' \alpha'^{-2} ({\bf \hat{u}}_j \cdot {\bf
\hat{r}}_{ij} ) - 2 \chi'^2 ({\bf \hat{u}}_i \cdot {\bf
\hat{r}}_{ij} )({\bf \hat{u}}_j \cdot {\bf
\hat{r}}_{ij} ) ({\bf \hat{u}}_i \cdot {\bf \hat{u}}_j)}{1 - \chi'^2
\left({\bf \hat{u}}_i \cdot {\bf \hat{u}}_j\right)^2} \right\},
\end {eqnarray*}
although many of the quantities and derivatives are identical with
those obtained for the range parameter. Ref. \onlinecite{Luckhurst90}
has a particularly good explanation of the choice of the Gay-Berne
parameters $\mu$ and $\nu$ for modeling liquid crystal molecules. An
excellent overview of the computational methods that can be used to
efficiently compute forces and torques for this potential can be found
in Ref. \onlinecite{Golubkov06}

The choices of parameters we have used in this study correspond to a
shape anisotropy of 3 for the chain portion of the molecule.  In
principle, this could be varied to allow for modeling of longer or
shorter chain lipid molecules. For these prolate ellipsoids, we have:
\begin{equation}
\begin{array}{rcl}
d & < & l \\
\epsilon^{r} & < & 1
\end{array}
\end{equation}
A sketch of the various structural elements of our molecular-scale
lipid / solvent model is shown in figure \ref{fig:lipidModel}.  The
actual parameters used in our simulations are given in table
\ref{tab:parameters}. 

\begin{figure}[htb]
\centering
\includegraphics[width=4in]{2lipidModel} 
\caption{The parameters defining the behavior of the lipid
models. $l / d$ is the ratio of the head group to body diameter.
Molecular bodies had a fixed aspect ratio of 3.0.  The solvent model
was a simplified 4-water bead ($\sigma_w \approx d$) that has been
used in other coarse-grained (DPD) simulations.  The dipolar strength
(and the temperature and pressure) were the only other parameters that
were varied systematically.\label{fig:lipidModel}}
\end{figure}

To take into account the permanent dipolar interactions of the
zwitterionic head groups, we place fixed dipole moments $\mu_{i}$ at
one end of the Gay-Berne particles.  The dipoles are oriented at an
angle $\theta = \pi / 2$ relative to the major axis.  These dipoles
are protected by a head ``bead'' with a range parameter which we have
varied between $1.20 d$ and $1.41 d$.  The head groups interact with
each other using a combination of Lennard-Jones,
\begin{eqnarray*}
V_{ij}(r_{ij}) = 4\epsilon_h \left[\left(\frac{\sigma_h}{r_{ij}}\right)^{12} -
\left(\frac{\sigma_h}{r_{ij}}\right)^6\right],
\end{eqnarray*}
and dipole-dipole,
\begin{eqnarray*}
V_{ij}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
\hat{r}}_{ij})) = \frac{|\mu|^2}{4 \pi \epsilon_0 r_{ij}^3}
\left[ \hat{\bf u}_i \cdot \hat{\bf u}_j - 3 (\hat{\bf u}_i \cdot
\hat{\bf r}_{ij})(\hat{\bf u}_j \cdot \hat{\bf r}_{ij}) \right]
\end{eqnarray*}
potentials.  
In these potentials, $\mathbf{\hat u}_i$ is the unit vector pointing
along dipole $i$ and $\mathbf{\hat r}_{ij}$ is the unit vector
pointing along the inter-dipole vector $\mathbf{r}_{ij}$.  

For the interaction between nonequivalent uniaxial ellipsoids (in this
case, between spheres and ellipsoids), the spheres are treated as
ellipsoids with an aspect ratio of 1 ($d = l$) and with an well depth
ratio ($\epsilon^r$) of 1 ($\epsilon^e = \epsilon^s$).  The form of
the Gay-Berne potential we are using was generalized by Cleaver {\it
et al.} and is appropriate for dissimilar uniaxial
ellipsoids.\cite{Cleaver96}

The solvent model in our simulations is identical to one used by
Marrink {\it et al.}  in their dissipative particle dynamics (DPD)
simulation of lipid bilayers.\cite{Marrink04} This solvent bead is a single
site that represents four water molecules (m = 72 amu) and has
comparable density and diffusive behavior to liquid water.  However,
since there are no electrostatic sites on these beads, this solvent
model cannot replicate the dielectric properties of water.
\begin{table*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{Potential parameters used for molecular-scale coarse-grained
lipid simulations}
\begin{tabular}{llccc}
\hline 
  & &  Head & Chain & Solvent \\
\hline 
$d$ (\AA) & & varied & 4.6  & 4.7 \\
$l$ (\AA) & & $= d$ & 13.8 & 4.7 \\
$\epsilon^s$ (kcal/mol) & & 0.185 & 1.0 & 0.8 \\
$\epsilon_r$ (well-depth aspect ratio)& & 1 & 0.2 &  1 \\
$m$ (amu) & & 196 & 760 & 72.06 \\
$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$) & & & & \\
\multicolumn{2}c{$I_{xx}$} & 1125 & 45000 & N/A \\
\multicolumn{2}c{$I_{yy}$} & 1125 & 45000 & N/A \\
\multicolumn{2}c{$I_{zz}$} &  0 &    9000 & N/A \\
$\mu$ (Debye) & & varied & 0 & 0 \\
\end{tabular}
\label{tab:parameters}
\end{center}
\end{minipage}
\end{table*}

A switching function has been applied to all potentials to smoothly
turn off the interactions between a range of $22$ and $25$ \AA.

The parameters that were systematically varied in this study were the
size of the head group ($\sigma_h$), the strength of the dipole moment
($\mu$), and the temperature of the system.  Values for $\sigma_h$
ranged from 5.5 \AA\  to 6.5 \AA\ .  If the width of the tails is
taken to be the unit of length, these head groups correspond to a
range from $1.2 d$ to $1.41 d$.  Since the solvent beads are nearly
identical in diameter to the tail ellipsoids, all distances that
follow will be measured relative to this unit of distance. 

\section{Experimental Methodology}
\label{sec:experiment}

To create unbiased bilayers, all simulations were started from two
perfectly flat monolayers separated by a 26 \AA\ gap between the
molecular bodies of the upper and lower leaves.  The separated
monolayers were evolved in a vaccum with $x-y$ anisotropic pressure
coupling. The length of $z$ axis of the simulations was fixed and a
constant surface tension was applied to enable real fluctuations of
the bilayer. Periodic boundary conditions were used, and $480-720$
lipid molecules were present in the simulations, depending on the size
of the head beads.  In all cases, the two monolayers spontaneously
collapsed into bilayer structures within 100 ps. Following this
collapse, all systems were equlibrated for $100$ ns at $300$ K.

The resulting bilayer structures were then solvated at a ratio of $6$
solvent beads (24 water molecules) per lipid. These configurations
were then equilibrated for another $30$ ns. All simulations utilizing
the solvent were carried out at constant pressure ($P=1$ atm) with
$3$D anisotropic coupling, and constant surface tension
($\gamma=0.015$ UNIT). Given the absence of fast degrees of freedom in
this model, a timestep of $50$ fs was utilized with excellent energy
conservation.  Data collection for structural properties of the
bilayers was carried out during a final 5 ns run following the solvent
equilibration.  All simulations were performed using the OOPSE
molecular modeling program.\cite{Meineke05}

\section{Results}
\label{sec:results}

Snapshots in Figure \ref{fig:phaseCartoon} show that the membrane is
more corrugated with increasing size of the head groups. The surface
is nearly flat when $\sigma_h=1.20 d$. With $\sigma_h=1.28 d$,
although the surface is still flat, the bilayer starts to splay
inward; the upper leaf of the bilayer is connected to the lower leaf
with an interdigitated line defect. Two periodicities with $100$ \AA\
wavelengths were observed in the simulation. This structure is very
similiar to the structure observed by de Vries and Lenz {\it et
al.}. The same basic structure is also observed when $\sigma_h=1.41
d$, but the wavelength of the surface corrugations depends sensitively
on the size of the ``head'' beads. From the undulation spectrum, the
corrugation is clearly non-thermal.
\begin{figure}[htb]
\centering
\includegraphics[width=4in]{phaseCartoon} 
\caption{A sketch to discribe the structure of the phases observed in
our simulations.\label{fig:phaseCartoon}}
\end{figure}

When $\sigma_h=1.35 d$, we observed another corrugated surface
morphology.  This structure is different from the asymmetric rippled
surface; there is no interdigitation between the upper and lower
leaves of the bilayer. Each leaf of the bilayer is broken into several
hemicylinderical sections, and opposite leaves are fitted together
much like roof tiles. Unlike the surface in which the upper
hemicylinder is always interdigitated on the leading or trailing edge
of lower hemicylinder, this ``symmetric'' ripple has no prefered
direction.  The corresponding structures are shown in Figure
\ref{fig:phaseCartoon} for elucidation of the detailed structures of
different phases.  The top panel in figure \ref{fig:phaseCartoon} is
the flat phase, the middle panel shows the asymmetric ripple phase
corresponding to $\sigma_h = 1.41 d$ and the lower panel shows the
symmetric ripple phase observed when $\sigma_h=1.35 d$. In the
symmetric ripple, the bilayer is continuous over the whole membrane,
however, in asymmetric ripple phase, the bilayer domains are connected
by thin interdigitated monolayers that share molecules between the
upper and lower leaves.
\begin{table*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{Phases, ripple wavelengths and amplitudes observed as a
function of the ratio between the head beads and the diameters of the
tails.  All lengths are normalized to the diameter of the tail
ellipsoids.}
\begin{tabular}{lccc}
\hline 
$\sigma_h / d$ & type of phase & $\lambda / d$ & $A / d$\\
\hline 
1.20 & flat & N/A & N/A \\
1.28 & asymmetric ripple or flat & 21.7 & N/A \\
1.35 & symmetric ripple & 17.2 & 2.2 \\
1.41 & asymmetric ripple & 15.4 & 1.5 \\
\end{tabular}
\label{tab:property}
\end{center}
\end{minipage}
\end{table*}

The membrane structures and the reduced wavelength $\lambda / d$,
reduced amplitude $A / d$ of the ripples are summarized in Table
\ref{tab:property}. The wavelength range is $15~21$ molecular bodies
and the amplitude is $1.5$ molecular bodies for asymmetric ripple and
$2.2$ for symmetric ripple. These values are consistent to the
experimental results.  Note, that given the lack of structural freedom
in the tails of our model lipids, the amplitudes observed from these
simulations are likely to underestimate of the true amplitudes.

\begin{figure}[htb]
\centering
\includegraphics[width=4in]{topDown} 
\caption{Top views of the flat (upper), asymmetric ripple (middle),
and symmetric ripple (lower) phases.  Note that the head-group dipoles
have formed head-to-tail chains in all three of these phases, but in
the two rippled phases, the dipolar chains are all aligned
{\it perpendicular} to the direction of the ripple.  The flat membrane
has multiple point defects in the dipolar orientational ordering, and
the dipolar ordering on the lower leaf of the bilayer can be in a
different direction from the upper leaf.\label{fig:topView}} 
\end{figure}

Figure \ref{fig:topView} shows snapshots of bird's-eye views of the
flat ($\sigma_h = 1.20 d$) and rippled ($\sigma_h = 1.41, 1.35 d$)
bilayers.  The directions of the dipoles on the head groups are
represented with two colored half spheres: blue (phosphate) and yellow
(amino).  For flat bilayers, the system exhibits signs of
orientational frustration; some disorder in the dipolar chains is
evident with kinks visible at the edges between different ordered
domains.  The lipids can also move independently of lipids in the
opposing leaf, so the ordering of the dipoles on one leaf is not
necessarily consistant with the ordering on the other.  These two
factors keep the total dipolar order parameter relatively low for the
flat phases.

With increasing head group size, the surface becomes corrugated, and
the dipoles cannot move as freely on the surface. Therefore, the
translational freedom of lipids in one layer is dependent upon the
position of lipids in the other layer.  As a result, the ordering of
the dipoles on head groups in one leaf is correlated with the ordering
in the other leaf.  Furthermore, as the membrane deforms due to the
corrugation, the symmetry of the allowed dipolar ordering on each leaf
is broken. The dipoles then self-assemble in a head-to-tail
configuration, and the dipolar order parameter increases dramatically.
However, the total polarization of the system is still close to zero.
This is strong evidence that the corrugated structure is an
antiferroelectric state.   It is also notable that the head-to-tail
arrangement of the dipoles is in a direction perpendicular to the wave
vector for the surface corrugation.  This is a similar finding to what
we observed in our earlier work on the elastic dipolar
membranes.\cite{Sun07} 

The $P_2$ order parameters (for both the molecular bodies and the head
group dipoles) have been calculated to quantify the ordering in these
phases.  $P_2 = 1$ implies a perfectly ordered structure, and $P_2 = 0$
implies complete orientational randomization. Figure \ref{fig:rP2}
shows the $P_2$ order parameter for the head-group dipoles increasing
with increasing head group size. When the heads of the lipid molecules
are small, the membrane is nearly flat. The dipolar ordering exhibits
frustrated orientational ordering in this circumstance. 

The ordering of the tails is essentially opposite to the ordering of
the dipoles on head group. The $P_2$ order parameter {\it decreases}
with increasing head size. This indicates that the surface is more
curved with larger head / tail size ratios. When the surface is flat,
all tails are pointing in the same direction (parallel to the normal
of the surface).  This simplified model appears to be exhibiting a
smectic A fluid phase, similar to the real $L_{\beta}$ phase.  We have
not observed a smectic C gel phase ($L_{\beta'}$) for this model
system.  Increasing the size of the heads, results in rapidly
decreasing $P_2$ ordering for the molecular bodies.

\begin{figure}[htb]
\centering
\includegraphics[width=\linewidth]{rP2} 
\caption{The $P_2$ order parameter as a function of the ratio of
$\sigma_h$ to $d$. \label{fig:rP2}}
\end{figure}

We studied the effects of the interactions between head groups on the
structure of lipid bilayer by changing the strength of the dipole.
Figure \ref{fig:sP2} shows how the $P_2$ order parameter changes with
increasing strength of the dipole. Generally the dipoles on the head
group are more ordered by increase in the strength of the interaction
between heads and are more disordered by decreasing the interaction
stength. When the interaction between the heads is weak enough, the
bilayer structure does not persist; all lipid molecules are solvated
directly in the water. The critial value of the strength of the dipole
depends on the head size. The perfectly flat surface melts at $5$
$0.03$ debye, the asymmetric rippled surfaces melt at $8$ $0.04$
$0.03$ debye, the symmetric rippled surfaces melt at $10$ $0.04$
debye. The ordering of the tails is the same as the ordering of the
dipoles except for the flat phase. Since the surface is already
perfect flat, the order parameter does not change much until the
strength of the dipole is $15$ debye. However, the order parameter
decreases quickly when the strength of the dipole is further
increased. The head groups of the lipid molecules are brought closer
by stronger interactions between them. For a flat surface, a large
amount of free volume between the head groups is available, but when
the head groups are brought closer, the tails will splay outward,
forming an inverse micelle. When $\sigma_h=1.28\sigma_0$, the $P_2$
order parameter decreases slightly after the strength of the dipole is
increased to $16$ debye. For rippled surfaces, there is less free
volume available between the head groups. Therefore there is little
effect on the structure of the membrane due to increasing dipolar
strength. However, the increase of the $P_2$ order parameter implies
the membranes are flatten by the increase of the strength of the
dipole. Unlike other systems that melt directly when the interaction
is weak enough, for $\sigma_h=1.41\sigma_0$, part of the membrane
melts into itself first. The upper leaf of the bilayer becomes totally
interdigitated with the lower leaf. This is different behavior than
what is exhibited with the interdigitated lines in the rippled phase
where only one interdigitated line connects the two leaves of bilayer.
\begin{figure}[htb]
\centering
\includegraphics[width=\linewidth]{sP2} 
\caption{The $P_2$ order parameter as a funtion of the strength of the
dipole.\label{fig:sP2}}
\end{figure}

Figure \ref{fig:tP2} shows the dependence of the order parameter on
temperature. The behavior of the $P_2$ order paramter is
straightforward. Systems are more ordered at low temperature, and more
disordered at high temperatures. When the temperature is high enough,
the membranes are instable. For flat surfaces ($\sigma_h=1.20\sigma_0$
and $\sigma_h=1.28\sigma_0$), when the temperature is increased to
$310$, the $P_2$ order parameter increases slightly instead of
decreases like ripple surface. This is an evidence of the frustration
of the dipolar ordering in each leaf of the lipid bilayer, at low
temperature, the systems are locked in a local minimum energy state,
with increase of the temperature, the system can jump out the local
energy well to find the lower energy state which is the longer range
orientational ordering. Like the dipolar ordering of the flat
surfaces, the ordering of the tails of the lipid molecules for ripple
membranes ($\sigma_h=1.35\sigma_0$ and $\sigma_h=1.41\sigma_0$) also
show some nonthermal characteristic. With increase of the temperature,
the $P_2$ order parameter decreases firstly, and increases afterward
when the temperature is greater than $290 K$. The increase of the
$P_2$ order parameter indicates a more ordered structure for the tails
of the lipid molecules which corresponds to a more flat surface. Since
our model lacks the detailed information on lipid tails, we can not
simulate the fluid phase with melted fatty acid chains. Moreover, the
formation of the tilted $L_{\beta'}$ phase also depends on the
organization of fatty groups on tails.
\begin{figure}[htb]
\centering
\includegraphics[width=\linewidth]{tP2} 
\caption{The $P_2$ order parameter as a funtion of
temperature.\label{fig:tP2}}
\end{figure}

\section{Discussion}
\label{sec:discussion}

\newpage
\bibliography{mdripple}
\end{document}
Revision:	3200
Committed:	Fri Jul 27 21:59:45 2007 UTC (16 years, 11 months ago) by gezelter
Content type:	application/x-tex
File size:	31166 byte(s)
Log Message:	friday edits