--- trunk/mdRipple/mdripple.tex	2007/07/19 22:24:21	3186
+++ trunk/mdRipple/mdripple.tex	2007/08/01 16:07:12	3202
@@ -1,304 +1,462 @@
 %\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 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
 
 \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}
 
-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 headgroups
+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 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 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 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.
 
-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 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}.
 
-\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]{lipidModels} 
+\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 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]
-\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.  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 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}
 
+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[height=4in]{lipidModel} 
+\includegraphics[width=4in]{2lipidModel} 
 \caption{The parameters defining the behavior of the lipid
-models. $\sigma_h / \sigma_0$ is the ratio of the head group to body
-diameter.  Molecular bodies all had an aspect ratio of 3.0.  The
-dipolar strength (and the temperature and pressure) wer the only other
-parameters that wer varied systematically.\label{fig:lipidModel}}
+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}
 
-\section{Experiment}
+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 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
+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 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}.
+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.
 
-\section{Results and Analysis}
+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 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.
+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=\linewidth]{phaseCartoon} 
+\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\sigma_0$, we observed another corrugated surface
-morphology. This structure is different from the asymmetric rippled
+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, the symmetric ripple has no prefered direction.
-The corresponding cartoons are shown in Figure
+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. 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.
+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{}
+\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/\sigma_0$ & type of phase & $\lambda/\sigma_0$ & $A/\sigma_0$\\
+$\sigma_h / d$ & type of phase & $\lambda / d$ & $A / d$\\
 \hline 
 1.20 & flat & N/A & N/A \\
-1.28 & asymmetric flat & 21.7 & 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}
@@ -307,130 +465,259 @@ The membrane structures and the reduced wavelength $\l
 \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~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.
 
-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), 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 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.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 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$
-$0.03$ debye, the asymmetric rippled surfaces melt at $8$ $0.04$
-$0.03$ debye, the symmetric rippled surfaces melt at $10$ $0.04$
-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. When $\sigma_h=1.28\sigma_0$, the $P_2$
-order parameter decreases slightly after the strength of the dipole is
-increased to $16$ debye. 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. However, the increase of the $P_2$ order parameter implies
-the membranes are flatten by the increase of the strength of the
-dipole. 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 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 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. For flat surfaces ($\sigma_h=1.20\sigma_0$
-and $\sigma_h=1.28\sigma_0$), when the temperature is increased to
-$310$, the $P_2$ order parameter increases slightly instead of
-decreases like ripple surface. This is an evidence of the frustration
-of the dipolar ordering in each leaf of the lipid bilayer, at low
-temperature, the systems are locked in a local minimum energy state,
-with increase of the temperature, the system can jump out the local
-energy well to find the lower energy state which is the longer range
-orientational ordering. Like the dipolar ordering of the flat
-surfaces, the ordering of the tails of the lipid molecules for ripple
-membranes ($\sigma_h=1.35\sigma_0$ and $\sigma_h=1.41\sigma_0$) also
-show some nonthermal characteristic. With increase of the temperature,
-the $P_2$ order parameter decreases firstly, and increases afterward
-when the temperature is greater than $290 K$. The increase of the
-$P_2$ order parameter indicates a more ordered structure for the tails
-of the lipid molecules which corresponds to a more flat surface. 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}
 
 \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}