--- trunk/xDissertation/md.tex	2008/01/30 16:01:02	3336
+++ trunk/xDissertation/md.tex	2008/03/05 23:34:04	3360
@@ -1 +1,897 @@
-\chapter{\label{chap:md}MOLECULAR DYNAMICS}
+\chapter{\label{chap:md}DIPOLAR ORDERING IN THE RIPPLE PHASES OF
+MOLECULAR-SCALE MODELS OF LIPID MEMBRANES}
+
+\section{Introduction}
+\label{mdsec: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} although Tenchov {\it et al.} have
+recently observed near-hexagonal packing in some phosphatidylcholine
+(PC) gel phases.\cite{Tenchov2001} 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}
+
+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} Recently, Kranenburg and Smit
+described the formation of symmetric ripple-like structures using a
+coarse grained solvent-head-tail bead model.\cite{Kranenburg2005}
+Their lipids consisted of a short chain of head beads tied to the two
+longer ``chains''. 
+
+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 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 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{sun:031602} 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{mdsec:method}.
+Details on the simulations are presented in section
+\ref{mdsec:experiment}, with results following in section
+\ref{mdsec:results}.  A final discussion of the role of dipolar heads in
+the ripple formation can be found in section
+\ref{mdsec:discussion}.
+
+\section{Computational Model}
+\label{mdsec:method}
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{./figures/mdLipidModels.pdf} 
+\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{mdfig:lipidModels}}
+\end{figure}
+
+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{mdfig: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{Gay1981} It can be thought of as a modification of the
+Gaussian overlap model originally described by Berne and
+Pechukas.\cite{Berne72} The potential is constructed in the familiar
+form of the Lennard-Jones function using orientation-dependent
+$\sigma$ and $\epsilon$ parameters,
+\begin{equation*}
+V_{ij}({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat
+r}_{ij}}) = 4\epsilon ({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
+{\mathbf{\hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i},
+{\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^{12}
+-\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
+{\mathbf{\hat r}_{ij}})+\sigma_0}\right)^6\right]
+\label{mdeq: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*}
+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 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{mdfig:lipidModel}.  The
+actual parameters used in our simulations are given in table
+\ref{mdtab:parameters}. 
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{./figures/md2LipidModel.pdf} 
+\caption{The parameters defining the behavior of the lipid
+models. $\sigma_h / 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 simulations.  The dipolar strength
+(and the temperature and pressure) were the only other parameters that
+were varied systematically.\label{mdfig: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}$.  
+
+Since the charge separation distance is so large in zwitterionic head
+groups (like the PC head groups), it would also be possible to use
+either point charges or a ``split dipole'' approximation,
+\begin{equation}
+V_{ij}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
+\hat{r}}_{ij})) = \frac{1}{{4\pi \epsilon_0 }}\left[ {\frac{{\mu _i  \cdot \mu _j }}{{R_{ij}^3 }} -
+\frac{{3\left( {\mu _i  \cdot r_{ij} } \right)\left( {\mu _i  \cdot
+r_{ij} } \right)}}{{R_{ij}^5 }}} \right]
+\end{equation}
+where $\mu _i$ and $\mu _j$ are the dipole moments of sites $i$ and
+$j$, $r_{ij}$ is vector between the two sites, and $R_{ij}$ is given
+by,
+\begin{equation}
+R_{ij}  = \sqrt {r_{ij}^2  + \frac{{d_i^2 }}{4} + \frac{{d_j^2
+}}{4}}.
+\end{equation}
+Here, $d_i$ and $d_j$ are charge separation distances associated with
+each of the two dipolar sites. This approximation to the multipole
+expansion maintains the fast fall-off of the multipole potentials but
+lacks the normal divergences when two polar groups get close to one
+another.
+
+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 similar to the one used by
+Marrink {\it et al.}  in their coarse grained simulations of lipid
+bilayers.\cite{Marrink2004} The 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.  Note that although we
+are using larger cutoff and switching radii than Marrink {\it et al.},
+our solvent density at 300 K remains at 0.944 g cm$^{-3}$, and the
+solvent diffuses at 0.43 \AA$^2 ps^{-1}$ (only twice as fast as liquid
+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{mdtab:parameters}
+\end{center}
+\end{minipage}
+\end{table*}
+
+\section{Experimental Methodology}
+\label{mdsec:experiment}
+
+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 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 small 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.  Orientational correlation functions and diffusion
+constants were computed from 30 ns simulations in the microcanonical
+(NVE) ensemble using the average volume from the end of the constant
+pressure and surface tension runs.  The timestep on these final
+molecular dynamics runs was 25 fs.  No appreciable changes in phase
+structure were noticed upon switching to a microcanonical ensemble.
+All simulations were performed using the {\sc oopse} molecular
+modeling program.\cite{Meineke2005}
+
+A switching function was applied to all potentials to smoothly turn
+off the interactions between a range of $22$ and $25$ \AA.  The
+switching function was the standard (cubic) function, 
+\begin{equation}
+s(r) = 
+        \begin{cases}
+        1 & \text{if $r \le r_{\text{sw}}$},\\
+        \frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2}
+        {(r_{\text{cut}} - r_{\text{sw}})^3} 
+        & \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\
+        0 & \text{if $r > r_{\text{cut}}$.}
+        \end{cases}
+\label{mdeq:dipoleSwitching}
+\end{equation}
+
+\section{Results}
+\label{mdsec:results}
+
+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}
+\centering
+\includegraphics[width=\linewidth]{./figures/mdPhaseCartoon.pdf} 
+\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{mdfig:phaseCartoon}} 
+\end{figure}
+
+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{mdfig: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{mdtab: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 \AA\ (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{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 / d$ & type of phase & bilayer spacing (\AA) & area per
+lipid (\AA$^2$) & $\lambda / d$ & $A / d$\\
+\hline 
+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{mdtab:property}
+\end{center}
+\end{minipage}
+\end{table*}
+
+The membrane structures and the reduced wavelength $\lambda / d$,
+reduced amplitude $A / d$ of the ripples are summarized in Table
+\ref{mdtab: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.
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{./figures/mdTopDown.pdf} 
+\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{mdfig:topView}}
+\end{figure}
+
+The principal method for observing orientational ordering in dipolar
+or liquid crystalline systems is the $P_2$ order parameter (defined
+as $1.5 \times \lambda_{max}$, where $\lambda_{max}$ is the largest
+eigenvalue of the matrix,
+\begin{equation}
+{\mathsf{S}} = \frac{1}{N} \sum_i \left(
+\begin{array}{ccc}
+        u^{x}_i u^{x}_i-\frac{1}{3} & u^{x}_i u^{y}_i & u^{x}_i u^{z}_i \\
+        u^{y}_i u^{x}_i  & u^{y}_i u^{y}_i -\frac{1}{3} & u^{y}_i u^{z}_i \\
+        u^{z}_i u^{x}_i & u^{z}_i u^{y}_i  & u^{z}_i u^{z}_i -\frac{1}{3} 
+\end{array} \right).
+\label{mdeq: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{mdfig: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{sun:031602}
+
+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{mdfig: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}
+\centering
+\includegraphics[width=\linewidth]{./figures/mdRP2.pdf} 
+\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{mdfig:rP2}}
+\end{figure}
+
+In addition to varying the size of the head groups, we studied the
+effects of the interactions between head groups on the structure of
+lipid bilayer by changing the strength of the dipoles.  Figure
+\ref{mdfig: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}
+\centering
+\includegraphics[width=\linewidth]{./figures/mdSP2.pdf} 
+\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{mdfig:sP2}}
+\end{figure}
+
+Figure \ref{mdfig: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}
+\centering
+\includegraphics[width=\linewidth]{./figures/mdTP2.pdf} 
+\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{mdfig:tP2}}
+\end{figure}
+
+Fig. \ref{mdfig: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}
+\centering
+\includegraphics[width=\linewidth]{./figures/mdPhaseDiagram.pdf} 
+\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{mdfig:phaseDiagram}}
+\end{figure}
+
+We have computed translational diffusion constants for lipid molecules
+from the mean-square displacement,
+\begin{equation}
+D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \rangle,
+\end{equation}
+of the lipid bodies. Translational diffusion constants for the
+different head-to-tail size ratios (all at 300 K) are shown in table
+\ref{mdtab:relaxation}.  We have also computed orientational correlation
+times for the head groups from fits of the second-order Legendre
+polynomial correlation function,
+\begin{equation}
+C_{\ell}(t)  =  \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf
+\mu}_{i}(0) \right) \rangle
+\end{equation}
+of the head group dipoles.  The orientational correlation functions
+appear to have multiple components in their decay: a fast ($12 \pm 2$
+ps) decay due to librational motion of the head groups, as well as
+moderate ($\tau^h_{\rm mid}$) and slow ($\tau^h_{\rm slow}$)
+components.  The fit values for the moderate and slow correlation
+times are listed in table \ref{mdtab:relaxation}.  Standard deviations
+in the fit time constants are quite large (on the order of the values
+themselves).
+
+Sparrman and Westlund used a multi-relaxation model for NMR lineshapes
+observed in gel, fluid, and ripple phases of DPPC and obtained
+estimates of a correlation time for water translational diffusion
+($\tau_c$) of 20 ns.\cite{Sparrman2003} This correlation time
+corresponds to water bound to small regions of the lipid membrane.
+They further assume that the lipids can explore only a single period
+of the ripple (essentially moving in a nearly one-dimensional path to
+do so), and the correlation time can therefore be used to estimate a
+value for the translational diffusion constant of $2.25 \times
+10^{-11} m^2 s^{-1}$.  The translational diffusion constants we obtain
+are in reasonable agreement with this experimentally determined
+value. However, the $T_2$ relaxation times obtained by Sparrman and
+Westlund are consistent with P-N vector reorientation timescales of
+2-5 ms.  This is substantially slower than even the slowest component
+we observe in the decay of our orientational correlation functions.
+Other than the dipole-dipole interactions, our head groups have no
+shape anisotropy which would force them to move as a unit with
+neighboring molecules.  This would naturally lead to P-N reorientation
+times that are too fast when compared with experimental measurements.
+
+\begin{table*}
+\begin{minipage}{\linewidth}
+\begin{center}
+\caption{Fit values for the rotational correlation times for the head
+groups ($\tau^h$) and molecular bodies ($\tau^b$) as well as the
+translational diffusion constants for the molecule as a function of
+the head-to-body width ratio.  All correlation functions and transport
+coefficients were computed from microcanonical simulations with an
+average temperture of 300 K.  In all of the phases, the head group
+correlation functions decay with an fast librational contribution ($12
+\pm 1$ ps).  There are additional moderate ($\tau^h_{\rm mid}$) and
+slow $\tau^h_{\rm slow}$ contributions to orientational decay that
+depend strongly on the phase exhibited by the lipids.  The symmetric
+ripple phase ($\sigma_h / d = 1.35$) appears to exhibit the slowest
+molecular reorientation.}
+\begin{tabular}{lcccc}
+\hline 
+$\sigma_h / d$ & $\tau^h_{\rm mid} (ns)$ & $\tau^h_{\rm
+slow} (\mu s)$ & $\tau^b (\mu s)$ & $D (\times 10^{-11} m^2 s^{-1})$ \\
+\hline 
+1.20 & $0.4$ &  $9.6$ & $9.5$ & $0.43(1)$ \\
+1.28 & $2.0$ & $13.5$ & $3.0$ & $5.91(3)$ \\
+1.35 & $3.2$ &  $4.0$ & $0.9$ & $3.42(1)$ \\
+1.41 & $0.3$ & $23.8$ & $6.9$ & $7.16(1)$ \\
+\end{tabular}
+\label{mdtab:relaxation}
+\end{center}
+\end{minipage}
+\end{table*}
+
+\section{Discussion}
+\label{mdsec: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{mdfig: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.
+
+The attractive forces holding the bilayer together could either be
+non-directional (as in the work of Kranenburg and
+Smit),\cite{Kranenburg2005} or directional (as we have utilized in
+these simulations).  The dipolar head groups are key for the
+maintaining the bilayer structures exhibited by this particular model;
+reducing the strength of the dipole has the tendency to make the
+rippled phase disappear.  The dipoles are likely to form attractive
+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}
+
+Experimental verification of our predictions of dipolar orientation
+correlating with the ripple direction would require knowing both the
+local orientation of a rippled region of the membrane (available via
+AFM studies of supported bilayers) as well as the local ordering of
+the membrane dipoles. Obtaining information about the local
+orientations of the membrane dipoles may be available from
+fluorescence detected linear dichroism (LD).  Benninger {\it et al.}
+have recently used axially-specific chromophores
+2-(4,4-difluoro-5,7-dimethyl-4-bora-3a,4a-diaza-s-indacene-3-pentanoyl)-1-hexadecanoyl-sn-glycero-3-phospocholine
+($\beta$-BODIPY FL C5-HPC or BODIPY-PC) and 3,3'
+dioctadecyloxacarbocyanine perchlorate (DiO) in their
+fluorescence-detected linear dichroism (LD) studies of plasma
+membranes of living cells.\cite{Benninger:2005qy} The DiO dye aligns
+its transition moment perpendicular to the membrane normal, while the
+BODIPY-PC transition dipole is parallel with the membrane normal.
+Without a doubt, using fluorescence detection of linear dichroism in
+concert with AFM surface scanning would be difficult experiments to
+carry out.  However, there is some hope of performing experiments to
+either verify or falsify the predictions of our simulations.
+
+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.