ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdRipple/mdripple.tex
(Generate patch)

Comparing trunk/mdRipple/mdripple.tex (file contents):
Revision 3198 by xsun, Wed Jul 25 22:42:56 2007 UTC vs.
Revision 3264 by gezelter, Thu Oct 18 21:07:03 2007 UTC

# Line 1 | Line 1
1   %\documentclass[aps,pre,twocolumn,amssymb,showpacs,floatfix]{revtex4}
2 < \documentclass[aps,pre,preprint,amssymb,showpacs]{revtex4}
2 > %\documentclass[aps,pre,preprint,amssymb]{revtex4}
3 > \documentclass[12pt]{article}
4 > \usepackage{times}
5 > \usepackage{mathptm}
6 > \usepackage{tabularx}
7 > \usepackage{setspace}
8 > \usepackage{amsmath}
9 > \usepackage{amssymb}
10   \usepackage{graphicx}
11 + \usepackage[ref]{overcite}
12 + \pagestyle{plain}
13 + \pagenumbering{arabic}
14 + \oddsidemargin 0.0cm \evensidemargin 0.0cm
15 + \topmargin -21pt \headsep 10pt
16 + \textheight 9.0in \textwidth 6.5in
17 + \brokenpenalty=10000
18 + \renewcommand{\baselinestretch}{1.2}
19 + \renewcommand\citemid{\ } % no comma in optional reference note
20  
21   \begin{document}
22 < \renewcommand{\thefootnote}{\fnsymbol{footnote}}
23 < \renewcommand{\theequation}{\arabic{section}.\arabic{equation}}
22 > %\renewcommand{\thefootnote}{\fnsymbol{footnote}}
23 > %\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}
24  
25 < %\bibliographystyle{aps}
25 > \bibliographystyle{achemso}
26  
27 < \title{Dipolar Ordering of the Ripple Phase in Lipid Membranes}
28 < \author{Xiuquan Sun and J. Daniel Gezelter}
29 < \email[E-mail:]{gezelter@nd.edu}
30 < \affiliation{Department of Chemistry and Biochemistry,\\
31 < University of Notre Dame, \\
27 > \title{Dipolar ordering in the ripple phases of molecular-scale models
28 > of lipid membranes}
29 > \author{Xiuquan Sun and J. Daniel Gezelter \\
30 > Department of Chemistry and Biochemistry,\\
31 > University of Notre Dame, \\
32   Notre Dame, Indiana 46556}
33  
34 + %\email[E-mail:]{gezelter@nd.edu}
35 +
36   \date{\today}
37  
38 + \maketitle
39 +
40   \begin{abstract}
41 < The ripple phase in phosphatidylcholine (PC) bilayers has never been
42 < completely explained.
41 > Symmetric and asymmetric ripple phases have been observed to form in
42 > molecular dynamics simulations of a simple molecular-scale lipid
43 > model. The lipid model consists of an dipolar head group and an
44 > ellipsoidal tail.  Within the limits of this model, an explanation for
45 > generalized membrane curvature is a simple mismatch in the size of the
46 > heads with the width of the molecular bodies.  The persistence of a
47 > {\it bilayer} structure requires strong attractive forces between the
48 > head groups.  One feature of this model is that an energetically
49 > favorable orientational ordering of the dipoles can be achieved by
50 > out-of-plane membrane corrugation.  The corrugation of the surface
51 > stabilizes the long range orientational ordering for the dipoles in the
52 > head groups which then adopt a bulk anti-ferroelectric state. We
53 > observe a common feature of the corrugated dipolar membranes: the wave
54 > vectors for the surface ripples are always found to be perpendicular
55 > to the dipole director axis.  
56   \end{abstract}
57  
58 < \pacs{}
59 < \maketitle
58 > %\maketitle
59 > \newpage
60  
61   \section{Introduction}
62   \label{sec:Int}
# Line 42 | Line 75 | within the gel phase.~\cite{Cevc87}
75   experimental results provide strong support for a 2-dimensional
76   hexagonal packing lattice of the lipid molecules within the ripple
77   phase.  This is a notable change from the observed lipid packing
78 < within the gel phase.~\cite{Cevc87}
78 > within the gel phase.~\cite{Cevc87} The X-ray diffraction work by
79 > Katsaras {\it et al.} showed that a rich phase diagram exhibiting both
80 > {\it asymmetric} and {\it symmetric} ripples is possible for lecithin
81 > bilayers.\cite{Katsaras00}
82  
83   A number of theoretical models have been presented to explain the
84   formation of the ripple phase. Marder {\it et al.} used a
85 < curvature-dependent Landau-de Gennes free-energy functional to predict
85 > curvature-dependent Landau-de~Gennes free-energy functional to predict
86   a rippled phase.~\cite{Marder84} This model and other related continuum
87   models predict higher fluidity in convex regions and that concave
88   portions of the membrane correspond to more solid-like regions.
# Line 72 | Line 108 | of lamellar stacks of hexagonal lattices to show that
108   regions.~\cite{Heimburg00} Kubica has suggested that a lattice model of
109   polar head groups could be valuable in trying to understand bilayer
110   phase formation.~\cite{Kubica02} Bannerjee used Monte Carlo simulations
111 < of lamellar stacks of hexagonal lattices to show that large headgroups
111 > of lamellar stacks of hexagonal lattices to show that large head groups
112   and molecular tilt with respect to the membrane normal vector can
113   cause bulk rippling.~\cite{Bannerjee02}
114  
115 < In contrast, few large-scale molecular modelling studies have been
115 > In contrast, few large-scale molecular modeling studies have been
116   done due to the large size of the resulting structures and the time
117   required for the phases of interest to develop.  With all-atom (and
118   even unified-atom) simulations, only one period of the ripple can be
119 < observed and only for timescales in the range of 10-100 ns.  One of
120 < the most interesting molecular simulations was carried out by De Vries
119 > observed and only for time scales in the range of 10-100 ns.  One of
120 > the most interesting molecular simulations was carried out by de~Vries
121   {\it et al.}~\cite{deVries05}. According to their simulation results,
122   the ripple consists of two domains, one resembling the gel bilayer,
123   while in the other, the two leaves of the bilayer are fully
124   interdigitated.  The mechanism for the formation of the ripple phase
125   suggested by their work is a packing competition between the head
126   groups and the tails of the lipid molecules.\cite{Carlson87} Recently,
127 < the ripple phase has also been studied by the XXX group using Monte
127 > the ripple phase has also been studied by Lenz and Schmid using Monte
128   Carlo simulations.\cite{Lenz07} Their structures are similar to the De
129   Vries {\it et al.} structures except that the connection between the
130   two leaves of the bilayer is a narrow interdigitated line instead of
# Line 103 | Line 139 | between headgroups and tails is strongly implicated as
139  
140   Although the organization of the tails of lipid molecules are
141   addressed by these molecular simulations and the packing competition
142 < between headgroups and tails is strongly implicated as the primary
142 > between head groups and tails is strongly implicated as the primary
143   driving force for ripple formation, questions about the ordering of
144 < the head groups in ripple phase has not been settled.
144 > the head groups in ripple phase have not been settled.
145  
146   In a recent paper, we presented a simple ``web of dipoles'' spin
147   lattice model which provides some physical insight into relationship
148   between dipolar ordering and membrane buckling.\cite{Sun2007} We found
149   that dipolar elastic membranes can spontaneously buckle, forming
150 < ripple-like topologies.  The driving force for the buckling in dipolar
151 < elastic membranes the antiferroelectric ordering of the dipoles, and
152 < this was evident in the ordering of the dipole director axis
153 < perpendicular to the wave vector of the surface ripples.  A similiar
150 > ripple-like topologies.  The driving force for the buckling of dipolar
151 > elastic membranes is the anti-ferroelectric ordering of the dipoles.
152 > This was evident in the ordering of the dipole director axis
153 > perpendicular to the wave vector of the surface ripples.  A similar
154   phenomenon has also been observed by Tsonchev {\it et al.} in their
155 < work on the spontaneous formation of dipolar molecules into curved
156 < nano-structures.\cite{Tsonchev04}
155 > work on the spontaneous formation of dipolar peptide chains into
156 > curved nano-structures.\cite{Tsonchev04,Tsonchev04II}
157  
158   In this paper, we construct a somewhat more realistic molecular-scale
159   lipid model than our previous ``web of dipoles'' and use molecular
# Line 133 | Line 169 | Our simple molecular-scale lipid model for studying th
169   \section{Computational Model}
170   \label{sec:method}
171  
172 + \begin{figure}[htb]
173 + \centering
174 + \includegraphics[width=4in]{lipidModels}
175 + \caption{Three different representations of DPPC lipid molecules,
176 + including the chemical structure, an atomistic model, and the
177 + head-body ellipsoidal coarse-grained model used in this
178 + work.\label{fig:lipidModels}}
179 + \end{figure}
180 +
181   Our simple molecular-scale lipid model for studying the ripple phase
182   is based on two facts: one is that the most essential feature of lipid
183   molecules is their amphiphilic structure with polar head groups and
184   non-polar tails. Another fact is that the majority of lipid molecules
185   in the ripple phase are relatively rigid (i.e. gel-like) which makes
186   some fraction of the details of the chain dynamics negligible.  Figure
187 < \ref{fig:lipidModels} shows the molecular strucure of a DPPC
187 > \ref{fig:lipidModels} shows the molecular structure of a DPPC
188   molecule, as well as atomistic and molecular-scale representations of
189   a DPPC molecule.  The hydrophilic character of the head group is
190   largely due to the separation of charge between the nitrogen and
# Line 155 | Line 200 | the point dipole rigidly in this orientation.  
200   nearly perpendicular to the tail, so we have fixed the direction of
201   the point dipole rigidly in this orientation.  
202  
158 \begin{figure}[htb]
159 \centering
160 \includegraphics[width=\linewidth]{lipidModels}
161 \caption{Three different representations of DPPC lipid molecules,
162 including the chemical structure, an atomistic model, and the
163 head-body ellipsoidal coarse-grained model used in this
164 work.\label{fig:lipidModels}}
165 \end{figure}
166
203   The ellipsoidal portions of the model interact via the Gay-Berne
204   potential which has seen widespread use in the liquid crystal
205 < community.  In its original form, the Gay-Berne potential was a single
206 < site model for the interactions of rigid ellipsoidal
205 > community.  Ayton and Voth have also used Gay-Berne ellipsoids for
206 > modeling large length-scale properties of lipid
207 > bilayers.\cite{Ayton01} In its original form, the Gay-Berne potential
208 > was a single site model for the interactions of rigid ellipsoidal
209   molecules.\cite{Gay81} It can be thought of as a modification of the
210   Gaussian overlap model originally described by Berne and
211   Pechukas.\cite{Berne72} The potential is constructed in the familiar
212   form of the Lennard-Jones function using orientation-dependent
213   $\sigma$ and $\epsilon$ parameters,
214 < \begin{eqnarray*}
214 > \begin{equation*}
215   V_{ij}({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat
216   r}_{ij}}) = 4\epsilon ({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
217   {\mathbf{\hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i},
# Line 181 | Line 219 | -\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_
219   -\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
220   {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^6\right]
221   \label{eq:gb}
222 < \end{eqnarray*}
222 > \end{equation*}
223  
186
187
224   The range $(\sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
225 < \hat{r}}))$, and strength $(\epsilon({\bf \hat{u}}_{i},{\bf
226 < \hat{u}}_{j},{\bf \hat{r}}))$ parameters
225 > \hat{r}}_{ij}))$, and strength $(\epsilon({\bf \hat{u}}_{i},{\bf
226 > \hat{u}}_{j},{\bf \hat{r}}_{ij}))$ parameters
227   are dependent on the relative orientations of the two molecules (${\bf
228   \hat{u}}_{i},{\bf \hat{u}}_{j}$) as well as the direction of the
229 < intermolecular separation (${\bf \hat{r}}$).  The functional forms for
230 < $\sigma({\bf
231 < \hat{u}}_{i},{\bf
232 < \hat{u}}_{j},{\bf \hat{r}})$ and $\epsilon({\bf \hat{u}}_{i},{\bf
233 < \hat{u}}_{j},{\bf \hat{r}}))$ are given elsewhere
234 < and will not be repeated here.  However, $\epsilon$ and $\sigma$ are
235 < governed by two anisotropy parameters,
236 < \begin {equation}
237 < \begin{array}{rcl}
238 < \chi & = & \frac
239 < {(\sigma_{e}/\sigma_{s})^2-1}{(\sigma_{e}/\sigma_{s})^2+1} \\ \\
240 < \chi\prime & = & \frac {1-(\epsilon_{e}/\epsilon_{s})^{1/\mu}}{1+(\epsilon_{e}/
241 < \epsilon_{s})^{1/\mu}}
242 < \end{array}
243 < \end{equation}
244 < In these equations, $\sigma$ and $\epsilon$ refer to the point of
245 < closest contact and the depth of the well in different orientations of
246 < the two molecules.  The subscript $s$ refers to the {\it side-by-side}
247 < configuration where $\sigma$ has it's smallest value,
248 < $\sigma_{s}$, and where the potential well is $\epsilon_{s}$ deep.
249 < The subscript $e$ refers to the {\it end-to-end} configuration where
250 < $\sigma$ is at it's largest value, $\sigma_{e}$ and where the well
251 < depth, $\epsilon_{e}$ is somewhat smaller than in the side-by-side
252 < configuration.  For the prolate ellipsoids we are using, we have
229 > intermolecular separation (${\bf \hat{r}}_{ij}$).  $\sigma$ and
230 > $\sigma_0$ are also governed by shape mixing and anisotropy variables,
231 > \begin {eqnarray*}
232 > \sigma_0 & = & \sqrt{d_i^2 + d_j^2} \\ \\
233 > \chi & = & \left[ \frac{\left( l_i^2 - d_i^2 \right) \left(l_j^2 -
234 > d_j^2 \right)}{\left( l_j^2 + d_i^2 \right) \left(l_i^2 +
235 > d_j^2 \right)}\right]^{1/2} \\ \\
236 > \alpha^2 & = & \left[ \frac{\left( l_i^2 - d_i^2 \right) \left(l_j^2 +
237 > d_i^2 \right)}{\left( l_j^2 - d_j^2 \right) \left(l_i^2 +
238 > d_j^2 \right)}\right]^{1/2},
239 > \end{eqnarray*}
240 > where $l$ and $d$ describe the length and width of each uniaxial
241 > ellipsoid.  These shape anisotropy parameters can then be used to
242 > calculate the range function,
243 > \begin{equation*}
244 > \sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf \hat{r}}_{ij}) = \sigma_{0}
245 > \left[ 1- \left\{ \frac{ \chi \alpha^2 ({\bf \hat{u}}_i \cdot {\bf
246 > \hat{r}}_{ij} ) + \chi \alpha^{-2} ({\bf \hat{u}}_j \cdot {\bf
247 > \hat{r}}_{ij} ) - 2 \chi^2 ({\bf \hat{u}}_i \cdot {\bf
248 > \hat{r}}_{ij} )({\bf \hat{u}}_j \cdot {\bf
249 > \hat{r}}_{ij} ) ({\bf \hat{u}}_i \cdot {\bf \hat{u}}_j)}{1 - \chi^2
250 > \left({\bf \hat{u}}_i \cdot {\bf \hat{u}}_j\right)^2} \right\}
251 > \right]^{-1/2}
252 > \end{equation*}
253 >
254 > Gay-Berne ellipsoids also have an energy scaling parameter,
255 > $\epsilon^s$, which describes the well depth for two identical
256 > ellipsoids in a {\it side-by-side} configuration.  Additionally, a well
257 > depth aspect ratio, $\epsilon^r = \epsilon^e / \epsilon^s$, describes
258 > the ratio between the well depths in the {\it end-to-end} and
259 > side-by-side configurations.  As in the range parameter, a set of
260 > mixing and anisotropy variables can be used to describe the well
261 > depths for dissimilar particles,
262 > \begin {eqnarray*}
263 > \epsilon_0 & = & \sqrt{\epsilon^s_i  * \epsilon^s_j} \\ \\
264 > \epsilon^r & = & \sqrt{\epsilon^r_i  * \epsilon^r_j} \\ \\
265 > \chi' & = & \frac{1 - (\epsilon^r)^{1/\mu}}{1 + (\epsilon^r)^{1/\mu}}
266 > \\ \\
267 > \alpha'^2 & = & \left[1 + (\epsilon^r)^{1/\mu}\right]^{-1}
268 > \end{eqnarray*}
269 > The form of the strength function is somewhat complicated,
270 > \begin {eqnarray*}
271 > \epsilon({\bf \hat{u}}_{i}, {\bf \hat{u}}_{j},{\bf \hat{r}}_{ij}) & = &
272 > \epsilon_{0}  \epsilon_{1}^{\nu}({\bf \hat{u}}_{i}.{\bf \hat{u}}_{j})
273 > \epsilon_{2}^{\mu}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
274 > \hat{r}}_{ij}) \\ \\
275 > \epsilon_{1}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j}) & = &
276 > \left[1-\chi^{2}({\bf \hat{u}}_{i}.{\bf
277 > \hat{u}}_{j})^{2}\right]^{-1/2} \\ \\
278 > \epsilon_{2}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf \hat{r}}_{ij}) &
279 > = &
280 > 1 - \left\{ \frac{ \chi' \alpha'^2 ({\bf \hat{u}}_i \cdot {\bf
281 > \hat{r}}_{ij} ) + \chi' \alpha'^{-2} ({\bf \hat{u}}_j \cdot {\bf
282 > \hat{r}}_{ij} ) - 2 \chi'^2 ({\bf \hat{u}}_i \cdot {\bf
283 > \hat{r}}_{ij} )({\bf \hat{u}}_j \cdot {\bf
284 > \hat{r}}_{ij} ) ({\bf \hat{u}}_i \cdot {\bf \hat{u}}_j)}{1 - \chi'^2
285 > \left({\bf \hat{u}}_i \cdot {\bf \hat{u}}_j\right)^2} \right\},
286 > \end {eqnarray*}
287 > although many of the quantities and derivatives are identical with
288 > those obtained for the range parameter. Ref. \citen{Luckhurst90}
289 > has a particularly good explanation of the choice of the Gay-Berne
290 > parameters $\mu$ and $\nu$ for modeling liquid crystal molecules. An
291 > excellent overview of the computational methods that can be used to
292 > efficiently compute forces and torques for this potential can be found
293 > in Ref. \citen{Golubkov06}
294 >
295 > The choices of parameters we have used in this study correspond to a
296 > shape anisotropy of 3 for the chain portion of the molecule.  In
297 > principle, this could be varied to allow for modeling of longer or
298 > shorter chain lipid molecules. For these prolate ellipsoids, we have:
299   \begin{equation}
300   \begin{array}{rcl}
301 < \sigma_{s} & < & \sigma_{e} \\
302 < \epsilon_{s} & > & \epsilon_{e}
301 > d & < & l \\
302 > \epsilon^{r} & < & 1
303   \end{array}
304   \end{equation}
305 < Ref. \onlinecite{Luckhurst90} has a particularly good explanation of the
306 < choice of the Gay-Berne parameters $\mu$ and $\nu$ for modeling liquid
307 < crystal molecules.
305 > A sketch of the various structural elements of our molecular-scale
306 > lipid / solvent model is shown in figure \ref{fig:lipidModel}.  The
307 > actual parameters used in our simulations are given in table
308 > \ref{tab:parameters}.
309  
310 < The breadth and length of tail are $\sigma_0$, $3\sigma_0$,
311 < corresponding to a shape anisotropy of 3 for the chain portion of the
312 < molecule.  In principle, this could be varied to allow for modeling of
313 < longer or shorter chain lipid molecules.
310 > \begin{figure}[htb]
311 > \centering
312 > \includegraphics[width=4in]{2lipidModel}
313 > \caption{The parameters defining the behavior of the lipid
314 > models. $l / d$ is the ratio of the head group to body diameter.
315 > Molecular bodies had a fixed aspect ratio of 3.0.  The solvent model
316 > was a simplified 4-water bead ($\sigma_w \approx d$) that has been
317 > used in other coarse-grained (DPD) simulations.  The dipolar strength
318 > (and the temperature and pressure) were the only other parameters that
319 > were varied systematically.\label{fig:lipidModel}}
320 > \end{figure}
321  
322   To take into account the permanent dipolar interactions of the
323 < zwitterionic head groups, we place fixed dipole moments $\mu_{i}$ at
324 < one end of the Gay-Berne particles.  The dipoles will be oriented at
325 < an angle $\theta = \pi / 2$ relative to the major axis.  These dipoles
326 < are protected by a head ``bead'' with a range parameter which we have
327 < varied between $1.20\sigma_0$ and $1.41\sigma_0$.  The head groups
328 < interact with each other using a combination of Lennard-Jones,
329 < \begin{eqnarray*}
330 < V_{ij} = 4\epsilon \left[\left(\frac{\sigma_h}{r_{ij}}\right)^{12} -
323 > zwitterionic head groups, we have placed fixed dipole moments $\mu_{i}$ at
324 > one end of the Gay-Berne particles.  The dipoles are oriented at an
325 > angle $\theta = \pi / 2$ relative to the major axis.  These dipoles
326 > are protected by a head ``bead'' with a range parameter ($\sigma_h$) which we have
327 > varied between $1.20 d$ and $1.41 d$.  The head groups interact with
328 > each other using a combination of Lennard-Jones,
329 > \begin{equation}
330 > V_{ij}(r_{ij}) = 4\epsilon_h \left[\left(\frac{\sigma_h}{r_{ij}}\right)^{12} -
331   \left(\frac{\sigma_h}{r_{ij}}\right)^6\right],
332 < \end{eqnarray*}
333 < and dipole,
334 < \begin{eqnarray*}
335 < V_{ij} = \frac{|\mu|^2}{4 \pi \epsilon_0 r_{ij}^3}
332 > \end{equation}
333 > and dipole-dipole,
334 > \begin{equation}
335 > V_{ij}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
336 > \hat{r}}_{ij})) = \frac{|\mu|^2}{4 \pi \epsilon_0 r_{ij}^3}
337   \left[ \hat{\bf u}_i \cdot \hat{\bf u}_j - 3 (\hat{\bf u}_i \cdot
338   \hat{\bf r}_{ij})(\hat{\bf u}_j \cdot \hat{\bf r}_{ij}) \right]
339 < \end{eqnarray*}
339 > \end{equation}
340   potentials.  
341   In these potentials, $\mathbf{\hat u}_i$ is the unit vector pointing
342   along dipole $i$ and $\mathbf{\hat r}_{ij}$ is the unit vector
343 < pointing along the inter-dipole vector $\mathbf{r}_{ij}$.
343 > pointing along the inter-dipole vector $\mathbf{r}_{ij}$.  
344  
345   For the interaction between nonequivalent uniaxial ellipsoids (in this
346 < case, between spheres and ellipsoids), the range parameter is
347 < generalized as\cite{Cleaver96}
348 < \begin{eqnarray*}
349 < \sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}}) =
350 < {\sigma_{0}} \left(1-\frac{1}{2}\chi\left[\frac{(\alpha{\mathbf{\hat r}_{ij}}
351 < \cdot {\mathbf{\hat u}_i} + \alpha^{-1}{\mathbf{\hat r}_{ij}} \cdot {\mathbf{\hat
261 < u}_j})^2}{1+\chi({\mathbf{\hat u}_i} \cdot {\mathbf{\hat u}_j})} +
262 < \frac{(\alpha{\mathbf{\hat r}_{ij}} \cdot {\mathbf{\hat u}_i} - \alpha^{-1}{\mathbf{\hat r}_{ij}}
263 < \cdot {\mathbf{\hat u}_j})^2}{1-\chi({\mathbf{\hat u}_i} \cdot
264 < {\mathbf{\hat u}_j})} \right] \right)^{-\frac{1}{2}}
265 < \end{eqnarray*}
266 < where $\alpha$ is given by
267 < \begin{eqnarray*}
268 < \alpha^2 =
269 < \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)}
270 < \right]^{\frac{1}{2}}
271 < \end{eqnarray*}
272 < the strength parameter has been adjusted as suggested by Cleaver {\it
273 < et al.}\cite{Cleaver96}  A switching function has been applied to all
274 < potentials to smoothly turn off the interactions between a range of  $22$  and $25$ \AA.
346 > case, between spheres and ellipsoids), the spheres are treated as
347 > ellipsoids with an aspect ratio of 1 ($d = l$) and with an well depth
348 > ratio ($\epsilon^r$) of 1 ($\epsilon^e = \epsilon^s$).  The form of
349 > the Gay-Berne potential we are using was generalized by Cleaver {\it
350 > et al.} and is appropriate for dissimilar uniaxial
351 > ellipsoids.\cite{Cleaver96}
352  
353 < The solvent model in our simulations is identical to one used by XXX
354 < in their dissipative particle dynamics (DPD) simulation of lipid
355 < bilayers.]cite{XXX} This solvent bead is a single site that represents
356 < four water molecules (m = 72 amu) and has comparable density and
357 < diffusive behavior to liquid water.  However, since there are no
358 < electrostatic sites on these beads, this solvent model cannot
359 < replicate the dielectric properties of water.
353 > The solvent model in our simulations is identical to one used by
354 > Marrink {\it et al.}  in their dissipative particle dynamics (DPD)
355 > simulation of lipid bilayers.\cite{Marrink04} This solvent bead is a
356 > single site that represents four water molecules (m = 72 amu) and has
357 > comparable density and diffusive behavior to liquid water.  However,
358 > since there are no electrostatic sites on these beads, this solvent
359 > model cannot replicate the dielectric properties of water.
360 >
361   \begin{table*}
362   \begin{minipage}{\linewidth}
363   \begin{center}
364 < \caption{}
365 < \begin{tabular}{lccc}
364 > \caption{Potential parameters used for molecular-scale coarse-grained
365 > lipid simulations}
366 > \begin{tabular}{llccc}
367   \hline
368 < N/A & Head & Chain & Solvent \\
368 >  & &  Head & Chain & Solvent \\
369   \hline
370 < $\sigma_0$ (\AA) & varied & 4.6 & 4.7 \\
371 < l (aspect ratio) & N/A & 3 & N/A \\
372 < $\epsilon_0$ (kcal/mol) & 0.185 & 1.0 & 0.8 \\
373 < $\epsilon_e$ (aspect ratio) & N/A & 0.2 & N/A \\
374 < M (amu) & 196 & 760 & 72.06112 \\
375 < $I_{xx}$, $I_{yy}$, $I_{zz}$ (amu \AA$^2$) & 1125, 1125, 0 & 45000, 45000, 9000 & N/A \\
376 < $\mu$ (Debye) & varied & N/A & N/A \\
370 > $d$ (\AA) & & varied & 4.6  & 4.7 \\
371 > $l$ (\AA) & & $= d$ & 13.8 & 4.7 \\
372 > $\epsilon^s$ (kcal/mol) & & 0.185 & 1.0 & 0.8 \\
373 > $\epsilon_r$ (well-depth aspect ratio)& & 1 & 0.2 &  1 \\
374 > $m$ (amu) & & 196 & 760 & 72.06 \\
375 > $\overleftrightarrow{\mathsf I}$ (amu \AA$^2$) & & & & \\
376 > \multicolumn{2}c{$I_{xx}$} & 1125 & 45000 & N/A \\
377 > \multicolumn{2}c{$I_{yy}$} & 1125 & 45000 & N/A \\
378 > \multicolumn{2}c{$I_{zz}$} &  0 &    9000 & N/A \\
379 > $\mu$ (Debye) & & varied & 0 & 0 \\
380   \end{tabular}
381   \label{tab:parameters}
382   \end{center}
383   \end{minipage}
384   \end{table*}
385  
304 \begin{figure}[htb]
305 \centering
306 \includegraphics[width=\linewidth]{2lipidModel}
307 \caption{The parameters defining the behavior of the lipid
308 models. $\sigma_h / \sigma_0$ is the ratio of the head group to body
309 diameter.  Molecular bodies had a fixed aspect ratio of 3.0.  The
310 solvent model was a simplified 4-water bead ($\sigma_w = 1.02
311 \sigma_0$) that has been used in other coarse-grained (DPD) simulations.
312 The dipolar strength (and the temperature and pressure) were the only
313 other parameters that were varied
314 systematically.\label{fig:lipidModel}}
315 \end{figure}
316
386   \section{Experimental Methodology}
387   \label{sec:experiment}
388  
389 + The parameters that were systematically varied in this study were the
390 + size of the head group ($\sigma_h$), the strength of the dipole moment
391 + ($\mu$), and the temperature of the system.  Values for $\sigma_h$
392 + ranged from 5.5 \AA\ to 6.5 \AA\ .  If the width of the tails is taken
393 + to be the unit of length, these head groups correspond to a range from
394 + $1.2 d$ to $1.41 d$.  Since the solvent beads are nearly identical in
395 + diameter to the tail ellipsoids, all distances that follow will be
396 + measured relative to this unit of distance.  Because the solvent we
397 + are using is non-polar and has a dielectric constant of 1, values for
398 + $\mu$ are sampled from a range that is somewhat smaller than the 20.6
399 + Debye dipole moment of the PC head groups.
400 +
401   To create unbiased bilayers, all simulations were started from two
402 < perfectly flat monolayers separated by a 20 \AA\ gap between the
402 > perfectly flat monolayers separated by a 26 \AA\ gap between the
403   molecular bodies of the upper and lower leaves.  The separated
404 < monolayers were evolved in a vaccum with $x-y$ anisotropic pressure
404 > monolayers were evolved in a vacuum with $x-y$ anisotropic pressure
405   coupling. The length of $z$ axis of the simulations was fixed and a
406   constant surface tension was applied to enable real fluctuations of
407 < the bilayer. Periodic boundaries were used, and $480-720$ lipid
408 < molecules were present in the simulations depending on the size of the
409 < head beads.  The two monolayers spontaneously collapse into bilayer
410 < structures within 100 ps, and following this collapse, all systems
411 < were equlibrated for $100$ ns at $300$ K.
407 > the bilayer. Periodic boundary conditions were used, and $480-720$
408 > lipid molecules were present in the simulations, depending on the size
409 > of the head beads.  In all cases, the two monolayers spontaneously
410 > collapsed into bilayer structures within 100 ps. Following this
411 > collapse, all systems were equilibrated for $100$ ns at $300$ K.
412  
413 < The resulting structures were then solvated at a ratio of $6$ DPD
413 > The resulting bilayer structures were then solvated at a ratio of $6$
414   solvent beads (24 water molecules) per lipid. These configurations
415 < were then equilibrated for another $30$ ns. All simulations with
416 < solvent were carried out at constant pressure ($P=1$ atm) by $3$D
417 < anisotropic coupling, and constant surface tension ($\gamma=0.015$
418 < UNIT). Given the absence of fast degrees of freedom in this model, a
419 < timestep of $50$ fs was utilized.  Data collection for structural
420 < properties of the bilayers was carried out during a final 5 ns run
421 < following the solvent equilibration.  All simulations were performed
422 < using the OOPSE molecular modeling program.\cite{Meineke05}
415 > were then equilibrated for another $30$ ns. All simulations utilizing
416 > the solvent were carried out at constant pressure ($P=1$ atm) with
417 > $3$D anisotropic coupling, and constant surface tension
418 > ($\gamma=0.015$ N/m). Given the absence of fast degrees of freedom in
419 > this model, a time step of $50$ fs was utilized with excellent energy
420 > conservation.  Data collection for structural properties of the
421 > bilayers was carried out during a final 5 ns run following the solvent
422 > equilibration.  All simulations were performed using the OOPSE
423 > molecular modeling program.\cite{Meineke05}
424  
425 + A switching function was applied to all potentials to smoothly turn
426 + off the interactions between a range of $22$ and $25$ \AA.
427 +
428   \section{Results}
429   \label{sec:results}
430  
431 < Snapshots in Figure \ref{fig:phaseCartoon} show that the membrane is
432 < more corrugated increasing size of the head groups. The surface is
433 < nearly flat when $\sigma_h=1.20\sigma_0$. With
434 < $\sigma_h=1.28\sigma_0$, although the surface is still flat, the
435 < bilayer starts to splay inward; the upper leaf of the bilayer is
436 < connected to the lower leaf with an interdigitated line defect. Two
437 < periodicities with $100$ \AA\ width were observed in the
438 < simulation. This structure is very similiar to the structure observed
439 < by de Vries and Lenz {\it et al.}. The same basic structure is also
440 < observed when $\sigma_h=1.41\sigma_0$, but the wavelength of the
441 < surface corrugations depends sensitively on the size of the ``head''
442 < beads. From the undulation spectrum, the corrugation is clearly
443 < non-thermal.
431 > The membranes in our simulations exhibit a number of interesting
432 > bilayer phases.  The surface topology of these phases depends most
433 > sensitively on the ratio of the size of the head groups to the width
434 > of the molecular bodies.  With heads only slightly larger than the
435 > bodies ($\sigma_h=1.20 d$) the membrane exhibits a flat bilayer.
436 >
437 > Increasing the head / body size ratio increases the local membrane
438 > curvature around each of the lipids.  With $\sigma_h=1.28 d$, the
439 > surface is still essentially flat, but the bilayer starts to exhibit
440 > signs of instability.  We have observed occasional defects where a
441 > line of lipid molecules on one leaf of the bilayer will dip down to
442 > interdigitate with the other leaf.  This gives each of the two bilayer
443 > leaves some local convexity near the line defect.  These structures,
444 > once developed in a simulation, are very stable and are spaced
445 > approximately 100 \AA\ away from each other.
446 >
447 > With larger heads ($\sigma_h = 1.35 d$) the membrane curvature
448 > resolves into a ``symmetric'' ripple phase.  Each leaf of the bilayer
449 > is broken into several convex, hemicylinderical sections, and opposite
450 > leaves are fitted together much like roof tiles.  There is no
451 > interdigitation between the upper and lower leaves of the bilayer.
452 >
453 > For the largest head / body ratios studied ($\sigma_h = 1.41 d$) the
454 > local curvature is substantially larger, and the resulting bilayer
455 > structure resolves into an asymmetric ripple phase.  This structure is
456 > very similar to the structures observed by both de~Vries {\it et al.}
457 > and Lenz {\it et al.}.  For a given ripple wave vector, there are two
458 > possible asymmetric ripples, which is not the case for the symmetric
459 > phase observed when $\sigma_h = 1.35 d$.
460 >
461   \begin{figure}[htb]
462   \centering
463 < \includegraphics[width=\linewidth]{phaseCartoon}
464 < \caption{A sketch to discribe the structure of the phases observed in
465 < our simulations.\label{fig:phaseCartoon}}
463 > \includegraphics[width=4in]{phaseCartoon}
464 > \caption{The role of the ratio between the head group size and the
465 > width of the molecular bodies is to increase the local membrane
466 > curvature.  With strong attractive interactions between the head
467 > groups, this local curvature can be maintained in bilayer structures
468 > through surface corrugation.  Shown above are three phases observed in
469 > these simulations.  With $\sigma_h = 1.20 d$, the bilayer maintains a
470 > flat topology.  For larger heads ($\sigma_h = 1.35 d$) the local
471 > curvature resolves into a symmetrically rippled phase with little or
472 > no interdigitation between the upper and lower leaves of the membrane.
473 > The largest heads studied ($\sigma_h = 1.41 d$) resolve into an
474 > asymmetric rippled phases with interdigitation between the two
475 > leaves.\label{fig:phaseCartoon}}
476   \end{figure}
477  
478 < When $\sigma_h=1.35\sigma_0$, we observed another corrugated surface
479 < morphology. This structure is different from the asymmetric rippled
480 < surface; there is no interdigitation between the upper and lower
481 < leaves of the bilayer. Each leaf of the bilayer is broken into several
482 < hemicylinderical sections, and opposite leaves are fitted together
483 < much like roof tiles. Unlike the surface in which the upper
484 < hemicylinder is always interdigitated on the leading or trailing edge
485 < of lower hemicylinder, the symmetric ripple has no prefered direction.
486 < The corresponding cartoons are shown in Figure
487 < \ref{fig:phaseCartoon} for elucidation of the detailed structures of
488 < different phases. Figure \ref{fig:phaseCartoon} (a) is the flat phase,
489 < (b) is the asymmetric ripple phase corresponding to the lipid
490 < organization when $\sigma_h=1.28\sigma_0$ and $\sigma_h=1.41\sigma_0$,
491 < and (c) is the symmetric ripple phase observed when
492 < $\sigma_h=1.35\sigma_0$. In symmetric ripple phase, the bilayer is
493 < continuous everywhere on the whole membrane, however, in asymmetric
494 < ripple phase, the bilayer is intermittent domains connected by thin
495 < interdigitated monolayer which consists of upper and lower leaves of
496 < the bilayer.
478 > Sample structures for the flat ($\sigma_h = 1.20 d$), symmetric
479 > ($\sigma_h = 1.35 d$, and asymmetric ($\sigma_h = 1.41 d$) ripple
480 > phases are shown in Figure \ref{fig:phaseCartoon}.  
481 >
482 > It is reasonable to ask how well the parameters we used can produce
483 > bilayer properties that match experimentally known values for real
484 > lipid bilayers.  Using a value of $l = 13.8$ \AA for the ellipsoidal
485 > tails and the fixed ellipsoidal aspect ratio of 3, our values for the
486 > area per lipid ($A$) and inter-layer spacing ($D_{HH}$) depend
487 > entirely on the size of the head bead relative to the molecular body.
488 > These values are tabulated in table \ref{tab:property}.  Kucera {\it
489 > et al.}  have measured values for the head group spacings for a number
490 > of PC lipid bilayers that range from 30.8 \AA (DLPC) to 37.8 (DPPC).
491 > They have also measured values for the area per lipid that range from
492 > 60.6
493 > \AA$^2$ (DMPC) to 64.2 \AA$^2$
494 > (DPPC).\cite{NorbertKucerka04012005,NorbertKucerka06012006} Only the
495 > largest of the head groups we modeled ($\sigma_h = 1.41 d$) produces
496 > bilayers (specifically the area per lipid) that resemble real PC
497 > bilayers.  The smaller head beads we used are perhaps better models
498 > for PE head groups.
499 >
500   \begin{table*}
501   \begin{minipage}{\linewidth}
502   \begin{center}
503 < \caption{}
504 < \begin{tabular}{lccc}
503 > \caption{Phase, bilayer spacing, area per lipid, ripple wavelength
504 > and amplitude observed as a function of the ratio between the head
505 > beads and the diameters of the tails.  Ripple wavelengths and
506 > amplitudes are normalized to the diameter of the tail ellipsoids.}
507 > \begin{tabular}{lccccc}
508   \hline
509 < $\sigma_h/\sigma_0$ & type of phase & $\lambda/\sigma_0$ & $A/\sigma_0$\\
509 > $\sigma_h / d$ & type of phase & bilayer spacing (\AA) & area per
510 > lipid (\AA$^2$) & $\lambda / d$ & $A / d$\\
511   \hline
512 < 1.20 & flat & N/A & N/A \\
513 < 1.28 & asymmetric flat & 21.7 & N/A \\
514 < 1.35 & symmetric ripple & 17.2 & 2.2 \\
515 < 1.41 & asymmetric ripple & 15.4 & 1.5 \\
512 > 1.20 & flat & 33.4 & 49.6 & N/A & N/A \\
513 > 1.28 & flat & 33.7 & 54.7 & N/A & N/A \\
514 > 1.35 & symmetric ripple & 42.9 & 51.7 & 17.2 & 2.2 \\
515 > 1.41 & asymmetric ripple & 37.1 & 63.1 & 15.4 & 1.5 \\
516   \end{tabular}
517   \label{tab:property}
518   \end{center}
519   \end{minipage}
520   \end{table*}
521  
522 < The membrane structures and the reduced wavelength $\lambda/\sigma_0$,
523 < reduced amplitude $A/\sigma_0$ of the ripples are summerized in Table
524 < \ref{tab:property}. The wavelength range is $15~21$ and the amplitude
525 < is $1.5$ for asymmetric ripple and $2.2$ for symmetric ripple. These
526 < values are consistent to the experimental results. Note, the
527 < amplitudes are underestimated without the melted tails in our
528 < simulations.
522 > The membrane structures and the reduced wavelength $\lambda / d$,
523 > reduced amplitude $A / d$ of the ripples are summarized in Table
524 > \ref{tab:property}. The wavelength range is $15 - 17$ molecular bodies
525 > and the amplitude is $1.5$ molecular bodies for asymmetric ripple and
526 > $2.2$ for symmetric ripple. These values are reasonably consistent
527 > with experimental measurements.\cite{Sun96,Katsaras00,Kaasgaard03}
528 > Note, that given the lack of structural freedom in the tails of our
529 > model lipids, the amplitudes observed from these simulations are
530 > likely to underestimate of the true amplitudes.
531  
532   \begin{figure}[htb]
533   \centering
534 < \includegraphics[width=\linewidth]{topDown}
535 < \caption{Top views of the flat (upper), asymmetric ripple (middle),
536 < and symmetric ripple (lower) phases.  Note that the head-group dipoles
537 < have formed head-to-tail chains in all three of these phases, but in
538 < the two rippled phases, the dipolar chains are all aligned
539 < {\it perpendicular} to the direction of the ripple.  The flat membrane
540 < has multiple point defects in the dipolar orientational ordering, and
541 < the dipolar ordering on the lower leaf of the bilayer can be in a
542 < different direction from the upper leaf.\label{fig:topView}}
534 > \includegraphics[width=4in]{topDown}
535 > \caption{Top views of the flat (upper), symmetric ripple (middle),
536 > and asymmetric ripple (lower) phases.  Note that the head-group
537 > dipoles have formed head-to-tail chains in all three of these phases,
538 > but in the two rippled phases, the dipolar chains are all aligned {\it
539 > perpendicular} to the direction of the ripple.  Note that the flat
540 > membrane has multiple vortex defects in the dipolar ordering, and the
541 > ordering on the lower leaf of the bilayer can be in an entirely
542 > different direction from the upper leaf.\label{fig:topView}}
543   \end{figure}
544  
545 < The $P_2$ order paramters (for molecular bodies and head group
546 < dipoles) have been calculated to clarify the ordering in these phases
547 < quantatively. $P_2=1$ means a perfectly ordered structure, and $P_2=0$
548 < implies orientational randomization. Figure \ref{fig:rP2} shows the
549 < $P_2$ order paramter of the dipoles on head group rising with
550 < increasing head group size. When the heads of the lipid molecules are
551 < small, the membrane is flat. The dipolar ordering is essentially
552 < frustrated on orientational ordering in this circumstance. Figure
553 < \ref{fig:topView} shows the snapshots of the top view for the flat system
554 < ($\sigma_h=1.20\sigma$) and rippled system
555 < ($\sigma_h=1.41\sigma$). The pointing direction of the dipoles on the
556 < head groups are represented by two colored half spheres from blue to
557 < yellow. For flat surfaces, the system obviously shows frustration on
558 < the dipolar ordering, there are kinks on the edge of defferent
559 < domains. Another reason is that the lipids can move independently in
560 < each monolayer, it is not nessasory for the direction of dipoles on
561 < one leaf is consistant to another layer, which makes total order
562 < parameter is relatively low. With increasing head group size, the
563 < surface is corrugated, and dipoles do not move as freely on the
564 < surface. Therefore, the translational freedom of lipids in one layer
565 < is dependent upon the position of lipids in another layer, as a
566 < result, the symmetry of the dipoles on head group in one layer is tied
567 < to the symmetry in the other layer. Furthermore, as the membrane
568 < deforms from two to three dimensions due to the corrugation, the
569 < symmetry of the ordering for the dipoles embedded on each leaf is
570 < broken. The dipoles then self-assemble in a head-tail configuration,
571 < and the order parameter increases dramaticaly. However, the total
451 < polarization of the system is still close to zero. This is strong
452 < evidence that the corrugated structure is an antiferroelectric
453 < state. From the snapshot in Figure \ref{}, the dipoles arrange as
454 < arrays along $Y$ axis and fall into head-to-tail configuration in each
455 < line, but every $3$ or $4$ lines of dipoles change their direction
456 < from neighbour lines. The system shows antiferroelectric
457 < charactoristic as a whole. The orientation of the dipolar is always
458 < perpendicular to the ripple wave vector. These results are consistent
459 < with our previous study on dipolar membranes.
545 > The principal method for observing orientational ordering in dipolar
546 > or liquid crystalline systems is the $P_2$ order parameter (defined
547 > as $1.5 \times \lambda_{max}$, where $\lambda_{max}$ is the largest
548 > eigenvalue of the matrix,
549 > \begin{equation}
550 > {\mathsf{S}} = \frac{1}{N} \sum_i \left(
551 > \begin{array}{ccc}
552 >        u^{x}_i u^{x}_i-\frac{1}{3} & u^{x}_i u^{y}_i & u^{x}_i u^{z}_i \\
553 >        u^{y}_i u^{x}_i  & u^{y}_i u^{y}_i -\frac{1}{3} & u^{y}_i u^{z}_i \\
554 >        u^{z}_i u^{x}_i & u^{z}_i u^{y}_i  & u^{z}_i u^{z}_i -\frac{1}{3}
555 > \end{array} \right).
556 > \label{eq:opmatrix}
557 > \end{equation}
558 > Here $u^{\alpha}_i$ is the $\alpha=x,y,z$ component of the unit vector
559 > for molecule $i$.  (Here, $\hat{\bf u}_i$ can refer either to the
560 > principal axis of the molecular body or to the dipole on the head
561 > group of the molecule.)  $P_2$ will be $1.0$ for a perfectly-ordered
562 > system and near $0$ for a randomized system.  Note that this order
563 > parameter is {\em not} equal to the polarization of the system.  For
564 > example, the polarization of a perfect anti-ferroelectric arrangement
565 > of point dipoles is $0$, but $P_2$ for the same system is $1$.  The
566 > eigenvector of $\mathsf{S}$ corresponding to the largest eigenvalue is
567 > familiar as the director axis, which can be used to determine a
568 > privileged axis for an orientationally-ordered system.  Since the
569 > molecular bodies are perpendicular to the head group dipoles, it is
570 > possible for the director axes for the molecular bodies and the head
571 > groups to be completely decoupled from each other.
572  
573 < The ordering of the tails is essentially opposite to the ordering of
574 < the dipoles on head group. The $P_2$ order parameter decreases with
575 < increasing head size. This indicates the surface is more curved with
576 < larger head groups. When the surface is flat, all tails are pointing
577 < in the same direction; in this case, all tails are parallel to the
578 < normal of the surface,(making this structure remindcent of the
579 < $L_{\beta}$ phase. Increasing the size of the heads, results in
573 > Figure \ref{fig:topView} shows snapshots of bird's-eye views of the
574 > flat ($\sigma_h = 1.20 d$) and rippled ($\sigma_h = 1.35, 1.41 d$)
575 > bilayers.  The directions of the dipoles on the head groups are
576 > represented with two colored half spheres: blue (phosphate) and yellow
577 > (amino).  For flat bilayers, the system exhibits signs of
578 > orientational frustration; some disorder in the dipolar head-to-tail
579 > chains is evident with kinks visible at the edges between differently
580 > ordered domains.  The lipids can also move independently of lipids in
581 > the opposing leaf, so the ordering of the dipoles on one leaf is not
582 > necessarily consistent with the ordering on the other.  These two
583 > factors keep the total dipolar order parameter relatively low for the
584 > flat phases.
585 >
586 > With increasing head group size, the surface becomes corrugated, and
587 > the dipoles cannot move as freely on the surface. Therefore, the
588 > translational freedom of lipids in one layer is dependent upon the
589 > position of the lipids in the other layer.  As a result, the ordering of
590 > the dipoles on head groups in one leaf is correlated with the ordering
591 > in the other leaf.  Furthermore, as the membrane deforms due to the
592 > corrugation, the symmetry of the allowed dipolar ordering on each leaf
593 > is broken. The dipoles then self-assemble in a head-to-tail
594 > configuration, and the dipolar order parameter increases dramatically.
595 > However, the total polarization of the system is still close to zero.
596 > This is strong evidence that the corrugated structure is an
597 > anti-ferroelectric state.  It is also notable that the head-to-tail
598 > arrangement of the dipoles is always observed in a direction
599 > perpendicular to the wave vector for the surface corrugation.  This is
600 > a similar finding to what we observed in our earlier work on the
601 > elastic dipolar membranes.\cite{Sun2007}
602 >
603 > The $P_2$ order parameters (for both the molecular bodies and the head
604 > group dipoles) have been calculated to quantify the ordering in these
605 > phases.  Figure \ref{fig:rP2} shows that the $P_2$ order parameter for
606 > the head-group dipoles increases with increasing head group size. When
607 > the heads of the lipid molecules are small, the membrane is nearly
608 > flat. Since the in-plane packing is essentially a close packing of the
609 > head groups, the head dipoles exhibit frustration in their
610 > orientational ordering.
611 >
612 > The ordering trends for the tails are essentially opposite to the
613 > ordering of the head group dipoles. The tail $P_2$ order parameter
614 > {\it decreases} with increasing head size. This indicates that the
615 > surface is more curved with larger head / tail size ratios. When the
616 > surface is flat, all tails are pointing in the same direction (normal
617 > to the bilayer surface).  This simplified model appears to be
618 > exhibiting a smectic A fluid phase, similar to the real $L_{\beta}$
619 > phase.  We have not observed a smectic C gel phase ($L_{\beta'}$) for
620 > this model system.  Increasing the size of the heads results in
621   rapidly decreasing $P_2$ ordering for the molecular bodies.
622 +
623   \begin{figure}[htb]
624   \centering
625   \includegraphics[width=\linewidth]{rP2}
626 < \caption{The $P_2$ order parameter as a funtion of the ratio of
627 < $\sigma_h$ to $\sigma_0$.\label{fig:rP2}}
626 > \caption{The $P_2$ order parameters for head groups (circles) and
627 > molecular bodies (squares) as a function of the ratio of head group
628 > size ($\sigma_h$) to the width of the molecular bodies ($d$). \label{fig:rP2}}
629   \end{figure}
630  
631 < We studied the effects of the interactions between head groups on the
632 < structure of lipid bilayer by changing the strength of the dipole.
633 < Figure \ref{fig:sP2} shows how the $P_2$ order parameter changes with
634 < increasing strength of the dipole. Generally the dipoles on the head
635 < group are more ordered by increase in the strength of the interaction
636 < between heads and are more disordered by decreasing the interaction
637 < stength. When the interaction between the heads is weak enough, the
638 < bilayer structure does not persist; all lipid molecules are solvated
639 < directly in the water. The critial value of the strength of the dipole
640 < depends on the head size. The perfectly flat surface melts at $5$
641 < $0.03$ debye, the asymmetric rippled surfaces melt at $8$ $0.04$
642 < $0.03$ debye, the symmetric rippled surfaces melt at $10$ $0.04$
643 < debye. The ordering of the tails is the same as the ordering of the
644 < dipoles except for the flat phase. Since the surface is already
645 < perfect flat, the order parameter does not change much until the
646 < strength of the dipole is $15$ debye. However, the order parameter
647 < decreases quickly when the strength of the dipole is further
648 < increased. The head groups of the lipid molecules are brought closer
649 < by stronger interactions between them. For a flat surface, a large
650 < amount of free volume between the head groups is available, but when
651 < the head groups are brought closer, the tails will splay outward,
652 < forming an inverse micelle. When $\sigma_h=1.28\sigma_0$, the $P_2$
653 < order parameter decreases slightly after the strength of the dipole is
654 < increased to $16$ debye. For rippled surfaces, there is less free
655 < volume available between the head groups. Therefore there is little
656 < effect on the structure of the membrane due to increasing dipolar
657 < strength. However, the increase of the $P_2$ order parameter implies
658 < the membranes are flatten by the increase of the strength of the
659 < dipole. Unlike other systems that melt directly when the interaction
660 < is weak enough, for $\sigma_h=1.41\sigma_0$, part of the membrane
661 < melts into itself first. The upper leaf of the bilayer becomes totally
662 < interdigitated with the lower leaf. This is different behavior than
663 < what is exhibited with the interdigitated lines in the rippled phase
664 < where only one interdigitated line connects the two leaves of bilayer.
631 > In addition to varying the size of the head groups, we studied the
632 > effects of the interactions between head groups on the structure of
633 > lipid bilayer by changing the strength of the dipoles.  Figure
634 > \ref{fig:sP2} shows how the $P_2$ order parameter changes with
635 > increasing strength of the dipole.  Generally, the dipoles on the head
636 > groups become more ordered as the strength of the interaction between
637 > heads is increased and become more disordered by decreasing the
638 > interaction strength.  When the interaction between the heads becomes
639 > too weak, the bilayer structure does not persist; all lipid molecules
640 > become dispersed in the solvent (which is non-polar in this
641 > molecular-scale model).  The critical value of the strength of the
642 > dipole depends on the size of the head groups.  The perfectly flat
643 > surface becomes unstable below $5$ Debye, while the  rippled
644 > surfaces melt at $8$ Debye (asymmetric) or $10$ Debye (symmetric).
645 >
646 > The ordering of the tails mirrors the ordering of the dipoles {\it
647 > except for the flat phase}. Since the surface is nearly flat in this
648 > phase, the order parameters are only weakly dependent on dipolar
649 > strength until it reaches $15$ Debye.  Once it reaches this value, the
650 > head group interactions are strong enough to pull the head groups
651 > close to each other and distort the bilayer structure. For a flat
652 > surface, a substantial amount of free volume between the head groups
653 > is normally available.  When the head groups are brought closer by
654 > dipolar interactions, the tails are forced to splay outward, first forming
655 > curved bilayers, and then inverted micelles.
656 >
657 > When $\sigma_h=1.28 d$, the $P_2$ order parameter decreases slightly
658 > when the strength of the dipole is increased above $16$ Debye. For
659 > rippled bilayers, there is less free volume available between the head
660 > groups. Therefore increasing dipolar strength weakly influences the
661 > structure of the membrane.  However, the increase in the body $P_2$
662 > order parameters implies that the membranes are being slightly
663 > flattened due to the effects of increasing head-group attraction.
664 >
665 > A very interesting behavior takes place when the head groups are very
666 > large relative to the molecular bodies ($\sigma_h = 1.41 d$) and the
667 > dipolar strength is relatively weak ($\mu < 9$ Debye). In this case,
668 > the two leaves of the bilayer become totally interdigitated with each
669 > other in large patches of the membrane.   With higher dipolar
670 > strength, the interdigitation is limited to single lines that run
671 > through the bilayer in a direction perpendicular to the ripple wave
672 > vector.
673 >
674   \begin{figure}[htb]
675   \centering
676   \includegraphics[width=\linewidth]{sP2}
677 < \caption{The $P_2$ order parameter as a funtion of the strength of the
678 < dipole.\label{fig:sP2}}
677 > \caption{The $P_2$ order parameters for head group dipoles (a) and
678 > molecular bodies (b) as a function of the strength of the dipoles.
679 > These order parameters are shown for four values of the head group /
680 > molecular width ratio ($\sigma_h / d$). \label{fig:sP2}}
681   \end{figure}
682  
683 < Figure \ref{fig:tP2} shows the dependence of the order parameter on
684 < temperature. The behavior of the $P_2$ order paramter is
685 < straightforward. Systems are more ordered at low temperature, and more
686 < disordered at high temperatures. When the temperature is high enough,
687 < the membranes are instable. For flat surfaces ($\sigma_h=1.20\sigma_0$
688 < and $\sigma_h=1.28\sigma_0$), when the temperature is increased to
689 < $310$, the $P_2$ order parameter increases slightly instead of
690 < decreases like ripple surface. This is an evidence of the frustration
691 < of the dipolar ordering in each leaf of the lipid bilayer, at low
692 < temperature, the systems are locked in a local minimum energy state,
693 < with increase of the temperature, the system can jump out the local
694 < energy well to find the lower energy state which is the longer range
695 < orientational ordering. Like the dipolar ordering of the flat
696 < surfaces, the ordering of the tails of the lipid molecules for ripple
697 < membranes ($\sigma_h=1.35\sigma_0$ and $\sigma_h=1.41\sigma_0$) also
698 < show some nonthermal characteristic. With increase of the temperature,
699 < the $P_2$ order parameter decreases firstly, and increases afterward
700 < when the temperature is greater than $290 K$. The increase of the
701 < $P_2$ order parameter indicates a more ordered structure for the tails
702 < of the lipid molecules which corresponds to a more flat surface. Since
703 < our model lacks the detailed information on lipid tails, we can not
704 < simulate the fluid phase with melted fatty acid chains. Moreover, the
705 < formation of the tilted $L_{\beta'}$ phase also depends on the
540 < organization of fatty groups on tails.
683 > Figure \ref{fig:tP2} shows the dependence of the order parameters on
684 > temperature.  As expected, systems are more ordered at low
685 > temperatures, and more disordered at high temperatures.  All of the
686 > bilayers we studied can become unstable if the temperature becomes
687 > high enough.  The only interesting feature of the temperature
688 > dependence is in the flat surfaces ($\sigma_h=1.20 d$ and
689 > $\sigma_h=1.28 d$).  Here, when the temperature is increased above
690 > $310$K, there is enough jostling of the head groups to allow the
691 > dipolar frustration to resolve into more ordered states.  This results
692 > in a slight increase in the $P_2$ order parameter above this
693 > temperature.
694 >
695 > For the rippled surfaces ($\sigma_h=1.35 d$ and $\sigma_h=1.41 d$),
696 > there is a slightly increased orientational ordering in the molecular
697 > bodies above $290$K.  Since our model lacks the detailed information
698 > about the behavior of the lipid tails, this is the closest the model
699 > can come to depicting the ripple ($P_{\beta'}$) to fluid
700 > ($L_{\alpha}$) phase transition.  What we are observing is a
701 > flattening of the rippled structures made possible by thermal
702 > expansion of the tightly-packed head groups.  The lack of detailed
703 > chain configurations also makes it impossible for this model to depict
704 > the ripple to gel ($L_{\beta'}$) phase transition.
705 >
706   \begin{figure}[htb]
707   \centering
708   \includegraphics[width=\linewidth]{tP2}
709 < \caption{The $P_2$ order parameter as a funtion of
710 < temperature.\label{fig:tP2}}
709 > \caption{The $P_2$ order parameters for head group dipoles (a) and
710 > molecular bodies (b) as a function of temperature.
711 > These order parameters are shown for four values of the head group /
712 > molecular width ratio ($\sigma_h / d$).\label{fig:tP2}}
713   \end{figure}
714  
715 + Fig. \ref{fig:phaseDiagram} shows a phase diagram for the model as a
716 + function of the head group / molecular width ratio ($\sigma_h / d$)
717 + and the strength of the head group dipole moment ($\mu$).  Note that
718 + the specific form of the bilayer phase is governed almost entirely by
719 + the head group / molecular width ratio, while the strength of the
720 + dipolar interactions between the head groups governs the stability of
721 + the bilayer phase.  Weaker dipoles result in unstable bilayer phases,
722 + while extremely strong dipoles can shift the equilibrium to an
723 + inverted micelle phase when the head groups are small.   Temperature
724 + has little effect on the actual bilayer phase observed, although higher
725 + temperatures can cause the unstable region to grow into the higher
726 + dipole region of this diagram.
727 +
728 + \begin{figure}[htb]
729 + \centering
730 + \includegraphics[width=\linewidth]{phaseDiagram}
731 + \caption{Phase diagram for the simple molecular model as a function
732 + of the head group / molecular width ratio ($\sigma_h / d$) and the
733 + strength of the head group dipole moment
734 + ($\mu$).\label{fig:phaseDiagram}}
735 + \end{figure}
736 +
737 +
738 + We have also computed orientational diffusion constants for the head
739 + groups from the relaxation of the second-order Legendre polynomial
740 + correlation function,
741 + \begin{eqnarray}
742 + C_{\ell}(t) & = & \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf
743 + \mu}_{i}(0) \right) \rangle  \\ \\
744 + & \approx & e^{-\ell(\ell + 1) \theta t},
745 + \end{eqnarray}
746 + of the head group dipoles.  In this last line, we have used a simple
747 + ``Debye''-like model for the relaxation of the correlation function,
748 + specifically in the case when $\ell = 2$.   The computed orientational
749 + diffusion constants are given in table \ref{tab:relaxation}.  The
750 + notable feature we observe is that the orientational diffusion
751 + constant for the head group exhibits an order of magnitude decrease
752 + upon entering the rippled phase.  Our orientational correlation times
753 + are substantially in excess of those provided by...
754 +
755 +
756 + \begin{table*}
757 + \begin{minipage}{\linewidth}
758 + \begin{center}
759 + \caption{Rotational diffusion constants for the head groups
760 + ($\theta_h$) and molecular bodies ($\theta_b$) as a function of the
761 + head-to-body width ratio.  The orientational mobility of the head
762 + groups experiences an {\it order of magnitude decrease} upon entering
763 + the rippled phase, which suggests that the rippling is tied to a
764 + freezing out of head group orientational freedom.  Uncertainties in
765 + the last digit are indicated by the values in parentheses.}
766 + \begin{tabular}{lcc}
767 + \hline
768 + $\sigma_h / d$ & $\theta_h (\mu s^{-1})$ & $\theta_b (1/fs)$ \\
769 + \hline
770 + 1.20 & $0.206(1) $ & $0.0175(5) $ \\
771 + 1.28 & $0.179(2) $ & $0.055(2)  $ \\
772 + 1.35 & $0.025(1) $ & $0.195(3)  $ \\
773 + 1.41 & $0.023(1) $ & $0.024(3)  $ \\
774 + \end{tabular}
775 + \label{tab:relaxation}
776 + \end{center}
777 + \end{minipage}
778 + \end{table*}
779 +
780   \section{Discussion}
781   \label{sec:discussion}
782  
783 + Symmetric and asymmetric ripple phases have been observed to form in
784 + our molecular dynamics simulations of a simple molecular-scale lipid
785 + model. The lipid model consists of an dipolar head group and an
786 + ellipsoidal tail.  Within the limits of this model, an explanation for
787 + generalized membrane curvature is a simple mismatch in the size of the
788 + heads with the width of the molecular bodies.  With heads
789 + substantially larger than the bodies of the molecule, this curvature
790 + should be convex nearly everywhere, a requirement which could be
791 + resolved either with micellar or cylindrical phases.
792 +
793 + The persistence of a {\it bilayer} structure therefore requires either
794 + strong attractive forces between the head groups or exclusionary
795 + forces from the solvent phase.  To have a persistent bilayer structure
796 + with the added requirement of convex membrane curvature appears to
797 + result in corrugated structures like the ones pictured in
798 + Fig. \ref{fig:phaseCartoon}.  In each of the sections of these
799 + corrugated phases, the local curvature near a most of the head groups
800 + is convex.  These structures are held together by the extremely strong
801 + and directional interactions between the head groups.
802 +
803 + Dipolar head groups are key for the maintaining the bilayer structures
804 + exhibited by this model.  The dipoles are likely to form head-to-tail
805 + configurations even in flat configurations, but the temperatures are
806 + high enough that vortex defects become prevalent in the flat phase.
807 + The flat phase we observed therefore appears to be substantially above
808 + the Kosterlitz-Thouless transition temperature for a planar system of
809 + dipoles with this set of parameters.  For this reason, it would be
810 + interesting to observe the thermal behavior of the flat phase at
811 + substantially lower temperatures.
812 +
813 + One feature of this model is that an energetically favorable
814 + orientational ordering of the dipoles can be achieved by forming
815 + ripples.  The corrugation of the surface breaks the symmetry of the
816 + plane, making vortex defects somewhat more expensive, and stabilizing
817 + the long range orientational ordering for the dipoles in the head
818 + groups.  Most of the rows of the head-to-tail dipoles are parallel to
819 + each other and the system adopts a bulk anti-ferroelectric state.  We
820 + believe that this is the first time the organization of the head
821 + groups in ripple phases has been addressed.
822 +
823 + Although the size-mismatch between the heads and molecular bodies
824 + appears to be the primary driving force for surface convexity, the
825 + persistence of the bilayer through the use of rippled structures is a
826 + function of the strong, attractive interactions between the heads.
827 + One important prediction we can make using the results from this
828 + simple model is that if the dipole-dipole interaction is the leading
829 + contributor to the head group attractions, the wave vectors for the
830 + ripples should always be found {\it perpendicular} to the dipole
831 + director axis.  This echoes the prediction we made earlier for simple
832 + elastic dipolar membranes, and may suggest experimental designs which
833 + will test whether this is really the case in the phosphatidylcholine
834 + $P_{\beta'}$ phases.  The dipole director axis should also be easily
835 + computable for the all-atom and coarse-grained simulations that have
836 + been published in the literature.\cite{deVries05}
837 +
838 + Although our model is simple, it exhibits some rich and unexpected
839 + behaviors.  It would clearly be a closer approximation to reality if
840 + we allowed bending motions between the dipoles and the molecular
841 + bodies, and if we replaced the rigid ellipsoids with ball-and-chain
842 + tails.  However, the advantages of this simple model (large system
843 + sizes, 50 fs time steps) allow us to rapidly explore the phase diagram
844 + for a wide range of parameters.  Our explanation of this rippling
845 + phenomenon will help us design more accurate molecular models for
846 + corrugated membranes and experiments to test whether or not
847 + dipole-dipole interactions exert an influence on membrane rippling.
848 + \newpage
849   \bibliography{mdripple}
850   \end{document}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines