trunk/mdRipple/mdripple.tex

%\documentclass[aps,pre,twocolumn,amssymb,showpacs,floatfix]{revtex4}
%\documentclass[aps,pre,preprint,amssymb]{revtex4}
\documentclass[12pt]{article}
\usepackage{times}
\usepackage{mathptm}
\usepackage{tabularx}
\usepackage{setspace}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage[ref]{overcite}
\pagestyle{plain}
\pagenumbering{arabic}
\oddsidemargin 0.0cm \evensidemargin 0.0cm
\topmargin -21pt \headsep 10pt
\textheight 9.0in \textwidth 6.5in
\brokenpenalty=10000
\renewcommand{\baselinestretch}{1.2}
\renewcommand\citemid{\ } % no comma in optional reference note

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

\bibliographystyle{achemso}

\title{Dipolar Ordering in Molecular-scale Models of the Ripple Phase
in Lipid Membranes} 
\author{Xiuquan Sun and J. Daniel Gezelter \\
Department of Chemistry and Biochemistry,\\
University of Notre Dame, \\ 
Notre Dame, Indiana 46556}

%\email[E-mail:]{gezelter@nd.edu}

\date{\today} 

\maketitle

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

%\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{equation*}
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{equation*}

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 {eqnarray*}
\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{eqnarray*}
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. \citen{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. \citen{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{equation}
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{equation}
and dipole-dipole,
\begin{equation}
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{equation}
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}

The principal method for observing orientational ordering in dipolar
or liquid crystalline systems is the $P_2$ order parameter (defined
as $1.5 \times \lambda_{max}$, where $\lambda_{max}$ is the largest
eigenvalue of the matrix,
\begin{equation}
{\mathsf{S}} = \frac{1}{N} \sum_i \left(
\begin{array}{ccc}
        u^{x}_i u^{x}_i-\frac{1}{3} & u^{x}_i u^{y}_i & u^{x}_i u^{z}_i \\
        u^{y}_i u^{x}_i  & u^{y}_i u^{y}_i -\frac{1}{3} & u^{y}_i u^{z}_i \\
        u^{z}_i u^{x}_i & u^{z}_i u^{y}_i  & u^{z}_i u^{z}_i -\frac{1}{3} 
\end{array} \right).
\label{eq:opmatrix}
\end{equation}
Here $u^{\alpha}_i$ is the $\alpha=x,y,z$ component of the unit vector
for molecule $i$.  (Here, $\hat{\bf u}_i$ can refer either to the
principal axis of the molecular body or to the dipole on the head
group of the molecule.)  $P_2$ will be $1.0$ for a perfectly-ordered
system and near $0$ for a randomized system.  Note that this order
parameter is {\em not} equal to the polarization of the system.  For
example, the polarization of a perfect anti-ferroelectric arrangement
of point dipoles is $0$, but $P_2$ for the same system is $1$.  The
eigenvector of $\mathsf{S}$ corresponding to the largest eigenvalue is
familiar as the director axis, which can be used to determine a
privileged axis for an orientationally-ordered system.  Since the
molecular bodies are perpendicular to the head group dipoles, it is
possible for the director axes for the molecular bodies and the head
groups to be completely decoupled from each other.

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 head-to-tail
chains is evident with kinks visible at the edges between differently
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 consistent 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 the 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 always observed 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{Sun2007}

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.  Figure \ref{fig:rP2} shows that the $P_2$ order parameter for
the head-group dipoles increases with increasing head group size. When
the heads of the lipid molecules are small, the membrane is nearly
flat. Since the in-plane packing is essentially a close packing of the
head groups, the head dipoles exhibit frustration in their 
orientational ordering.

The ordering trends for the tails are essentially opposite to the
ordering of the head group dipoles. The tail $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 (normal
to the bilayer 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 parameters for head groups (circles) and
molecular bodies (squares) as a function of the ratio of head group
size ($\sigma_h$) to the width of the molecular bodies ($d$). \label{fig:rP2}}
\end{figure}

In addition to varying the size of the head groups, we studied the
effects of the interactions between head groups on the structure of
lipid bilayer by changing the strength of the dipoles.  Figure
\ref{fig:sP2} shows how the $P_2$ order parameter changes with
increasing strength of the dipole.  Generally, the dipoles on the head
groups become more ordered as the strength of the interaction between
heads is increased and become more disordered by decreasing the
interaction stength.  When the interaction between the heads becomes
too weak, the bilayer structure does not persist; all lipid molecules
become dispersed in the solvent (which is non-polar in this
molecular-scale model).  The critial value of the strength of the
dipole depends on the size of the head groups.  The perfectly flat
surface becomes unstable below $5$ Debye, while the  rippled
surfaces melt at $8$ Debye (asymmetric) or $10$ Debye (symmetric). 

The ordering of the tails mirrors the ordering of the dipoles {\it
except for the flat phase}. Since the surface is nearly flat in this
phase, the order parameters are only weakly dependent on dipolar
strength until it reaches $15$ Debye.  Once it reaches this value, the
head group interactions are strong enough to pull the head groups
close to each other and distort the bilayer structure. For a flat
surface, a substantial amount of free volume between the head groups
is normally available.  When the head groups are brought closer by
dipolar interactions, the tails are forced to splay outward, forming
first curved bilayers, and then inverted micelles. 

When $\sigma_h=1.28 d$, the $P_2$ order parameter decreases slightly
when the strength of the dipole is increased above $16$ debye. For
rippled bilayers, there is less free volume available between the head
groups. Therefore increasing dipolar strength weakly influences the
structure of the membrane.  However, the increase in the body $P_2$
order parameters implies that the membranes are being slightly
flattened due to the effects of increasing head-group attraction.

A very interesting behavior takes place when the head groups are very
large relative to the molecular bodies ($\sigma_h = 1.41 d$) and the
dipolar strength is relatively weak ($\mu < 9$ Debye). In this case,
the two leaves of the bilayer become totally interdigitated with each
other in large patches of the membrane.   With higher dipolar
strength, the interdigitation is limited to single lines that run
through the bilayer in a direction perpendicular to the ripple wave
vector. 

\begin{figure}[htb]
\centering
\includegraphics[width=\linewidth]{sP2} 
\caption{The $P_2$ order parameters for head group dipoles (a) and
molecular bodies (b) as a function of the strength of the dipoles.
These order parameters are shown for four values of the head group /
molecular width ratio ($\sigma_h / d$). \label{fig:sP2}}
\end{figure}

Figure \ref{fig:tP2} shows the dependence of the order parameters on
temperature.  As expected, systems are more ordered at low
temperatures, and more disordered at high temperatures.  All of the
bilayers we studied can become unstable if the temperature becomes
high enough.  The only interesting feature of the temperature
dependence is in the flat surfaces ($\sigma_h=1.20 d$ and
$\sigma_h=1.28 d$).  Here, when the temperature is increased above
$310$K, there is enough jostling of the head groups to allow the
dipolar frustration to resolve into more ordered states.  This results
in a slight increase in the $P_2$ order parameter above this
temperature.

For the rippled surfaces ($\sigma_h=1.35 d$ and $\sigma_h=1.41 d$),
there is a slightly increased orientational ordering in the molecular
bodies above $290$K.  Since our model lacks the detailed information
about the behavior of the lipid tails, this is the closest the model
can come to depicting the ripple ($P_{\beta'}$) to fluid
($L_{\alpha}$) phase transition.  What we are observing is a
flattening of the rippled structures made possible by thermal
expansion of the tightly-packed head groups.  The lack of detailed
chain configurations also makes it impossible for this model to depict
the ripple to gel ($L_{\beta'}$) phase transition.

\begin{figure}[htb]
\centering
\includegraphics[width=\linewidth]{tP2} 
\caption{The $P_2$ order parameters for head group dipoles (a) and
molecular bodies (b) as a function of temperature.
These order parameters are shown for four values of the head group /
molecular width ratio ($\sigma_h / d$).\label{fig:tP2}}
\end{figure}

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

The ripple phases have been observed in our molecular dynamic
simulations using a simple molecular lipid model. The lipid model
consists of an anisotropic interacting dipolar head group and an
ellipsoid shape tail. According to our simulations, the explanation of
the formation for the ripples are originated in the size mismatch
between the head groups and the tails. The ripple phases are only
observed in the studies using larger head group lipid models. However,
there is a mismatch betweent the size of the head groups and the size
of the tails in the simulations of the flat surface. This indicates
the competition between the anisotropic dipolar interaction and the
packing of the tails also plays a major role for formation of the
ripple phase. The larger head groups provide more free volume for the
tails, while these hydrophobic ellipsoids trying to be close to each
other, this gives the origin of the spontanous curvature of the
surface, which is believed as the beginning of the ripple phases. The
lager head groups cause the spontanous curvature inward for both of
leaves of the bilayer. This results in a steric strain when the tails
of two leaves too close to each other. The membrane has to be broken
to release this strain. There are two ways to arrange these broken
curvatures: symmetric and asymmetric ripples. Both of the ripple
phases have been observed in our studies. The difference between these
two ripples is that the bilayer is continuum in the symmetric ripple
phase and is disrupt in the asymmetric ripple phase.

Dipolar head groups are the key elements for the maintaining of the
bilayer structure. The lipids are solvated in water when lowering the
the strength of the dipole on the head groups. The long range
orientational ordering of the dipoles can be achieved by forming the
ripples, although the dipoles are likely to form head-to-tail
configurations even in flat surface, the frustration prevents the
formation of the long range orientational ordering for dipoles. The
corrugation of the surface breaks the frustration and stablizes the
long range oreintational ordering for the dipoles in the head groups
of the lipid molecules. Many rows of the head-to-tail dipoles are
parallel to each other and adopt the antiferroelectric state as a
whole. This is the first time the organization of the head groups in
ripple phases of the lipid bilayer has been addressed.

The most important prediction we can make using the results from this
simple model is that if dipolar ordering is driving the surface
corrugation, the wave vectors for the ripples should always found to
be {\it perpendicular} to the dipole director axis.  This prediction
should suggest experimental designs which test whether this is really
true in the phosphatidylcholine $P_{\beta'}$ phases.  The dipole
director axis should also be easily computable for the all-atom and
coarse-grained simulations that have been published in the literature.

Although our model is simple, it exhibits some rich and unexpected
behaviors.  It would clearly be a closer approximation to the reality
if we allowed greater translational freedom to the dipoles and
replaced the somewhat artificial lattice packing and the harmonic
elastic tension with more realistic molecular modeling potentials.
What we have done is to present a simple model which exhibits bulk
non-thermal corrugation, and our explanation of this rippling
phenomenon will help us design more accurate molecular models for
corrugated membranes and experiments to test whether rippling is
dipole-driven or not.

\newpage
\bibliography{mdripple}
\end{document}
Revision:	3202
Committed:	Wed Aug 1 16:07:12 2007 UTC (18 years, 3 months ago) by gezelter
Content type:	application/x-tex
File size:	37013 byte(s)
Log Message:	nearing completion