| 1 |
< |
\chapter{\label{chap:mc}MONTE CARLO} |
| 1 |
> |
\chapter{\label{chap:mc}SPONTANEOUS CORRUGATION OF DIPOLAR MEMBRANES} |
| 2 |
> |
|
| 3 |
> |
\section{Introduction} |
| 4 |
> |
\label{mc:sec:Int} |
| 5 |
> |
|
| 6 |
> |
The properties of polymeric membranes are known to depend sensitively |
| 7 |
> |
on the details of the internal interactions between the constituent |
| 8 |
> |
monomers. A flexible membrane will always have a competition between |
| 9 |
> |
the energy of curvature and the in-plane stretching energy and will be |
| 10 |
> |
able to buckle in certain limits of surface tension and |
| 11 |
> |
temperature.\cite{Safran94} The buckling can be non-specific and |
| 12 |
> |
centered at dislocation~\cite{Seung1988} or grain-boundary |
| 13 |
> |
defects,\cite{Carraro1993} or it can be directional and cause long |
| 14 |
> |
``roof-tile'' or tube-like structures to appear in |
| 15 |
> |
partially-polymerized phospholipid vesicles.\cite{Mutz1991} |
| 16 |
> |
|
| 17 |
> |
One would expect that anisotropic local interactions could lead to |
| 18 |
> |
interesting properties of the buckled membrane. We report here on the |
| 19 |
> |
buckling behavior of a membrane composed of harmonically-bound, but |
| 20 |
> |
freely-rotating electrostatic dipoles. The dipoles have strongly |
| 21 |
> |
anisotropic local interactions and the membrane exhibits coupling |
| 22 |
> |
between the buckling and the long-range ordering of the dipoles. |
| 23 |
> |
|
| 24 |
> |
Buckling behavior in liquid crystalline and biological membranes is a |
| 25 |
> |
well-known phenomenon. Relatively pure phosphatidylcholine (PC) |
| 26 |
> |
bilayers form a corrugated or ``rippled'' phase ($P_{\beta'}$) which |
| 27 |
> |
appears as an intermediate phase between the gel ($L_\beta$) and fluid |
| 28 |
> |
($L_{\alpha}$) phases. The $P_{\beta'}$ phase has attracted |
| 29 |
> |
substantial experimental interest over the past 30 |
| 30 |
> |
years,~\cite{Sun96,Katsaras00,Copeland80,Meyer96,Kaasgaard03} and |
| 31 |
> |
there have been a number of theoretical |
| 32 |
> |
approaches~\cite{Marder84,Carlson87,Goldstein88,McCullough90,Lubensky93,Misbah98,Heimburg00,Kubica02,Bannerjee02} |
| 33 |
> |
(and some heroic |
| 34 |
> |
simulations~\cite{Ayton02,Jiang04,Brannigan04a,deVries05,deJoannis06}) |
| 35 |
> |
undertaken to try to explain this phase, but to date, none have looked |
| 36 |
> |
specifically at the contribution of the dipolar character of the lipid |
| 37 |
> |
head groups towards this corrugation. Lipid chain interdigitation |
| 38 |
> |
certainly plays a major role, and the structures of the ripple phase |
| 39 |
> |
are highly ordered. The model we investigate here lacks chain |
| 40 |
> |
interdigitation (as well as the chains themselves!) and will not be |
| 41 |
> |
detailed enough to rule in favor of (or against) any of these |
| 42 |
> |
explanations for the $P_{\beta'}$ phase. |
| 43 |
> |
|
| 44 |
> |
Membranes containing electrostatic dipoles can also exhibit the |
| 45 |
> |
flexoelectric effect,\cite{Todorova2004,Harden2006,Petrov2006} which |
| 46 |
> |
is the ability of mechanical deformations to result in electrostatic |
| 47 |
> |
organization of the membrane. This phenomenon is a curvature-induced |
| 48 |
> |
membrane polarization which can lead to potential differences across a |
| 49 |
> |
membrane. Reverse flexoelectric behavior (in which applied currents |
| 50 |
> |
effect membrane curvature) has also been observed. Explanations of |
| 51 |
> |
the details of these effects have typically utilized membrane |
| 52 |
> |
polarization perpendicular to the face of the |
| 53 |
> |
membrane,\cite{Petrov2006} and the effect has been observed in both |
| 54 |
> |
biological,\cite{Raphael2000} bent-core liquid |
| 55 |
> |
crystalline,\cite{Harden2006} and polymer-dispersed liquid crystalline |
| 56 |
> |
membranes.\cite{Todorova2004} |
| 57 |
> |
|
| 58 |
> |
The problem with using atomistic and even coarse-grained approaches to |
| 59 |
> |
study membrane buckling phenomena is that only a relatively small |
| 60 |
> |
number of periods of the corrugation (i.e. one or two) can be |
| 61 |
> |
realistically simulated given current technology. Also, simulations |
| 62 |
> |
of lipid bilayers are traditionally carried out with periodic boundary |
| 63 |
> |
conditions in two or three dimensions and these have the potential to |
| 64 |
> |
enhance the periodicity of the system at that wavelength. To avoid |
| 65 |
> |
this pitfall, we are using a model which allows us to have |
| 66 |
> |
sufficiently large systems so that we are not causing artificial |
| 67 |
> |
corrugation through the use of periodic boundary conditions. |
| 68 |
> |
|
| 69 |
> |
The simplest dipolar membrane is one in which the dipoles are located |
| 70 |
> |
on fixed lattice sites. Ferroelectric states (with long-range dipolar |
| 71 |
> |
order) can be observed in dipolar systems with non-triangular |
| 72 |
> |
packings. However, {\em triangularly}-packed 2-D dipolar systems are |
| 73 |
> |
inherently frustrated and one would expect a dipolar-disordered phase |
| 74 |
> |
to be the lowest free energy |
| 75 |
> |
configuration.\cite{Toulouse1977,Marland1979} Dipolar lattices already |
| 76 |
> |
have rich phase behavior, but in order to allow the membrane to |
| 77 |
> |
buckle, a single degree of freedom (translation normal to the membrane |
| 78 |
> |
face) must be added to each of the dipoles. It would also be possible |
| 79 |
> |
to allow complete translational freedom. This approach |
| 80 |
> |
is similar in character to a number of elastic Ising models that have |
| 81 |
> |
been developed to explain interesting mechanical properties in |
| 82 |
> |
magnetic alloys.\cite{Renard1966,Zhu2005,Zhu2006,Jiang2006} |
| 83 |
> |
|
| 84 |
> |
What we present here is an attempt to find the simplest dipolar model |
| 85 |
> |
which will exhibit buckling behavior. We are using a modified XYZ |
| 86 |
> |
lattice model; details of the model can be found in section |
| 87 |
> |
\ref{mc:sec:model}, results of Monte Carlo simulations using this model |
| 88 |
> |
are presented in section |
| 89 |
> |
\ref{mc:sec:results}, and section \ref{mc:sec:discussion} contains our conclusions. |
| 90 |
> |
|
| 91 |
> |
\section{2-D Dipolar Membrane} |
| 92 |
> |
\label{mc:sec:model} |
| 93 |
> |
|
| 94 |
> |
The point of developing this model was to arrive at the simplest |
| 95 |
> |
possible theoretical model which could exhibit spontaneous corrugation |
| 96 |
> |
of a two-dimensional dipolar medium. Since molecules in polymerized |
| 97 |
> |
membranes and in the $P_{\beta'}$ ripple phase have limited |
| 98 |
> |
translational freedom, we have chosen a lattice to support the dipoles |
| 99 |
> |
in the x-y plane. The lattice may be either triangular (lattice |
| 100 |
> |
constants $a/b = |
| 101 |
> |
\sqrt{3}$) or distorted. However, each dipole has 3 degrees of |
| 102 |
> |
freedom. They may move freely {\em out} of the x-y plane (along the |
| 103 |
> |
$z$ axis), and they have complete orientational freedom ($0 <= \theta |
| 104 |
> |
<= \pi$, $0 <= \phi < 2 |
| 105 |
> |
\pi$). This is essentially a modified X-Y-Z model with translational |
| 106 |
> |
freedom along the z-axis. |
| 107 |
> |
|
| 108 |
> |
The potential energy of the system, |
| 109 |
> |
\begin{equation} |
| 110 |
> |
\begin{split} |
| 111 |
> |
V = \sum_i &\left( \sum_{j>i} \frac{|\mu|^2}{4\pi \epsilon_0 r_{ij}^3} \left[ |
| 112 |
> |
{\mathbf{\hat u}_i} \cdot {\mathbf{\hat u}_j} - |
| 113 |
> |
3({\mathbf{\hat u}_i} \cdot {\mathbf{\hat |
| 114 |
> |
r}_{ij}})({\mathbf{\hat u}_j} \cdot {\mathbf{\hat r}_{ij}})\right] \right. \\ |
| 115 |
> |
& \left. + \sum_{j \in NN_i}^6 \frac{k_r}{2}\left( |
| 116 |
> |
r_{ij}-\sigma \right)^2 \right) |
| 117 |
> |
\end{split} |
| 118 |
> |
\label{mceq:pot} |
| 119 |
> |
\end{equation} |
| 120 |
> |
|
| 121 |
> |
In this potential, $\mathbf{\hat u}_i$ is the unit vector pointing |
| 122 |
> |
along dipole $i$ and $\mathbf{\hat r}_{ij}$ is the unit vector |
| 123 |
> |
pointing along the inter-dipole vector $\mathbf{r}_{ij}$. The entire |
| 124 |
> |
potential is governed by three parameters, the dipolar strength |
| 125 |
> |
($\mu$), the harmonic spring constant ($k_r$) and the preferred |
| 126 |
> |
intermolecular spacing ($\sigma$). In practice, we set the value of |
| 127 |
> |
$\sigma$ to the average inter-molecular spacing from the planar |
| 128 |
> |
lattice, yielding a potential model that has only two parameters for a |
| 129 |
> |
particular choice of lattice constants $a$ (along the $x$-axis) and |
| 130 |
> |
$b$ (along the $y$-axis). We also define a set of reduced parameters |
| 131 |
> |
based on the length scale ($\sigma$) and the energy of the harmonic |
| 132 |
> |
potential at a deformation of 2 $\sigma$ ($\epsilon = k_r \sigma^2 / |
| 133 |
> |
2$). Using these two constants, we perform our calculations using |
| 134 |
> |
reduced distances, ($r^{*} = r / \sigma$), temperatures ($T^{*} = 2 |
| 135 |
> |
k_B T / k_r \sigma^2$), densities ($\rho^{*} = N \sigma^2 / L_x L_y$), |
| 136 |
> |
and dipole moments ($\mu^{*} = \mu / \sqrt{4 \pi \epsilon_0 \sigma^5 |
| 137 |
> |
k_r / 2}$). It should be noted that the density ($\rho^{*}$) depends |
| 138 |
> |
only on the mean particle spacing in the $x-y$ plane; the lattice is |
| 139 |
> |
fully populated. |
| 140 |
> |
|
| 141 |
> |
To investigate the phase behavior of this model, we have performed a |
| 142 |
> |
series of Metropolis Monte Carlo simulations of moderately-sized (34.3 |
| 143 |
> |
$\sigma$ on a side) patches of membrane hosted on both triangular |
| 144 |
> |
($\gamma = a/b = \sqrt{3}$) and distorted ($\gamma \neq \sqrt{3}$) |
| 145 |
> |
lattices. The linear extent of one edge of the monolayer was $20 a$ |
| 146 |
> |
and the system was kept roughly square. The average distance that |
| 147 |
> |
coplanar dipoles were positioned from their six nearest neighbors was |
| 148 |
> |
1 $\sigma$ (on both triangular and distorted lattices). Typical |
| 149 |
> |
system sizes were 1360 dipoles for the triangular lattices and |
| 150 |
> |
840-2800 dipoles for the distorted lattices. Two-dimensional periodic |
| 151 |
> |
boundary conditions were used, and the cutoff for the dipole-dipole |
| 152 |
> |
interaction was set to 4.3 $\sigma$. This cutoff is roughly 2.5 times |
| 153 |
> |
the typical real-space electrostatic cutoff for molecular systems. |
| 154 |
> |
Since dipole-dipole interactions decay rapidly with distance, and |
| 155 |
> |
since the intrinsic three-dimensional periodicity of the Ewald sum can |
| 156 |
> |
give artifacts in 2-d systems, we have chosen not to use it in these |
| 157 |
> |
calculations. Although the Ewald sum has been reformulated to handle |
| 158 |
> |
2-D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89} these |
| 159 |
> |
methods are computationally expensive,\cite{Spohr97,Yeh99} and are not |
| 160 |
> |
necessary in this case. All parameters ($T^{*}$, $\mu^{*}$, and |
| 161 |
> |
$\gamma$) were varied systematically to study the effects of these |
| 162 |
> |
parameters on the formation of ripple-like phases. |
| 163 |
> |
|
| 164 |
> |
\section{Results and Analysis} |
| 165 |
> |
\label{mc:sec:results} |
| 166 |
> |
|
| 167 |
> |
\subsection{Dipolar Ordering and Coexistence Temperatures} |
| 168 |
> |
The principal method for observing the orientational ordering |
| 169 |
> |
transition in dipolar or liquid crystalline systems is the $P_2$ order |
| 170 |
> |
parameter (defined as $1.5 \times \lambda_{max}$, where |
| 171 |
> |
$\lambda_{max}$ is the largest eigenvalue of the matrix, |
| 172 |
> |
\begin{equation} |
| 173 |
> |
{\mathsf{S}} = \frac{1}{N} \sum_i \left( |
| 174 |
> |
\begin{array}{ccc} |
| 175 |
> |
u^{x}_i u^{x}_i-\frac{1}{3} & u^{x}_i u^{y}_i & u^{x}_i u^{z}_i \\ |
| 176 |
> |
u^{y}_i u^{x}_i & u^{y}_i u^{y}_i -\frac{1}{3} & u^{y}_i u^{z}_i \\ |
| 177 |
> |
u^{z}_i u^{x}_i & u^{z}_i u^{y}_i & u^{z}_i u^{z}_i -\frac{1}{3} |
| 178 |
> |
\end{array} \right). |
| 179 |
> |
\label{mceq:opmatrix} |
| 180 |
> |
\end{equation} |
| 181 |
> |
Here $u^{\alpha}_i$ is the $\alpha=x,y,z$ component of the unit vector |
| 182 |
> |
for dipole $i$. $P_2$ will be $1.0$ for a perfectly-ordered system |
| 183 |
> |
and near $0$ for a randomized system. Note that this order parameter |
| 184 |
> |
is {\em not} equal to the polarization of the system. For example, |
| 185 |
> |
the polarization of the perfect anti-ferroelectric system is $0$, but |
| 186 |
> |
$P_2$ for an anti-ferroelectric system is $1$. The eigenvector of |
| 187 |
> |
$\mathsf{S}$ corresponding to the largest eigenvalue is familiar as |
| 188 |
> |
the director axis, which can be used to determine a privileged dipolar |
| 189 |
> |
axis for dipole-ordered systems. The top panel in Fig. \ref{mcfig:phase} |
| 190 |
> |
shows the values of $P_2$ as a function of temperature for both |
| 191 |
> |
triangular ($\gamma = 1.732$) and distorted ($\gamma=1.875$) |
| 192 |
> |
lattices. |
| 193 |
> |
|
| 194 |
> |
\begin{figure} |
| 195 |
> |
\includegraphics[width=\linewidth]{./figures/mcPhase.pdf} |
| 196 |
> |
\caption{\label{mcfig:phase} Top panel: The $P_2$ dipolar order parameter as |
| 197 |
> |
a function of temperature for both triangular ($\gamma = 1.732$) and |
| 198 |
> |
distorted ($\gamma = 1.875$) lattices. Bottom Panel: The phase |
| 199 |
> |
diagram for the dipolar membrane model. The line denotes the division |
| 200 |
> |
between the dipolar ordered (anti-ferroelectric) and disordered phases. |
| 201 |
> |
An enlarged view near the triangular lattice is shown inset.} |
| 202 |
> |
\end{figure} |
| 203 |
> |
|
| 204 |
> |
There is a clear order-disorder transition in evidence from this data. |
| 205 |
> |
Both the triangular and distorted lattices have dipolar-ordered |
| 206 |
> |
low-temperature phases, and orientationally-disordered high |
| 207 |
> |
temperature phases. The coexistence temperature for the triangular |
| 208 |
> |
lattice is significantly lower than for the distorted lattices, and |
| 209 |
> |
the bulk polarization is approximately $0$ for both dipolar ordered |
| 210 |
> |
and disordered phases. This gives strong evidence that the dipolar |
| 211 |
> |
ordered phase is anti-ferroelectric. We have verified that this |
| 212 |
> |
dipolar ordering transition is not a function of system size by |
| 213 |
> |
performing identical calculations with systems twice as large. The |
| 214 |
> |
transition is equally smooth at all system sizes that were studied. |
| 215 |
> |
Additionally, we have repeated the Monte Carlo simulations over a wide |
| 216 |
> |
range of lattice ratios ($\gamma$) to generate a dipolar |
| 217 |
> |
order/disorder phase diagram. The bottom panel in Fig. \ref{mcfig:phase} |
| 218 |
> |
shows that the triangular lattice is a low-temperature cusp in the |
| 219 |
> |
$T^{*}-\gamma$ phase diagram. |
| 220 |
> |
|
| 221 |
> |
This phase diagram is remarkable in that it shows an |
| 222 |
> |
anti-ferroelectric phase near $\gamma=1.732$ where one would expect |
| 223 |
> |
lattice frustration to result in disordered phases at all |
| 224 |
> |
temperatures. Observations of the configurations in this phase show |
| 225 |
> |
clearly that the system has accomplished dipolar ordering by forming |
| 226 |
> |
large ripple-like structures. We have observed anti-ferroelectric |
| 227 |
> |
ordering in all three of the equivalent directions on the triangular |
| 228 |
> |
lattice, and the dipoles have been observed to organize perpendicular |
| 229 |
> |
to the membrane normal (in the plane of the membrane). It is |
| 230 |
> |
particularly interesting to note that the ripple-like structures have |
| 231 |
> |
also been observed to propagate in the three equivalent directions on |
| 232 |
> |
the lattice, but the {\em direction of ripple propagation is always |
| 233 |
> |
perpendicular to the dipole director axis}. A snapshot of a typical |
| 234 |
> |
anti-ferroelectric rippled structure is shown in |
| 235 |
> |
Fig. \ref{mcfig:snapshot}. |
| 236 |
> |
|
| 237 |
> |
\begin{figure} |
| 238 |
> |
\includegraphics[width=\linewidth]{./figures/mcSnapshot.pdf} |
| 239 |
> |
\caption{\label{mcfig:snapshot} Top and Side views of a representative |
| 240 |
> |
configuration for the dipolar ordered phase supported on the |
| 241 |
> |
triangular lattice. Note the anti-ferroelectric ordering and the long |
| 242 |
> |
wavelength buckling of the membrane. Dipolar ordering has been |
| 243 |
> |
observed in all three equivalent directions on the triangular lattice, |
| 244 |
> |
and the ripple direction is always perpendicular to the director axis |
| 245 |
> |
for the dipoles.} |
| 246 |
> |
\end{figure} |
| 247 |
> |
|
| 248 |
> |
Although the snapshot in Fig. \ref{mcfig:snapshot} gives the appearance |
| 249 |
> |
of three-row stair-like structures, these appear to be transient. On |
| 250 |
> |
average, the corrugation of the membrane is a relatively smooth, |
| 251 |
> |
long-wavelength phenomenon, with occasional steep drops between |
| 252 |
> |
adjacent lines of anti-aligned dipoles. |
| 253 |
> |
|
| 254 |
> |
The height-dipole correlation function ($C_{\textrm{hd}}(r, \cos |
| 255 |
> |
\theta)$) makes the connection between dipolar ordering and the wave |
| 256 |
> |
vector of the ripple even more explicit. $C_{\textrm{hd}}(r, \cos |
| 257 |
> |
\theta)$ is an angle-dependent pair distribution function. The angle |
| 258 |
> |
($\theta$) is the angle between the intermolecular vector |
| 259 |
> |
$\vec{r}_{ij}$ and direction of dipole $i$, |
| 260 |
> |
\begin{equation} |
| 261 |
> |
C_{\textrm{hd}} = \frac{\langle \frac{1}{n(r)} \sum_{i}\sum_{j>i} |
| 262 |
> |
h_i \cdot h_j \delta(r - r_{ij}) \delta(\cos \theta_{ij} - |
| 263 |
> |
\cos \theta)\rangle} {\langle h^2 \rangle} |
| 264 |
> |
\end{equation} |
| 265 |
> |
where $\cos \theta_{ij} = \hat{\mu}_{i} \cdot \hat{r}_{ij}$ and |
| 266 |
> |
$\hat{r}_{ij} = \vec{r}_{ij} / r_{ij}$. $n(r)$ is the number of |
| 267 |
> |
dipoles found in a cylindrical shell between $r$ and $r+\delta r$ of |
| 268 |
> |
the central particle. Fig. \ref{mcfig:CrossCorrelation} shows contours |
| 269 |
> |
of this correlation function for both anti-ferroelectric, rippled |
| 270 |
> |
membranes as well as for the dipole-disordered portion of the phase |
| 271 |
> |
diagram. |
| 272 |
> |
|
| 273 |
> |
\begin{figure} |
| 274 |
> |
\includegraphics[width=\linewidth]{./figures/mcHdc.pdf} |
| 275 |
> |
\caption{\label{mcfig:CrossCorrelation} Contours of the height-dipole |
| 276 |
> |
correlation function as a function of the dot product between the |
| 277 |
> |
dipole ($\hat{\mu}$) and inter-dipole separation vector ($\hat{r}$) |
| 278 |
> |
and the distance ($r$) between the dipoles. Perfect height |
| 279 |
> |
correlation (contours approaching 1) are present in the ordered phase |
| 280 |
> |
when the two dipoles are in the same head-to-tail line. |
| 281 |
> |
Anti-correlation (contours below 0) is only seen when the inter-dipole |
| 282 |
> |
vector is perpendicular to the dipoles. In the dipole-disordered |
| 283 |
> |
portion of the phase diagram, there is only weak correlation in the |
| 284 |
> |
dipole direction and this correlation decays rapidly to zero for |
| 285 |
> |
intermolecular vectors that are not dipole-aligned.} |
| 286 |
> |
\end{figure} |
| 287 |
> |
|
| 288 |
> |
The height-dipole correlation function gives a map of how the topology |
| 289 |
> |
of the membrane surface varies with angular deviation around a given |
| 290 |
> |
dipole. The upper panel of Fig. \ref{mcfig:CrossCorrelation} shows that |
| 291 |
> |
in the anti-ferroelectric phase, the dipole heights are strongly |
| 292 |
> |
correlated for dipoles in head-to-tail arrangements, and this |
| 293 |
> |
correlation persists for very long distances (up to 15 $\sigma$). For |
| 294 |
> |
portions of the membrane located perpendicular to a given dipole, the |
| 295 |
> |
membrane height becomes anti-correlated at distances of 10 $\sigma$. |
| 296 |
> |
The correlation function is relatively smooth; there are no steep |
| 297 |
> |
jumps or steps, so the stair-like structures in |
| 298 |
> |
Fig. \ref{mcfig:snapshot} are indeed transient and disappear when |
| 299 |
> |
averaged over many configurations. In the dipole-disordered phase, |
| 300 |
> |
the height-dipole correlation function is relatively flat (and hovers |
| 301 |
> |
near zero). The only significant height correlations are for axial |
| 302 |
> |
dipoles at very short distances ($r \approx |
| 303 |
> |
\sigma$). |
| 304 |
> |
|
| 305 |
> |
\subsection{Discriminating Ripples from Thermal Undulations} |
| 306 |
> |
|
| 307 |
> |
In order to be sure that the structures we have observed are actually |
| 308 |
> |
a rippled phase and not simply thermal undulations, we have computed |
| 309 |
> |
the undulation spectrum, |
| 310 |
> |
\begin{equation} |
| 311 |
> |
h(\vec{q}) = A^{-1/2} \int d\vec{r} |
| 312 |
> |
h(\vec{r}) e^{-i \vec{q}\cdot\vec{r}} |
| 313 |
> |
\end{equation} |
| 314 |
> |
where $h(\vec{r})$ is the height of the membrane at location $\vec{r} |
| 315 |
> |
= (x,y)$.~\cite{Safran94,Seifert97} In simple (and more complicated) |
| 316 |
> |
elastic continuum models, it can shown that in the $NVT$ ensemble, the |
| 317 |
> |
absolute value of the undulation spectrum can be written, |
| 318 |
> |
\begin{equation} |
| 319 |
> |
\langle | h(q) |^2 \rangle_{NVT} = \frac{k_B T}{k_c q^4 + |
| 320 |
> |
\gamma q^2}, |
| 321 |
> |
\label{mceq:fit} |
| 322 |
> |
\end{equation} |
| 323 |
> |
where $k_c$ is the bending modulus for the membrane, and $\gamma$ is |
| 324 |
> |
the mechanical surface tension.~\cite{Safran94} The systems studied in |
| 325 |
> |
this paper have essentially zero bending moduli ($k_c$) and relatively |
| 326 |
> |
large mechanical surface tensions ($\gamma$), so a much simpler form |
| 327 |
> |
can be written, |
| 328 |
> |
\begin{equation} |
| 329 |
> |
\langle | h(q) |^2 \rangle_{NVT} = \frac{k_B T}{\gamma q^2}. |
| 330 |
> |
\label{mceq:fit2} |
| 331 |
> |
\end{equation} |
| 332 |
> |
|
| 333 |
> |
The undulation spectrum is computed by superimposing a rectangular |
| 334 |
> |
grid on top of the membrane, and by assigning height ($h(\vec{r})$) |
| 335 |
> |
values to the grid from the average of all dipoles that fall within a |
| 336 |
> |
given $\vec{r}+d\vec{r}$ grid area. Empty grid pixels are assigned |
| 337 |
> |
height values by interpolation from the nearest neighbor pixels. A |
| 338 |
> |
standard 2-d Fourier transform is then used to obtain $\langle | |
| 339 |
> |
h(q)|^2 \rangle$. Alternatively, since the dipoles sit on a Bravais |
| 340 |
> |
lattice, one could use the heights of the lattice points themselves as |
| 341 |
> |
the grid for the Fourier transform (without interpolating to a square |
| 342 |
> |
grid). However, if lateral translational freedom is added to this |
| 343 |
> |
model (a likely extension), an interpolated grid method for computing |
| 344 |
> |
undulation spectra will be required. |
| 345 |
> |
|
| 346 |
> |
As mentioned above, the best fits to our undulation spectra are |
| 347 |
> |
obtained by setting the value of $k_c$ to 0. In Fig. \ref{mcfig:fit} we |
| 348 |
> |
show typical undulation spectra for two different regions of the phase |
| 349 |
> |
diagram along with their fits from the Landau free energy approach |
| 350 |
> |
(Eq. \ref{mceq:fit2}). In the high-temperature disordered phase, the |
| 351 |
> |
Landau fits can be nearly perfect, and from these fits we can estimate |
| 352 |
> |
the tension in the surface. In reduced units, typical values of |
| 353 |
> |
$\gamma^{*} = \gamma / \epsilon = 2500$ are obtained for the |
| 354 |
> |
disordered phase ($\gamma^{*} = 2551.7$ in the top panel of |
| 355 |
> |
Fig. \ref{mcfig:fit}). |
| 356 |
> |
|
| 357 |
> |
Typical values of $\gamma^{*}$ in the dipolar-ordered phase are much |
| 358 |
> |
higher than in the dipolar-disordered phase ($\gamma^{*} = 73,538$ in |
| 359 |
> |
the lower panel of Fig. \ref{mcfig:fit}). For the dipolar-ordered |
| 360 |
> |
triangular lattice near the coexistence temperature, we also observe |
| 361 |
> |
long wavelength undulations that are far outliers to the fits. That |
| 362 |
> |
is, the Landau free energy fits are well within error bars for most of |
| 363 |
> |
the other points, but can be off by {\em orders of magnitude} for a |
| 364 |
> |
few low frequency components. |
| 365 |
> |
|
| 366 |
> |
We interpret these outliers as evidence that these low frequency modes |
| 367 |
> |
are {\em non-thermal undulations}. We take this as evidence that we |
| 368 |
> |
are actually seeing a rippled phase developing in this model system. |
| 369 |
> |
|
| 370 |
> |
\begin{figure} |
| 371 |
> |
\includegraphics[width=\linewidth]{./figures/mcLogFit.pdf} |
| 372 |
> |
\caption{\label{mcfig:fit} Evidence that the observed ripples are {\em |
| 373 |
> |
not} thermal undulations is obtained from the 2-d Fourier transform |
| 374 |
> |
$\langle |h^{*}(\vec{q})|^2 \rangle$ of the height profile ($\langle |
| 375 |
> |
h^{*}(x,y) \rangle$). Rippled samples show low-wavelength peaks that |
| 376 |
> |
are outliers on the Landau free energy fits by an order of magnitude. |
| 377 |
> |
Samples exhibiting only thermal undulations fit Eq. \ref{mceq:fit} |
| 378 |
> |
remarkably well.} |
| 379 |
> |
\end{figure} |
| 380 |
> |
|
| 381 |
> |
\subsection{Effects of Potential Parameters on Amplitude and Wavelength} |
| 382 |
> |
|
| 383 |
> |
We have used two different methods to estimate the amplitude and |
| 384 |
> |
periodicity of the ripples. The first method requires projection of |
| 385 |
> |
the ripples onto a one dimensional rippling axis. Since the rippling |
| 386 |
> |
is always perpendicular to the dipole director axis, we can define a |
| 387 |
> |
ripple vector as follows. The largest eigenvector, $s_1$, of the |
| 388 |
> |
$\mathsf{S}$ matrix in Eq. \ref{mceq:opmatrix} is projected onto a |
| 389 |
> |
planar director axis, |
| 390 |
> |
\begin{equation} |
| 391 |
> |
\vec{d} = \left(\begin{array}{c} |
| 392 |
> |
\vec{s}_1 \cdot \hat{i} \\ |
| 393 |
> |
\vec{s}_1 \cdot \hat{j} \\ |
| 394 |
> |
0 |
| 395 |
> |
\end{array} \right). |
| 396 |
> |
\end{equation} |
| 397 |
> |
($\hat{i}$, $\hat{j}$ and $\hat{k}$ are unit vectors along the $x$, |
| 398 |
> |
$y$, and $z$ axes, respectively.) The rippling axis is in the plane of |
| 399 |
> |
the membrane and is perpendicular to the planar director axis, |
| 400 |
> |
\begin{equation} |
| 401 |
> |
\vec{q}_{\mathrm{rip}} = \vec{d} \times \hat{k} |
| 402 |
> |
\end{equation} |
| 403 |
> |
We can then find the height profile of the membrane along the ripple |
| 404 |
> |
axis by projecting heights of the dipoles to obtain a one-dimensional |
| 405 |
> |
height profile, $h(q_{\mathrm{rip}})$. Ripple wavelengths can be |
| 406 |
> |
estimated from the largest non-thermal low-frequency component in the |
| 407 |
> |
Fourier transform of $h(q_{\mathrm{rip}})$. Amplitudes can be |
| 408 |
> |
estimated by measuring peak-to-trough distances in |
| 409 |
> |
$h(q_{\mathrm{rip}})$ itself. |
| 410 |
> |
|
| 411 |
> |
A second, more accurate, and simpler method for estimating ripple |
| 412 |
> |
shape is to extract the wavelength and height information directly |
| 413 |
> |
from the largest non-thermal peak in the undulation spectrum. For |
| 414 |
> |
large-amplitude ripples, the two methods give similar results. The |
| 415 |
> |
one-dimensional projection method is more prone to noise (particularly |
| 416 |
> |
in the amplitude estimates for the distorted lattices). We report |
| 417 |
> |
amplitudes and wavelengths taken directly from the undulation spectrum |
| 418 |
> |
below. |
| 419 |
> |
|
| 420 |
> |
In the triangular lattice ($\gamma = \sqrt{3}$), the rippling is |
| 421 |
> |
observed for temperatures ($T^{*}$) from $61-122$. The wavelength of |
| 422 |
> |
the ripples is remarkably stable at 21.4~$\sigma$ for all but the |
| 423 |
> |
temperatures closest to the order-disorder transition. At $T^{*} = |
| 424 |
> |
122$, the wavelength drops to 17.1~$\sigma$. |
| 425 |
> |
|
| 426 |
> |
The dependence of the amplitude on temperature is shown in the top |
| 427 |
> |
panel of Fig. \ref{mcfig:Amplitude}. The rippled structures shrink |
| 428 |
> |
smoothly as the temperature rises towards the order-disorder |
| 429 |
> |
transition. The wavelengths and amplitudes we observe are |
| 430 |
> |
surprisingly close to the $\Lambda / 2$ phase observed by Kaasgaard |
| 431 |
> |
{\it et al.} in their work on PC-based lipids.\cite{Kaasgaard03} |
| 432 |
> |
However, this is coincidental agreement based on a choice of 7~\AA~as |
| 433 |
> |
the mean spacing between lipids. |
| 434 |
> |
|
| 435 |
> |
\begin{figure} |
| 436 |
> |
\includegraphics[width=\linewidth]{./figures/mcProperties_sq.pdf} |
| 437 |
> |
\caption{\label{mcfig:Amplitude} a) The amplitude $A^{*}$ of the ripples |
| 438 |
> |
vs. temperature for a triangular lattice. b) The amplitude $A^{*}$ of |
| 439 |
> |
the ripples vs. dipole strength ($\mu^{*}$) for both the triangular |
| 440 |
> |
lattice (circles) and distorted lattice (squares). The reduced |
| 441 |
> |
temperatures were kept fixed at $T^{*} = 94$ for the triangular |
| 442 |
> |
lattice and $T^{*} = 106$ for the distorted lattice (approximately 2/3 |
| 443 |
> |
of the order-disorder transition temperature for each lattice).} |
| 444 |
> |
\end{figure} |
| 445 |
> |
|
| 446 |
> |
The ripples can be made to disappear by increasing the internal |
| 447 |
> |
elastic tension (i.e. by increasing $k_r$ or equivalently, reducing |
| 448 |
> |
the dipole moment). The amplitude of the ripples depends critically |
| 449 |
> |
on the strength of the dipole moments ($\mu^{*}$) in Eq. \ref{mceq:pot}. |
| 450 |
> |
If the dipoles are weakened substantially (below $\mu^{*}$ = 20) at a |
| 451 |
> |
fixed temperature of 94, the membrane loses dipolar ordering |
| 452 |
> |
and the ripple structures. The ripples reach a peak amplitude of |
| 453 |
> |
3.7~$\sigma$ at a dipolar strength of 25. We show the dependence |
| 454 |
> |
of ripple amplitude on the dipolar strength in |
| 455 |
> |
Fig. \ref{mcfig:Amplitude}. |
| 456 |
> |
|
| 457 |
> |
\subsection{Distorted lattices} |
| 458 |
> |
|
| 459 |
> |
We have also investigated the effect of the lattice geometry by |
| 460 |
> |
changing the ratio of lattice constants ($\gamma$) while keeping the |
| 461 |
> |
average nearest-neighbor spacing constant. The anti-ferroelectric state |
| 462 |
> |
is accessible for all $\gamma$ values we have used, although the |
| 463 |
> |
distorted triangular lattices prefer a particular director axis due to |
| 464 |
> |
the anisotropy of the lattice. |
| 465 |
> |
|
| 466 |
> |
Our observation of rippling behavior was not limited to the triangular |
| 467 |
> |
lattices. In distorted lattices the anti-ferroelectric phase can |
| 468 |
> |
develop nearly instantaneously in the Monte Carlo simulations, and |
| 469 |
> |
these dipolar-ordered phases tend to be remarkably flat. Whenever |
| 470 |
> |
rippling has been observed in these distorted lattices |
| 471 |
> |
(e.g. $\gamma = 1.875$), we see relatively short ripple wavelengths |
| 472 |
> |
(14 $\sigma$) and amplitudes of 2.4~$\sigma$. These ripples are |
| 473 |
> |
weakly dependent on dipolar strength (see Fig. \ref{mcfig:Amplitude}), |
| 474 |
> |
although below a dipolar strength of $\mu^{*} = 20$, the membrane |
| 475 |
> |
loses dipolar ordering and displays only thermal undulations. |
| 476 |
> |
|
| 477 |
> |
The ripple phase does {\em not} appear at all values of $\gamma$. We |
| 478 |
> |
have only observed non-thermal undulations in the range $1.625 < |
| 479 |
> |
\gamma < 1.875$. Outside this range, the order-disorder transition in |
| 480 |
> |
the dipoles remains, but the ordered dipolar phase has only thermal |
| 481 |
> |
undulations. This is one of our strongest pieces of evidence that |
| 482 |
> |
rippling is a symmetry-breaking phenomenon for triangular and |
| 483 |
> |
nearly-triangular lattices. |
| 484 |
> |
|
| 485 |
> |
\subsection{Effects of System Size} |
| 486 |
> |
To evaluate the effect of finite system size, we have performed a |
| 487 |
> |
series of simulations on the triangular lattice at a reduced |
| 488 |
> |
temperature of 122, which is just below the order-disorder transition |
| 489 |
> |
temperature ($T^{*} = 139$). These conditions are in the |
| 490 |
> |
dipole-ordered and rippled portion of the phase diagram. These are |
| 491 |
> |
also the conditions that should be most susceptible to system size |
| 492 |
> |
effects. |
| 493 |
> |
|
| 494 |
> |
\begin{figure} |
| 495 |
> |
\includegraphics[width=\linewidth]{./figures/mcSystemSize.pdf} |
| 496 |
> |
\caption{\label{mcfig:systemsize} The ripple wavelength (top) and |
| 497 |
> |
amplitude (bottom) as a function of system size for a triangular |
| 498 |
> |
lattice ($\gamma=1.732$) at $T^{*} = 122$.} |
| 499 |
> |
\end{figure} |
| 500 |
> |
|
| 501 |
> |
There is substantial dependence on system size for small (less than |
| 502 |
> |
29~$\sigma$) periodic boxes. Notably, there are resonances apparent |
| 503 |
> |
in the ripple amplitudes at box lengths of 17.3 and 29.5 $\sigma$. |
| 504 |
> |
For larger systems, the behavior of the ripples appears to have |
| 505 |
> |
stabilized and is on a trend to slightly smaller amplitudes (and |
| 506 |
> |
slightly longer wavelengths) than were observed from the 34.3 $\sigma$ |
| 507 |
> |
box sizes that were used for most of the calculations. |
| 508 |
> |
|
| 509 |
> |
It is interesting to note that system sizes which are multiples of the |
| 510 |
> |
default ripple wavelength can enhance the amplitude of the observed |
| 511 |
> |
ripples, but appears to have only a minor effect on the observed |
| 512 |
> |
wavelength. It would, of course, be better to use system sizes that |
| 513 |
> |
were many multiples of the ripple wavelength to be sure that the |
| 514 |
> |
periodic box is not driving the phenomenon, but at the largest system |
| 515 |
> |
size studied (70 $\sigma$ $\times$ 70 $\sigma$), the number of dipoles |
| 516 |
> |
(5440) made long Monte Carlo simulations prohibitively expensive. |
| 517 |
> |
|
| 518 |
> |
\section{Discussion} |
| 519 |
> |
\label{mc:sec:discussion} |
| 520 |
> |
|
| 521 |
> |
We have been able to show that a simple dipolar lattice model which |
| 522 |
> |
contains only molecular packing (from the lattice), anisotropy (in the |
| 523 |
> |
form of electrostatic dipoles) and a weak elastic tension (in the form |
| 524 |
> |
of a nearest-neighbor harmonic potential) is capable of exhibiting |
| 525 |
> |
stable long-wavelength non-thermal surface corrugations. The best |
| 526 |
> |
explanation for this behavior is that the ability of the dipoles to |
| 527 |
> |
translate out of the plane of the membrane is enough to break the |
| 528 |
> |
symmetry of the triangular lattice and allow the energetic benefit |
| 529 |
> |
from the formation of a bulk anti-ferroelectric phase. Were the weak |
| 530 |
> |
elastic tension absent from our model, it would be possible for the |
| 531 |
> |
entire lattice to ``tilt'' using $z$-translation. Tilting the lattice |
| 532 |
> |
in this way would yield an effectively non-triangular lattice which |
| 533 |
> |
would avoid dipolar frustration altogether. With the elastic tension |
| 534 |
> |
in place, bulk tilt causes a large strain, and the least costly way to |
| 535 |
> |
release this strain is between two rows of anti-aligned dipoles. |
| 536 |
> |
These ``breaks'' will result in rippled or sawtooth patterns in the |
| 537 |
> |
membrane, and allow small stripes of membrane to form |
| 538 |
> |
anti-ferroelectric regions that are tilted relative to the averaged |
| 539 |
> |
membrane normal. |
| 540 |
> |
|
| 541 |
> |
Although the dipole-dipole interaction is the major driving force for |
| 542 |
> |
the long range orientational ordered state, the formation of the |
| 543 |
> |
stable, smooth ripples is a result of the competition between the |
| 544 |
> |
elastic tension and the dipole-dipole interactions. This statement is |
| 545 |
> |
supported by the variation in $\mu^{*}$. Substantially weaker dipoles |
| 546 |
> |
relative to the surface tension can cause the corrugated phase to |
| 547 |
> |
disappear. |
| 548 |
> |
|
| 549 |
> |
The packing of the dipoles into a nearly-triangular lattice is clearly |
| 550 |
> |
an important piece of the puzzle. The dipolar head groups of lipid |
| 551 |
> |
molecules are sterically (as well as electrostatically) anisotropic, |
| 552 |
> |
and would not pack in triangular arrangements without the steric |
| 553 |
> |
interference of adjacent molecular bodies. Since we only see rippled |
| 554 |
> |
phases in the neighborhood of $\gamma=\sqrt{3}$, this implies that |
| 555 |
> |
even if this dipolar mechanism is the correct explanation for the |
| 556 |
> |
ripple phase in realistic bilayers, there would still be a role played |
| 557 |
> |
by the lipid chains in the in-plane organization of the triangularly |
| 558 |
> |
ordered phases which could support ripples. The present model is |
| 559 |
> |
certainly not detailed enough to answer exactly what drives the |
| 560 |
> |
formation of the $P_{\beta'}$ phase in real lipids, but suggests some |
| 561 |
> |
avenues for further experiments. |
| 562 |
> |
|
| 563 |
> |
The most important prediction we can make using the results from this |
| 564 |
> |
simple model is that if dipolar ordering is driving the surface |
| 565 |
> |
corrugation, the wave vectors for the ripples should always found to |
| 566 |
> |
be {\it perpendicular} to the dipole director axis. This prediction |
| 567 |
> |
should suggest experimental designs which test whether this is really |
| 568 |
> |
true in the phosphatidylcholine $P_{\beta'}$ phases. The dipole |
| 569 |
> |
director axis should also be easily computable for the all-atom and |
| 570 |
> |
coarse-grained simulations that have been published in the literature. |
| 571 |
> |
|
| 572 |
> |
Our other observation about the ripple and dipolar directionality is |
| 573 |
> |
that the dipole director axis can be found to be parallel to any of |
| 574 |
> |
the three equivalent lattice vectors in the triangular lattice. |
| 575 |
> |
Defects in the ordering of the dipoles can cause the dipole director |
| 576 |
> |
(and consequently the surface corrugation) of small regions to be |
| 577 |
> |
rotated relative to each other by 120$^{\circ}$. This is a similar |
| 578 |
> |
behavior to the domain rotation seen in the AFM studies of Kaasgaard |
| 579 |
> |
{\it et al.}\cite{Kaasgaard03} |
| 580 |
> |
|
| 581 |
> |
Although our model is simple, it exhibits some rich and unexpected |
| 582 |
> |
behaviors. It would clearly be a closer approximation to the reality |
| 583 |
> |
if we allowed greater translational freedom to the dipoles and |
| 584 |
> |
replaced the somewhat artificial lattice packing and the harmonic |
| 585 |
> |
elastic tension with more realistic molecular modeling potentials. |
| 586 |
> |
What we have done is to present a simple model which exhibits bulk |
| 587 |
> |
non-thermal corrugation, and our explanation of this rippling |
| 588 |
> |
phenomenon will help us design more accurate molecular models for |
| 589 |
> |
corrugated membranes and experiments to test whether rippling is |
| 590 |
> |
dipole-driven or not. |
