trunk/ripplePaper/ripple.tex

\documentclass[aps,endfloats*,preprint,amssymb]{revtex4}
\usepackage{epsfig}
\usepackage{times}
\usepackage{mathptm}

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

\bibliographystyle{pccp}

\title{Symmetry breaking and the $P_{\beta'}$ Ripple phase}
\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. We present a simple (XYZ) spin-lattice model
that allows spins to vary their elevation as well as their
orientation. The extra degree of freedom allows hexagonal lattices of
spins to find states that break out of the normally frustrated
randomized states and are stabilized by long-range antiferroelectric
ordering. To break out of the frustrated states, the spins must form
``rippled'' phases that make the lattices effectively
non-hexagonal. Our XYZ model contains a hydrophobic interaction and
dipole-dipole interactions to describe the interaction potential for
model lipid molecules.  We find non-thermal ripple phases and note
that the wave vectors for the ripples are always perpendicular to the
director axis for the dipoles. Non-hexagonal lattices of dipoles are
not inherently frustrated, and are therefore less likely to form
ripple phases because they can easily form low-energy
antiferroelectric states.  We see that the dipolar order-disorder
transition is substantially lower for hexagonal lattices and that the
ordered phase for this lattice is clearly rippled.
\end{abstract}

\maketitle


\section{Introduction}
\label{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}

Large-scale molecular dynamics simulations have also been performed on
rippled phases using united atom as well as molecular scale
models. De~Vries {\it et al.} studied the structure of lecithin ripple
phases via molecular dynamics and their simulations seem to support
the coexistence models (i.e. fluid-like chain dynamics was observed in
the kink regions).~\cite{deVries05} 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} Brannigan,
Tamboli and Brown have used a molecular scale model to elucidate the
role of molecular shape on membrane phase behavior and
elasticity.~\cite{Brannigan04b} They have also observed a buckled hexatic
phase with strong tail and moderate alignment attractions.~\cite{Brannigan04a}

Ferroelectric states (with long-range dipolar order) can be observed
in dipolar systems with non-hexagonal packings.  However, {\em
hexagonally}-packed 2-D dipolar systems are inherently frustrated and
one would expect a dipolar-disordered phase to be the lowest free
energy configuration.  Concomitantly, it would seem unlikely that a
frustrated lattice in a dipolar-disordered state could exhibit the
long-range periodicity in the range of 100-600 \AA (as exhibited in
the ripple phases studied by Kaasgard {\it et al.}).~\cite{Kaasgaard03}

The various theoretical models have attributed membrane rippling to
various causes which appear contradictory.  We are left with a number
of open questions: 1) Are inter-layer interactions required to explain
the ripple, or can a single bilayer (or even a single leaf) exhibit
the rippling?  2) To what degree is the dipolar anisotropy of the head
group important in determining the rippling?  3) Is chain fluidity
required? (i.e. are the coexistence models necessary to explain the
ripple phenomenon?) 4) How could a state with long-range order be
formed using a substrate consisting of 2-D hexagonally-packed dipolar
molecules?  What we present here is an attempt to find the simplest
model which will exhibit this phenomenon.  We are using a very simple
modified XYZ lattice model; details of the model can be found in
section \ref{sec:model}, results of Monte Carlo simulations using this
model are presented in section
\ref{sec:results}, and section \ref{sec:discussion} contains our conclusions.

\section{The Web-of-Dipoles Model}
\label{sec:model}
The model used in our simulations is shown schematically in
Figs. \ref{fmod1} and \ref{fmod2}.

\begin{figure}[ht]
\centering
\caption{The modified X-Y-Z model used in our simulations. Point dipoles are
represented as arrows. Dipoles are locked to the lattice points in x-y
plane and connect to their nearest neighbors with harmonic
potentials. The lattice parameters $a$ and $b$ are indicated above.}
\includegraphics[width=\linewidth]{picture/WebOfDipoles.eps}
\label{fmod1}
\end{figure}

\begin{figure}[ht]
\centering
\caption{The 6 coordinates describing the state of a 2-dipole system
in our extended X-Y-Z model. $z_i$ is the height of dipole $i$ from 
an arbitrary x-y plane, $\theta_i$ is the angle that the dipole makes
with the laboratory z-axis and $\phi_i$ is the angle between
the projection of the dipole on x-y plane with the x axis.}
\includegraphics[width=\linewidth]{picture/xyz.eps}
\label{fmod2}
\end{figure}

In this model, lipid molecules are represented by point-dipoles (which
is a reasonable approximation to the zwitterionic head groups of the
phosphatidylcholine head groups).  The dipoles are locked in place to
their original lattice sites on the x-y plane.  The original lattice
may be either hexagonal ($a/b = \sqrt{3}$) or non-hexagonal.  However,
each dipole has 3 degrees of freedom.  They may move freely {\em out} of the
x-y plane (along the $z$ axis), and they have complete orientational
freedom ($0 <= \theta <= \pi$, $0 <= \phi < 2 \pi$).  This is a
modified X-Y-Z model with translational freedom along the z-axis.

The potential energy of the system,
\begin{equation}
V = \sum_i \left[ \sum_{j>i} V^{\mathrm{dd}}_{ij} + \frac{1}{2}\sum_{j
\in NN_i}^6 V^{\mathrm{harm}}_{ij} \right]
\end{equation}
The dipolar head groups interact via a traditional point-dipolar
electrostatic potential,
\begin{equation}
V^{\mathrm{dd}}_{ij} = \frac{|\mu|^2}{4\pi \epsilon_0 r_{ij}^3} \left[
{\mathbf{\hat u}_i} \cdot {\mathbf{\hat u}_j} -
3({\mathbf{\hat u}_i} \cdot {\mathbf{\hat
r}_{ij}})({\mathbf{\hat u}_j} \cdot {\mathbf{\hat r}_{ij}})
\right],
\label{eq:vdd}
\end{equation}
and the hydrophobic interactions are approximated with a nearest
neighbor sum of harmonic interactions,
\begin{equation}
V^{\mathrm{harm}}_{ij} = \frac{k_r}{2} \left(r_{ij}-r_0\right)^2
\end{equation}
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}$.  The entire
potential is governed by three parameters, the dipolar strength
($\mu$), the harmonic spring constant ($k_r$) and the preferred
intermolecular spacing ($r_0$).  In practice, we set the value of
$r_0$ to the average inter-molecular spacing from the planar lattice,
yielding a potential model that has only two parameters for a
particular choice of lattice constants $a$ (along the $x$-axis) and $b$
(along the $y$-axis). 

To investigate the phase behavior of this model, we have performed a
series of Metropolis Monte Carlo simulations of moderately-sized (24
nm on a side) patches of membrane hosted on both hexagonal ($\gamma =
a/b = \sqrt{3}$) and non-hexagonal ($\gamma \neq \sqrt{3}$) lattices.
The linear extent of one edge of the monolayer was $20 a$ and the
system was kept roughly square. The average distance that coplanar
dipoles were positioned from their six nearest neighbors was $7$ \AA.
Typical system sizes were 1360 lipids for the hexagonal lattices and
840-2800 lipids for the non-hexagonal lattices.  Periodic boundary
conditions were used, and the cutoff for the dipole-dipole interaction
was set to 30 \AA.  All parameters ($T$, $\mu$, $k_r$, $\gamma$) were
varied systematically to study the effects of these parameters on the
formation of ripple-like phases.

\section{Results and Analysis}
\label{sec:results}
 
\subsection{Dipolar Ordering and Coexistence Temperatures}
The principal method for observing the orientational ordering
transition in dipolar 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 dipole $i$.  $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 the perfect antiferroelectric system is $0$, but
$P_2$ for an antiferroelectric 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 priveleged dipolar
axis for dipole-ordered systems.  Fig. \ref{t-op} shows the values of
$P_2$ as a function of temperature for both hexagonal ($\gamma =
1.732$) and non-hexagonal ($\gamma=1.875$) lattices.

\begin{figure}[ht]
\centering
\caption{The $P_2$ dipolar order parameter as a function of
temperature for both hexagonal ($\gamma = 1.732$) and non-hexagonal
($\gamma = 1.875$) lattices}
\includegraphics[width=\linewidth]{picture/t-orderpara.eps}
\label{t-op}
\end{figure}

There is a clear order-disorder transition in evidence from this data.
Both the hexagonal and non-hexagonal lattices have dipolar-ordered
low-temperature phases, and orientationally-disordered high
temperature phases.  The coexistence temperature for the hexagonal
lattice is significantly lower than for the non-hexagonal lattices,
and the bulk polarization is approximately $0$ for both dipolar
ordered and disordered phases.  This gives strong evidence that the
dipolar ordered phase is antiferroelectric.  We have repeated the
Monte Carlo simulations over a wide range of lattice ratios ($\gamma$)
to generate a dipolar order/disorder phase diagram.  Fig. \ref{phase}
shows that the hexagonal lattice is a low-temperature cusp in the
$T-\gamma$ phase diagram.

\begin{figure}[ht]
\centering
\caption{The phase diagram for the web-of-dipoles model.  The line
denotes the division between the dipolar ordered (antiferroelectric)
and disordered phases.  An enlarged view near the hexagonal lattice is
shown inset.}
\includegraphics[width=\linewidth]{picture/phase.eps}
\label{phase}
\end{figure}

This phase diagram is remarkable in that it shows an antiferroelectric
phase near $\gamma=1.732$ where one would expect lattice frustration
to result in disordered phases at all temperatures.  Observations of
the configurations in this phase show clearly that the system has
accomplished dipolar orderering by forming large ripple-like
structures.  We have observed antiferroelectric ordering in all three
of the equivalent directions on the hexagonal lattice, and the dipoles
have been observed to organize perpendicular to the membrane normal
(in the plane of the membrane).  It is particularly interesting to
note that the ripple-like structures have also been observed to
propagate in the three equivalent directions on the lattice, but the
{\em direction of ripple propagation is always perpendicular to the
dipole director axis}.  A snapshot of a typical antiferroelectric
rippled structure is shown in Fig. \ref{fig:snapshot}.

\begin{figure}[ht]
\centering
\caption{Top and Side views of a representative configuration for the
dipolar ordered phase supported on the hexagonal lattice. Note the
antiferroelectric ordering and the long wavelength buckling of the
membrane.  Dipolar ordering has been observed in all three equivalent
directions on the hexagonal lattice, and the ripple direction is
always perpendicular to the director axis for the dipoles.}
\includegraphics[width=\linewidth]{picture/snapshot.eps}
\label{fig:snapshot}
\end{figure}

\subsection{Discriminating Ripples from Thermal Undulations}

In order to be sure that the structures we have observed are actually
a rippled phase and not merely thermal undulations, we have computed
the undulation spectrum,
\begin{equation}
h(\vec{q}) = A^{-1/2} \int d\vec{r}
h(\vec{r}) e^{-i \vec{q}\cdot\vec{r}}
\end{equation}
where $h(\vec{r})$ is the height of the membrane at location $\vec{r}
= (x,y)$.~\cite{Safran94} In simple (and more complicated) elastic
continuum models, Brannigan {\it et al.} have shown that in the $NVT$
ensemble, the absolute value of the undulation spectrum can be
written,
\begin{equation}
\langle | h(q)|^2 \rangle_{NVT} = \frac{k_B T}{k_c |\vec{q}|^4 +
\tilde{\gamma}|\vec{q}|^2},
\label{eq:fit}
\end{equation}
where $k_c$ is the bending modulus for the membrane, and
$\tilde{\gamma}$ is the mechanical surface
tension.~\cite{Brannigan04b}

The undulation spectrum is computed by superimposing a rectangular
grid on top of the membrane, and by assigning height ($h(\vec{r})$)
values to the grid from the average of all dipoles that fall within a
given $\vec{r}+d\vec{r}$ grid area.  Empty grid pixels are assigned
height values by interpolation from the nearest neighbor pixels.  A
standard 2-d Fourier transform is then used to obtain $\langle |
h(q)|^2 \rangle$.

The systems studied in this paper have relatively small bending moduli
($k_c$) and relatively large mechanical surface tensions
($\tilde{\gamma}$).  In practice, the best fits to our undulation
spectra are obtained by approximating the value of $k_c$ to 0.  In
Fig. \ref{fig:fit} we show typical undulation spectra for two
different regions of the phase diagram along with their fits from the
Landau free energy approach (Eq. \ref{eq:fit}).  In the
high-temperature disordered phase, the Landau fits can be nearly
perfect, and from these fits we can estimate the bending modulus and
the mechanical surface tension.

For the dipolar-ordered hexagonal lattice near the coexistence
temperature, however, we observe long wavelength undulations that are
far outliers to the fits.  That is, the Landau free energy fits are
well within error bars for all other points, but can be off by {\em
orders of magnitude} for a few (but not all) low frequency
components.  

We interpret these outliers as evidence that these low frequency modes
are {\em non-thermal undulations} which is clear evidence that we are
actually seeing a rippled phase developing in this model system.

\begin{figure}[ht]
\centering
\caption{Evidence that the observed ripples are {\em not} thermal
undulations is obtained from the 2-d fourier transform $\langle
|h(\vec{q})|^2 \rangle$ of the height profile ($\langle h(x,y)
\rangle$). Rippled samples show low-wavelength peaks that are
outliers on the Landau free energy fits.  Samples exhibiting only
thermal undulations fit Eq. \ref{eq:fit} remarkably well.}
\includegraphics[width=\linewidth]{picture/fit.eps}
\label{fig:fit}
\end{figure}

\subsection{Effects of Parameters on Ripple Amplitude and Wavelength}

We have used two different methods to estimate the amplitude and
periodicity of the ripples.  The first method requires projection of
the ripples onto a one dimensional rippling axis. Since the rippling
is always perpendicular to the dipole director axis, we can define a
ripple vector as follows.  The largest eigenvector, $s_1$, of the
$\mathsf{S}$ matrix in Eq. \ref{eq:opmatrix} is projected onto a
planar director axis,
\begin{equation}
\vec{d} = \left(\begin{array}{c}
\vec{s}_1 \cdot \hat{i} \\
\vec{s}_1 \cdot \hat{j} \\
0
\end{array} \right).
\end{equation}
($\hat{i}$, $\hat{j}$ and $\hat{k}$ are unit vectors along the $x$,
$y$, and $z$ axes, respectively.)  The rippling axis is in the plane of
the membrane and is perpendicular to the planar director axis,
\begin{equation}
\vec{q}_{\mathrm{rip}} = \vec{d} \times \hat{k}
\end{equation}
We can then find the height profile of the membrane along the ripple
axis by projecting heights of the dipoles to obtain a one-dimensional
height profile, $h(q_{\mathrm{rip}})$. Ripple wavelengths can be
estimated from the largest non-thermal low-frequency component in the
fourier transform of $h(q_{\mathrm{rip}})$.  Amplitudes can be
estimated by measuring peak-to-trough distances in
$h(q_{\mathrm{rip}})$ itself.

A second, more accurate, and simpler method for estimating ripple
shape is to extract the wavelength and height information directly
from the largest non-thermal peak in the undulation spectrum.  For
large-amplitude ripples, the two methods give similar results.  The
one-dimensional projection method is more prone to noise (particularly
in the amplitude estimates for the non-hexagonal lattices).  We report
amplitudes and wavelengths taken directly from the undulation spectrum
below.

In the hexagonal lattice ($\gamma = \sqrt{3}$), the rippling is
observed from $150-300$ K.  The wavelength of the ripples is
remarkably stable at 150~\AA~for all but the temperatures closest to
the order-disorder transition.  At 300 K, the wavelength drops to 120
\AA.  

The dependence of the amplitude on temperature is shown in
Fig. \ref{fig:t-a}.  The rippled structures shrink smoothly as the
temperature rises towards the order-disorder transition.  The
wavelengths and amplitudes we observe are surprisingly close to the
$\Lambda / 2$ phase observed by Kaasgaard {\it et al.} in their work
on PC-based lipids,\cite{Kaasgaard03} although this may be
coincidental agreement given our choice of parameters.

\begin{figure}[ht]
\centering
\caption{ The amplitude $A$ of the ripples vs. temperature for a 
hexagonal lattice.}
\includegraphics[width=\linewidth]{picture/t-a-error.eps}
\label{fig:t-a}
\end{figure}

The ripples can be made to disappear by increasing the internal
surface tension (i.e. by increasing $k_r$).  In Fig. \ref{fig:kr-a}
we show the ripple amplitude as a function of the internal spring
constant for non-dipolar part of the lipid interaction potential.
Weaker ``hydrophobic'' interactions allow the lipid structure to be
dominated by the dipoles, and stronger ``hydrophobic'' interactions
result in much flatter membranes.  Section \ref{sec:discussion}
contains further discussion of this effect.

\begin{figure}[ht]
\centering
\caption{The amplitude $A$ of the ripples vs. the harmonic binding 
constant $k_r$ for both the hexagonal lattice (circles) and
non-hexagonal lattice (squares).  In both simulations the dipole
strength ($\mu$) was 7 Debye and the temperature was 260K.}
\includegraphics[width=\linewidth]{picture/k-a-error.eps}
\label{fig:kr-a}
\end{figure}

The amplitude of the ripples depends critically on the strength of the
dipole moments ($\mu$) in Eq. \ref{eq:vdd}.  If the dipoles are
weakened substantially (below $\mu$ = 5 Debye) at a fixed temperature
of 230 K, the membrane loses dipolar ordering and the ripple
structures. The ripples reach a peak amplitude
of 26~\AA~at a dipolar strength of 9 Debye.  We show the dependence of
ripple amplitude on the dipolar strength in Fig. \ref{fig:s-a}.

\begin{figure}[ht]
\centering
\caption{The amplitude $A$ of the ripples vs. dipole strength ($\mu$)
for both the hexagonal lattice (circles) and non-hexagonal lattice
(squares).  In both simulations the dipole 
strength ($k_r$) was kept constant at a value of $1.987 \times
10^{-4}$ kcal mol$^{-1}$ \AA$^{-2}$.  The temperatures were also kept
fixed at 230K for the hexagonal lattice and 260K for the non-hexagonal
lattice (approximately 2/3 of the order-disorder transition
temperature for each lattice).}
\includegraphics[width=\linewidth]{picture/A-s.eps}
\label{fig:s-a}
\end{figure}

\subsection{Non-hexagonal lattices}

We have also investigated the effect of the lattice geometry by
changing the ratio of lattice constants ($\gamma$) while keeping the
average nearest-neighbor spacing constant. The antiferroelectric state
is accessible for all $\gamma$ values we have used, although the
distorted hexagonal lattices prefer a particular director axis due to
the anisotropy of the lattice.

Our observation of rippling behavior was not limited to the hexagonal
lattices.  In non-hexagonal lattices the antiferroelectric phase can
develop nearly instantaneously in the Monte Carlo simulations, and
these dipolar-ordered phases tend to be remarkably flat.  Whenever
rippling has been observed in these non-hexagonal lattices
(e.g. $\gamma = 1.875$), we see relatively short ripple wavelengths
(98 \AA) and amplitudes of 17 \AA.  These ripples are weakly dependent
on dipolar strength (see Fig. \ref{fig:s-a}), although below a dipolar
strength of 5.5 Debye, the membrane loses dipolar ordering and
displays only thermal undulations.

The rippling in non-hexagonal lattices also shows a strong dependence
on the internal surface tension ($k_r$).  It is possible to make the
ripples disappear by increasing the internal tension.  The low-tension
limit appears to result in somewhat smaller ripples than in the
hexagonal lattice (see Fig. \ref{fig:kr-a}).

The ripple phase does {\em not} appear at all values of $\gamma$.  We
have only observed non-thermal undulations in the range $1.625 <
\gamma < 1.875$.  Outside this range, the order-disorder transition in
the dipoles remains, but the ordered dipolar phase has only thermal
undulations.  This is one of our strongest pieces of evidence that
rippling is a symmetry-breaking phenomenon for hexagonal and
nearly-hexagonal lattices.

\subsection{Effects of System Size}
To evaluate the effect of finite system size, we have performed a 
series of simulations on the hexagonal lattice at a temperature of 300K,
which is just below the order-disorder transition temperature (340K).
These conditions are in the dipole-ordered and rippled portion of the phase
diagram.  These are also the conditions that should be most susceptible to
system size effects.  The wavelength and amplitude of the observed
ripples as a function of system size are shown in Fig. \ref{fig:systemsize}.

\begin{figure}[ht]
\centering
\caption{The ripple wavelength (top) and amplitude (bottom) as a function of
system size for a hexagonal lattice ($\gamma=1.732$) at 300K.}
\includegraphics[width=\linewidth]{picture/SystemSize.eps}
\label{fig:systemsize}
\end{figure}

There is substantial dependence on system size for small (less than 200 \AA) 
periodic boxes.  Notably, there are resonances apparent in the ripple 
amplitudes at box lengths of 121 \AA and 206 \AA.  For larger systems, 
the behavior of the ripples appears to have stabilized and is on a trend to 
slightly smaller amplitudes (and slightly longer wavelenghts) than were 
observed from the 240 \AA box sizes that were used for most of the calculations. 

It is interesting to note that system sizes which are multiples of the
default ripple wavelength can enhance the amplitude of the observed ripples,
but appears to have only a minor effect on the observed wavelength.  It would,
of course, be better to use system sizes that were many multiples of the ripple
wavelength to be sure that the periodic box is not driving the phenomenon, but at
the largest system size studied (485 \AA $\times$ 485 \AA), the number of 
molecules (5440) made long Monte Carlo simulations prohibitively expensive.
We recognize this as a possible flaw of our model for bilayer rippling, but
it is a flaw that will plague any molecular-scale computational model for 
this phenomenon.

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

We have been able to show that a simple lattice model for membranes
which contains only molecular packing (from the lattice), head-group
anisotropy (in the form of electrostatic dipoles) and ``hydrophobic''
interactions (in the form of a nearest-neighbor harmonic potential) is
capable of exhibiting stable long-wavelength non-thermal ripple
structures.  The best explanation for this behavior is that the
ability of the molecules to translate out of the plane of the membrane
is enough to break the symmetry of the hexagonal lattice and allow the
enormous energetic benefit from the formation of a bulk
antiferroelectric phase.  Were the hydrophobic interactions absent
from our model, it would be possible for the entire lattice to
``tilt'' using $z$-translation.  Tilting the lattice in this way would
yield an effectively non-hexagonal lattice which would avoid dipolar
frustration altogether.  With the hydrophobic interactions, bulk tilt
would cause a large strain, and the simplest way to release this
strain is along line defects.  Line defects will result in rippled or
sawtooth patterns in the membrane, and allow small ``stripes'' of
membrane to form antiferroelectric regions that are tilted relative to
the averaged membrane normal.

Although the dipole-dipole interaction is the major driving force for
the long range orientational ordered state, the formation of the
stable, smooth ripples is a result of the competition between the
hydrophobic and dipole-dipole interactions.  This statement is
supported by the variations in both $\mu$ and $k_r$.  Substantially
weaker dipoles or stronger hydrophobic forces can both cause the
ripple phase to disappear.

Molecular packing also plays a role in the formation of the ripple
phase.  It would be surprising if strongly anisotropic head groups
would be able to pack in hexagonal lattices without the underlying
steric interactions between the rest of the molecular bodies.  Since
we only see rippled phases in the neighborhood of $\gamma=\sqrt{3}$,
this implies that there is a role played by the lipid chains in the
organization of the hexagonally ordered phases which support ripples.

Our simple model would clearly be a closer approximation to reality if
we allowed greater translational freedom to the dipoles and replaced
the somewhat artificial lattice packing and the harmonic mimic of the
hydrophobic interaction with more realistic molecular modelling
potentials.  What we have done is to present an extremely simple model
which exhibits bulk non-thermal rippling, and our explanation of the
rippling phenomenon will help us design more accurate molecular models
for the rippling phenomenon.

\clearpage

\bibliography{ripple}

\clearpage

\end{document}
Revision:	3006
Committed:	Tue Sep 19 14:45:15 2006 UTC (18 years, 10 months ago) by gezelter
Content type:	application/x-tex
File size:	29570 byte(s)
Log Message:	latest version