--- trunk/mdRipple/mdripple.tex	2007/06/25 21:16:17	3147
+++ trunk/mdRipple/mdripple.tex	2007/08/01 16:07:12	3202
@@ -1,176 +1,723 @@
 %\documentclass[aps,pre,twocolumn,amssymb,showpacs,floatfix]{revtex4}
-\documentclass[aps,pre,preprint,amssymb,showpacs]{revtex4}
+%\documentclass[aps,pre,preprint,amssymb]{revtex4}
+\documentclass[12pt]{article}
+\usepackage{times}
+\usepackage{mathptm}
+\usepackage{tabularx}
+\usepackage{setspace}
+\usepackage{amsmath}
+\usepackage{amssymb}
 \usepackage{graphicx}
+\usepackage[ref]{overcite}
+\pagestyle{plain}
+\pagenumbering{arabic}
+\oddsidemargin 0.0cm \evensidemargin 0.0cm
+\topmargin -21pt \headsep 10pt
+\textheight 9.0in \textwidth 6.5in
+\brokenpenalty=10000
+\renewcommand{\baselinestretch}{1.2}
+\renewcommand\citemid{\ } % no comma in optional reference note
 
 \begin{document}
-\renewcommand{\thefootnote}{\fnsymbol{footnote}}
-\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}
+%\renewcommand{\thefootnote}{\fnsymbol{footnote}}
+%\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}
 
-%\bibliographystyle{aps}
+\bibliographystyle{achemso}
 
-\title{}
-\author{Xiuquan Sun and J. Daniel Gezelter}
-\email[E-mail:]{gezelter@nd.edu}
-\affiliation{Department of Chemistry and Biochemistry,\\
-University of Notre Dame, \\
+\title{Dipolar Ordering in Molecular-scale Models of the Ripple Phase
+in Lipid Membranes} 
+\author{Xiuquan Sun and J. Daniel Gezelter \\
+Department of Chemistry and Biochemistry,\\
+University of Notre Dame, \\ 
 Notre Dame, Indiana 46556}
 
+%\email[E-mail:]{gezelter@nd.edu}
+
 \date{\today} 
 
-\begin{abstract}
+\maketitle
 
+\begin{abstract}
+The ripple phase in phosphatidylcholine (PC) bilayers has never been
+completely explained. 
 \end{abstract}
 
-\pacs{}
-\maketitle
+%\maketitle
 
-Our idea for developing a simple and reasonable lipid model to study
-the ripple pahse of lipid bilayers 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 dominant numbers of lipid molecules are very rigid in ripple
-phase which allows the details of the lipid molecules neglectable. In
-our model, lipid molecules are represented by rigid bodies made of one
-head sphere with a point dipole sitting on it and one ellipsoid tail,
-the direction of the dipole is fixed to be perpendicular to the
-tail. The breadth and length of tail are $\sigma_0$, $3\sigma_0$. The
-diameter of heads varies from $1.20\sigma_0$ to $1.41\sigma_0$.  The
-model of the solvent in our simulations is inspired by the idea of
-``DPD'' water. Every four water molecules are reprsented by one
-sphere.
+\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}
 
-Spheres interact each other with Lennard-Jones potential, ellipsoids
-interact each other with Gay-Berne potential, dipoles interact each
-other with typical dipole potential, spheres interact ellipsoids with
-LJ-GB potential. All potentials are truncated at $25$ \AA and shifted
-at $22$ \AA.
+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.
 
-To make the simulations less expensive and to observe long-time range
-behavior of the lipid membranes, all simulaitons were started from two
-sepetated monolayers in the vaccum with $x-y$ anisotropic pressure
-coupling, length of $z$ axis of the simulations was fixed to prevent
-the shrinkage of the simulation boxes due to the free volume outside
-of the bilayer, and a constant surface tension was applied to enable
-the fluctuation of the surface. Periodic boundaries were used. There
-were $480-720$ lipid molecules in simulations according to different
-size of the heads. All the simulations were stablized for $100$ ns at
-$300$ K. The resulted structures were solvated in the water (about
-$6$ DPD water/lipid molecule) as the initial configurations for another
-$30$ ns relaxation. All simulations with water were carried out at
-constant pressure ($P=1$bar) by $3$D anisotropic coupling, and
-constant surface tension ($\gamma=0.015$). Time step was
-$50$ fs. Simulations were performed by using OOPSE package.
+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}.
 
-Snap shots show that the membrane is more corrugated with increasing
-the size of the head groups. The surface is nearly perfect flat when
-$\sigma_h$ is $1.20\sigma_0$. At $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 a interdigitated line
-defect. Two periodicities with $100$\AA width were observed in the
-simulation. This structure is very similiar to OTHER PAPER. The same
-structure was also observed when $\sigma_h=1.41\sigma_0$. However, the
-surface of the membrane is corrugated, and the periodicity of the
-connection between upper and lower leaf membrane is shorter. From the
-undulation spectrum of the surface (the exact form is in OUR PREVIOUS
-PAPER), the corrugation is non-thermal fluctuation, and we are
-confident to identify it as the ripple phase. The width of one ripple
-is about $71$ \AA, and amplitude is about $7$ \AA. When
-$\sigma_h=1.35\sigma_0$, we observed another corrugated surface with
-$79$ \AA width and $10$ \AA amplitude. This structure is different to
-the previous rippled surface, there is no connection between upper and
-lower leaf of the bilayer. Each leaf of the bilayer is broken to
-several curved pieces, the broken position is mounted into the center
-of opposite piece in another leaf. Unlike another corrugated surface
-in which the upper leaf of the surface is always connected to the
-lower leaf from one direction, this ripple of this surface is
-isotropic. Therefore, we claim this is a symmetric ripple phase.
+\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}
 
-The $P_2$ order paramter is calculated to understand the phase
-behavior quantatively. $P_2=1$ means a perfect ordered structure, and
-$P_2=0$ means a random structure. The method can be found in OUR
-PAPER. Fig. shows $P_2$ order paramter of the dipoles on head group
-raises with increasing the size of the head group. When head of lipid
-molecule is small, the membrane is flat and shows strong two
-dimensional characters, dipoles are frustrated on orientational
-ordering in this circumstance. 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 the
-size of head group, the surface is being more corrugated, dipoles are
-not allowed to move freely on the surface, they are
-localized. 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
-consistant to the symmetry in another layer. Furthermore, the membrane
-tranlates from a two dimensional system to a three dimensional system
-by the corrugation, the symmetry of the ordering for the two
-dimensional dipoles on the head group of lipid molecules is broken, on
-a distorted lattice, dipoles are ordered on a head to tail energy
-state, the order parameter is increased dramaticly. However, the total
-polarization of the system is close to zero, which is a strong
-evidence it is a antiferroelectric state. The orientation of the
-dipole ordering is alway perpendicular to the ripple vector. These
-results are consistant to our previous study on similar system. The
-ordering of the tails are opposite to the ordering of the dipoles on
-head group, the $P_2$ order parameter decreases with increasing the
-size of head. This indicates the surface is more curved with larger
-head. When surface is flat, all tails are pointing to the same
-direction, in this case, all tails are parallal to the normal of the
-surface, which shares the same structure with $L_{\beta}$ phase. For the
-size of head being $1.28\sigma_0$, the surface starts to splay inward,
-however, the surface is still flat, therefore, although the order
-parameter is lower, it still indicates a very flat surface. Further
-increasing the size of the head, the order parameter drops dramaticly,
-the surface is rippled.
+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{equation*}
+V_{ij}({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat
+r}_{ij}}) = 4\epsilon ({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
+{\mathbf{\hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i},
+{\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^{12}
+-\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
+{\mathbf{\hat r}_{ij}})+\sigma_0}\right)^6\right]
+\label{eq:gb}
+\end{equation*}
 
-We studied the effects of interaction between head groups on the
-structure of lipid bilayer by changing the strength of the dipole. The
-fig. shows the $P_2$ order parameter changing with strength of the
-dipole. Generally the dipoles on the head group are more ordered with
-increasing the interaction between heads and more disordered with
-decreasing the interaction between heads. When the interaction between
-heads is weak enough, the bilayer structure is not persisted any more,
-all lipid molecules are melted in the water. The critial value of the
-strength of the dipole is various for different system. The perfect
-flat surface melts at $5$ debye, the asymmetric rippled surfaces melt
-at $8$ debye, the symmetric rippled surfaces melt at $10$ debye. This
-indicates that the flat phase is the most stable state, the asymmetric
-ripple phase is second stalbe state, and the symmetric ripple phase is
-the most unstable state. 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 group of the lipid molecules are brought closer by
-strenger interaction between them. For a flat surface, a mount of free
-volume between head groups is available, when the head groups are
-brought closer, the surface will splay outward to be a inverse
-micelle. For rippled surfaces, there is few free volume available on
-between the head groups, they can be closer, therefore there are
-little effect on the structure of the membrane. Another interesting
-fact, unlike other systems melting directly when the interaction is
-weak enough, for $\sigma_h$ is $1.41\sigma_0$, part of the membrane
-melts into itself first, the upper leaf of the bilayer is totally
-interdigitated with the lower leaf, this is different with the
-interdigitated lines in rippled phase where only one interdigitated
-line connects the two leaves of bilayer.
+The range $(\sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
+\hat{r}}_{ij}))$, and strength $(\epsilon({\bf \hat{u}}_{i},{\bf
+\hat{u}}_{j},{\bf \hat{r}}_{ij}))$ parameters 
+are dependent on the relative orientations of the two molecules (${\bf
+\hat{u}}_{i},{\bf \hat{u}}_{j}$) as well as the direction of the
+intermolecular separation (${\bf \hat{r}}_{ij}$).  $\sigma$ and
+$\sigma_0$ are also governed by shape mixing and anisotropy variables,
+\begin {eqnarray*}
+\sigma_0 & = & \sqrt{d_i^2 + d_j^2} \\ \\
+\chi & = & \left[ \frac{\left( l_i^2 - d_i^2 \right) \left(l_j^2 -
+d_j^2 \right)}{\left( l_j^2 + d_i^2 \right) \left(l_i^2 +
+d_j^2 \right)}\right]^{1/2} \\ \\
+\alpha^2 & = & \left[ \frac{\left( l_i^2 - d_i^2 \right) \left(l_j^2 +
+d_i^2 \right)}{\left( l_j^2 - d_j^2 \right) \left(l_i^2 +
+d_j^2 \right)}\right]^{1/2},
+\end{eqnarray*}
+where $l$ and $d$ describe the length and width of each uniaxial
+ellipsoid.  These shape anisotropy parameters can then be used to
+calculate the range function,
+\begin{equation*}
+\sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf \hat{r}}_{ij}) = \sigma_{0}
+ \left[ 1- \left\{ \frac{ \chi \alpha^2 ({\bf \hat{u}}_i \cdot {\bf
+\hat{r}}_{ij} ) + \chi \alpha^{-2} ({\bf \hat{u}}_j \cdot {\bf
+\hat{r}}_{ij} ) - 2 \chi^2 ({\bf \hat{u}}_i \cdot {\bf
+\hat{r}}_{ij} )({\bf \hat{u}}_j \cdot {\bf
+\hat{r}}_{ij} ) ({\bf \hat{u}}_i \cdot {\bf \hat{u}}_j)}{1 - \chi^2
+\left({\bf \hat{u}}_i \cdot {\bf \hat{u}}_j\right)^2} \right\}
+\right]^{-1/2}
+\end{equation*}
 
+Gay-Berne ellipsoids also have an energy scaling parameter,
+$\epsilon^s$, which describes the well depth for two identical
+ellipsoids in a {\it side-by-side} configuration.  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. \citen{Luckhurst90}
+has a particularly good explanation of the choice of the Gay-Berne
+parameters $\mu$ and $\nu$ for modeling liquid crystal molecules. An
+excellent overview of the computational methods that can be used to
+efficiently compute forces and torques for this potential can be found
+in Ref. \citen{Golubkov06}
 
-Fig. shows the changing of the order parameter with temperature. The
-behavior of the $P_2$ orderparamter is straightforword. Systems are
-more ordered at low temperature, and more disordered at high
-temperature. When the temperature is high enough, the membranes are
-discontinuted. The structures are stable during the changing of the
-temperature. Since our model lacks the detail information for tails of
-lipid molecules, we did not simulate the fluid phase with a melted
-fatty chains. Moreover, the formation of the tilted $L_{\beta'}$ phase
-also depends on the organization of fatty groups on tails, we did not
-observe it either.
+The choices of parameters we have used in this study correspond to a
+shape anisotropy of 3 for the chain portion of the molecule.  In
+principle, this could be varied to allow for modeling of longer or
+shorter chain lipid molecules. For these prolate ellipsoids, we have:
+\begin{equation}
+\begin{array}{rcl}
+d & < & l \\
+\epsilon^{r} & < & 1
+\end{array}
+\end{equation}
+A sketch of the various structural elements of our molecular-scale
+lipid / solvent model is shown in figure \ref{fig:lipidModel}.  The
+actual parameters used in our simulations are given in table
+\ref{tab:parameters}. 
 
+\begin{figure}[htb]
+\centering
+\includegraphics[width=4in]{2lipidModel} 
+\caption{The parameters defining the behavior of the lipid
+models. $l / d$ is the ratio of the head group to body diameter.
+Molecular bodies had a fixed aspect ratio of 3.0.  The solvent model
+was a simplified 4-water bead ($\sigma_w \approx d$) that has been
+used in other coarse-grained (DPD) simulations.  The dipolar strength
+(and the temperature and pressure) were the only other parameters that
+were varied systematically.\label{fig:lipidModel}}
+\end{figure}
+
+To take into account the permanent dipolar interactions of the
+zwitterionic head groups, we 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{equation}
+V_{ij}(r_{ij}) = 4\epsilon_h \left[\left(\frac{\sigma_h}{r_{ij}}\right)^{12} -
+\left(\frac{\sigma_h}{r_{ij}}\right)^6\right],
+\end{equation}
+and dipole-dipole,
+\begin{equation}
+V_{ij}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
+\hat{r}}_{ij})) = \frac{|\mu|^2}{4 \pi \epsilon_0 r_{ij}^3}
+\left[ \hat{\bf u}_i \cdot \hat{\bf u}_j - 3 (\hat{\bf u}_i \cdot
+\hat{\bf r}_{ij})(\hat{\bf u}_j \cdot \hat{\bf r}_{ij}) \right]
+\end{equation}
+potentials.  
+In these potentials, $\mathbf{\hat u}_i$ is the unit vector pointing
+along dipole $i$ and $\mathbf{\hat r}_{ij}$ is the unit vector
+pointing along the inter-dipole vector $\mathbf{r}_{ij}$.  
+
+For the interaction between nonequivalent uniaxial ellipsoids (in this
+case, between spheres and ellipsoids), the spheres are treated as
+ellipsoids with an aspect ratio of 1 ($d = l$) and with an well depth
+ratio ($\epsilon^r$) of 1 ($\epsilon^e = \epsilon^s$).  The form of
+the Gay-Berne potential we are using was generalized by Cleaver {\it
+et al.} and is appropriate for dissimilar uniaxial
+ellipsoids.\cite{Cleaver96}
+
+The solvent model in our simulations is 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) & & $= d$ & 13.8 & 4.7 \\
+$\epsilon^s$ (kcal/mol) & & 0.185 & 1.0 & 0.8 \\
+$\epsilon_r$ (well-depth aspect ratio)& & 1 & 0.2 &  1 \\
+$m$ (amu) & & 196 & 760 & 72.06 \\
+$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$) & & & & \\
+\multicolumn{2}c{$I_{xx}$} & 1125 & 45000 & N/A \\
+\multicolumn{2}c{$I_{yy}$} & 1125 & 45000 & N/A \\
+\multicolumn{2}c{$I_{zz}$} &  0 &    9000 & N/A \\
+$\mu$ (Debye) & & varied & 0 & 0 \\
+\end{tabular}
+\label{tab:parameters}
+\end{center}
+\end{minipage}
+\end{table*}
+
+A switching function has been applied to all potentials to smoothly
+turn off the interactions between a range of $22$ and $25$ \AA.
+
+The parameters that were systematically varied in this study were the
+size of the head group ($\sigma_h$), the strength of the dipole moment
+($\mu$), and the temperature of the system.  Values for $\sigma_h$
+ranged from 5.5 \AA\  to 6.5 \AA\ .  If the width of the tails is
+taken to be the unit of length, these head groups correspond to a
+range from $1.2 d$ to $1.41 d$.  Since the solvent beads are nearly
+identical in diameter to the tail ellipsoids, all distances that
+follow will be measured relative to this unit of distance. 
+
+\section{Experimental Methodology}
+\label{sec:experiment}
+
+To create unbiased bilayers, all simulations were started from two
+perfectly flat monolayers separated by a 26 \AA\ gap between the
+molecular bodies of the upper and lower leaves.  The separated
+monolayers were evolved in a 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 boundary conditions were used, and $480-720$
+lipid molecules were present in the simulations, depending on the size
+of the head beads.  In all cases, the two monolayers spontaneously
+collapsed into bilayer structures within 100 ps. Following this
+collapse, all systems were equlibrated for $100$ ns at $300$ K.
+
+The resulting bilayer structures were then solvated at a ratio of $6$
+solvent beads (24 water molecules) per lipid. These configurations
+were then equilibrated for another $30$ ns. All simulations utilizing
+the solvent were carried out at constant pressure ($P=1$ atm) with
+$3$D anisotropic coupling, and 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 with excellent energy
+conservation.  Data collection for structural properties of the
+bilayers was carried out during a final 5 ns run following the solvent
+equilibration.  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 with increasing size of the head groups. The surface
+is nearly flat when $\sigma_h=1.20 d$. With $\sigma_h=1.28 d$,
+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\
+wavelengths 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
+d$, 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 d$, 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, this ``symmetric'' ripple has no prefered
+direction.  The corresponding structures are shown in Figure
+\ref{fig:phaseCartoon} for elucidation of the detailed structures of
+different phases.  The top panel in figure \ref{fig:phaseCartoon} is
+the flat phase, the middle panel shows the asymmetric ripple phase
+corresponding to $\sigma_h = 1.41 d$ and the lower panel shows the
+symmetric ripple phase observed when $\sigma_h=1.35 d$. In the
+symmetric ripple, the bilayer is continuous over the whole membrane,
+however, in asymmetric ripple phase, the bilayer domains are connected
+by thin interdigitated monolayers that share molecules between the
+upper and lower leaves.
+\begin{table*}
+\begin{minipage}{\linewidth}
+\begin{center}
+\caption{Phases, ripple wavelengths and amplitudes observed as a
+function of the ratio between the head beads and the diameters of the
+tails.  All lengths are normalized to the diameter of the tail
+ellipsoids.}
+\begin{tabular}{lccc}
+\hline 
+$\sigma_h / d$ & type of phase & $\lambda / d$ & $A / d$\\
+\hline 
+1.20 & flat & N/A & N/A \\
+1.28 & asymmetric ripple or 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 / d$,
+reduced amplitude $A / d$ of the ripples are summarized in Table
+\ref{tab:property}. The wavelength range is $15~21$ molecular bodies
+and the amplitude is $1.5$ molecular bodies for asymmetric ripple and
+$2.2$ for symmetric ripple. These values are consistent to the
+experimental results.  Note, that given the lack of structural freedom
+in the tails of our model lipids, the amplitudes observed from these
+simulations are likely to underestimate of the true amplitudes.
+
+\begin{figure}[htb]
+\centering
+\includegraphics[width=4in]{topDown} 
+\caption{Top views of the flat (upper), 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 principal method for observing orientational ordering in dipolar
+or liquid crystalline systems is the $P_2$ order parameter (defined
+as $1.5 \times \lambda_{max}$, where $\lambda_{max}$ is the largest
+eigenvalue of the matrix,
+\begin{equation}
+{\mathsf{S}} = \frac{1}{N} \sum_i \left(
+\begin{array}{ccc}
+        u^{x}_i u^{x}_i-\frac{1}{3} & u^{x}_i u^{y}_i & u^{x}_i u^{z}_i \\
+        u^{y}_i u^{x}_i  & u^{y}_i u^{y}_i -\frac{1}{3} & u^{y}_i u^{z}_i \\
+        u^{z}_i u^{x}_i & u^{z}_i u^{y}_i  & u^{z}_i u^{z}_i -\frac{1}{3} 
+\end{array} \right).
+\label{eq:opmatrix}
+\end{equation}
+Here $u^{\alpha}_i$ is the $\alpha=x,y,z$ component of the unit vector
+for molecule $i$.  (Here, $\hat{\bf u}_i$ can refer either to the
+principal axis of the molecular body or to the dipole on the head
+group of the molecule.)  $P_2$ will be $1.0$ for a perfectly-ordered
+system and near $0$ for a randomized system.  Note that this order
+parameter is {\em not} equal to the polarization of the system.  For
+example, the polarization of a perfect anti-ferroelectric arrangement
+of point dipoles is $0$, but $P_2$ for the same system is $1$.  The
+eigenvector of $\mathsf{S}$ corresponding to the largest eigenvalue is
+familiar as the director axis, which can be used to determine a
+privileged axis for an orientationally-ordered system.  Since the
+molecular bodies are perpendicular to the head group dipoles, it is
+possible for the director axes for the molecular bodies and the head
+groups to be completely decoupled from each other.
+
+Figure \ref{fig:topView} shows snapshots of bird's-eye views of the
+flat ($\sigma_h = 1.20 d$) and rippled ($\sigma_h = 1.41, 1.35 d$)
+bilayers.  The directions of the dipoles on the head groups are
+represented with two colored half spheres: blue (phosphate) and yellow
+(amino).  For flat bilayers, the system exhibits signs of
+orientational frustration; some disorder in the dipolar head-to-tail
+chains is evident with kinks visible at the edges between differently
+ordered domains.  The lipids can also move independently of lipids in
+the opposing leaf, so the ordering of the dipoles on one leaf is not
+necessarily consistent with the ordering on the other.  These two
+factors keep the total dipolar order parameter relatively low for the
+flat phases.
+
+With increasing head group size, the surface becomes corrugated, and
+the dipoles cannot move as freely on the surface. Therefore, the
+translational freedom of lipids in one layer is dependent upon the
+position of the lipids in the other layer.  As a result, the ordering of
+the dipoles on head groups in one leaf is correlated with the ordering
+in the other leaf.  Furthermore, as the membrane deforms due to the
+corrugation, the symmetry of the allowed dipolar ordering on each leaf
+is broken. The dipoles then self-assemble in a head-to-tail
+configuration, and the dipolar order parameter increases dramatically.
+However, the total polarization of the system is still close to zero.
+This is strong evidence that the corrugated structure is an
+antiferroelectric state.  It is also notable that the head-to-tail
+arrangement of the dipoles is always observed in a direction
+perpendicular to the wave vector for the surface corrugation.  This is
+a similar finding to what we observed in our earlier work on the
+elastic dipolar membranes.\cite{Sun2007}
+
+The $P_2$ order parameters (for both the molecular bodies and the head
+group dipoles) have been calculated to quantify the ordering in these
+phases.  Figure \ref{fig:rP2} shows that the $P_2$ order parameter for
+the head-group dipoles increases with increasing head group size. When
+the heads of the lipid molecules are small, the membrane is nearly
+flat. Since the in-plane packing is essentially a close packing of the
+head groups, the head dipoles exhibit frustration in their 
+orientational ordering.
+
+The ordering trends for the tails are essentially opposite to the
+ordering of the head group dipoles. The tail $P_2$ order parameter
+{\it decreases} with increasing head size. This indicates that the
+surface is more curved with larger head / tail size ratios. When the
+surface is flat, all tails are pointing in the same direction (normal
+to the bilayer surface).  This simplified model appears to be
+exhibiting a smectic A fluid phase, similar to the real $L_{\beta}$
+phase.  We have not observed a smectic C gel phase ($L_{\beta'}$) for
+this model system.  Increasing the size of the heads results in
+rapidly decreasing $P_2$ ordering for the molecular bodies.
+
+\begin{figure}[htb]
+\centering
+\includegraphics[width=\linewidth]{rP2} 
+\caption{The $P_2$ order parameters for head groups (circles) and
+molecular bodies (squares) as a function of the ratio of head group
+size ($\sigma_h$) to the width of the molecular bodies ($d$). \label{fig:rP2}}
+\end{figure}
+
+In addition to varying the size of the head groups, we studied the
+effects of the interactions between head groups on the structure of
+lipid bilayer by changing the strength of the dipoles.  Figure
+\ref{fig:sP2} shows how the $P_2$ order parameter changes with
+increasing strength of the dipole.  Generally, the dipoles on the head
+groups become more ordered as the strength of the interaction between
+heads is increased and become more disordered by decreasing the
+interaction stength.  When the interaction between the heads becomes
+too weak, the bilayer structure does not persist; all lipid molecules
+become dispersed in the solvent (which is non-polar in this
+molecular-scale model).  The critial value of the strength of the
+dipole depends on the size of the head groups.  The perfectly flat
+surface becomes unstable below $5$ Debye, while the  rippled
+surfaces melt at $8$ Debye (asymmetric) or $10$ Debye (symmetric). 
+
+The ordering of the tails mirrors the ordering of the dipoles {\it
+except for the flat phase}. Since the surface is nearly flat in this
+phase, the order parameters are only weakly dependent on dipolar
+strength until it reaches $15$ Debye.  Once it reaches this value, the
+head group interactions are strong enough to pull the head groups
+close to each other and distort the bilayer structure. For a flat
+surface, a substantial amount of free volume between the head groups
+is normally available.  When the head groups are brought closer by
+dipolar interactions, the tails are forced to splay outward, forming
+first curved bilayers, and then inverted micelles. 
+
+When $\sigma_h=1.28 d$, the $P_2$ order parameter decreases slightly
+when the strength of the dipole is increased above $16$ debye. For
+rippled bilayers, there is less free volume available between the head
+groups. Therefore increasing dipolar strength weakly influences the
+structure of the membrane.  However, the increase in the body $P_2$
+order parameters implies that the membranes are being slightly
+flattened due to the effects of increasing head-group attraction.
+
+A very interesting behavior takes place when the head groups are very
+large relative to the molecular bodies ($\sigma_h = 1.41 d$) and the
+dipolar strength is relatively weak ($\mu < 9$ Debye). In this case,
+the two leaves of the bilayer become totally interdigitated with each
+other in large patches of the membrane.   With higher dipolar
+strength, the interdigitation is limited to single lines that run
+through the bilayer in a direction perpendicular to the ripple wave
+vector. 
+
+\begin{figure}[htb]
+\centering
+\includegraphics[width=\linewidth]{sP2} 
+\caption{The $P_2$ order parameters for head group dipoles (a) and
+molecular bodies (b) as a function of the strength of the dipoles.
+These order parameters are shown for four values of the head group /
+molecular width ratio ($\sigma_h / d$). \label{fig:sP2}}
+\end{figure}
+
+Figure \ref{fig:tP2} shows the dependence of the order parameters on
+temperature.  As expected, systems are more ordered at low
+temperatures, and more disordered at high temperatures.  All of the
+bilayers we studied can become unstable if the temperature becomes
+high enough.  The only interesting feature of the temperature
+dependence is in the flat surfaces ($\sigma_h=1.20 d$ and
+$\sigma_h=1.28 d$).  Here, when the temperature is increased above
+$310$K, there is enough jostling of the head groups to allow the
+dipolar frustration to resolve into more ordered states.  This results
+in a slight increase in the $P_2$ order parameter above this
+temperature.
+
+For the rippled surfaces ($\sigma_h=1.35 d$ and $\sigma_h=1.41 d$),
+there is a slightly increased orientational ordering in the molecular
+bodies above $290$K.  Since our model lacks the detailed information
+about the behavior of the lipid tails, this is the closest the model
+can come to depicting the ripple ($P_{\beta'}$) to fluid
+($L_{\alpha}$) phase transition.  What we are observing is a
+flattening of the rippled structures made possible by thermal
+expansion of the tightly-packed head groups.  The lack of detailed
+chain configurations also makes it impossible for this model to depict
+the ripple to gel ($L_{\beta'}$) phase transition.
+
+\begin{figure}[htb]
+\centering
+\includegraphics[width=\linewidth]{tP2} 
+\caption{The $P_2$ order parameters for head group dipoles (a) and
+molecular bodies (b) as a function of temperature.
+These order parameters are shown for four values of the head group /
+molecular width ratio ($\sigma_h / d$).\label{fig:tP2}}
+\end{figure}
+
+\section{Discussion}
+\label{sec:discussion}
+
+The ripple phases have been observed in our molecular dynamic
+simulations using a simple molecular lipid model. The lipid model
+consists of an anisotropic interacting dipolar head group and an
+ellipsoid shape tail. According to our simulations, the explanation of
+the formation for the ripples are originated in the size mismatch
+between the head groups and the tails. The ripple phases are only
+observed in the studies using larger head group lipid models. However,
+there is a mismatch betweent the size of the head groups and the size
+of the tails in the simulations of the flat surface. This indicates
+the competition between the anisotropic dipolar interaction and the
+packing of the tails also plays a major role for formation of the
+ripple phase. The larger head groups provide more free volume for the
+tails, while these hydrophobic ellipsoids trying to be close to each
+other, this gives the origin of the spontanous curvature of the
+surface, which is believed as the beginning of the ripple phases. The
+lager head groups cause the spontanous curvature inward for both of
+leaves of the bilayer. This results in a steric strain when the tails
+of two leaves too close to each other. The membrane has to be broken
+to release this strain. There are two ways to arrange these broken
+curvatures: symmetric and asymmetric ripples. Both of the ripple
+phases have been observed in our studies. The difference between these
+two ripples is that the bilayer is continuum in the symmetric ripple
+phase and is disrupt in the asymmetric ripple phase.
+
+Dipolar head groups are the key elements for the maintaining of the
+bilayer structure. The lipids are solvated in water when lowering the
+the strength of the dipole on the head groups. The long range
+orientational ordering of the dipoles can be achieved by forming the
+ripples, although the dipoles are likely to form head-to-tail
+configurations even in flat surface, the frustration prevents the
+formation of the long range orientational ordering for dipoles. The
+corrugation of the surface breaks the frustration and stablizes the
+long range oreintational ordering for the dipoles in the head groups
+of the lipid molecules. Many rows of the head-to-tail dipoles are
+parallel to each other and adopt the antiferroelectric state as a
+whole. This is the first time the organization of the head groups in
+ripple phases of the lipid bilayer has been addressed.
+
+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.
+
+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.
+
+\newpage
 \bibliography{mdripple}
 \end{document}