trunk/xDissertation/mc.tex

\chapter{\label{chap:mc}SPONTANEOUS CORRUGATION OF DIPOLAR MEMBRANES}

\section{Introduction}
\label{mc:sec:Int}

The properties of polymeric membranes are known to depend sensitively
on the details of the internal interactions between the constituent
monomers.  A flexible membrane will always have a competition between
the energy of curvature and the in-plane stretching energy and will be
able to buckle in certain limits of surface tension and
temperature.\cite{Safran94} The buckling can be non-specific and
centered at dislocation~\cite{Seung1988} or grain-boundary
defects,\cite{Carraro1993} or it can be directional and cause long
``roof-tile'' or tube-like structures to appear in
partially-polymerized phospholipid vesicles.\cite{Mutz1991}

One would expect that anisotropic local interactions could lead to
interesting properties of the buckled membrane.  We report here on the
buckling behavior of a membrane composed of harmonically-bound, but
freely-rotating electrostatic dipoles.  The dipoles have strongly
anisotropic local interactions and the membrane exhibits coupling
between the buckling and the long-range ordering of the dipoles.

Buckling behavior in liquid crystalline and biological membranes is a
well-known phenomenon.  Relatively pure phosphatidylcholine (PC)
bilayers form a corrugated or ``rippled'' phase ($P_{\beta'}$) which
appears as an intermediate phase between the gel ($L_\beta$) and fluid
($L_{\alpha}$) phases.  The $P_{\beta'}$ phase has attracted
substantial experimental interest over the past 30
years,~\cite{Sun96,Katsaras00,Copeland80,Meyer96,Kaasgaard03} and
there have been a number of theoretical
approaches~\cite{Marder84,Carlson87,Goldstein88,McCullough90,Lubensky93,Misbah98,Heimburg00,Kubica02,Bannerjee02}
(and some heroic
simulations~\cite{Ayton02,Jiang04,Brannigan04a,deVries05,deJoannis06})
undertaken to try to explain this phase, but to date, none have looked
specifically at the contribution of the dipolar character of the lipid
head groups towards this corrugation.  Lipid chain interdigitation
certainly plays a major role, and the structures of the ripple phase
are highly ordered.  The model we investigate here lacks chain
interdigitation (as well as the chains themselves!) and will not be
detailed enough to rule in favor of (or against) any of these
explanations for the $P_{\beta'}$ phase.

Membranes containing electrostatic dipoles can also exhibit the
flexoelectric effect,\cite{Todorova2004,Harden2006,Petrov2006} which
is the ability of mechanical deformations to result in electrostatic
organization of the membrane.  This phenomenon is a curvature-induced
membrane polarization which can lead to potential differences across a
membrane.  Reverse flexoelectric behavior (in which applied currents
effect membrane curvature) has also been observed.  Explanations of
the details of these effects have typically utilized membrane
polarization perpendicular to the face of the
membrane,\cite{Petrov2006} and the effect has been observed in both
biological,\cite{Raphael2000} bent-core liquid
crystalline,\cite{Harden2006} and polymer-dispersed liquid crystalline
membranes.\cite{Todorova2004}

The problem with using atomistic and even coarse-grained approaches to
study membrane buckling phenomena is that only a relatively small
number of periods of the corrugation (i.e. one or two) can be
realistically simulated given current technology.  Also, simulations
of lipid bilayers are traditionally carried out with periodic boundary
conditions in two or three dimensions and these have the potential to
enhance the periodicity of the system at that wavelength.  To avoid
this pitfall, we are using a model which allows us to have
sufficiently large systems so that we are not causing artificial
corrugation through the use of periodic boundary conditions.

The simplest dipolar membrane is one in which the dipoles are located
on fixed lattice sites. Ferroelectric states (with long-range dipolar
order) can be observed in dipolar systems with non-triangular
packings.  However, {\em triangularly}-packed 2-D dipolar systems are
inherently frustrated and one would expect a dipolar-disordered phase
to be the lowest free energy
configuration.\cite{Toulouse1977,Marland1979} Dipolar lattices already
have rich phase behavior, but in order to allow the membrane to
buckle, a single degree of freedom (translation normal to the membrane
face) must be added to each of the dipoles.  It would also be possible
to allow complete translational freedom.  This approach
is similar in character to a number of elastic Ising models that have
been developed to explain interesting mechanical properties in
magnetic alloys.\cite{Renard1966,Zhu2005,Zhu2006,Jiang2006}

What we present here is an attempt to find the simplest dipolar model
which will exhibit buckling behavior.  We are using a modified XYZ
lattice model; details of the model can be found in section
\ref{mc:sec:model}, results of Monte Carlo simulations using this model
are presented in section
\ref{mc:sec:results}, and section \ref{mc:sec:discussion} contains our conclusions.

\section{2-D Dipolar Membrane}
\label{mc:sec:model}

The point of developing this model was to arrive at the simplest
possible theoretical model which could exhibit spontaneous corrugation
of a two-dimensional dipolar medium.  Since molecules in polymerized
membranes and in the $P_{\beta'}$ ripple phase have limited
translational freedom, we have chosen a lattice to support the dipoles
in the x-y plane.  The lattice may be either triangular (lattice
constants $a/b =
\sqrt{3}$) or distorted.  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 essentially a modified X-Y-Z model with translational
freedom along the z-axis.

The potential energy of the system,
\begin{equation}
\begin{split}
V = \sum_i  &\left( \sum_{j>i} \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] \right. \\
  & \left. + \sum_{j \in NN_i}^6 \frac{k_r}{2}\left(
r_{ij}-\sigma \right)^2 \right)
\end{split}
\label{mceq:pot}
\end{equation}

In this potential, $\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 ($\sigma$).  In practice, we set the value of
$\sigma$ 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).  We also define a set of reduced parameters
based on the length scale ($\sigma$) and the energy of the harmonic
potential at a deformation of 2 $\sigma$ ($\epsilon = k_r \sigma^2 /
2$).  Using these two constants, we perform our calculations using
reduced distances, ($r^{*} = r / \sigma$), temperatures ($T^{*} = 2
k_B T / k_r \sigma^2$), densities ($\rho^{*} = N \sigma^2 / L_x L_y$),
and dipole moments ($\mu^{*} = \mu / \sqrt{4 \pi \epsilon_0 \sigma^5
k_r / 2}$).  It should be noted that the density ($\rho^{*}$) depends
only on the mean particle spacing in the $x-y$ plane; the lattice is
fully populated.

To investigate the phase behavior of this model, we have performed a
series of Me\-trop\-o\-lis Monte Carlo simulations of moderately-sized
(34.3 $\sigma$ on a side) patches of membrane hosted on both
triangular ($\gamma = a/b = \sqrt{3}$) and distorted ($\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 1 $\sigma$ (on both triangular and distorted lattices).
Typical system sizes were 1360 dipoles for the triangular lattices and
840-2800 dipoles for the distorted lattices.  Two-dimensional periodic
boundary conditions were used, and the cutoff for the dipole-dipole
interaction was set to 4.3 $\sigma$.  This cutoff is roughly 2.5 times
the typical real-space electrostatic cutoff for molecular systems.
Since dipole-dipole interactions decay rapidly with distance, and
since the intrinsic three-dimensional periodicity of the Ewald sum can
give artifacts in 2-d systems, we have chosen not to use it in these
calculations.  Although the Ewald sum has been reformulated to handle
2-D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89} these
methods are computationally expensive,\cite{Spohr97,Yeh99} and are not
necessary in this case.  All parameters ($T^{*}$, $\mu^{*}$, and
$\gamma$) were varied systematically to study the effects of these
parameters on the formation of ripple-like phases. The error bars in
our results are one $\sigma$ on each side of the average values, where
$\sigma$ is the standard deviation obtained from repeated observations
of many configurations.

\section{Results and Analysis}
\label{mc:sec:results}
 
\subsection{Dipolar Ordering and Coexistence Temperatures}
The principal method for observing the orientational ordering
transition 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{mceq: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 anti-ferroelectric system is $0$, but
$P_2$ for an anti-ferroelectric 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 dipolar
axis for dipole-ordered systems.  The top panel in Fig. \ref{mcfig:phase}
shows the values of $P_2$ as a function of temperature for both
triangular ($\gamma = 1.732$) and distorted ($\gamma=1.875$)
lattices.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/mcPhase.pdf}
\caption[ The $P_2$ dipolar order parameter as
a function of temperature and the phase diagram for the dipolar
membrane model]{\label{mcfig:phase} Top panel: The $P_2$ dipolar order
parameter as a function of temperature for both triangular ($\gamma =
1.732$) and distorted ($\gamma = 1.875$) lattices.  Bottom Panel: The
phase diagram for the dipolar membrane model.  The line denotes the
division between the dipolar ordered (anti-ferroelectric) and
disordered phases.  An enlarged view near the triangular lattice is
shown inset.}
\end{figure}

There is a clear order-disorder transition in evidence from this data.
Both the triangular and distorted lattices have dipolar-ordered
low-temperature phases, and ori\-en\-ta\-tion\-al\-ly-disordered high
temperature phases.  The coexistence temperature for the triangular
lattice is significantly lower than for the distorted 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 anti-ferroelectric.  We have verified that this
dipolar ordering transition is not a function of system size by
performing identical calculations with systems twice as large.  The
transition is equally smooth at all system sizes that were studied.
Additionally, we have repeated the Monte Carlo simulations over a wide
range of lattice ratios ($\gamma$) to generate a dipolar
order/disorder phase diagram.  The bottom panel in
Fig. \ref{mcfig:phase} shows that the triangular lattice is a
low-temperature cusp in the $T^{*}-\gamma$ phase diagram.

This phase diagram is remarkable in that it shows an
anti-ferroelectric 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 ordering by forming
large ripple-like structures.  We have observed anti-ferroelectric
ordering in all three of the equivalent directions on the triangular
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
anti-ferroelectric rippled structure is shown in
Fig. \ref{mcfig:snapshot}.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/mcSnapshot.pdf}
\caption[ Top and Side views of a representative
configuration for the dipolar ordered phase supported on the
triangular lattice]{\label{mcfig:snapshot} Top and Side views of a
representative configuration for the dipolar ordered phase supported
on the triangular lattice. Note the anti-ferroelectric ordering and
the long wavelength buckling of the membrane.  Dipolar ordering has
been observed in all three equivalent directions on the triangular
lattice, and the ripple direction is always perpendicular to the
director axis for the dipoles.}
\end{figure}

Although the snapshot in Fig. \ref{mcfig:snapshot} gives the appearance
of three-row stair-like structures, these appear to be transient.  On
average, the corrugation of the membrane is a relatively smooth,
long-wavelength phenomenon, with occasional steep drops between
adjacent lines of anti-aligned dipoles.

The height-dipole correlation function ($C_{\textrm{hd}}(r, \cos
\theta)$) makes the connection between dipolar ordering and the wave
vector of the ripple even more explicit.  $C_{\textrm{hd}}(r, \cos
\theta)$ is an angle-dependent pair distribution function. The angle
($\theta$) is the angle between the intermolecular vector
$\vec{r}_{ij}$ and direction of dipole $i$,
\begin{equation}
C_{\textrm{hd}} = \frac{\langle \frac{1}{n(r)} \sum_{i}\sum_{j>i}
 h_i \cdot h_j \delta(r - r_{ij}) \delta(\cos \theta_{ij} -
\cos \theta)\rangle} {\langle h^2 \rangle}
\end{equation}
where $\cos \theta_{ij} = \hat{\mu}_{i} \cdot \hat{r}_{ij}$ and
$\hat{r}_{ij} = \vec{r}_{ij} / r_{ij}$.  $n(r)$ is the number of
dipoles found in a cylindrical shell between $r$ and $r+\delta r$ of
the central particle. Fig. \ref{mcfig:CrossCorrelation} shows contours
of this correlation function for both anti-ferroelectric, rippled
membranes as well as for the dipole-disordered portion of the phase
diagram.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/mcHdc.pdf}
\caption[Contours of the height-dipole
correlation function]{\label{mcfig:CrossCorrelation} Contours of the
height-dipole correlation function as a function of the dot product
between the dipole ($\hat{\mu}$) and inter-dipole separation vector
($\hat{r}$) and the distance ($r$) between the dipoles.  Perfect
height correlation (contours approaching 1) are present in the ordered
phase when the two dipoles are in the same head-to-tail line.
Anti-correlation (contours below 0) is only seen when the inter-dipole
vector is perpendicular to the dipoles.  In the dipole-disordered
portion of the phase diagram, there is only weak correlation in the
dipole direction and this correlation decays rapidly to zero for
intermolecular vectors that are not dipole-aligned.}
\end{figure}

The height-dipole correlation function gives a map of how the topology
of the membrane surface varies with angular deviation around a given
dipole.  The upper panel of Fig. \ref{mcfig:CrossCorrelation} shows that
in the anti-ferroelectric phase, the dipole heights are strongly
correlated for dipoles in head-to-tail arrangements, and this
correlation persists for very long distances (up to 15 $\sigma$).  For
portions of the membrane located perpendicular to a given dipole, the
membrane height becomes anti-correlated at distances of 10 $\sigma$.
The correlation function is relatively smooth; there are no steep
jumps or steps, so the stair-like structures in
Fig. \ref{mcfig:snapshot} are indeed transient and disappear when
averaged over many configurations.  In the dipole-disordered phase,
the height-dipole correlation function is relatively flat (and hovers
near zero).  The only significant height correlations are for axial
dipoles at very short distances ($r \approx
\sigma$).

\subsection{Discriminating Ripples from Thermal Undulations}

In order to be sure that the structures we have observed are actually
a rippled phase and not simply 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,Seifert97} In simple (and more complicated)
elastic continuum models, it can 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 q^4 +
\gamma q^2},
\label{mceq:fit}
\end{equation}
where $k_c$ is the bending modulus for the membrane, and $\gamma$ is
the mechanical surface tension.~\cite{Safran94} The systems studied in
this paper have essentially zero bending moduli ($k_c$) and relatively
large mechanical surface tensions ($\gamma$), so a much simpler form
can be written,
\begin{equation}
\langle | h(q) |^2 \rangle_{NVT} = \frac{k_B T}{\gamma q^2}.
\label{mceq:fit2}
\end{equation}

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$.  Alternatively, since the dipoles sit on a Bravais
lattice, one could use the heights of the lattice points themselves as
the grid for the Fourier transform (without interpolating to a square
grid).  However, if lateral translational freedom is added to this
model (a likely extension), an interpolated grid method for computing
undulation spectra will be required.

As mentioned above, the best fits to our undulation spectra are
obtained by setting the value of $k_c$ to 0.  In Fig. \ref{mcfig: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{mceq:fit2}).  In the high-temperature disordered phase, the
Landau fits can be nearly perfect, and from these fits we can estimate
the tension in the surface.  In reduced units, typical values of
$\gamma^{*} = \gamma / \epsilon = 2500$ are obtained for the
disordered phase ($\gamma^{*} = 2551.7$ in the top panel of
Fig. \ref{mcfig:fit}).

Typical values of $\gamma^{*}$ in the dipolar-ordered phase are much
higher than in the dipolar-disordered phase ($\gamma^{*} = 73,538$ in
the lower panel of Fig. \ref{mcfig:fit}).  For the dipolar-ordered
triangular lattice near the coexistence temperature, we also observe
long wavelength undulations that are far outliers to the fits.  That
is, the Landau free energy fits are well within error bars for most of
the other points, but can be off by {\em orders of magnitude} for a
few low frequency components.

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

\begin{figure}
\includegraphics[width=\linewidth]{./figures/mcLogFit.pdf}
\caption[Evidence that the observed ripples are {\em not} thermal
undulations]{\label{mcfig:fit} 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 by an order of
magnitude.  Samples exhibiting only thermal undulations fit
Eq. \ref{mceq:fit} remarkably well.}
\end{figure}

\subsection{Effects of Potential Parameters on 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{mceq: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 distorted lattices).  We report
amplitudes and wavelengths taken directly from the undulation spectrum
below.

In the triangular lattice ($\gamma = \sqrt{3}$), the rippling is
observed for temperatures ($T^{*}$) from $61-122$.  The wavelength of
the ripples is remarkably stable at 21.4~$\sigma$ for all but the
temperatures closest to the order-disorder transition.  At $T^{*} =
122$, the wavelength drops to 17.1~$\sigma$.

The dependence of the amplitude on temperature is shown in the top
panel of Fig. \ref{mcfig:Amplitude}.  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}
However, this is coincidental agreement based on a choice of 7~\AA~as
the mean spacing between lipids. 

\begin{figure}
\includegraphics[width=\linewidth]{./figures/mcProperties_sq.pdf}
\caption[ The amplitude $A^{*}$ of the ripples
vs. temperature and dipole strength
($\mu^{*}$)]{\label{mcfig:Amplitude} a) The amplitude $A^{*}$ of the
ripples vs. temperature for a triangular lattice. b) The amplitude
$A^{*}$ of the ripples vs. dipole strength ($\mu^{*}$) for both the
triangular lattice (circles) and distorted lattice (squares).  The
reduced temperatures were kept fixed at $T^{*} = 94$ for the
triangular lattice and $T^{*} = 106$ for the distorted lattice
(approximately 2/3 of the order-disorder transition temperature for
each lattice).}
\end{figure}

The ripples can be made to disappear by increasing the internal
elastic tension (i.e. by increasing $k_r$ or equivalently, reducing
the dipole moment).  The amplitude of the ripples depends critically
on the strength of the dipole moments ($\mu^{*}$) in Eq. \ref{mceq:pot}.
If the dipoles are weakened substantially (below $\mu^{*}$ = 20) at a
fixed temperature of 94, the membrane loses dipolar ordering
and the ripple structures. The ripples reach a peak amplitude of
3.7~$\sigma$ at a dipolar strength of 25.  We show the dependence
of ripple amplitude on the dipolar strength in
Fig. \ref{mcfig:Amplitude}.

\subsection{Distorted 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 anti-ferroelectric state
is accessible for all $\gamma$ values we have used, although the
distorted triangular lattices prefer a particular director axis due to
the anisotropy of the lattice.

Our observation of rippling behavior was not limited to the triangular
lattices.  In distorted lattices the anti-ferroelectric 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 distorted lattices
(e.g. $\gamma = 1.875$), we see relatively short ripple wavelengths
(14 $\sigma$) and amplitudes of 2.4~$\sigma$.  These ripples are
weakly dependent on dipolar strength (see Fig. \ref{mcfig:Amplitude}),
although below a dipolar strength of $\mu^{*} = 20$, the membrane
loses dipolar ordering and displays only thermal undulations.

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 triangular and
nearly-triangular lattices.

\subsection{Effects of System Size}
To evaluate the effect of finite system size, we have performed a
series of simulations on the triangular lattice at a reduced
temperature of 122, which is just below the order-disorder transition
temperature ($T^{*} = 139$).  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.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/mcSystemSize.pdf}
\caption[The ripple wavelength and amplitude as a function of system
size]{\label{mcfig:systemsize} The ripple wavelength (top) and
amplitude (bottom) as a function of system size for a triangular
lattice ($\gamma=1.732$) at $T^{*} = 122$.}
\end{figure}

There is substantial dependence on system size for small (less than
29~$\sigma$) periodic boxes.  Notably, there are resonances apparent
in the ripple amplitudes at box lengths of 17.3 and 29.5 $\sigma$.
For larger systems, the behavior of the ripples appears to have
stabilized and is on a trend to slightly smaller amplitudes (and
slightly longer wavelengths) than were observed from the 34.3 $\sigma$
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 (70 $\sigma$ $\times$ 70 $\sigma$), the number of dipoles
(5440) made long Monte Carlo simulations prohibitively expensive. 

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

We have been able to show that a simple dipolar lattice model which
contains only molecular packing (from the lattice), anisotropy (in the
form of electrostatic dipoles) and a weak elastic tension (in the form
of a nearest-neighbor harmonic potential) is capable of exhibiting
stable long-wavelength non-thermal surface corrugations.  The best
explanation for this behavior is that the ability of the dipoles to
translate out of the plane of the membrane is enough to break the
symmetry of the triangular lattice and allow the energetic benefit
from the formation of a bulk anti-ferroelectric phase.  Were the weak
elastic tension 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-triangular lattice which
would avoid dipolar frustration altogether.  With the elastic tension
in place, bulk tilt causes a large strain, and the least costly way to
release this strain is between two rows of anti-aligned dipoles.
These ``breaks'' will result in rippled or sawtooth patterns in the
membrane, and allow small stripes of membrane to form
anti-ferroelectric 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
elastic tension and the dipole-dipole interactions.  This statement is
supported by the variation in $\mu^{*}$.  Substantially weaker dipoles
relative to the surface tension can cause the corrugated phase to
disappear.

The packing of the dipoles into a nearly-triangular lattice is clearly
an important piece of the puzzle.  The dipolar head groups of lipid
molecules are sterically (as well as electrostatically) anisotropic,
and would not pack in triangular arrangements without the steric
interference of adjacent molecular bodies.  Since we only see rippled
phases in the neighborhood of $\gamma=\sqrt{3}$, this implies that
even if this dipolar mechanism is the correct explanation for the
ripple phase in realistic bilayers, there would still be a role played
by the lipid chains in the in-plane organization of the triangularly
ordered phases which could support ripples.  The present model is
certainly not detailed enough to answer exactly what drives the
formation of the $P_{\beta'}$ phase in real lipids, but suggests some
avenues for further experiments.

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

Our other observation about the ripple and dipolar directionality is
that the dipole director axis can be found to be parallel to any of
the three equivalent lattice vectors in the triangular lattice.
Defects in the ordering of the dipoles can cause the dipole director
(and consequently the surface corrugation) of small regions to be
rotated relative to each other by 120$^{\circ}$.  This is a similar
behavior to the domain rotation seen in the AFM studies of Kaasgaard
{\it et al.}\cite{Kaasgaard03}   

Although our model is simple, it exhibits some rich and unexpected
behaviors.  It would clearly be a closer approximation to the reality
if we allowed greater translational freedom to the dipoles and
replaced the somewhat artificial lattice packing and the harmonic
elastic tension with more realistic molecular modeling potentials.
What we have done is to present a simple model which exhibits bulk
non-thermal corrugation, and our explanation of this rippling
phenomenon will help us design more accurate molecular models for
corrugated membranes and experiments to test whether rippling is
dipole-driven or not.
Revision:	3383
Committed:	Wed Apr 16 21:56:34 2008 UTC (16 years, 5 months ago) by xsun
Content type:	application/x-tex
File size:	32364 byte(s)
Log Message:	formatting check