release_0_1/ripplePaper/ripple.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb} 
\usepackage{endfloat}
\usepackage{epsf}
\usepackage{berkeley}
\usepackage{graphicx}
\usepackage[ref]{overcite}
\usepackage{tabularx}
\pagestyle{plain}
\pagenumbering{arabic}
\oddsidemargin 0.0cm \evensidemargin 0.0cm
\topmargin -21pt \headsep 10pt
\textheight 9.0in \textwidth 6.5in
\brokenpenalty=10000
\renewcommand{\baselinestretch}{1.2}
\renewcommand\citemid{\ } % no comma in optional reference note

\begin{document}

\title{Ripple Phase of the Lipid Bilayers: A Monte Carlo Simulation}
\author{Xiuquan Sun and J. Daniel Gezelter\footnote{Corresponding author. Email: gezelter@nd.edu} \\ 
Department of Chemistry and Biochemistry \\ 
University of Notre Dame \\ 
Notre Dame, Indiana 46556}

\maketitle

\begin{abstract}
The molecular explanation for the origin and properties of the ripple
phase is addressed in this paper. A model which contains the surface
tension and dipole-dipole interactions is used to describe the
potential for a monolayer of simple point dipoles. The simulations are
carried out using Monte Carlo method. It is shown asymmetry of the
translational freedom of the dipoles breaks the symmetry of the
hexagonal lattice and allow antiferroelectric ordering of the
dipoles. The existence of the ripples only depends on the dipolar
property of the system. The structure of the ripples is affected by
surface tension. Only close to the hexagonal lattice, can the ripple
phase be reached. Surface has the lowest transition temperature on
hexagonal lattice elucidates the reason of the existence of the ripple
phase in organism. A mechanism for the phase transition of the lipid
bilayer is proposed.
\end{abstract}

\section{Introduction}
\label{Int}
\indent
Fully hydrated lipids will aggregate spontaneously to form bilayers
which exhibit a variety of phases according to temperature and
composition. Among these phases, a periodic rippled
phase---($P_{\beta'}$) phase is found as an intermediate phase during
the phase transition. This ripple phase can be obtained through either
cooling the lipids from fluid ($L_{\beta'}$) phase or heating from gel
($L_\beta$) phase.  The ripple phase attracts lots of researches from
chemists in the past 30 years. Most structural information of the
ripple phase was obtained by the X-ray diffraction and freeze-fracture
electron microscopy
(FFEM)\cite{Copeland80,Meyer96,Sun96,Katsaras00}. Recently, atomic
force microscopy (AFM) is used as one of these
tools\cite{Kaasgaard03}. All these experimental results strongly
support a 2-Dimensional hexagonal packing lattice for the ripple phase
which is different to the gel phase.  Numerous models were built to
explain the formation of the ripple
phase\cite{Goldstein88,McCullough90,Lubensky93,Tieleman96,Misbah98,Heimburg00,Kubica02,Banerjee02}. However,
the origin of the ripple phase is still on debate.  The behavior of
the dipolar materials in the bulk attracts lots of
interests\cite{Luttinger46,Weis92,Ayton95,Ayton97}. The
ferroelectric state is observed for this kind of system, however, the
frustrated state is found in the 2-D hexagonal lattice of the dipolar
materials, the long range orientational ordered state can not be
formed in this situation. The experimental results show that the
periodicity of the ripples is in the range of 100-600 \AA
\cite{Kaasgaard03}, it is a pretty long range ordered state. So, we
may ask ourselves: {\it ``How could this long range ordered state be
formed in a hexagonal lattice surface?''} We addressed this problem
for a dipolar monolayer using Monte Carlo (MC) simulation.

\section{Model and calculation method}
\label{Mod}

The model used in our simulations is shown in Fig. \ref{fmod1} and Fig. \ref{fmod2}. 
\begin{figure}
\centering
\includegraphics[width=\linewidth]{picture/lattice.eps}
\caption{The modified X-Y-Z model in the simulations. The dipoles are
represented by the arrows. Dipoles are locked to the lattice points
in x-y plane and connect to their nearest neighbors with harmonic
potentials.}
\label{fmod1}
\end{figure}
\begin{figure}
\includegraphics[width=\linewidth]{picture/xyz.eps}
\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
the initial x-y plane, $\theta_i$ is the angle that the dipole is away
from the z axis and $\phi_i$ is the angle between the projection of
the dipole on x-y plane with the x axis.}
\label{fmod2}
\end{figure}
The lipids are represented by the simple point-dipole. During the
simulations, dipoles are locked (in the x-y plane) to lattice points
of hexagonal (or distorted) lattice. Each dipole can move freely out
of the plane and has complete orientational freedom. This is a
modified X-Y-Z model with translational freedom along the z-axis. The
potential of the system
\begin{equation}
V = \sum _i {\sum _{j\in NN_i}^6 {{\frac{k_r}{2}} (r_{ij}-r_0)^2}} +
V_{\text
{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
\label{tp}
\end{equation}
where 
\[ \sum _i {\sum _{j\in NN_i}^6 {{\frac{k_r}{2}} (r_{ij}-r_0)^2}} \]
and 
\[ V_{\text {dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) =  \sum _i {\sum _{j>i} {{\frac{|\mu_i||\mu_j|}{4\pi \epsilon_0 r_{ij}^3}} \biggl[ {\boldsymbol{\hat u}_i} \cdot {\boldsymbol{\hat u}_j} - 3({\boldsymbol{\hat u}_i} \cdot {\mathbf{\hat r}_{ij}})({\boldsymbol{\hat u}_j} \cdot {\mathbf{\hat r}_{ij}}) \biggr]}} \]
are the surface tension and the dipole-dipole interactions. In our
simulation, the surface tension for every dipole is represented by the
harmonic potential with its six nearest neighbors. $r_{ij}$ is the
distance between dipole $i$ and dipole $j$, $r_0$ is the lattice
distance in the x-y plane between dipole $i$ and $j$, $k_r$ is the
surface energy and corresponds to $k_BT$, $k_B$ is the Bolzmann's
constant. For the dipole-dipole interaction part, $\mathbf{r}_{ij}$ is
the vector starting at atom $i$ pointing towards $j$, and
$\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ are the
orientational degrees of freedom for atoms $i$ and $j$
respectively. The magnitude of the dipole moment of atom $i$ is
$|\mu_i|$ which is referred as the strength of the dipole $s$,
$\boldsymbol{\hat{u}}_i$ is the standard unit orientation vector of
$\boldsymbol{\Omega}_i$, and $\mathbf{\hat{r}}_{ij}$ is the unit
vector pointing along $\mathbf{r}_{ij}$
($\mathbf{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$). The unit
of the temperature ($T$) is $kelvin$, the strength of the dipole ($s$)
is $Debye$, the surface energy ($k_r$) is $k_B$---Bolzmann's
constant. For convenience, we will omit the units in the following
discussion.  The order parameter $P_2$ is defined as $1.5 \times
\lambda_{max}$, where $\lambda_{max}$ is the largest eigenvalue of the
matrix $\mathsf S$
\begin{equation}
{\mathsf{S}} = 
        \begin{pmatrix} 
        u_{x}u_{x}-\frac{1}{3} & u_{x}u_{y} & u_{x}u_{z} \\
        u_{y}u_{x} & u_{y}u_{y}-\frac{1}{3} & u_{y}u_{z} \\
        u_{z}u_{x} & u_{z}u_{y} & u_{z}u_{z}-\frac{1}{3} 
        \end{pmatrix},
\label{opmatrix}
\end{equation}
and $u_{\alpha}$ is the $\alpha$ element of the dipole moment averaged
over all particles and configurations. $P_2$ will be $1.0$ for a
perfect ordered system or $0$ for a random one. Note this order
parameter is not equal to the polarization of the system, for example,
the polarization of the perfect antiferroelectric system is $0$, but
$P_2$ is $1.0$. The eigenvector of this matrix is the direction axis
which can detect the direction of the dipoles. The periodicity and
amplitude of the ripples is given by the fast Fourier transform (FFT)
of the perpendicular axis of the direction axis.  To detect the
lattice of the system, $\gamma = {aLat}/{bLat}$ is defined, where
$aLat$, $bLat$ are the lattice distance in X and Y direction
respectively. $\gamma = \sqrt 3$ for the hexagonal lattice. The length
of the monolayer in X axis is $20 \times aLat$ and the system is
roughly square. The average distance that dipoles are from their six
nearest neighbors is $7$ \AA. So, for the hexagonal lattice, the size
of the monolayer is about $250$ \AA $\times$ $250$ \AA \ which is
large enough for the formation of some types of the ripples. In all
simulations, $10^8$ Monte Carlo moves are attempted, the results are
judged by standard Metropolis algorithm. Periodic boundary condition
are used. The cutoff for the long range dipole-dipole interactions is
set to 30 \AA.
%The $P_2$ order parameter allows us to measure the amount of
%directional ordering that exists in the bodies of the molecules making
%up the bilayer. Each lipid molecule can be thought of as a cylindrical
%rod with the head group at the top. If all of the rods are perfectly
%aligned, the $P_2$ order parameter will be $1.0$. If the rods are
%completely disordered, the $P_2$ order parameter will be 0. For a
%collection of unit vectors pointing along the principal axes of the
%rods, the $P_2$ order parameter can be solved via the following
%method.\cite{zannoni94}
%
%Define an ordering tensor $\overleftrightarrow{\mathsf{Q}}$, such that,
%
%where the $u_{i\alpha}$ is the $\alpha$ element of the unit vector
%$\mathbf{\hat{u}}_i$, and the sum over $i$ averages over the whole
%collection of unit vectors. This allows the tensor to be written:
%\begin{equation}
%\overleftrightarrow{\mathsf{Q}} = \frac{1}{N}\sum_i^N \biggl[
%       \mathbf{\hat{u}}_i \otimes \mathbf{\hat{u}}_i 
%       - \frac{1}{3} \cdot \mathsf{1} \biggr].
%\label{lipidEq:po2}
%\end{equation}
%
%After constructing the tensor, diagonalizing
%$\overleftrightarrow{\mathsf{Q}}$ yields three eigenvalues and
%eigenvectors. The eigenvector associated with the largest eigenvalue,
%$\lambda_{\text{max}}$, is the director axis  for the system of unit
%vectors. The director axis is the average direction all of the unit vectors
%are pointing. The $P_2$ order parameter is then simply
%\begin{equation}
%\langle P_2 \rangle = \frac{3}{2}\lambda_{\text{max}}.
%\label{lipidEq:po3}
%\end{equation}
%
%\begin{figure}
%\begin{center}
%\includegraphics[scale=0.3]{/home/maul/gezelter/xsun/Documents/ripple/picture/lattice.eps}
%\caption{ The lattice\label{lat}}
%\end{center}
%\end{figure}

\section{Results and discussion}
\label{Res}

\subsection{Hexagonal}
\label{Hex}
%Fig. \ref{frip} shows the typical simulation results for the hexagonal system when $T = 300$, $s = 7$, $k_r = 0.1$.
%\begin{figure}
%\centering
%\epsfbox{/home/maul/gezelter/xsun/Documents/ripple/picture/rippletop.eps}
%\epsfbox{/home/maul/gezelter/xsun/Documents/ripple/picture/rippleside.eps}
%\caption{A snapshot of our simulation results. The filled circle indicates the position of the dipole, the tail attached on it points out the direction of the dipole. (a)Top view of the monolayer. (b)Side view of the monolayer}
%\label{frip}
%\end{figure}
From the results of the simulation at $T = 300$ for hexagonal lattice, the system is in an antiferroelectric state. Every $3$ or $4$ arrows of the dipoles form a strip whose direction is opposite to its neighbors. $P_2$ is about $0.7$. The ripple is formed clearly. The simulation results shows the ripple has equal opportunity to be formed along different directions, this is due to the isotropic property of the hexagonal lattice.
We use the last configuration of this simulation as the initial condition to increase the system to $T = 400$ every $10\ kelvin$, at the same time decrease the temperature to $T = 100$ every $10\ kelvin$. The trend that the order parameter varies with temperature is plotted in Fig. \ref{t-op}.
\begin{figure}
\begin{center}
\includegraphics[width=\linewidth]{picture/hexorderpara.eps}
\caption{ The orderparameter $P_2$ vs temperature T at hexagonal lattice.\label{t-op}}
\end{center}
\end{figure}
The $P_2 \approx 0.9$ for $T = 100$ implies that the system is in a
highly ordered state. As the temperature increases, the order
parameter is decreasing gradually before $T = 300$, from $T = 310$ the
order parameter drops dramatically, get to nearly $0$ at $T =
400$. This means the system reaches a random state from an ordered
state. The phase transition occurs at $T \approx 340$. At the
temperature range the ripples formed, the structure is fairly stable
with the temperature changing, we can say this structure is in one of
the energy minimum of the energy surface. The amplitude of the ripples
is around $15$ \AA. With the temperature changing, the amplitude of
the ripples is stable also. This is contrast with our general
knowledge that ripples will increase with thermal energy of the system
increasing.  To understand the origin and property of the ripples, we
need look at the potential of our system, which is $V = V_{\text
{surface tension}} + V_{\text {dipole}}$. There are two parts of
it. The intense of the $V_{\text {surface tension}}$ is controlled by
$k_r$ which is the surface energy, and the intense of the $ V_{\text
{dipole}}$ is controlled by $s$ which is the strength of the
dipoles. So, according to adjusting these two parameters, we can get
the further insight into this problem.  At first, we fixed the value
of $s = 7$, and vary $k_r$, the results are shown in
Fig. \ref{kr-a-hf}.
\begin{figure}
\begin{center}
\includegraphics[width=\linewidth]{picture/kr_amplitude.eps}
\caption{ The amplitude $A$ of the ripples at $T = 300$ vs $k_r$ for hexagonal lattice. The height fluctuation $h_f$ vs $k_r$ is shown inset for the same situation.\label{kr-a-hf}}
\end{center}
\end{figure}
When $k_r < 0.1$, due to the small surface tension part, the dipoles
can go far away from their neighbors, lots of noise make the ripples
undiscernable. So, we start from $k_r = 0.1$. However, even when $k_r
= 0$, which means the surface tension is turned off, the
antiferroelectric state still can be reached. This strongly supports
the dipole-dipole interaction is the major driving force to form the
long range orietational ordered state. From Fig. \ref{kr-a-hf}, the
amplitude decreases as the $k_r$ increasing, actually, when
$k_r > 0.7$, although the FFT results still show the values of amplitudes, 
the ripples disappear. From the inset of the
Fig. \ref{kr-a-hf}, the trend of the fluctuation of height of the dipoles---$h_f$
with $k_r$ is similar to the amplitude.
Here $h_f = < h^2 > - {< h >}^2$, $h$ is the $z$
coordinate of the dipoles, $<>$ means $h$ averaged by all dipoles and
configurations. The decreasing of the height fluctuation is due to the
increasing of the surface tension with increasing the $k_r$.
No ripple is observed
when $k_r > 0.7$. When $k_r > 0.7$, the surface tension part of the total
potential of the system dominate the structure of the monolayer, the
dipoles will be kept as near as possible with their neighbors, the
whole system is fairly flat under this situation, and the ripples
disappear.  Then we investigate the role of the dipole-dipole
interactions by fixing the $k_r$ to be $0.1$. This long range
orientational ordered state is very sensitive to the value of $s$ for
hexagonal lattice. For $s = 6$, only local orientational ordering
occurs, when $s$ is even smaller, the system is on a random state. For
$s \geq 9$, the system enters a frustrated state, the amplitude is
hard to tell, however, from observation, the amplitude does not change
too much. We will fully discuss this problem using a distorted
hexagonal lattice.  In brief, asymmetry of the translational freedom
of the dipoles breaks the symmetry of the hexagonal lattice and allow
antiferroelectric ordering of the dipoles.  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 surface tension and
dipole-dipole interaction.

\subsection{Non-hexagonal}
\label{Nhe}
We also investigate the effect of lattice type by changing
$\gamma$. The antiferroelectric state is accessible for all $\gamma$
we use, and will melt with temperature increasing, unlike hexagonal
lattice, the distorted hexagonal lattices prefer a particular director
axis due to their anisotropic property. The phase diagram for this
system is shown in Fig. \ref{phase}.
\begin{figure}
\begin{center}
\includegraphics[width=\linewidth]{picture/phase.eps}
\caption{ The phase diagram with temperature $T$ and lattice variable $\gamma$. The enlarged view near the hexagonal lattice is shown inset.\label{phase}}
\end{center}
\end{figure}
$T_c$ is the transition temperature. The hexagonal lattice has the
lowest $T_c$, and $T_c$ goes up with lattice being more
distorted. There is only two phases in our diagram. When we do
annealing for all the system, the antiferroelectric phase is fairly
stable, although the spin glass is accessible for $\gamma \leq
\sqrt{3}$ if the simulations is started from the random initial
configuration. So, we consider the antiferroelectric phase as a local
minimum energy state even at low temperature. From the inset of
Fig. \ref{phase}, at the hexagonal lattice, $T_c$ changes
quickly. $T_c$ increases more quickly for $\gamma$ getting larger than
$\gamma$ getting smaller. The reason is that: although the average
distance between dipole and its neighbors is same for all types of
lattices, $V_\text{dipole} \propto 1/r_{ij}^3$ in our model, the
change of the lattice spacing in one direction is more effective than
another in this range of $\gamma$. There is another type of
antiferroelectric state when the lattice is far away from the
hexagonal one. Unlike the antiferroelectric state of the hexagonal
lattice which is composed of the strips that have $3$ or $4$ rows of
same direction dipoles, the strips in this type of antiferroelectric
state have $1$, $2$ or $3$ rows of same direction dipoles. In our
phase diagram, this difference is not shown.  However, only when
$\gamma$ is close to $\sqrt{3}$, the long range spatial
ordering---ripple is still maintained. The surface is flat when
$\gamma \ll \sqrt{3}$, and randomly fluctuate due to the appearance of
another type antiferroelectric state when $\gamma \gg \sqrt{3}$. The
change of the lattice type changes the contribution of the surface
tension and the dipole-dipole interaction for the total potential of
the system. For $\gamma \ll \sqrt{3}$, the total potential is
dominated by the surface tension part, so, the surface is flat. For
$\gamma \gg \sqrt{3}$, the total potential is dominated by the
dipole-dipole interaction part, it is very easy to introduce too much
noise to make the ripples undiscernable. In our simulations, the
amplitude of the ripples for distorted hexagonal lattice is larger
than that for hexagonal lattice in the small range around the
hexagonal lattice. The reason is still not clear. A possible
explanation is that the distribution of the dipole-dipole interaction
through the surface is anisotropic in the distorted hexagonal
lattice. Another possibility is that the hexagonal lattice has many
translational local minimum, it has not entered the more rippled state
for our reasonable simulation period.  We investigate the effect of
the strength of the dipole $s$ to the amplitude of the ripples for
$\gamma = 1.875$, $k_r = 0.1$, $T = 260$. Under this situation, the
system reaches the equilibrium very quickly, and the ripples are
fairly stable. The results are shown in Fig. \ref{samplitude}.
\begin{figure}
\begin{center}
\includegraphics[width=\linewidth]{picture/samplitude.eps}
\caption{ The amplitude of ripples $A$ at $T = 260$ vs strength of dipole $s$ for $\gamma = 1.875$. The orderparameter $P_2$ vs $s$ is shown inset at the same situation.\label{samplitude}} 
\end{center}
\end{figure}
For small $s$, there is no long range ordering in the system, so, we
start from $s = 7$, and we use the rippled state as the initial
configuration for all the simulations to reduce the noise. There is no
considerable change of the amplitude in our simulations. At first, the
system is under the competition of the surface tension and
dipole-dipole interactions, increasing $s$ will make the dipole-dipole
interactions more contribute to the total potential and the amplitude
of the ripples is increased a little bit. After the total potential is
totally dominated by the dipole-dipole interactions, the amplitude
does not change too much. This result indicates that the ripples are
the natural property of the dipolar system, the existence of the
ripples does not depend on the surface tension. The orderparameter
increases with increasing the strength of the dipole.

\section{Conclusion}
\label{Con}
In conclusion, the molecular explanation of the origin of the long
range ordering of the hexagonal lattice is given by our
simulations. Asymmetry of the translational freedom of the dipoles
breaks the symmetry of the hexagonal lattice and allow
antiferroelectric ordering of the dipoles.  The simulation results
demonstrate that the dipole-dipole interaction is the major driving
force for the long range orientational ordered state.  According to
the study of the effect of the surface tension and the dipole-dipole
interaction, we find ripples are the natural property of the dipolar
system. Its existence does not depend on the surface tension, however,
a stable, smooth ripple phase is a result of the competition between
surface tension and dipole-dipole interaction, and when surface
tension is large enough to dominate the total potential, the amplitude
of the ripples can be determined by it.  The ripple phase can only be
reached near the hexagonal lattice. Under same condition, the
amplitude of the ripples for hexagonal lattice is smaller than that
for distorted hexagonal lattice. The reason is not clear, however, we
think it is a result of the anisotropic distribution of the
dipole-dipole interaction through the surface in the distorted
hexagonal lattice.  From the phase diagram, the reason of the
existence of the ripple phase in organism is elucidated. To melt at
the body temperature and perform its bio-function, the lipid bilayer
must have a relative low transition temperature which can be realized
near the hexagonal lattice, and the ripple phase is a natural phase
for dipolar system at the hexagonal lattice. So, with the temperature
increasing, the lipid bilayer undergoes a translational adjustment to
enter the ripple phase to lower the transition temperature for the
gel-liquid phase transition, then it can enter the liquid phase even
at a low temperature.

\newpage
\bibliographystyle{jcp}
\bibliography{reference.bib}
\end{document}
Revision:	2139
Committed:	Mon Mar 28 20:49:45 2005 UTC (19 years, 3 months ago)
Content type:	application/x-tex
File size:	22501 byte(s)
Log Message:	This commit was manufactured by cvs2svn to create tag 'release_0_1'.