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 3147 by xsun, Mon Jun 25 21:16:17 2007 UTC vs.
Revision 3265 by xsun, Fri Oct 19 18:40:38 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{}
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 < \begin{abstract}
38 > \maketitle
39  
40 + \begin{abstract}
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 < Our idea for developing a simple and reasonable lipid model to study
62 < the ripple pahse of lipid bilayers is based on two facts: one is that
63 < the most essential feature of lipid molecules is their amphiphilic
64 < structure with polar head groups and non-polar tails. Another fact is
65 < that dominant numbers of lipid molecules are very rigid in ripple
66 < phase which allows the details of the lipid molecules neglectable. In
67 < our model, lipid molecules are represented by rigid bodies made of one
68 < head sphere with a point dipole sitting on it and one ellipsoid tail,
69 < the direction of the dipole is fixed to be perpendicular to the
70 < tail. The breadth and length of tail are $\sigma_0$, $3\sigma_0$. The
71 < diameter of heads varies from $1.20\sigma_0$ to $1.41\sigma_0$.  The
72 < model of the solvent in our simulations is inspired by the idea of
73 < ``DPD'' water. Every four water molecules are reprsented by one
74 < sphere.
61 > \section{Introduction}
62 > \label{sec:Int}
63 > Fully hydrated lipids will aggregate spontaneously to form bilayers
64 > which exhibit a variety of phases depending on their temperatures and
65 > compositions. Among these phases, a periodic rippled phase
66 > ($P_{\beta'}$) appears as an intermediate phase between the gel
67 > ($L_\beta$) and fluid ($L_{\alpha}$) phases for relatively pure
68 > phosphatidylcholine (PC) bilayers.  The ripple phase has attracted
69 > substantial experimental interest over the past 30 years. Most
70 > structural information of the ripple phase has been obtained by the
71 > X-ray diffraction~\cite{Sun96,Katsaras00} and freeze-fracture electron
72 > microscopy (FFEM).~\cite{Copeland80,Meyer96} Recently, Kaasgaard {\it
73 > et al.} used atomic force microscopy (AFM) to observe ripple phase
74 > morphology in bilayers supported on mica.~\cite{Kaasgaard03} The
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} 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
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.
89 + Carlson and Sethna used a packing-competition model (in which head
90 + groups and chains have competing packing energetics) to predict the
91 + formation of a ripple-like phase.  Their model predicted that the
92 + high-curvature portions have lower-chain packing and correspond to
93 + more fluid-like regions.  Goldstein and Leibler used a mean-field
94 + approach with a planar model for {\em inter-lamellar} interactions to
95 + predict rippling in multilamellar phases.~\cite{Goldstein88} McCullough
96 + and Scott proposed that the {\em anisotropy of the nearest-neighbor
97 + interactions} coupled to hydrophobic constraining forces which
98 + restrict height differences between nearest neighbors is the origin of
99 + the ripple phase.~\cite{McCullough90} Lubensky and MacKintosh
100 + introduced a Landau theory for tilt order and curvature of a single
101 + membrane and concluded that {\em coupling of molecular tilt to membrane
102 + curvature} is responsible for the production of
103 + ripples.~\cite{Lubensky93} Misbah, Duplat and Houchmandzadeh proposed
104 + that {\em inter-layer dipolar interactions} can lead to ripple
105 + instabilities.~\cite{Misbah98} Heimburg presented a {\em coexistence
106 + model} for ripple formation in which he postulates that fluid-phase
107 + line defects cause sharp curvature between relatively flat gel-phase
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 head groups
112 + and molecular tilt with respect to the membrane normal vector can
113 + cause bulk rippling.~\cite{Bannerjee02}
114  
115 < Spheres interact each other with Lennard-Jones potential, ellipsoids
116 < interact each other with Gay-Berne potential, dipoles interact each
117 < other with typical dipole potential, spheres interact ellipsoids with
118 < LJ-GB potential. All potentials are truncated at $25$ \AA and shifted
119 < at $22$ \AA.
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 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 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
131 > the fully interdigitated domain.  The symmetric ripple phase was also
132 > observed by Lenz {\it et al.}, and their work supports other claims
133 > that the mismatch between the size of the head group and tail of the
134 > lipid molecules is the driving force for the formation of the ripple
135 > phase. Ayton and Voth have found significant undulations in
136 > zero-surface-tension states of membranes simulated via dissipative
137 > particle dynamics, but their results are consistent with purely
138 > thermal undulations.~\cite{Ayton02}
139  
140 + Although the organization of the tails of lipid molecules are
141 + addressed by these molecular simulations and the packing competition
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 have not been settled.
145  
146 < To make the simulations less expensive and to observe long-time range
147 < behavior of the lipid membranes, all simulaitons were started from two
148 < sepetated monolayers in the vaccum with $x-y$ anisotropic pressure
149 < coupling, length of $z$ axis of the simulations was fixed to prevent
150 < the shrinkage of the simulation boxes due to the free volume outside
151 < of the bilayer, and a constant surface tension was applied to enable
152 < the fluctuation of the surface. Periodic boundaries were used. There
153 < were $480-720$ lipid molecules in simulations according to different
154 < size of the heads. All the simulations were stablized for $100$ ns at
155 < $300$ K. The resulted structures were solvated in the water (about
156 < $6$ DPD water/lipid molecule) as the initial configurations for another
61 < $30$ ns relaxation. All simulations with water were carried out at
62 < constant pressure ($P=1$bar) by $3$D anisotropic coupling, and
63 < constant surface tension ($\gamma=0.015$). Time step was
64 < $50$ fs. Simulations were performed by using OOPSE package.
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 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 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
160 + dynamics simulations to elucidate the role of the head group dipoles
161 + in the formation and morphology of the ripple phase.  We describe our
162 + model and computational methodology in section \ref{sec:method}.
163 + Details on the simulations are presented in section
164 + \ref{sec:experiment}, with results following in section
165 + \ref{sec:results}.  A final discussion of the role of dipolar heads in
166 + the ripple formation can be found in section
167 + \ref{sec:discussion}.
168  
169 < Snap shots show that the membrane is more corrugated with increasing
170 < the size of the head groups. The surface is nearly perfect flat when
69 < $\sigma_h$ is $1.20\sigma_0$. At $1.28\sigma_0$, although the surface
70 < is still flat, the bilayer starts to splay inward, the upper leaf of
71 < the bilayer is connected to the lower leaf with a interdigitated line
72 < defect. Two periodicities with $100$\AA width were observed in the
73 < simulation. This structure is very similiar to OTHER PAPER. The same
74 < structure was also observed when $\sigma_h=1.41\sigma_0$. However, the
75 < surface of the membrane is corrugated, and the periodicity of the
76 < connection between upper and lower leaf membrane is shorter. From the
77 < undulation spectrum of the surface (the exact form is in OUR PREVIOUS
78 < PAPER), the corrugation is non-thermal fluctuation, and we are
79 < confident to identify it as the ripple phase. The width of one ripple
80 < is about $71$ \AA, and amplitude is about $7$ \AA. When
81 < $\sigma_h=1.35\sigma_0$, we observed another corrugated surface with
82 < $79$ \AA width and $10$ \AA amplitude. This structure is different to
83 < the previous rippled surface, there is no connection between upper and
84 < lower leaf of the bilayer. Each leaf of the bilayer is broken to
85 < several curved pieces, the broken position is mounted into the center
86 < of opposite piece in another leaf. Unlike another corrugated surface
87 < in which the upper leaf of the surface is always connected to the
88 < lower leaf from one direction, this ripple of this surface is
89 < isotropic. Therefore, we claim this is a symmetric ripple phase.
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 < The $P_2$ order paramter is calculated to understand the phase
182 < behavior quantatively. $P_2=1$ means a perfect ordered structure, and
183 < $P_2=0$ means a random structure. The method can be found in OUR
184 < PAPER. Fig. shows $P_2$ order paramter of the dipoles on head group
185 < raises with increasing the size of the head group. When head of lipid
186 < molecule is small, the membrane is flat and shows strong two
187 < dimensional characters, dipoles are frustrated on orientational
188 < ordering in this circumstance. Another reason is that the lipids can
189 < move independently in each monolayer, it is not nessasory for the
190 < direction of dipoles on one leaf is consistant to another layer, which
191 < makes total order parameter is relatively low. With increasing the
192 < size of head group, the surface is being more corrugated, dipoles are
193 < not allowed to move freely on the surface, they are
194 < localized. Therefore, the translational freedom of lipids in one layer
195 < is dependent upon the position of lipids in another layer, as a
196 < result, the symmetry of the dipoles on head group in one layer is
197 < consistant to the symmetry in another layer. Furthermore, the membrane
198 < tranlates from a two dimensional system to a three dimensional system
199 < by the corrugation, the symmetry of the ordering for the two
200 < dimensional dipoles on the head group of lipid molecules is broken, on
201 < a distorted lattice, dipoles are ordered on a head to tail energy
113 < state, the order parameter is increased dramaticly. However, the total
114 < polarization of the system is close to zero, which is a strong
115 < evidence it is a antiferroelectric state. The orientation of the
116 < dipole ordering is alway perpendicular to the ripple vector. These
117 < results are consistant to our previous study on similar system. The
118 < ordering of the tails are opposite to the ordering of the dipoles on
119 < head group, the $P_2$ order parameter decreases with increasing the
120 < size of head. This indicates the surface is more curved with larger
121 < head. When surface is flat, all tails are pointing to the same
122 < direction, in this case, all tails are parallal to the normal of the
123 < surface, which shares the same structure with $L_{\beta}$ phase. For the
124 < size of head being $1.28\sigma_0$, the surface starts to splay inward,
125 < however, the surface is still flat, therefore, although the order
126 < parameter is lower, it still indicates a very flat surface. Further
127 < increasing the size of the head, the order parameter drops dramaticly,
128 < the surface is rippled.
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 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
191 > phosphate groups.  The zwitterionic nature of the PC headgroups leads
192 > to abnormally large dipole moments (as high as 20.6 D), and this
193 > strongly polar head group interacts strongly with the solvating water
194 > layers immediately surrounding the membrane.  The hydrophobic tail
195 > consists of fatty acid chains.  In our molecular scale model, lipid
196 > molecules have been reduced to these essential features; the fatty
197 > acid chains are represented by an ellipsoid with a dipolar ball
198 > perched on one end to represent the effects of the charge-separated
199 > head group.  In real PC lipids, the direction of the dipole is
200 > nearly perpendicular to the tail, so we have fixed the direction of
201 > the point dipole rigidly in this orientation.  
202  
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.  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{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},
218 + {\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^{12}
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{equation*}
223  
224 < We studied the effects of interaction between head groups on the
225 < structure of lipid bilayer by changing the strength of the dipole. The
226 < fig. shows the $P_2$ order parameter changing with strength of the
227 < dipole. Generally the dipoles on the head group are more ordered with
228 < increasing the interaction between heads and more disordered with
229 < decreasing the interaction between heads. When the interaction between
230 < heads is weak enough, the bilayer structure is not persisted any more,
231 < all lipid molecules are melted in the water. The critial value of the
232 < strength of the dipole is various for different system. The perfect
233 < flat surface melts at $5$ debye, the asymmetric rippled surfaces melt
234 < at $8$ debye, the symmetric rippled surfaces melt at $10$ debye. This
235 < indicates that the flat phase is the most stable state, the asymmetric
236 < ripple phase is second stalbe state, and the symmetric ripple phase is
237 < the most unstable state. The ordering of the tails is the same as the
238 < ordering of the dipoles except for the flat phase. Since the surface
239 < is already perfect flat, the order parameter does not change much
240 < until the strength of the dipole is $15$ debye. However, the order
241 < parameter decreases quickly when the strength of the dipole is further
242 < increased. The head group of the lipid molecules are brought closer by
243 < strenger interaction between them. For a flat surface, a mount of free
244 < volume between head groups is available, when the head groups are
245 < brought closer, the surface will splay outward to be a inverse
246 < micelle. For rippled surfaces, there is few free volume available on
247 < between the head groups, they can be closer, therefore there are
248 < little effect on the structure of the membrane. Another interesting
249 < fact, unlike other systems melting directly when the interaction is
250 < weak enough, for $\sigma_h$ is $1.41\sigma_0$, part of the membrane
251 < melts into itself first, the upper leaf of the bilayer is totally
252 < interdigitated with the lower leaf, this is different with the
160 < interdigitated lines in rippled phase where only one interdigitated
161 < line connects the two leaves of bilayer.
224 > The range $(\sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
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}}_{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 < Fig. shows the changing of the order parameter with temperature. The
296 < behavior of the $P_2$ orderparamter is straightforword. Systems are
297 < more ordered at low temperature, and more disordered at high
298 < temperature. When the temperature is high enough, the membranes are
299 < discontinuted. The structures are stable during the changing of the
300 < temperature. Since our model lacks the detail information for tails of
301 < lipid molecules, we did not simulate the fluid phase with a melted
302 < fatty chains. Moreover, the formation of the tilted $L_{\beta'}$ phase
303 < also depends on the organization of fatty groups on tails, we did not
304 < observe it either.
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 > d & < & l \\
302 > \epsilon^{r} & < & 1
303 > \end{array}
304 > \end{equation}
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 + \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 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{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{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}$.  
344 +
345 + For the interaction between nonequivalent uniaxial ellipsoids (in this
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
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{Potential parameters used for molecular-scale coarse-grained
365 + lipid simulations}
366 + \begin{tabular}{llccc}
367 + \hline
368 +  & &  Head & Chain & Solvent \\
369 + \hline
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 +
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 26 \AA\ gap between the
403 + molecular bodies of the upper and lower leaves.  The separated
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 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 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 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 + 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=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 + 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{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 / d$ & type of phase & bilayer spacing (\AA) & area per
510 + lipid (\AA$^2$) & $\lambda / d$ & $A / d$\\
511 + \hline
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 / 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=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 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 + 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 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 + 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 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 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 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 + We have computed translational diffusion coefficients for lipid
738 + molecules from the mean square displacement,
739 + \begin{eqnarray}
740 + \langle {|\left({\bf r}_{i}(t) - {\bt r}_{i}(0) \right)|}^2 \rangle \\ \\
741 + & = & 6Dt
742 + \end{eqnarray}
743 + of the lipid bodies. The values of the translational diffusion
744 + coefficient for different head-to-tail size ratio are shown in table
745 + \ref{tab:relaxation}.
746 +
747 + We have also computed orientational diffusion constants for the head
748 + groups from the relaxation of the second-order Legendre polynomial
749 + correlation function,
750 + \begin{eqnarray}
751 + C_{\ell}(t) & = & \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf
752 + \mu}_{i}(0) \right) \rangle  \\ \\
753 + & \approx & e^{-\ell(\ell + 1) \theta t},
754 + \end{eqnarray}
755 + of the head group dipoles.  In this last line, we have used a simple
756 + ``Debye''-like model for the relaxation of the correlation function,
757 + specifically in the case when $\ell = 2$.   The computed orientational
758 + diffusion constants are given in table \ref{tab:relaxation}.  The
759 + notable feature we observe is that the orientational diffusion
760 + constant for the head group exhibits an order of magnitude decrease
761 + upon entering the rippled phase.  Our orientational correlation times
762 + are substantially in excess of those provided by...
763 +
764 +
765 + \begin{table*}
766 + \begin{minipage}{\linewidth}
767 + \begin{center}
768 + \caption{Rotational diffusion constants for the head groups
769 + ($\theta_h$) and molecular bodies ($\theta_b$) as well as the
770 + translational diffusion coefficients for the molecule as a function of
771 + the head-to-body width ratio.  The orientational mobility of the head
772 + groups experiences an {\it order of magnitude decrease} upon entering
773 + the rippled phase, which suggests that the rippling is tied to a
774 + freezing out of head group orientational freedom.  Uncertainties in
775 + the last digit are indicated by the values in parentheses.}
776 + \begin{tabular}{lccc}
777 + \hline
778 + $\sigma_h / d$ & $\theta_h (\mu s^{-1})$ & $\theta_b (1/fs)$ & $D (
779 + \times 10^{-11} m^2 s^{-1} \\
780 + \hline
781 + 1.20 & $0.206(1) $ & $0.0175(5) $ & $0.43(1)$ \\
782 + 1.28 & $0.179(2) $ & $0.055(2)  $ & $5.91(3)$ \\
783 + 1.35 & $0.025(1) $ & $0.195(3)  $ & $3.42(1)$ \\
784 + 1.41 & $0.023(1) $ & $0.024(3)  $ & $7.16(1)$ \\
785 + \end{tabular}
786 + \label{tab:relaxation}
787 + \end{center}
788 + \end{minipage}
789 + \end{table*}
790 +
791 + \section{Discussion}
792 + \label{sec:discussion}
793 +
794 + Symmetric and asymmetric ripple phases have been observed to form in
795 + our molecular dynamics simulations of a simple molecular-scale lipid
796 + model. The lipid model consists of an dipolar head group and an
797 + ellipsoidal tail.  Within the limits of this model, an explanation for
798 + generalized membrane curvature is a simple mismatch in the size of the
799 + heads with the width of the molecular bodies.  With heads
800 + substantially larger than the bodies of the molecule, this curvature
801 + should be convex nearly everywhere, a requirement which could be
802 + resolved either with micellar or cylindrical phases.
803 +
804 + The persistence of a {\it bilayer} structure therefore requires either
805 + strong attractive forces between the head groups or exclusionary
806 + forces from the solvent phase.  To have a persistent bilayer structure
807 + with the added requirement of convex membrane curvature appears to
808 + result in corrugated structures like the ones pictured in
809 + Fig. \ref{fig:phaseCartoon}.  In each of the sections of these
810 + corrugated phases, the local curvature near a most of the head groups
811 + is convex.  These structures are held together by the extremely strong
812 + and directional interactions between the head groups.
813 +
814 + Dipolar head groups are key for the maintaining the bilayer structures
815 + exhibited by this model.  The dipoles are likely to form head-to-tail
816 + configurations even in flat configurations, but the temperatures are
817 + high enough that vortex defects become prevalent in the flat phase.
818 + The flat phase we observed therefore appears to be substantially above
819 + the Kosterlitz-Thouless transition temperature for a planar system of
820 + dipoles with this set of parameters.  For this reason, it would be
821 + interesting to observe the thermal behavior of the flat phase at
822 + substantially lower temperatures.
823 +
824 + One feature of this model is that an energetically favorable
825 + orientational ordering of the dipoles can be achieved by forming
826 + ripples.  The corrugation of the surface breaks the symmetry of the
827 + plane, making vortex defects somewhat more expensive, and stabilizing
828 + the long range orientational ordering for the dipoles in the head
829 + groups.  Most of the rows of the head-to-tail dipoles are parallel to
830 + each other and the system adopts a bulk anti-ferroelectric state.  We
831 + believe that this is the first time the organization of the head
832 + groups in ripple phases has been addressed.
833 +
834 + Although the size-mismatch between the heads and molecular bodies
835 + appears to be the primary driving force for surface convexity, the
836 + persistence of the bilayer through the use of rippled structures is a
837 + function of the strong, attractive interactions between the heads.
838 + One important prediction we can make using the results from this
839 + simple model is that if the dipole-dipole interaction is the leading
840 + contributor to the head group attractions, the wave vectors for the
841 + ripples should always be found {\it perpendicular} to the dipole
842 + director axis.  This echoes the prediction we made earlier for simple
843 + elastic dipolar membranes, and may suggest experimental designs which
844 + will test whether this is really the case in the phosphatidylcholine
845 + $P_{\beta'}$ phases.  The dipole director axis should also be easily
846 + computable for the all-atom and coarse-grained simulations that have
847 + been published in the literature.\cite{deVries05}
848 +
849 + Although our model is simple, it exhibits some rich and unexpected
850 + behaviors.  It would clearly be a closer approximation to reality if
851 + we allowed bending motions between the dipoles and the molecular
852 + bodies, and if we replaced the rigid ellipsoids with ball-and-chain
853 + tails.  However, the advantages of this simple model (large system
854 + sizes, 50 fs time steps) allow us to rapidly explore the phase diagram
855 + for a wide range of parameters.  Our explanation of this rippling
856 + phenomenon will help us design more accurate molecular models for
857 + corrugated membranes and experiments to test whether or not
858 + dipole-dipole interactions exert an influence on membrane rippling.
859 + \newpage
860   \bibliography{mdripple}
861   \end{document}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines