--- trunk/mdRipple/mdripple.tex	2007/07/13 22:01:52	3174
+++ trunk/mdRipple/mdripple.tex	2007/10/19 18:40:38	3265
@@ -1,405 +1,861 @@
 %\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{Dipolar Ordering of the Ripple Phase in Lipid Membranes}
-\author{Xiuquan Sun and J. Daniel Gezelter}
-\email[E-mail:]{gezelter@nd.edu}
-\affiliation{Department of Chemistry and Biochemistry,\\
-University of Notre Dame, \\
+\title{Dipolar ordering in the ripple phases of molecular-scale models
+of 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}
+Symmetric and asymmetric ripple phases have been observed to form in
+molecular dynamics simulations of a simple molecular-scale lipid
+model. The lipid model consists of an dipolar head group and an
+ellipsoidal tail.  Within the limits of this model, an explanation for
+generalized membrane curvature is a simple mismatch in the size of the
+heads with the width of the molecular bodies.  The persistence of a
+{\it bilayer} structure requires strong attractive forces between the
+head groups.  One feature of this model is that an energetically
+favorable orientational ordering of the dipoles can be achieved by
+out-of-plane membrane corrugation.  The corrugation of the surface
+stabilizes the long range orientational ordering for the dipoles in the
+head groups which then adopt a bulk anti-ferroelectric state. We
+observe a common feature of the corrugated dipolar membranes: the wave
+vectors for the surface ripples are always found to be perpendicular
+to the dipole director axis.  
 \end{abstract}
 
-\pacs{}
-\maketitle
+%\maketitle
+\newpage
 
 \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} The X-ray diffraction work by
+Katsaras {\it et al.} showed that a rich phase diagram exhibiting both
+{\it asymmetric} and {\it symmetric} ripples is possible for lecithin
+bilayers.\cite{Katsaras00}
 
-As one of the most important components in the formation of the
-biomembrane, lipid molecules attracted numerous studies in the past
-several decades. Due to their amphiphilic structure, when dispersed in
-water, lipids can self-assemble to construct a bilayer structure. The
-phase behavior of lipid membrane is well understood. The gel-fluid
-phase transition is known as main phase transition. However, there is
-an intermediate phase between gel and fluid phase for some lipid (like
-phosphatidycholine (PC)) membranes. This intermediate phase
-distinguish itself from other phases by its corrugated membrane
-surface, therefore this phase is termed ``ripple'' ($P_{\beta '}$)
-phase. The phase transition between gel-fluid and ripple phase is
-called pretransition. Since the pretransition usually occurs in room
-temperature, there might be some important biofuntions carried by the
-ripple phase for the living organism.
+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 head groups
+and molecular tilt with respect to the membrane normal vector can
+cause bulk rippling.~\cite{Bannerjee02}
 
-The ripple phase is observed experimentally by x-ray diffraction
-~\cite{Sun96,Katsaras00}, freeze-fracture electron microscopy
-(FFEM)~\cite{Copeland80,Meyer96}, and atomic force microscopy (AFM)
-recently~\cite{Kaasgaard03}. The experimental studies suggest two
-kinds of ripple structures: asymmetric (sawtooth like) and symmetric
-(sinusoidal like) ripple phases. Substantial number of theoretical
-explaination applied on the formation of the ripple
-phases~\cite{Marder84,Carlson87,Goldstein88,McCullough90,Lubensky93,Misbah98,Heimburg00,Kubica02,Bannerjee02}.
-In contrast, few molecular modelling have been done due to the large
-size of the resulting structures and the time required for the phases
-of interest to develop. One of the interesting molecular simulations
-was carried out by De Vries and Marrink {\it et
-al.}~\cite{deVries05}. According to their dynamic simulation results,
-the ripple consists of two domains, one is gel bilayer, and in the
-other domain, the upper and lower leaves of the bilayer are fully
-interdigitated. The mechanism of the formation of the ripple phase in
-their work suggests the theory that the packing competition between
-head group and tail of lipid molecules is the driving force for the
-formation of the ripple phase~\cite{Carlson87}. Recently, the ripple
-phase is also studied by using monte carlo simulation~\cite{Lenz07},
-the ripple structure is similar to the results of Marrink except that
-the connection of the upper and lower leaves of the bilayer is an
-interdigitated line instead of the fully interdigitated
-domain. Furthermore, the symmetric ripple phase was also observed in
-their work. They claimed 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.
+In contrast, few large-scale molecular modeling 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 time scales 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 organizations of the tails of lipid molecules are
-addressed by these molecular simulations, the ordering of the head
-group in ripple phase is still not settlement. We developed a simple
-``web of dipoles'' spin lattice model which provides some physical
-insight in our previous studies~\cite{Sun2007}, we found the dipoles
-on head groups of the lipid molecules are ordered in an
-antiferroelectric state. The similiar phenomenon is also observed by
-Tsonchev {\it et al.} when they studied the formation of the
-nanotube\cite{Tsonchev04}.
+Although the organization of the tails of lipid molecules are
+addressed by these molecular simulations and the packing competition
+between head groups and tails is strongly implicated as the primary
+driving force for ripple formation, questions about the ordering of
+the head groups in ripple phase have not been settled.
 
-In this paper, we made a more realistic coarse-grained lipid model to
-understand the primary driving force for membrane corrugation and to
-elucidate the organization of the anisotropic interacting head group
-via molecular dynamics simulation. We will talk about our model and
-methodology in section \ref{sec:method}, and details of the simulation
-in section \ref{sec:experiment}. The results are shown in section
-\ref{sec:results}. At last, we will discuss the results in section
+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 of dipolar
+elastic membranes is the anti-ferroelectric ordering of the dipoles.
+This was evident in the ordering of the dipole director axis
+perpendicular to the wave vector of the surface ripples.  A similar
+phenomenon has also been observed by Tsonchev {\it et al.} in their
+work on the spontaneous formation of dipolar peptide chains into
+curved nano-structures.\cite{Tsonchev04,Tsonchev04II}
+
+In this paper, we construct a somewhat more realistic molecular-scale
+lipid model than our previous ``web of dipoles'' and use molecular
+dynamics simulations to elucidate the role of the head group dipoles
+in the formation and morphology of the ripple phase.  We describe our
+model and computational methodology in section \ref{sec:method}.
+Details on the simulations are presented in section
+\ref{sec:experiment}, with results following in section
+\ref{sec:results}.  A final discussion of the role of dipolar heads in
+the ripple formation can be found in section
 \ref{sec:discussion}.
 
-\section{Methodology and Model}
+\section{Computational Model}
 \label{sec:method}
 
-Our idea for developing a simple and reasonable lipid model to study
-the ripple phase 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. The
-lipid model is shown in Figure \ref{fig:lipidMM}. Figure
-\ref{fig:lipidMM}a shows the molecular strucure of a DPPC molecule. The
-hydrophilic character of the head group is the effect of the strong
-dipole composed by a positive charge sitting on the nitrogen and a
-negative charge on the phosphate. The hydrophobic tail consists of
-fatty acid chains. In our model shown in Figure \ref{fig:lipidMM}b,
-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.
-
 \begin{figure}[htb]
 \centering
-\includegraphics[width=\linewidth]{lipidMM} 
-\caption{The molecular structure of a DPPC molecule and the
-coars-grained model for PC molecules.\label{fig:lipidMM}}
+\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}
 
-Spheres interact each other with Lennard-Jones potential
-\begin{eqnarray*}
-V_{ij} = 4\epsilon \left[\left(\frac{\sigma_0}{r_{ij}}\right)^{12} -
-\left(\frac{\sigma_0}{r_{ij}}\right)^6\right]
-\end{eqnarray*}
-here, $\sigma_0$ is diameter of the spherical particle. $r_{ij}$ is
-the distance between two spheres. $\epsilon$ is the well depth.
-Dipoles interact each other with typical dipole potential
-\begin{eqnarray*}
-V_{ij} = \frac{|\mu|^2}{4 \pi \epsilon_0 r_{ij}^3}
-\left[ \hat{\bf u}_i \cdot \hat{\bf u}_j - 3 (\hat{\bf u}_i \cdot
-\hat{\bf r}_{ij})(\hat{\bf u}_j \cdot \hat{\bf r}_{ij}) \right]
-\end{eqnarray*}
-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}$. Identical
-ellipsoids interact each other with Gay-Berne potential.
-\begin{eqnarray*}
+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 structure 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
+modeling 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]
-\end{eqnarray*} 
-where $\sigma_0 = \sqrt 2\sigma_s$, and orietation-dependent range
-parameter is given by
-\begin{eqnarray*}
-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}}) = 
-{\sigma_{0}} \left(1-\frac{1}{2}\chi\left[\frac{({\mathbf{\hat r}_{ij}}
-\cdot {\mathbf{\hat u}_i} + {\mathbf{\hat r}_{ij}} \cdot {\mathbf{\hat
-u}_j})^2}{1+\chi({\mathbf{\hat u}_i} \cdot {\mathbf{\hat u}_j})} +
-\frac{({\mathbf{\hat r}_{ij}} \cdot {\mathbf{\hat u}_i} - {\mathbf{\hat r}_{ij}}
-\cdot {\mathbf{\hat u}_j})^2}{1-\chi({\mathbf{\hat u}_i} \cdot
-{\mathbf{\hat u}_j})} \right] \right)^{-\frac{1}{2}}
-\end{eqnarray*} 
-and the strength anisotropy function is,
-\begin{eqnarray*}
-\epsilon ({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}}) = 
-{\epsilon^\nu}({\mathbf{\hat u}_i}, {\mathbf{\hat
-u}_j}){\epsilon'^\mu}({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
-{\mathbf{\hat r}_{ij}})
-\end{eqnarray*} 
-with $\nu$ and $\mu$ being adjustable exponent, and
-$\epsilon({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j})$,
-$\epsilon'({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat
-r}_{ij}})$ defined as
-\begin{eqnarray*} 
-\epsilon({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}) =
-\epsilon_0\left[1-\chi^2({\mathbf{\hat u}_i} \cdot {\mathbf{\hat
-u}_j})^2\right]^{-\frac{1}{2}} 
+\label{eq:gb}
+\end{equation*}
+
+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*}
-\begin{eqnarray*}
-\epsilon'({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}}) = 
-1-\frac{1}{2}\chi'\left[\frac{({\mathbf{\hat r}_{ij}} \cdot {\mathbf{\hat
-u}_i} + {\mathbf{\hat r}_{ij}} \cdot {\mathbf{\hat
-u}_j})^2}{1+\chi'({\mathbf{\hat u}_i} \cdot {\mathbf{\hat u}_j})} +
-\frac{({\mathbf{\hat r}_{ij}} \cdot {\mathbf{\hat u}_i} - {\mathbf{\hat r}_{ij}}
-\cdot {\mathbf{\hat u}_j})^2}{1-\chi'({\mathbf{\hat u}_i} \cdot
-{\mathbf{\hat u}_j})} \right]
+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.  Additionally, 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 diameter dependent parameter $\chi$ is given by
-\begin{eqnarray*} 
-\chi = \frac{({\sigma_s}^2 -
-{\sigma_e}^2)}{({\sigma_s}^2 + {\sigma_e}^2)} 
-\end{eqnarray*}
-and the well depth dependent parameter $\chi'$ is given by
-\begin{eqnarray*} 
-\chi' = \frac{({\epsilon_s}^{\frac{1}{\mu}} -
-{\epsilon_e}^{\frac{1}{\mu}})}{({\epsilon_s}^{\frac{1}{\mu}} +
-{\epsilon_e}^{\frac{1}{\mu}})} 
-\end{eqnarray*}
-$\sigma_s$ is the side-by-side width, and $\sigma_e$ is the end-to-end
-length. $\epsilon_s$ is the side-by-side well depth, and $\epsilon_e$
-is the end-to-end well depth. For the interaction between
-nonequivalent uniaxial ellipsoids (in this case, between spheres and
-ellipsoids), the range parameter is generalized as\cite{Cleaver96}
-\begin{eqnarray*}
-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}}) = 
-{\sigma_{0}} \left(1-\frac{1}{2}\chi\left[\frac{(\alpha{\mathbf{\hat r}_{ij}}
-\cdot {\mathbf{\hat u}_i} + \alpha^{-1}{\mathbf{\hat r}_{ij}} \cdot {\mathbf{\hat
-u}_j})^2}{1+\chi({\mathbf{\hat u}_i} \cdot {\mathbf{\hat u}_j})} +
-\frac{(\alpha{\mathbf{\hat r}_{ij}} \cdot {\mathbf{\hat u}_i} - \alpha^{-1}{\mathbf{\hat r}_{ij}}
-\cdot {\mathbf{\hat u}_j})^2}{1-\chi({\mathbf{\hat u}_i} \cdot
-{\mathbf{\hat u}_j})} \right] \right)^{-\frac{1}{2}}
-\end{eqnarray*}
-where $\alpha$ is given by
-\begin{eqnarray*}
-\alpha^2 =
-\left[\frac{\left({\sigma_e}_i^2-{\sigma_s}_i^2\right)\left({\sigma_e}_j^2+{\sigma_s}_i^2\right)}{\left({\sigma_e}_j^2-{\sigma_s}_j^2\right)\left({\sigma_e}_i^2+{\sigma_s}_j^2\right)}
-\right]^{\frac{1}{2}}
-\end{eqnarray*}
-the strength parameter is adjusted by the suggestion of
-\cite{Cleaver96}. All potentials are truncated at $25$ \AA\ and
-shifted at $22$ \AA.
+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}
 
-\section{Experiment}
+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 have placed 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 ($\sigma_h$) 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*}
+
+\section{Experimental Methodology}
 \label{sec:experiment}
 
-To make the simulations less expensive and to observe long-time
-behavior of the lipid membranes, all simulations were started from two
-separate monolayers in the vaccum with $x-y$ anisotropic pressure
+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.  Because the solvent we
+are using is non-polar and has a dielectric constant of 1, values for
+$\mu$ are sampled from a range that is somewhat smaller than the 20.6
+Debye dipole moment of the PC head groups.
+
+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 vacuum with $x-y$ anisotropic pressure
 coupling. The length of $z$ axis of the simulations was fixed and a
 constant surface tension was applied to enable real fluctuations of
-the bilayer. Periodic boundaries were used. There were $480-720$ lipid
-molecules in the simulations depending on the size of the head
-beads. All the simulations were equlibrated for $100$ ns at $300$
-K. The resulting structures were solvated in water ($6$ DPD
-water/lipid molecule). These configurations were relaxed for another
-$30$ ns relaxation. All simulations with water were carried out at
-constant pressure ($P=1$ atm) by $3$D anisotropic coupling, and
-constant surface tension ($\gamma=0.015$). Given the absence of fast
-degrees of freedom in this model, a timestep of $50$ fs was
-utilized. Simulations were performed by using OOPSE
-package\cite{Meineke05}.
-
-\section{Results and Analysis}
+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 equilibrated 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$ N/m). Given the absence of fast degrees of freedom in
+this model, a time step 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}
+
+A switching function was applied to all potentials to smoothly turn
+off the interactions between a range of $22$ and $25$ \AA.
+
+\section{Results}
 \label{sec:results}
 
-Snapshots in Figure \ref{fig:phaseCartoon} show that the membrane is
-more corrugated increasing size of the head groups. The surface is
-nearly flat when $\sigma_h=1.20\sigma_0$. With
-$\sigma_h=1.28\sigma_0$, although the surface is still flat, the
-bilayer starts to splay inward; the upper leaf of the bilayer is
-connected to the lower leaf with an interdigitated line defect. Two
-periodicities with $100$ \AA\ width were observed in the
-simulation. This structure is very similiar to the structure observed
-by de Vries and Lenz {\it et al.}. The same basic structure is also
-observed when $\sigma_h=1.41\sigma_0$, but the wavelength of the
-surface corrugations depends sensitively on the size of the ``head''
-beads. From the undulation spectrum, the corrugation is clearly
-non-thermal.
+The membranes in our simulations exhibit a number of interesting
+bilayer phases.  The surface topology of these phases depends most
+sensitively on the ratio of the size of the head groups to the width
+of the molecular bodies.  With heads only slightly larger than the
+bodies ($\sigma_h=1.20 d$) the membrane exhibits a flat bilayer. 
+
+Increasing the head / body size ratio increases the local membrane
+curvature around each of the lipids.  With $\sigma_h=1.28 d$, the
+surface is still essentially flat, but the bilayer starts to exhibit
+signs of instability.  We have observed occasional defects where a
+line of lipid molecules on one leaf of the bilayer will dip down to
+interdigitate with the other leaf.  This gives each of the two bilayer
+leaves some local convexity near the line defect.  These structures,
+once developed in a simulation, are very stable and are spaced
+approximately 100 \AA\ away from each other.
+
+With larger heads ($\sigma_h = 1.35 d$) the membrane curvature
+resolves into a ``symmetric'' ripple phase.  Each leaf of the bilayer
+is broken into several convex, hemicylinderical sections, and opposite
+leaves are fitted together much like roof tiles.  There is no
+interdigitation between the upper and lower leaves of the bilayer.
+
+For the largest head / body ratios studied ($\sigma_h = 1.41 d$) the
+local curvature is substantially larger, and the resulting bilayer
+structure resolves into an asymmetric ripple phase.  This structure is
+very similar to the structures observed by both de~Vries {\it et al.}
+and Lenz {\it et al.}.  For a given ripple wave vector, there are two
+possible asymmetric ripples, which is not the case for the symmetric
+phase observed when $\sigma_h = 1.35 d$.
+
 \begin{figure}[htb]
 \centering
-\includegraphics[width=\linewidth]{phaseCartoon} 
-\caption{A sketch to discribe the structure of the phases observed in
-our simulations.\label{fig:phaseCartoon}}
+\includegraphics[width=4in]{phaseCartoon} 
+\caption{The role of the ratio between the head group size and the
+width of the molecular bodies is to increase the local membrane
+curvature.  With strong attractive interactions between the head
+groups, this local curvature can be maintained in bilayer structures
+through surface corrugation.  Shown above are three phases observed in
+these simulations.  With $\sigma_h = 1.20 d$, the bilayer maintains a
+flat topology.  For larger heads ($\sigma_h = 1.35 d$) the local
+curvature resolves into a symmetrically rippled phase with little or
+no interdigitation between the upper and lower leaves of the membrane.
+The largest heads studied ($\sigma_h = 1.41 d$) resolve into an
+asymmetric rippled phases with interdigitation between the two
+leaves.\label{fig:phaseCartoon}} 
 \end{figure}
 
-When $\sigma_h=1.35\sigma_0$, we observed another corrugated surface
-morphology. This structure is different from the asymmetric rippled
-surface; there is no interdigitation between the upper and lower
-leaves of the bilayer. Each leaf of the bilayer is broken into several
-hemicylinderical sections, and opposite leaves are fitted together
-much like roof tiles. Unlike the surface in which the upper
-hemicylinder is always interdigitated on the leading or trailing edge
-of lower hemicylinder, the symmetric ripple has no prefered direction.
-The corresponding cartoons are shown in Figure
-\ref{fig:phaseCartoon} for elucidation of the detailed structures of
-different phases. Figure \ref{fig:phaseCartoon} (a) is the flat phase,
-(b) is the asymmetric ripple phase corresponding to the lipid
-organization when $\sigma_h=1.28\sigma_0$ and $\sigma_h=1.41\sigma_0$,
-and (c) is the symmetric ripple phase observed when
-$\sigma_h=1.35\sigma_0$. In symmetric ripple phase, the bilayer is
-continuous everywhere on the whole membrane, however, in asymmetric
-ripple phase, the bilayer is intermittent domains connected by thin
-interdigitated monolayer which consists of upper and lower leaves of
-the bilayer.
+Sample structures for the flat ($\sigma_h = 1.20 d$), symmetric
+($\sigma_h = 1.35 d$, and asymmetric ($\sigma_h = 1.41 d$) ripple
+phases are shown in Figure \ref{fig:phaseCartoon}.  
+
+It is reasonable to ask how well the parameters we used can produce
+bilayer properties that match experimentally known values for real
+lipid bilayers.  Using a value of $l = 13.8$ \AA for the ellipsoidal
+tails and the fixed ellipsoidal aspect ratio of 3, our values for the
+area per lipid ($A$) and inter-layer spacing ($D_{HH}$) depend
+entirely on the size of the head bead relative to the molecular body.
+These values are tabulated in table \ref{tab:property}.  Kucera {\it
+et al.}  have measured values for the head group spacings for a number
+of PC lipid bilayers that range from 30.8 \AA (DLPC) to 37.8 (DPPC).
+They have also measured values for the area per lipid that range from
+60.6
+\AA$^2$ (DMPC) to 64.2 \AA$^2$
+(DPPC).\cite{NorbertKucerka04012005,NorbertKucerka06012006} Only the
+largest of the head groups we modeled ($\sigma_h = 1.41 d$) produces
+bilayers (specifically the area per lipid) that resemble real PC
+bilayers.  The smaller head beads we used are perhaps better models
+for PE head groups.
+
 \begin{table*}
 \begin{minipage}{\linewidth}
 \begin{center}
-\caption{}
-\begin{tabular}{lccc}
+\caption{Phase, bilayer spacing, area per lipid, ripple wavelength
+and amplitude observed as a function of the ratio between the head
+beads and the diameters of the tails.  Ripple wavelengths and
+amplitudes are normalized to the diameter of the tail ellipsoids.}
+\begin{tabular}{lccccc}
 \hline 
-$\sigma_h/\sigma_0$ & type of phase & $\lambda/\sigma_0$ & $A/\sigma_0$\\
+$\sigma_h / d$ & type of phase & bilayer spacing (\AA) & area per
+lipid (\AA$^2$) & $\lambda / d$ & $A / d$\\
 \hline 
-1.20 & flat & N/A & N/A \\
-1.28 & asymmetric flat & 21.7 & N/A \\
-1.35 & symmetric ripple & 17.2 & 2.2 \\
-1.41 & asymmetric ripple & 15.4 & 1.5 \\
+1.20 & flat & 33.4 & 49.6 & N/A & N/A \\
+1.28 & flat & 33.7 & 54.7 & N/A & N/A \\
+1.35 & symmetric ripple & 42.9 & 51.7 & 17.2 & 2.2 \\
+1.41 & asymmetric ripple & 37.1 & 63.1 & 15.4 & 1.5 \\
 \end{tabular}
 \label{tab:property}
 \end{center}
 \end{minipage}
 \end{table*}
 
-The membrane structures and the reduced wavelength $\lambda/\sigma_0$,
-reduced amplitude $A/\sigma_0$ of the ripples are summerized in Table
-\ref{tab:property}. The wavelength range is $15~21$ and the amplitude
-is $1.5$ for asymmetric ripple and $2.2$ for symmetric ripple. These
-values are consistent to the experimental results. Note, the
-amplitudes are underestimated without the melted tails in our
-simulations.
+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 - 17$ molecular bodies
+and the amplitude is $1.5$ molecular bodies for asymmetric ripple and
+$2.2$ for symmetric ripple. These values are reasonably consistent
+with experimental measurements.\cite{Sun96,Katsaras00,Kaasgaard03}
+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.
 
-The $P_2$ order paramters (for molecular bodies and head group
-dipoles) have been calculated to clarify the ordering in these phases
-quantatively. $P_2=1$ means a perfectly ordered structure, and $P_2=0$
-implies orientational randomization. Figure \ref{fig:rP2} shows the
-$P_2$ order paramter of the dipoles on head group rising with
-increasing head group size. When the heads of the lipid molecules are
-small, the membrane is flat. The dipolar ordering is essentially
-frustrated on orientational ordering in this circumstance. Another
-reason is that the lipids can move independently in each monolayer, it
-is not nessasory for the direction of dipoles on one leaf is
-consistant to another layer, which makes total order parameter is
-relatively low. With increasing head group size, the surface is
-corrugated, and dipoles do not move as freely on the
-surface. Therefore, the translational freedom of lipids in one layer
-is dependent upon the position of lipids in another layer, as a
-result, the symmetry of the dipoles on head group in one layer is tied
-to the symmetry in the other layer. Furthermore, as the membrane
-deforms from two to three dimensions due to the corrugation, the
-symmetry of the ordering for the dipoles embedded on each leaf is
-broken. The dipoles then self-assemble in a head-tail configuration,
-and the order parameter increases dramaticaly. However, the total
-polarization of the system is still close to zero. This is strong
-evidence that the corrugated structure is an antiferroelectric
-state. The orientation of the dipolar is always perpendicular to the
-ripple wave vector. These results are consistent with our previous
-study on dipolar membranes.
+\begin{figure}[htb]
+\centering
+\includegraphics[width=4in]{topDown} 
+\caption{Top views of the flat (upper), symmetric ripple (middle),
+and asymmetric 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.  Note that the flat
+membrane has multiple vortex defects in the dipolar ordering, and the
+ordering on the lower leaf of the bilayer can be in an entirely
+different direction from the upper leaf.\label{fig:topView}}
+\end{figure}
 
-The ordering of the tails is essentially opposite to the ordering of
-the dipoles on head group. The $P_2$ order parameter decreases with
-increasing head size. This indicates the surface is more curved with
-larger head groups. When the surface is flat, all tails are pointing
-in the same direction; in this case, all tails are parallel to the
-normal of the surface,(making this structure remindcent of the
-$L_{\beta}$ phase. Increasing the size of the heads, results in
+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.35, 1.41 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
+anti-ferroelectric 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 parameter as a funtion of the ratio of
-$\sigma_h$ to $\sigma_0$.\label{fig: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}
 
-We studied the effects of the interactions between head groups on the
-structure of lipid bilayer by changing the strength of the dipole.
-Figure \ref{fig:sP2} shows how the $P_2$ order parameter changes with
-increasing strength of the dipole. Generally the dipoles on the head
-group are more ordered by increase in the strength of the interaction
-between heads and are more disordered by decreasing the interaction
-stength. When the interaction between the heads is weak enough, the
-bilayer structure does not persist; all lipid molecules are solvated
-directly in the water. The critial value of the strength of the dipole
-depends on the head size. The perfectly flat surface melts at $5$
-debye, the asymmetric rippled surfaces melt at $8$ debye, the
-symmetric rippled surfaces melt at $10$ debye. The ordering of the
-tails is the same as the ordering of the dipoles except for the flat
-phase. Since the surface is already perfect flat, the order parameter
-does not change much until the strength of the dipole is $15$
-debye. However, the order parameter decreases quickly when the
-strength of the dipole is further increased. The head groups of the
-lipid molecules are brought closer by stronger interactions between
-them. For a flat surface, a large amount of free volume between the
-head groups is available, but when the head groups are brought closer,
-the tails will splay outward, forming an inverse micelle. For rippled
-surfaces, there is less free volume available between the head
-groups. Therefore there is little effect on the structure of the
-membrane due to increasing dipolar strength. Unlike other systems that
-melt directly when the interaction is weak enough, for
-$\sigma_h=1.41\sigma_0$, part of the membrane melts into itself
-first. The upper leaf of the bilayer becomes totally interdigitated
-with the lower leaf. This is different behavior than what is exhibited
-with the interdigitated lines in the rippled phase where only one
-interdigitated line connects the two leaves of bilayer.
+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 strength.  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 critical 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, first forming
+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 parameter as a funtion of the strength of the
-dipole.\label{fig: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 parameter on
-temperature. The behavior of the $P_2$ order paramter is
-straightforward. Systems are more ordered at low temperature, and more
-disordered at high temperatures. When the temperature is high enough,
-the membranes are instable. Since our model lacks the detailed
-information on lipid tails, we can not simulate the fluid phase with
-melted fatty acid chains. Moreover, the formation of the tilted
-$L_{\beta'}$ phase also depends on the organization of fatty groups on
-tails.
+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 parameter as a funtion of
-temperature.\label{fig: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}
 
+Fig. \ref{fig:phaseDiagram} shows a phase diagram for the model as a
+function of the head group / molecular width ratio ($\sigma_h / d$)
+and the strength of the head group dipole moment ($\mu$).  Note that
+the specific form of the bilayer phase is governed almost entirely by
+the head group / molecular width ratio, while the strength of the
+dipolar interactions between the head groups governs the stability of
+the bilayer phase.  Weaker dipoles result in unstable bilayer phases,
+while extremely strong dipoles can shift the equilibrium to an
+inverted micelle phase when the head groups are small.   Temperature
+has little effect on the actual bilayer phase observed, although higher
+temperatures can cause the unstable region to grow into the higher
+dipole region of this diagram.
+
+\begin{figure}[htb]
+\centering
+\includegraphics[width=\linewidth]{phaseDiagram} 
+\caption{Phase diagram for the simple molecular model as a function
+of the head group / molecular width ratio ($\sigma_h / d$) and the
+strength of the head group dipole moment
+($\mu$).\label{fig:phaseDiagram}}
+\end{figure}
+
+We have computed translational diffusion coefficients for lipid
+molecules from the mean square displacement,
+\begin{eqnarray}
+\langle {|\left({\bf r}_{i}(t) - {\bt r}_{i}(0) \right)|}^2 \rangle \\ \\
+& = & 6Dt
+\end{eqnarray}
+of the lipid bodies. The values of the translational diffusion
+coefficient for different head-to-tail size ratio are shown in table
+\ref{tab:relaxation}.
+
+We have also computed orientational diffusion constants for the head
+groups from the relaxation of the second-order Legendre polynomial
+correlation function,
+\begin{eqnarray}
+C_{\ell}(t) & = & \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf
+\mu}_{i}(0) \right) \rangle  \\ \\
+& \approx & e^{-\ell(\ell + 1) \theta t},
+\end{eqnarray}
+of the head group dipoles.  In this last line, we have used a simple
+``Debye''-like model for the relaxation of the correlation function,
+specifically in the case when $\ell = 2$.   The computed orientational
+diffusion constants are given in table \ref{tab:relaxation}.  The
+notable feature we observe is that the orientational diffusion
+constant for the head group exhibits an order of magnitude decrease
+upon entering the rippled phase.  Our orientational correlation times
+are substantially in excess of those provided by...
+
+
+\begin{table*}
+\begin{minipage}{\linewidth}
+\begin{center}
+\caption{Rotational diffusion constants for the head groups
+($\theta_h$) and molecular bodies ($\theta_b$) as well as the
+translational diffusion coefficients for the molecule as a function of
+the head-to-body width ratio.  The orientational mobility of the head
+groups experiences an {\it order of magnitude decrease} upon entering
+the rippled phase, which suggests that the rippling is tied to a
+freezing out of head group orientational freedom.  Uncertainties in
+the last digit are indicated by the values in parentheses.}
+\begin{tabular}{lccc}
+\hline 
+$\sigma_h / d$ & $\theta_h (\mu s^{-1})$ & $\theta_b (1/fs)$ & $D (
+\times 10^{-11} m^2 s^{-1} \\
+\hline 
+1.20 & $0.206(1) $ & $0.0175(5) $ & $0.43(1)$ \\
+1.28 & $0.179(2) $ & $0.055(2)  $ & $5.91(3)$ \\
+1.35 & $0.025(1) $ & $0.195(3)  $ & $3.42(1)$ \\
+1.41 & $0.023(1) $ & $0.024(3)  $ & $7.16(1)$ \\
+\end{tabular}
+\label{tab:relaxation}
+\end{center}
+\end{minipage}
+\end{table*}
+
 \section{Discussion}
 \label{sec:discussion}
 
+Symmetric and asymmetric ripple phases have been observed to form in
+our molecular dynamics simulations of a simple molecular-scale lipid
+model. The lipid model consists of an dipolar head group and an
+ellipsoidal tail.  Within the limits of this model, an explanation for
+generalized membrane curvature is a simple mismatch in the size of the
+heads with the width of the molecular bodies.  With heads
+substantially larger than the bodies of the molecule, this curvature
+should be convex nearly everywhere, a requirement which could be
+resolved either with micellar or cylindrical phases.
+
+The persistence of a {\it bilayer} structure therefore requires either
+strong attractive forces between the head groups or exclusionary
+forces from the solvent phase.  To have a persistent bilayer structure
+with the added requirement of convex membrane curvature appears to
+result in corrugated structures like the ones pictured in
+Fig. \ref{fig:phaseCartoon}.  In each of the sections of these
+corrugated phases, the local curvature near a most of the head groups
+is convex.  These structures are held together by the extremely strong
+and directional interactions between the head groups.
+
+Dipolar head groups are key for the maintaining the bilayer structures
+exhibited by this model.  The dipoles are likely to form head-to-tail
+configurations even in flat configurations, but the temperatures are
+high enough that vortex defects become prevalent in the flat phase.
+The flat phase we observed therefore appears to be substantially above
+the Kosterlitz-Thouless transition temperature for a planar system of
+dipoles with this set of parameters.  For this reason, it would be
+interesting to observe the thermal behavior of the flat phase at
+substantially lower temperatures.
+
+One feature of this model is that an energetically favorable
+orientational ordering of the dipoles can be achieved by forming
+ripples.  The corrugation of the surface breaks the symmetry of the
+plane, making vortex defects somewhat more expensive, and stabilizing
+the long range orientational ordering for the dipoles in the head
+groups.  Most of the rows of the head-to-tail dipoles are parallel to
+each other and the system adopts a bulk anti-ferroelectric state.  We
+believe that this is the first time the organization of the head
+groups in ripple phases has been addressed.
+
+Although the size-mismatch between the heads and molecular bodies
+appears to be the primary driving force for surface convexity, the
+persistence of the bilayer through the use of rippled structures is a
+function of the strong, attractive interactions between the heads.
+One important prediction we can make using the results from this
+simple model is that if the dipole-dipole interaction is the leading
+contributor to the head group attractions, the wave vectors for the
+ripples should always be found {\it perpendicular} to the dipole
+director axis.  This echoes the prediction we made earlier for simple
+elastic dipolar membranes, and may suggest experimental designs which
+will test whether this is really the case 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.\cite{deVries05}
+
+Although our model is simple, it exhibits some rich and unexpected
+behaviors.  It would clearly be a closer approximation to reality if
+we allowed bending motions between the dipoles and the molecular
+bodies, and if we replaced the rigid ellipsoids with ball-and-chain
+tails.  However, the advantages of this simple model (large system
+sizes, 50 fs time steps) allow us to rapidly explore the phase diagram
+for a wide range of parameters.  Our explanation of this rippling
+phenomenon will help us design more accurate molecular models for
+corrugated membranes and experiments to test whether or not
+dipole-dipole interactions exert an influence on membrane rippling.
+\newpage
 \bibliography{mdripple}
 \end{document}