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 the ripple phases of molecular-scale models
of 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}
Symmetric and asymmetric ripple phases have been observed to form in
molecular dynamics simulations of a simple molecular-scale lipid
model. The lipid model consists of an dipolar head group and an
ellipsoidal tail.  Within the limits of this model, an explanation for
generalized membrane curvature is a simple mismatch in the size of the
heads with the width of the molecular bodies.  The persistence of a
{\it bilayer} structure requires strong attractive forces between the
head groups.  One feature of this model is that an energetically
favorable orientational ordering of the dipoles can be achieved by
out-of-plane membrane corrugation.  The corrugation of the surface
stabilizes the long range orientational ordering for the dipoles in the
head groups which then adopt a bulk anti-ferroelectric state. We
observe a common feature of the corrugated dipolar membranes: the wave
vectors for the surface ripples are always found to be perpendicular
to the dipole director axis.  
\end{abstract}

%\maketitle
\newpage

\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} although Tenchov {\it et al.} have
recently observed near-hexagonal packing in some phosphatidylcholine
(PC) gel phases.\cite{Tenchov2001} The X-ray diffraction work by
Katsaras {\it et al.} showed that a rich phase diagram exhibiting both
{\it asymmetric} and {\it symmetric} ripples is possible for lecithin
bilayers.\cite{Katsaras00}

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 head groups
and molecular tilt with respect to the membrane normal vector can
cause bulk rippling.~\cite{Bannerjee02} Recently, Kranenburg and Smit
described the formation of symmetric ripple-like structures using a
coarse grained solvent-head-tail bead model.\cite{Kranenburg2005}
Their lipids consisted of a short chain of head beads tied to the two
longer ``chains''. 

In contrast, few large-scale molecular modeling 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 time scales 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 head groups and tails is strongly implicated as the primary
driving force for ripple formation, questions about the ordering of
the head groups in ripple phase have 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 of dipolar
elastic membranes is the anti-ferroelectric ordering of the dipoles.
This was evident in the ordering of the dipole director axis
perpendicular to the wave vector of the surface ripples.  A similar
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 structure 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
modeling 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.  Additionally, 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 have placed 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 ($\sigma_h$) 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}$.  

Since the charge separation distance is so large in zwitterionic head
groups (like the PC head groups), it would also be possible to use
either point charges or a ``split dipole'' approximation,
\begin{equation}
V_{ij}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
\hat{r}}_{ij})) = \frac{1}{{4\pi \epsilon_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]
\end{equation}
where $\mu _i$ and $\mu _j$ are the dipole moments of sites $i$ and
$j$, $r_{ij}$ is vector between the two sites, and $R_{ij}$ is given
by,
\begin{equation}
R_{ij}  = \sqrt {r_{ij}^2  + \frac{{d_i^2 }}{4} + \frac{{d_j^2
}}{4}}.
\end{equation}
Here, $d_i$ and $d_j$ are effect charge separation distances
associated with each of the two dipolar sites. 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.

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 similar to the one used by
Marrink {\it et al.}  in their coarse grained simulations of lipid
bilayers.\cite{Marrink04} The 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.  Note that although we
are using larger cutoff and switching radii than Marrink {\it et al.},
our solvent density at 300 K remains at 0.944 g cm$^{-3}$, and the
solvent diffuses at 0.43 $\AA^2 ps^{-1}$ (roughly twice as fast as
liquid 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*}

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

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.  Because the solvent we
are using is non-polar and has a dielectric constant of 1, values for
$\mu$ are sampled from a range that is somewhat smaller than the 20.6
Debye dipole moment of the PC head groups.

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 vacuum 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 equilibrated 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 small constant surface tension
($\gamma=0.015$ N/m). Given the absence of fast degrees of freedom in
this model, a time step 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.  Orientational correlation functions and diffusion
constants were computed from 30 ns simulations in the microcanonical
(NVE) ensemble using the average volume from the end of the constant
pressure and surface tension runs.  The timestep on these final
molecular dynamics runs was 25 fs.  No appreciable changes in phase
structure were noticed upon switching to a microcanonical ensemble.
All simulations were performed using the {\sc oopse} molecular
modeling program.\cite{Meineke05}

A switching function was applied to all potentials to smoothly turn
off the interactions between a range of $22$ and $25$ \AA.  The
switching function was the standard (cubic) function, 
\begin{equation}
s(r) = 
        \begin{cases}
        1 & \text{if $r \le r_{\text{sw}}$},\\
        \frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2}
        {(r_{\text{cut}} - r_{\text{sw}})^3} 
        & \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\
        0 & \text{if $r > r_{\text{cut}}$.}
        \end{cases}
\label{eq:dipoleSwitching}
\end{equation}

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

The membranes in our simulations exhibit a number of interesting
bilayer phases.  The surface topology of these phases depends most
sensitively on the ratio of the size of the head groups to the width
of the molecular bodies.  With heads only slightly larger than the
bodies ($\sigma_h=1.20 d$) the membrane exhibits a flat bilayer. 

Increasing the head / body size ratio increases the local membrane
curvature around each of the lipids.  With $\sigma_h=1.28 d$, the
surface is still essentially flat, but the bilayer starts to exhibit
signs of instability.  We have observed occasional defects where a
line of lipid molecules on one leaf of the bilayer will dip down to
interdigitate with the other leaf.  This gives each of the two bilayer
leaves some local convexity near the line defect.  These structures,
once developed in a simulation, are very stable and are spaced
approximately 100 \AA\ away from each other.

With larger heads ($\sigma_h = 1.35 d$) the membrane curvature
resolves into a ``symmetric'' ripple phase.  Each leaf of the bilayer
is broken into several convex, hemicylinderical sections, and opposite
leaves are fitted together much like roof tiles.  There is no
interdigitation between the upper and lower leaves of the bilayer.

For the largest head / body ratios studied ($\sigma_h = 1.41 d$) the
local curvature is substantially larger, and the resulting bilayer
structure resolves into an asymmetric ripple phase.  This structure is
very similar to the structures observed by both de~Vries {\it et al.}
and Lenz {\it et al.}.  For a given ripple wave vector, there are two
possible asymmetric ripples, which is not the case for the symmetric
phase observed when $\sigma_h = 1.35 d$.

\begin{figure}[htb]
\centering
\includegraphics[width=4in]{phaseCartoon} 
\caption{The role of the ratio between the head group size and the
width of the molecular bodies is to increase the local membrane
curvature.  With strong attractive interactions between the head
groups, this local curvature can be maintained in bilayer structures
through surface corrugation.  Shown above are three phases observed in
these simulations.  With $\sigma_h = 1.20 d$, the bilayer maintains a
flat topology.  For larger heads ($\sigma_h = 1.35 d$) the local
curvature resolves into a symmetrically rippled phase with little or
no interdigitation between the upper and lower leaves of the membrane.
The largest heads studied ($\sigma_h = 1.41 d$) resolve into an
asymmetric rippled phases with interdigitation between the two
leaves.\label{fig:phaseCartoon}} 
\end{figure}

Sample structures for the flat ($\sigma_h = 1.20 d$), symmetric
($\sigma_h = 1.35 d$, and asymmetric ($\sigma_h = 1.41 d$) ripple
phases are shown in Figure \ref{fig:phaseCartoon}.  

It is reasonable to ask how well the parameters we used can produce
bilayer properties that match experimentally known values for real
lipid bilayers.  Using a value of $l = 13.8$ \AA for the ellipsoidal
tails and the fixed ellipsoidal aspect ratio of 3, our values for the
area per lipid ($A$) and inter-layer spacing ($D_{HH}$) depend
entirely on the size of the head bead relative to the molecular body.
These values are tabulated in table \ref{tab:property}.  Kucera {\it
et al.}  have measured values for the head group spacings for a number
of PC lipid bilayers that range from 30.8 \AA (DLPC) to 37.8 (DPPC).
They have also measured values for the area per lipid that range from
60.6
\AA$^2$ (DMPC) to 64.2 \AA$^2$
(DPPC).\cite{NorbertKucerka04012005,NorbertKucerka06012006} Only the
largest of the head groups we modeled ($\sigma_h = 1.41 d$) produces
bilayers (specifically the area per lipid) that resemble real PC
bilayers.  The smaller head beads we used are perhaps better models
for PE head groups.

\begin{table*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{Phase, bilayer spacing, area per lipid, ripple wavelength
and amplitude observed as a function of the ratio between the head
beads and the diameters of the tails.  Ripple wavelengths and
amplitudes are normalized to the diameter of the tail ellipsoids.}
\begin{tabular}{lccccc}
\hline 
$\sigma_h / d$ & type of phase & bilayer spacing (\AA) & area per
lipid (\AA$^2$) & $\lambda / d$ & $A / d$\\
\hline 
1.20 & flat & 33.4 & 49.6 & N/A & N/A \\
1.28 & flat & 33.7 & 54.7 & N/A & N/A \\
1.35 & symmetric ripple & 42.9 & 51.7 & 17.2 & 2.2 \\
1.41 & asymmetric ripple & 37.1 & 63.1 & 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 - 17$ molecular bodies
and the amplitude is $1.5$ molecular bodies for asymmetric ripple and
$2.2$ for symmetric ripple. These values are reasonably consistent
with experimental measurements.\cite{Sun96,Katsaras00,Kaasgaard03}
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), symmetric ripple (middle),
and asymmetric 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.  Note that the flat
membrane has multiple vortex defects in the dipolar ordering, and the
ordering on the lower leaf of the bilayer can be in an entirely
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.35, 1.41 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
anti-ferroelectric 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 strength.  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 critical 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, first forming
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}

Fig. \ref{fig:phaseDiagram} shows a phase diagram for the model as a
function of the head group / molecular width ratio ($\sigma_h / d$)
and the strength of the head group dipole moment ($\mu$).  Note that
the specific form of the bilayer phase is governed almost entirely by
the head group / molecular width ratio, while the strength of the
dipolar interactions between the head groups governs the stability of
the bilayer phase.  Weaker dipoles result in unstable bilayer phases,
while extremely strong dipoles can shift the equilibrium to an
inverted micelle phase when the head groups are small.   Temperature
has little effect on the actual bilayer phase observed, although higher
temperatures can cause the unstable region to grow into the higher
dipole region of this diagram.

\begin{figure}[htb]
\centering
\includegraphics[width=\linewidth]{phaseDiagram} 
\caption{Phase diagram for the simple molecular model as a function
of the head group / molecular width ratio ($\sigma_h / d$) and the
strength of the head group dipole moment
($\mu$).\label{fig:phaseDiagram}}
\end{figure}

We have computed translational diffusion constants for lipid molecules
from the mean-square displacement,
\begin{equation}
D = \lim_{t\rightarrow \infty} \frac{1}{6 t} \langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \rangle,
\end{equation}
of the lipid bodies. Translational diffusion constants for the
different head-to-tail size ratios (all at 300 K) are shown in table
\ref{tab:relaxation}.  We have also computed orientational correlation
times for the head groups from fits of the second-order Legendre
polynomial correlation function,
\begin{equation}
C_{\ell}(t)  =  \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf
\mu}_{i}(0) \right) 
\end{equation}
of the head group dipoles.  The orientational correlation functions
appear to have multiple components in their decay: a fast ($12 \pm 2$
ps) decay due to librational motion of the head groups, as well as
moderate ($\tau^h_{\rm mid}$) and slow ($\tau^h_{\rm slow}$)
components.  The fit values for the moderate and slow correlation
times are listed in table \ref{tab:relaxation}.  Standard deviations
in the fit time constants are quite large (on the order of the values
themselves).

Sparrman and Westlund used a multi-relaxation model for NMR lineshapes
observed in gel, fluid, and ripple phases of DPPC and obtained
estimates of a correlation time for water translational diffusion
($\tau_c$) of 20 ns.\cite{Sparrman2003} This correlation time
corresponds to water bound to small regions of the lipid membrane.
They further assume that the lipids can explore only a single period
of the ripple (essentially moving in a nearly one-dimensional path to
do so), and the correlation time can therefore be used to estimate a
value for the translational diffusion constant of $2.25 \times
10^{-11} m^2 s^{-1}$.  The translational diffusion constants we obtain
are in reasonable agreement with this experimentally determined
value. However, the $T_2$ relaxation times obtained by Sparrman and
Westlund are consistent with P-N vector reorientation timescales of
2-5 ms.  This is substantially slower than even the slowest component
we observe in the decay of our orientational correlation functions.
Other than the dipole-dipole interactions, our head groups have no
shape anisotropy which would force them to move as a unit with
neighboring molecules.  This would naturally lead to P-N reorientation
times that are too fast when compared with experimental measurements.

\begin{table*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{Fit values for the rotational correlation times for the head
groups ($\tau^h$) and molecular bodies ($\tau^b$) as well as the
translational diffusion constants for the molecule as a function of
the head-to-body width ratio (all at 300 K).  In all of the phases,
the head group correlation functions decay with an fast librational
contribution ($12 \pm 1$ ps).  There are additional moderate
($\tau^h_{\rm mid}$) and slow $\tau^h_{\rm slow}$ contributions to
orientational decay that depend strongly on the phase exhibited by the
lipids.  The symmetric ripple phase ($\sigma_h / d = 1.35$) appears to
exhibit the slowest molecular reorientation.}
\begin{tabular}{lcccc}
\hline 
$\sigma_h / d$ & $\tau^h_{\rm mid} (ns)$ & $\tau^h_{\rm
slow} (\mu s)$ & $\tau_b (\mu s)$ & $D (\times 10^{-11} m^2 s^{-1})$ \\
\hline 
1.20 & $0.4$ &  $9.6$ & $9.5$ & $0.43(1)$ \\
1.28 & $2.0$ & $13.5$ & $3.0$ & $5.91(3)$ \\
1.35 & $3.2$ &  $4.0$ & $0.9$ & $3.42(1)$ \\
1.41 & $0.3$ & $23.8$ & $6.9$ & $7.16(1)$ \\
\end{tabular}
\label{tab:relaxation}
\end{center}
\end{minipage}
\end{table*}

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

Symmetric and asymmetric ripple phases have been observed to form in
our molecular dynamics simulations of a simple molecular-scale lipid
model. The lipid model consists of an dipolar head group and an
ellipsoidal tail.  Within the limits of this model, an explanation for
generalized membrane curvature is a simple mismatch in the size of the
heads with the width of the molecular bodies.  With heads
substantially larger than the bodies of the molecule, this curvature
should be convex nearly everywhere, a requirement which could be
resolved either with micellar or cylindrical phases.

The persistence of a {\it bilayer} structure therefore requires either
strong attractive forces between the head groups or exclusionary
forces from the solvent phase.  To have a persistent bilayer structure
with the added requirement of convex membrane curvature appears to
result in corrugated structures like the ones pictured in
Fig. \ref{fig:phaseCartoon}.  In each of the sections of these
corrugated phases, the local curvature near a most of the head groups
is convex.  These structures are held together by the extremely strong
and directional interactions between the head groups.

The attractive forces holding the bilayer together could either be
non-directional (as in the work of Kranenburg and
Smit),\cite{Kranenburg2005} or directional (as we have utilized in
these simulations).  The dipolar head groups are key for the
maintaining the bilayer structures exhibited by this particular model;
reducing the strength of the dipole has the tendency to make the
rippled phase disappear.  The dipoles are likely to form attractive
head-to-tail configurations even in flat configurations, but the
temperatures are high enough that vortex defects become prevalent in
the flat phase.  The flat phase we observed therefore appears to be
substantially above the Kosterlitz-Thouless transition temperature for
a planar system of dipoles with this set of parameters.  For this
reason, it would be interesting to observe the thermal behavior of the
flat phase at substantially lower temperatures.

One feature of this model is that an energetically favorable
orientational ordering of the dipoles can be achieved by forming
ripples.  The corrugation of the surface breaks the symmetry of the
plane, making vortex defects somewhat more expensive, and stabilizing
the long range orientational ordering for the dipoles in the head
groups.  Most of the rows of the head-to-tail dipoles are parallel to
each other and the system adopts a bulk anti-ferroelectric state.  We
believe that this is the first time the organization of the head
groups in ripple phases has been addressed.

Although the size-mismatch between the heads and molecular bodies
appears to be the primary driving force for surface convexity, the
persistence of the bilayer through the use of rippled structures is a
function of the strong, attractive interactions between the heads.
One important prediction we can make using the results from this
simple model is that if the dipole-dipole interaction is the leading
contributor to the head group attractions, the wave vectors for the
ripples should always be found {\it perpendicular} to the dipole
director axis.  This echoes the prediction we made earlier for simple
elastic dipolar membranes, and may suggest experimental designs which
will test whether this is really the case 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.\cite{deVries05}

Although our model is simple, it exhibits some rich and unexpected
behaviors.  It would clearly be a closer approximation to reality if
we allowed bending motions between the dipoles and the molecular
bodies, and if we replaced the rigid ellipsoids with ball-and-chain
tails.  However, the advantages of this simple model (large system
sizes, 50 fs time steps) allow us to rapidly explore the phase diagram
for a wide range of parameters.  Our explanation of this rippling
phenomenon will help us design more accurate molecular models for
corrugated membranes and experiments to test whether or not
dipole-dipole interactions exert an influence on membrane rippling.
\newpage
\bibliography{mdripple}
\end{document}
Revision:	3270
Committed:	Fri Oct 26 22:19:44 2007 UTC (17 years, 9 months ago) by xsun
Content type:	application/x-tex
File size:	48015 byte(s)
Log Message:	fix citation.