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}

\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} = 4\epsilon \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} = \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 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) & & 1 & 3 & 1 \\
$\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.06112 \\
$\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.

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

To create unbiased bilayers, all simulations were started from two
perfectly flat monolayers separated by a 20 \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 boundaries were used, and $480-720$ lipid
molecules were present in the simulations depending on the size of the
head beads.  The two monolayers spontaneously collapse into bilayer
structures within 100 ps, and following this collapse, all systems
were equlibrated for $100$ ns at $300$ K. 

The resulting structures were then solvated at a ratio of $6$ DPD
solvent beads (24 water molecules) per lipid. These configurations
were then equilibrated for another $30$ ns. All simulations with
solvent were carried out at constant pressure ($P=1$ atm) by $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.  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 increasing size of the head groups. The surface is
nearly flat when $\sigma_h=1.20\sigma_0$. With
$\sigma_h=1.28\sigma_0$, 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\ width 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\sigma_0$, 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\sigma_0$, 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, the symmetric ripple has no prefered direction.
The corresponding cartoons are shown in Figure
\ref{fig:phaseCartoon} for elucidation of the detailed structures of
different phases. Figure \ref{fig:phaseCartoon} (a) is the flat phase,
(b) is the asymmetric ripple phase corresponding to the lipid
organization when $\sigma_h=1.28\sigma_0$ and $\sigma_h=1.41\sigma_0$,
and (c) is the symmetric ripple phase observed when
$\sigma_h=1.35\sigma_0$. In symmetric ripple phase, the bilayer is
continuous everywhere on the whole membrane, however, in asymmetric
ripple phase, the bilayer is intermittent domains connected by thin
interdigitated monolayer which consists of upper and lower leaves of
the bilayer.
\begin{table*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{}
\begin{tabular}{lccc}
\hline 
$\sigma_h/\sigma_0$ & type of phase & $\lambda/\sigma_0$ & $A/\sigma_0$\\
\hline 
1.20 & flat & N/A & N/A \\
1.28 & asymmetric 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/\sigma_0$,
reduced amplitude $A/\sigma_0$ of the ripples are summerized in Table
\ref{tab:property}. The wavelength range is $15~21$ and the amplitude
is $1.5$ for asymmetric ripple and $2.2$ for symmetric ripple. These
values are consistent to the experimental results. Note, the
amplitudes are underestimated without the melted tails in our
simulations.

\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 $P_2$ order paramters (for molecular bodies and head group
dipoles) have been calculated to clarify the ordering in these phases
quantatively. $P_2=1$ means a perfectly ordered structure, and $P_2=0$
implies orientational randomization. Figure \ref{fig:rP2} shows the
$P_2$ order paramter of the dipoles on head group rising with
increasing head group size. When the heads of the lipid molecules are
small, the membrane is flat. The dipolar ordering is essentially
frustrated on orientational ordering in this circumstance. Figure
\ref{fig:topView} shows the snapshots of the top view for the flat system
($\sigma_h=1.20\sigma$) and rippled system
($\sigma_h=1.41\sigma$). The pointing direction of the dipoles on the
head groups are represented by two colored half spheres from blue to
yellow. For flat surfaces, the system obviously shows frustration on
the dipolar ordering, there are kinks on the edge of defferent
domains. Another reason is that the lipids can move independently in
each monolayer, it is not nessasory for the direction of dipoles on
one leaf is consistant to another layer, which makes total order
parameter is relatively low. With increasing head group size, the
surface is corrugated, and dipoles do not move as freely on the
surface. Therefore, the translational freedom of lipids in one layer
is dependent upon the position of lipids in another layer, as a
result, the symmetry of the dipoles on head group in one layer is tied
to the symmetry in the other layer. Furthermore, as the membrane
deforms from two to three dimensions due to the corrugation, the
symmetry of the ordering for the dipoles embedded on each leaf is
broken. The dipoles then self-assemble in a head-tail configuration,
and the order parameter increases dramaticaly. However, the total
polarization of the system is still close to zero. This is strong
evidence that the corrugated structure is an antiferroelectric
state. From the snapshot in Figure \ref{}, the dipoles arrange as
arrays along $Y$ axis and fall into head-to-tail configuration in each
line, but every $3$ or $4$ lines of dipoles change their direction
from neighbour lines. The system shows antiferroelectric
charactoristic as a whole. The orientation of the dipolar is always
perpendicular to the ripple wave vector. These results are consistent
with our previous study on dipolar membranes.

The ordering of the tails is essentially opposite to the ordering of
the dipoles on head group. The $P_2$ order parameter decreases with
increasing head size. This indicates the surface is more curved with
larger head groups. When the surface is flat, all tails are pointing
in the same direction; in this case, all tails are parallel to the
normal of the surface,(making this structure remindcent of the
$L_{\beta}$ phase. 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 funtion of the ratio of
$\sigma_h$ to $\sigma_0$.\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:	3199
Committed:	Thu Jul 26 19:47:06 2007 UTC (18 years ago) by gezelter
Content type:	application/x-tex
File size:	29904 byte(s)
Log Message:	Edits