--- trunk/mdRipple/mdripple.tex 2007/08/01 16:07:12 3202 +++ trunk/mdRipple/mdripple.tex 2008/02/29 22:02:20 3351 @@ -1,13 +1,11 @@ -%\documentclass[aps,pre,twocolumn,amssymb,showpacs,floatfix]{revtex4} -%\documentclass[aps,pre,preprint,amssymb]{revtex4} \documentclass[12pt]{article} +\usepackage{graphicx} \usepackage{times} \usepackage{mathptm} \usepackage{tabularx} \usepackage{setspace} \usepackage{amsmath} \usepackage{amssymb} -\usepackage{graphicx} \usepackage[ref]{overcite} \pagestyle{plain} \pagenumbering{arabic} @@ -18,14 +16,15 @@ \renewcommand{\baselinestretch}{1.2} \renewcommand\citemid{\ } % no comma in optional reference note + \begin{document} %\renewcommand{\thefootnote}{\fnsymbol{footnote}} %\renewcommand{\theequation}{\arabic{section}.\arabic{equation}} \bibliographystyle{achemso} -\title{Dipolar Ordering in Molecular-scale Models of the Ripple Phase -in Lipid Membranes} +\title{Dipolar ordering in the ripple phases of molecular-scale models +of lipid membranes} \author{Xiuquan Sun and J. Daniel Gezelter \\ Department of Chemistry and Biochemistry,\\ University of Notre Dame, \\ @@ -38,11 +37,25 @@ The ripple phase in phosphatidylcholine (PC) bilayers \maketitle \begin{abstract} -The ripple phase in phosphatidylcholine (PC) bilayers has never been -completely explained. +Symmetric and asymmetric ripple phases have been observed to form in +molecular dynamics simulations of a simple molecular-scale lipid +model. The lipid model consists of an dipolar head group and an +ellipsoidal tail. Within the limits of this model, an explanation for +generalized membrane curvature is a simple mismatch in the size of the +heads with the width of the molecular bodies. The persistence of a +{\it bilayer} structure requires strong attractive forces between the +head groups. One feature of this model is that an energetically +favorable orientational ordering of the dipoles can be achieved by +out-of-plane membrane corrugation. The corrugation of the surface +stabilizes the long range orientational ordering for the dipoles in the +head groups which then adopt a bulk anti-ferroelectric state. We +observe a common feature of the corrugated dipolar membranes: the wave +vectors for the surface ripples are always found to be perpendicular +to the dipole director axis. \end{abstract} %\maketitle +\newpage \section{Introduction} \label{sec:Int} @@ -61,46 +74,56 @@ within the gel phase.~\cite{Cevc87} experimental results provide strong support for a 2-dimensional hexagonal packing lattice of the lipid molecules within the ripple phase. This is a notable change from the observed lipid packing -within the gel phase.~\cite{Cevc87} +within the gel phase,~\cite{Cevc87} although Tenchov {\it et al.} have +recently observed near-hexagonal packing in some phosphatidylcholine +(PC) gel phases.\cite{Tenchov2001} The X-ray diffraction work by +Katsaras {\it et al.} showed that a rich phase diagram exhibiting both +{\it asymmetric} and {\it symmetric} ripples is possible for lecithin +bilayers.\cite{Katsaras00} A number of theoretical models have been presented to explain the formation of the ripple phase. Marder {\it et al.} used a -curvature-dependent Landau-de Gennes free-energy functional to predict -a rippled phase.~\cite{Marder84} This model and other related continuum -models predict higher fluidity in convex regions and that concave -portions of the membrane correspond to more solid-like regions. -Carlson and Sethna used a packing-competition model (in which head -groups and chains have competing packing energetics) to predict the -formation of a ripple-like phase. Their model predicted that the -high-curvature portions have lower-chain packing and correspond to -more fluid-like regions. Goldstein and Leibler used a mean-field -approach with a planar model for {\em inter-lamellar} interactions to -predict rippling in multilamellar phases.~\cite{Goldstein88} McCullough -and Scott proposed that the {\em anisotropy of the nearest-neighbor -interactions} coupled to hydrophobic constraining forces which -restrict height differences between nearest neighbors is the origin of -the ripple phase.~\cite{McCullough90} Lubensky and MacKintosh -introduced a Landau theory for tilt order and curvature of a single -membrane and concluded that {\em coupling of molecular tilt to membrane -curvature} is responsible for the production of -ripples.~\cite{Lubensky93} Misbah, Duplat and Houchmandzadeh proposed -that {\em inter-layer dipolar interactions} can lead to ripple -instabilities.~\cite{Misbah98} Heimburg presented a {\em coexistence -model} for ripple formation in which he postulates that fluid-phase -line defects cause sharp curvature between relatively flat gel-phase -regions.~\cite{Heimburg00} Kubica has suggested that a lattice model of -polar head groups could be valuable in trying to understand bilayer -phase formation.~\cite{Kubica02} Bannerjee used Monte Carlo simulations -of lamellar stacks of hexagonal lattices to show that large headgroups -and molecular tilt with respect to the membrane normal vector can -cause bulk rippling.~\cite{Bannerjee02} +curvature-dependent Landau-de~Gennes free-energy functional to predict +a rippled phase.~\cite{Marder84} This model and other related +continuum models predict higher fluidity in convex regions and that +concave portions of the membrane correspond to more solid-like +regions. Carlson and Sethna used a packing-competition model (in +which head groups and chains have competing packing energetics) to +predict the formation of a ripple-like phase. Their model predicted +that the high-curvature portions have lower-chain packing and +correspond to more fluid-like regions. Goldstein and Leibler used a +mean-field approach with a planar model for {\em inter-lamellar} +interactions to predict rippling in multilamellar +phases.~\cite{Goldstein88} McCullough and Scott proposed that the {\em +anisotropy of the nearest-neighbor interactions} coupled to +hydrophobic constraining forces which restrict height differences +between nearest neighbors is the origin of the ripple +phase.~\cite{McCullough90} Lubensky and MacKintosh introduced a Landau +theory for tilt order and curvature of a single membrane and concluded +that {\em coupling of molecular tilt to membrane curvature} is +responsible for the production of ripples.~\cite{Lubensky93} Misbah, +Duplat and Houchmandzadeh proposed that {\em inter-layer dipolar +interactions} can lead to ripple instabilities.~\cite{Misbah98} +Heimburg presented a {\em coexistence model} for ripple formation in +which he postulates that fluid-phase line defects cause sharp +curvature between relatively flat gel-phase regions.~\cite{Heimburg00} +Kubica has suggested that a lattice model of polar head groups could +be valuable in trying to understand bilayer phase +formation.~\cite{Kubica02} Bannerjee used Monte Carlo simulations of +lamellar stacks of hexagonal lattices to show that large head groups +and molecular tilt with respect to the membrane normal vector can +cause bulk rippling.~\cite{Bannerjee02} Recently, Kranenburg and Smit +described the formation of symmetric ripple-like structures using a +coarse grained solvent-head-tail bead model.\cite{Kranenburg2005} +Their lipids consisted of a short chain of head beads tied to the two +longer ``chains''. -In contrast, few large-scale molecular modelling studies have been +In contrast, few large-scale molecular modeling studies have been done due to the large size of the resulting structures and the time required for the phases of interest to develop. With all-atom (and even unified-atom) simulations, only one period of the ripple can be -observed and only for timescales in the range of 10-100 ns. One of -the most interesting molecular simulations was carried out by De Vries +observed and only for time scales in the range of 10-100 ns. One of +the most interesting molecular simulations was carried out by de~Vries {\it et al.}~\cite{deVries05}. According to their simulation results, the ripple consists of two domains, one resembling the gel bilayer, while in the other, the two leaves of the bilayer are fully @@ -122,18 +145,18 @@ between headgroups and tails is strongly implicated as Although the organization of the tails of lipid molecules are addressed by these molecular simulations and the packing competition -between headgroups and tails is strongly implicated as the primary +between head groups and tails is strongly implicated as the primary driving force for ripple formation, questions about the ordering of -the head groups in ripple phase has not been settled. +the head groups in ripple phase have not been settled. In a recent paper, we presented a simple ``web of dipoles'' spin lattice model which provides some physical insight into relationship between dipolar ordering and membrane buckling.\cite{Sun2007} We found that dipolar elastic membranes can spontaneously buckle, forming -ripple-like topologies. The driving force for the buckling in dipolar -elastic membranes the antiferroelectric ordering of the dipoles, and -this was evident in the ordering of the dipole director axis -perpendicular to the wave vector of the surface ripples. A similiar +ripple-like topologies. The driving force for the buckling of dipolar +elastic membranes is the anti-ferroelectric ordering of the dipoles. +This was evident in the ordering of the dipole director axis +perpendicular to the wave vector of the surface ripples. A similar phenomenon has also been observed by Tsonchev {\it et al.} in their work on the spontaneous formation of dipolar peptide chains into curved nano-structures.\cite{Tsonchev04,Tsonchev04II} @@ -167,7 +190,7 @@ some fraction of the details of the chain dynamics neg non-polar tails. Another fact is that the majority of lipid molecules in the ripple phase are relatively rigid (i.e. gel-like) which makes some fraction of the details of the chain dynamics negligible. Figure -\ref{fig:lipidModels} shows the molecular strucure of a DPPC +\ref{fig:lipidModels} shows the molecular structure of a DPPC molecule, as well as atomistic and molecular-scale representations of a DPPC molecule. The hydrophilic character of the head group is largely due to the separation of charge between the nitrogen and @@ -186,7 +209,7 @@ modelling large length-scale properties of lipid The ellipsoidal portions of the model interact via the Gay-Berne potential which has seen widespread use in the liquid crystal community. Ayton and Voth have also used Gay-Berne ellipsoids for -modelling large length-scale properties of lipid +modeling large length-scale properties of lipid bilayers.\cite{Ayton01} In its original form, the Gay-Berne potential was a single site model for the interactions of rigid ellipsoidal molecules.\cite{Gay81} It can be thought of as a modification of the @@ -236,7 +259,7 @@ ellipsoids in a {\it side-by-side} configuration. Add Gay-Berne ellipsoids also have an energy scaling parameter, $\epsilon^s$, which describes the well depth for two identical -ellipsoids in a {\it side-by-side} configuration. Additionaly, a well +ellipsoids in a {\it side-by-side} configuration. Additionally, a well depth aspect ratio, $\epsilon^r = \epsilon^e / \epsilon^s$, describes the ratio between the well depths in the {\it end-to-end} and side-by-side configurations. As in the range parameter, a set of @@ -294,19 +317,19 @@ models. $l / d$ is the ratio of the head group to body \centering \includegraphics[width=4in]{2lipidModel} \caption{The parameters defining the behavior of the lipid -models. $l / d$ is the ratio of the head group to body diameter. +models. $\sigma_h / d$ is the ratio of the head group to body diameter. Molecular bodies had a fixed aspect ratio of 3.0. The solvent model was a simplified 4-water bead ($\sigma_w \approx d$) that has been -used in other coarse-grained (DPD) simulations. The dipolar strength +used in other coarse-grained simulations. The dipolar strength (and the temperature and pressure) were the only other parameters that were varied systematically.\label{fig:lipidModel}} \end{figure} To take into account the permanent dipolar interactions of the -zwitterionic head groups, we place fixed dipole moments $\mu_{i}$ at +zwitterionic head groups, we have placed fixed dipole moments $\mu_{i}$ at one end of the Gay-Berne particles. The dipoles are oriented at an angle $\theta = \pi / 2$ relative to the major axis. These dipoles -are protected by a head ``bead'' with a range parameter which we have +are protected by a head ``bead'' with a range parameter ($\sigma_h$) which we have varied between $1.20 d$ and $1.41 d$. The head groups interact with each other using a combination of Lennard-Jones, \begin{equation} @@ -325,6 +348,28 @@ For the interaction between nonequivalent uniaxial ell along dipole $i$ and $\mathbf{\hat r}_{ij}$ is the unit vector pointing along the inter-dipole vector $\mathbf{r}_{ij}$. +Since the charge separation distance is so large in zwitterionic head +groups (like the PC head groups), it would also be possible to use +either point charges or a ``split dipole'' approximation, +\begin{equation} +V_{ij}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf +\hat{r}}_{ij})) = \frac{1}{{4\pi \epsilon_0 }}\left[ {\frac{{\mu _i \cdot \mu _j }}{{R_{ij}^3 }} - +\frac{{3\left( {\mu _i \cdot r_{ij} } \right)\left( {\mu _i \cdot +r_{ij} } \right)}}{{R_{ij}^5 }}} \right] +\end{equation} +where $\mu _i$ and $\mu _j$ are the dipole moments of sites $i$ and +$j$, $r_{ij}$ is vector between the two sites, and $R_{ij}$ is given +by, +\begin{equation} +R_{ij} = \sqrt {r_{ij}^2 + \frac{{d_i^2 }}{4} + \frac{{d_j^2 +}}{4}}. +\end{equation} +Here, $d_i$ and $d_j$ are charge separation distances associated with +each of the two dipolar sites. This approximation to the multipole +expansion maintains the fast fall-off of the multipole potentials but +lacks the normal divergences when two polar groups get close to one +another. + For the interaction between nonequivalent uniaxial ellipsoids (in this case, between spheres and ellipsoids), the spheres are treated as ellipsoids with an aspect ratio of 1 ($d = l$) and with an well depth @@ -333,13 +378,18 @@ The solvent model in our simulations is identical to o et al.} and is appropriate for dissimilar uniaxial ellipsoids.\cite{Cleaver96} -The solvent model in our simulations is identical to one used by -Marrink {\it et al.} in their dissipative particle dynamics (DPD) -simulation of lipid bilayers.\cite{Marrink04} This solvent bead is a single -site that represents four water molecules (m = 72 amu) and has -comparable density and diffusive behavior to liquid water. However, -since there are no electrostatic sites on these beads, this solvent -model cannot replicate the dielectric properties of water. +The solvent model in our simulations is similar to the one used by +Marrink {\it et al.} in their coarse grained simulations of lipid +bilayers.\cite{Marrink04} The solvent bead is a single site that +represents four water molecules (m = 72 amu) and has comparable +density and diffusive behavior to liquid water. However, since there +are no electrostatic sites on these beads, this solvent model cannot +replicate the dielectric properties of water. Note that although we +are using larger cutoff and switching radii than Marrink {\it et al.}, +our solvent density at 300 K remains at 0.944 g cm$^{-3}$, and the +solvent diffuses at 0.43 $\AA^2 ps^{-1}$ (only twice as fast as liquid +water). + \begin{table*} \begin{minipage}{\linewidth} \begin{center} @@ -365,100 +415,154 @@ A switching function has been applied to all potential \end{minipage} \end{table*} -A switching function has been applied to all potentials to smoothly -turn off the interactions between a range of $22$ and $25$ \AA. +\section{Experimental Methodology} +\label{sec:experiment} The parameters that were systematically varied in this study were the size of the head group ($\sigma_h$), the strength of the dipole moment ($\mu$), and the temperature of the system. Values for $\sigma_h$ -ranged from 5.5 \AA\ to 6.5 \AA\ . If the width of the tails is -taken to be the unit of length, these head groups correspond to a -range from $1.2 d$ to $1.41 d$. Since the solvent beads are nearly -identical in diameter to the tail ellipsoids, all distances that -follow will be measured relative to this unit of distance. +ranged from 5.5 \AA\ to 6.5 \AA. If the width of the tails is taken +to be the unit of length, these head groups correspond to a range from +$1.2 d$ to $1.41 d$. Since the solvent beads are nearly identical in +diameter to the tail ellipsoids, all distances that follow will be +measured relative to this unit of distance. Because the solvent we +are using is non-polar and has a dielectric constant of 1, values for +$\mu$ are sampled from a range that is somewhat smaller than the 20.6 +Debye dipole moment of the PC head groups. -\section{Experimental Methodology} -\label{sec:experiment} - To create unbiased bilayers, all simulations were started from two perfectly flat monolayers separated by a 26 \AA\ gap between the molecular bodies of the upper and lower leaves. The separated -monolayers were evolved in a vaccum with $x-y$ anisotropic pressure +monolayers were evolved in a vacuum with $x-y$ anisotropic pressure coupling. The length of $z$ axis of the simulations was fixed and a constant surface tension was applied to enable real fluctuations of the bilayer. Periodic boundary conditions were used, and $480-720$ lipid molecules were present in the simulations, depending on the size of the head beads. In all cases, the two monolayers spontaneously collapsed into bilayer structures within 100 ps. Following this -collapse, all systems were equlibrated for $100$ ns at $300$ K. +collapse, all systems were equilibrated for $100$ ns at $300$ K. The resulting bilayer structures were then solvated at a ratio of $6$ solvent beads (24 water molecules) per lipid. These configurations were then equilibrated for another $30$ ns. All simulations utilizing the solvent were carried out at constant pressure ($P=1$ atm) with -$3$D anisotropic coupling, and constant surface tension -($\gamma=0.015$ UNIT). Given the absence of fast degrees of freedom in -this model, a timestep of $50$ fs was utilized with excellent energy +$3$D anisotropic coupling, and small constant surface tension +($\gamma=0.015$ N/m). Given the absence of fast degrees of freedom in +this model, a time step of $50$ fs was utilized with excellent energy conservation. Data collection for structural properties of the bilayers was carried out during a final 5 ns run following the solvent -equilibration. All simulations were performed using the OOPSE -molecular modeling program.\cite{Meineke05} +equilibration. Orientational correlation functions and diffusion +constants were computed from 30 ns simulations in the microcanonical +(NVE) ensemble using the average volume from the end of the constant +pressure and surface tension runs. The timestep on these final +molecular dynamics runs was 25 fs. No appreciable changes in phase +structure were noticed upon switching to a microcanonical ensemble. +All simulations were performed using the {\sc oopse} molecular +modeling program.\cite{Meineke05} +A switching function was applied to all potentials to smoothly turn +off the interactions between a range of $22$ and $25$ \AA. The +switching function was the standard (cubic) function, +\begin{equation} +s(r) = + \begin{cases} + 1 & \text{if $r \le r_{\text{sw}}$},\\ + \frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2} + {(r_{\text{cut}} - r_{\text{sw}})^3} + & \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\ + 0 & \text{if $r > r_{\text{cut}}$.} + \end{cases} +\label{eq:dipoleSwitching} +\end{equation} + \section{Results} \label{sec:results} -Snapshots in Figure \ref{fig:phaseCartoon} show that the membrane is -more corrugated with increasing size of the head groups. The surface -is nearly flat when $\sigma_h=1.20 d$. With $\sigma_h=1.28 d$, -although the surface is still flat, the bilayer starts to splay -inward; the upper leaf of the bilayer is connected to the lower leaf -with an interdigitated line defect. Two periodicities with $100$ \AA\ -wavelengths were observed in the simulation. This structure is very -similiar to the structure observed by de Vries and Lenz {\it et -al.}. The same basic structure is also observed when $\sigma_h=1.41 -d$, but the wavelength of the surface corrugations depends sensitively -on the size of the ``head'' beads. From the undulation spectrum, the -corrugation is clearly non-thermal. +The membranes in our simulations exhibit a number of interesting +bilayer phases. The surface topology of these phases depends most +sensitively on the ratio of the size of the head groups to the width +of the molecular bodies. With heads only slightly larger than the +bodies ($\sigma_h=1.20 d$) the membrane exhibits a flat bilayer. + +Increasing the head / body size ratio increases the local membrane +curvature around each of the lipids. With $\sigma_h=1.28 d$, the +surface is still essentially flat, but the bilayer starts to exhibit +signs of instability. We have observed occasional defects where a +line of lipid molecules on one leaf of the bilayer will dip down to +interdigitate with the other leaf. This gives each of the two bilayer +leaves some local convexity near the line defect. These structures, +once developed in a simulation, are very stable and are spaced +approximately 100 \AA\ away from each other. + +With larger heads ($\sigma_h = 1.35 d$) the membrane curvature +resolves into a ``symmetric'' ripple phase. Each leaf of the bilayer +is broken into several convex, hemicylinderical sections, and opposite +leaves are fitted together much like roof tiles. There is no +interdigitation between the upper and lower leaves of the bilayer. + +For the largest head / body ratios studied ($\sigma_h = 1.41 d$) the +local curvature is substantially larger, and the resulting bilayer +structure resolves into an asymmetric ripple phase. This structure is +very similar to the structures observed by both de~Vries {\it et al.} +and Lenz {\it et al.}. For a given ripple wave vector, there are two +possible asymmetric ripples, which is not the case for the symmetric +phase observed when $\sigma_h = 1.35 d$. + \begin{figure}[htb] \centering \includegraphics[width=4in]{phaseCartoon} -\caption{A sketch to discribe the structure of the phases observed in -our simulations.\label{fig:phaseCartoon}} +\caption{The role of the ratio between the head group size and the +width of the molecular bodies is to increase the local membrane +curvature. With strong attractive interactions between the head +groups, this local curvature can be maintained in bilayer structures +through surface corrugation. Shown above are three phases observed in +these simulations. With $\sigma_h = 1.20 d$, the bilayer maintains a +flat topology. For larger heads ($\sigma_h = 1.35 d$) the local +curvature resolves into a symmetrically rippled phase with little or +no interdigitation between the upper and lower leaves of the membrane. +The largest heads studied ($\sigma_h = 1.41 d$) resolve into an +asymmetric rippled phases with interdigitation between the two +leaves.\label{fig:phaseCartoon}} \end{figure} -When $\sigma_h=1.35 d$, we observed another corrugated surface -morphology. This structure is different from the asymmetric rippled -surface; there is no interdigitation between the upper and lower -leaves of the bilayer. Each leaf of the bilayer is broken into several -hemicylinderical sections, and opposite leaves are fitted together -much like roof tiles. Unlike the surface in which the upper -hemicylinder is always interdigitated on the leading or trailing edge -of lower hemicylinder, this ``symmetric'' ripple has no prefered -direction. The corresponding structures are shown in Figure -\ref{fig:phaseCartoon} for elucidation of the detailed structures of -different phases. The top panel in figure \ref{fig:phaseCartoon} is -the flat phase, the middle panel shows the asymmetric ripple phase -corresponding to $\sigma_h = 1.41 d$ and the lower panel shows the -symmetric ripple phase observed when $\sigma_h=1.35 d$. In the -symmetric ripple, the bilayer is continuous over the whole membrane, -however, in asymmetric ripple phase, the bilayer domains are connected -by thin interdigitated monolayers that share molecules between the -upper and lower leaves. +Sample structures for the flat ($\sigma_h = 1.20 d$), symmetric +($\sigma_h = 1.35 d$, and asymmetric ($\sigma_h = 1.41 d$) ripple +phases are shown in Figure \ref{fig:phaseCartoon}. + +It is reasonable to ask how well the parameters we used can produce +bilayer properties that match experimentally known values for real +lipid bilayers. Using a value of $l = 13.8$ \AA for the ellipsoidal +tails and the fixed ellipsoidal aspect ratio of 3, our values for the +area per lipid ($A$) and inter-layer spacing ($D_{HH}$) depend +entirely on the size of the head bead relative to the molecular body. +These values are tabulated in table \ref{tab:property}. Kucera {\it +et al.} have measured values for the head group spacings for a number +of PC lipid bilayers that range from 30.8 \AA\ (DLPC) to 37.8 \AA\ (DPPC). +They have also measured values for the area per lipid that range from +60.6 +\AA$^2$ (DMPC) to 64.2 \AA$^2$ +(DPPC).\cite{NorbertKucerka04012005,NorbertKucerka06012006} Only the +largest of the head groups we modeled ($\sigma_h = 1.41 d$) produces +bilayers (specifically the area per lipid) that resemble real PC +bilayers. The smaller head beads we used are perhaps better models +for PE head groups. + \begin{table*} \begin{minipage}{\linewidth} \begin{center} -\caption{Phases, ripple wavelengths and amplitudes observed as a -function of the ratio between the head beads and the diameters of the -tails. All lengths are normalized to the diameter of the tail -ellipsoids.} -\begin{tabular}{lccc} +\caption{Phase, bilayer spacing, area per lipid, ripple wavelength +and amplitude observed as a function of the ratio between the head +beads and the diameters of the tails. Ripple wavelengths and +amplitudes are normalized to the diameter of the tail ellipsoids.} +\begin{tabular}{lccccc} \hline -$\sigma_h / d$ & type of phase & $\lambda / d$ & $A / d$\\ +$\sigma_h / d$ & type of phase & bilayer spacing (\AA) & area per +lipid (\AA$^2$) & $\lambda / d$ & $A / d$\\ \hline -1.20 & flat & N/A & N/A \\ -1.28 & asymmetric ripple or flat & 21.7 & N/A \\ -1.35 & symmetric ripple & 17.2 & 2.2 \\ -1.41 & asymmetric ripple & 15.4 & 1.5 \\ +1.20 & flat & 33.4 & 49.6 & N/A & N/A \\ +1.28 & flat & 33.7 & 54.7 & N/A & N/A \\ +1.35 & symmetric ripple & 42.9 & 51.7 & 17.2 & 2.2 \\ +1.41 & asymmetric ripple & 37.1 & 63.1 & 15.4 & 1.5 \\ \end{tabular} \label{tab:property} \end{center} @@ -467,24 +571,25 @@ reduced amplitude $A / d$ of the ripples are summarize The membrane structures and the reduced wavelength $\lambda / d$, reduced amplitude $A / d$ of the ripples are summarized in Table -\ref{tab:property}. The wavelength range is $15~21$ molecular bodies +\ref{tab:property}. The wavelength range is $15 - 17$ molecular bodies and the amplitude is $1.5$ molecular bodies for asymmetric ripple and -$2.2$ for symmetric ripple. These values are consistent to the -experimental results. Note, that given the lack of structural freedom -in the tails of our model lipids, the amplitudes observed from these -simulations are likely to underestimate of the true amplitudes. +$2.2$ for symmetric ripple. These values are reasonably consistent +with experimental measurements.\cite{Sun96,Katsaras00,Kaasgaard03} +Note, that given the lack of structural freedom in the tails of our +model lipids, the amplitudes observed from these simulations are +likely to underestimate of the true amplitudes. \begin{figure}[htb] \centering \includegraphics[width=4in]{topDown} -\caption{Top views of the flat (upper), asymmetric ripple (middle), -and symmetric ripple (lower) phases. Note that the head-group dipoles -have formed head-to-tail chains in all three of these phases, but in -the two rippled phases, the dipolar chains are all aligned -{\it perpendicular} to the direction of the ripple. The flat membrane -has multiple point defects in the dipolar orientational ordering, and -the dipolar ordering on the lower leaf of the bilayer can be in a -different direction from the upper leaf.\label{fig:topView}} +\caption{Top views of the flat (upper), symmetric ripple (middle), +and asymmetric ripple (lower) phases. Note that the head-group +dipoles have formed head-to-tail chains in all three of these phases, +but in the two rippled phases, the dipolar chains are all aligned {\it +perpendicular} to the direction of the ripple. Note that the flat +membrane has multiple vortex defects in the dipolar ordering, and the +ordering on the lower leaf of the bilayer can be in an entirely +different direction from the upper leaf.\label{fig:topView}} \end{figure} The principal method for observing orientational ordering in dipolar @@ -516,7 +621,7 @@ flat ($\sigma_h = 1.20 d$) and rippled ($\sigma_h = 1. groups to be completely decoupled from each other. Figure \ref{fig:topView} shows snapshots of bird's-eye views of the -flat ($\sigma_h = 1.20 d$) and rippled ($\sigma_h = 1.41, 1.35 d$) +flat ($\sigma_h = 1.20 d$) and rippled ($\sigma_h = 1.35, 1.41 d$) bilayers. The directions of the dipoles on the head groups are represented with two colored half spheres: blue (phosphate) and yellow (amino). For flat bilayers, the system exhibits signs of @@ -539,7 +644,7 @@ antiferroelectric state. It is also notable that the configuration, and the dipolar order parameter increases dramatically. However, the total polarization of the system is still close to zero. This is strong evidence that the corrugated structure is an -antiferroelectric state. It is also notable that the head-to-tail +anti-ferroelectric state. It is also notable that the head-to-tail arrangement of the dipoles is always observed in a direction perpendicular to the wave vector for the surface corrugation. This is a similar finding to what we observed in our earlier work on the @@ -580,10 +685,10 @@ interaction stength. When the interaction between the increasing strength of the dipole. Generally, the dipoles on the head groups become more ordered as the strength of the interaction between heads is increased and become more disordered by decreasing the -interaction stength. When the interaction between the heads becomes +interaction strength. When the interaction between the heads becomes too weak, the bilayer structure does not persist; all lipid molecules become dispersed in the solvent (which is non-polar in this -molecular-scale model). The critial value of the strength of the +molecular-scale model). The critical value of the strength of the dipole depends on the size of the head groups. The perfectly flat surface becomes unstable below $5$ Debye, while the rippled surfaces melt at $8$ Debye (asymmetric) or $10$ Debye (symmetric). @@ -596,11 +701,11 @@ dipolar interactions, the tails are forced to splay ou close to each other and distort the bilayer structure. For a flat surface, a substantial amount of free volume between the head groups is normally available. When the head groups are brought closer by -dipolar interactions, the tails are forced to splay outward, forming -first curved bilayers, and then inverted micelles. +dipolar interactions, the tails are forced to splay outward, first forming +curved bilayers, and then inverted micelles. When $\sigma_h=1.28 d$, the $P_2$ order parameter decreases slightly -when the strength of the dipole is increased above $16$ debye. For +when the strength of the dipole is increased above $16$ Debye. For rippled bilayers, there is less free volume available between the head groups. Therefore increasing dipolar strength weakly influences the structure of the membrane. However, the increase in the body $P_2$ @@ -657,67 +762,194 @@ molecular width ratio ($\sigma_h / d$).\label{fig:tP2} molecular width ratio ($\sigma_h / d$).\label{fig:tP2}} \end{figure} -\section{Discussion} -\label{sec:discussion} +Fig. \ref{fig:phaseDiagram} shows a phase diagram for the model as a +function of the head group / molecular width ratio ($\sigma_h / d$) +and the strength of the head group dipole moment ($\mu$). Note that +the specific form of the bilayer phase is governed almost entirely by +the head group / molecular width ratio, while the strength of the +dipolar interactions between the head groups governs the stability of +the bilayer phase. Weaker dipoles result in unstable bilayer phases, +while extremely strong dipoles can shift the equilibrium to an +inverted micelle phase when the head groups are small. Temperature +has little effect on the actual bilayer phase observed, although higher +temperatures can cause the unstable region to grow into the higher +dipole region of this diagram. -The ripple phases have been observed in our molecular dynamic -simulations using a simple molecular lipid model. The lipid model -consists of an anisotropic interacting dipolar head group and an -ellipsoid shape tail. According to our simulations, the explanation of -the formation for the ripples are originated in the size mismatch -between the head groups and the tails. The ripple phases are only -observed in the studies using larger head group lipid models. However, -there is a mismatch betweent the size of the head groups and the size -of the tails in the simulations of the flat surface. This indicates -the competition between the anisotropic dipolar interaction and the -packing of the tails also plays a major role for formation of the -ripple phase. The larger head groups provide more free volume for the -tails, while these hydrophobic ellipsoids trying to be close to each -other, this gives the origin of the spontanous curvature of the -surface, which is believed as the beginning of the ripple phases. The -lager head groups cause the spontanous curvature inward for both of -leaves of the bilayer. This results in a steric strain when the tails -of two leaves too close to each other. The membrane has to be broken -to release this strain. There are two ways to arrange these broken -curvatures: symmetric and asymmetric ripples. Both of the ripple -phases have been observed in our studies. The difference between these -two ripples is that the bilayer is continuum in the symmetric ripple -phase and is disrupt in the asymmetric ripple phase. +\begin{figure}[htb] +\centering +\includegraphics[width=\linewidth]{phaseDiagram} +\caption{Phase diagram for the simple molecular model as a function +of the head group / molecular width ratio ($\sigma_h / d$) and the +strength of the head group dipole moment +($\mu$).\label{fig:phaseDiagram}} +\end{figure} -Dipolar head groups are the key elements for the maintaining of the -bilayer structure. The lipids are solvated in water when lowering the -the strength of the dipole on the head groups. The long range -orientational ordering of the dipoles can be achieved by forming the -ripples, although the dipoles are likely to form head-to-tail -configurations even in flat surface, the frustration prevents the -formation of the long range orientational ordering for dipoles. The -corrugation of the surface breaks the frustration and stablizes the -long range oreintational ordering for the dipoles in the head groups -of the lipid molecules. Many rows of the head-to-tail dipoles are -parallel to each other and adopt the antiferroelectric state as a -whole. This is the first time the organization of the head groups in -ripple phases of the lipid bilayer has been addressed. +We have computed translational diffusion constants for lipid molecules +from the mean-square displacement, +\begin{equation} +D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \rangle, +\end{equation} +of the lipid bodies. Translational diffusion constants for the +different head-to-tail size ratios (all at 300 K) are shown in table +\ref{tab:relaxation}. We have also computed orientational correlation +times for the head groups from fits of the second-order Legendre +polynomial correlation function, +\begin{equation} +C_{\ell}(t) = \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf +\mu}_{i}(0) \right) \rangle +\end{equation} +of the head group dipoles. The orientational correlation functions +appear to have multiple components in their decay: a fast ($12 \pm 2$ +ps) decay due to librational motion of the head groups, as well as +moderate ($\tau^h_{\rm mid}$) and slow ($\tau^h_{\rm slow}$) +components. The fit values for the moderate and slow correlation +times are listed in table \ref{tab:relaxation}. Standard deviations +in the fit time constants are quite large (on the order of the values +themselves). -The most important prediction we can make using the results from this -simple model is that if dipolar ordering is driving the surface -corrugation, the wave vectors for the ripples should always found to -be {\it perpendicular} to the dipole director axis. This prediction -should suggest experimental designs which test whether this is really -true in the phosphatidylcholine $P_{\beta'}$ phases. The dipole -director axis should also be easily computable for the all-atom and -coarse-grained simulations that have been published in the literature. +Sparrman and Westlund used a multi-relaxation model for NMR lineshapes +observed in gel, fluid, and ripple phases of DPPC and obtained +estimates of a correlation time for water translational diffusion +($\tau_c$) of 20 ns.\cite{Sparrman2003} This correlation time +corresponds to water bound to small regions of the lipid membrane. +They further assume that the lipids can explore only a single period +of the ripple (essentially moving in a nearly one-dimensional path to +do so), and the correlation time can therefore be used to estimate a +value for the translational diffusion constant of $2.25 \times +10^{-11} m^2 s^{-1}$. The translational diffusion constants we obtain +are in reasonable agreement with this experimentally determined +value. However, the $T_2$ relaxation times obtained by Sparrman and +Westlund are consistent with P-N vector reorientation timescales of +2-5 ms. This is substantially slower than even the slowest component +we observe in the decay of our orientational correlation functions. +Other than the dipole-dipole interactions, our head groups have no +shape anisotropy which would force them to move as a unit with +neighboring molecules. This would naturally lead to P-N reorientation +times that are too fast when compared with experimental measurements. +\begin{table*} +\begin{minipage}{\linewidth} +\begin{center} +\caption{Fit values for the rotational correlation times for the head +groups ($\tau^h$) and molecular bodies ($\tau^b$) as well as the +translational diffusion constants for the molecule as a function of +the head-to-body width ratio. All correlation functions and transport +coefficients were computed from microcanonical simulations with an +average temperture of 300 K. In all of the phases, the head group +correlation functions decay with an fast librational contribution ($12 +\pm 1$ ps). There are additional moderate ($\tau^h_{\rm mid}$) and +slow $\tau^h_{\rm slow}$ contributions to orientational decay that +depend strongly on the phase exhibited by the lipids. The symmetric +ripple phase ($\sigma_h / d = 1.35$) appears to exhibit the slowest +molecular reorientation.} +\begin{tabular}{lcccc} +\hline +$\sigma_h / d$ & $\tau^h_{\rm mid} (ns)$ & $\tau^h_{\rm +slow} (\mu s)$ & $\tau^b (\mu s)$ & $D (\times 10^{-11} m^2 s^{-1})$ \\ +\hline +1.20 & $0.4$ & $9.6$ & $9.5$ & $0.43(1)$ \\ +1.28 & $2.0$ & $13.5$ & $3.0$ & $5.91(3)$ \\ +1.35 & $3.2$ & $4.0$ & $0.9$ & $3.42(1)$ \\ +1.41 & $0.3$ & $23.8$ & $6.9$ & $7.16(1)$ \\ +\end{tabular} +\label{tab:relaxation} +\end{center} +\end{minipage} +\end{table*} + +\section{Discussion} +\label{sec:discussion} + +Symmetric and asymmetric ripple phases have been observed to form in +our molecular dynamics simulations of a simple molecular-scale lipid +model. The lipid model consists of an dipolar head group and an +ellipsoidal tail. Within the limits of this model, an explanation for +generalized membrane curvature is a simple mismatch in the size of the +heads with the width of the molecular bodies. With heads +substantially larger than the bodies of the molecule, this curvature +should be convex nearly everywhere, a requirement which could be +resolved either with micellar or cylindrical phases. + +The persistence of a {\it bilayer} structure therefore requires either +strong attractive forces between the head groups or exclusionary +forces from the solvent phase. To have a persistent bilayer structure +with the added requirement of convex membrane curvature appears to +result in corrugated structures like the ones pictured in +Fig. \ref{fig:phaseCartoon}. In each of the sections of these +corrugated phases, the local curvature near a most of the head groups +is convex. These structures are held together by the extremely strong +and directional interactions between the head groups. + +The attractive forces holding the bilayer together could either be +non-directional (as in the work of Kranenburg and +Smit),\cite{Kranenburg2005} or directional (as we have utilized in +these simulations). The dipolar head groups are key for the +maintaining the bilayer structures exhibited by this particular model; +reducing the strength of the dipole has the tendency to make the +rippled phase disappear. The dipoles are likely to form attractive +head-to-tail configurations even in flat configurations, but the +temperatures are high enough that vortex defects become prevalent in +the flat phase. The flat phase we observed therefore appears to be +substantially above the Kosterlitz-Thouless transition temperature for +a planar system of dipoles with this set of parameters. For this +reason, it would be interesting to observe the thermal behavior of the +flat phase at substantially lower temperatures. + +One feature of this model is that an energetically favorable +orientational ordering of the dipoles can be achieved by forming +ripples. The corrugation of the surface breaks the symmetry of the +plane, making vortex defects somewhat more expensive, and stabilizing +the long range orientational ordering for the dipoles in the head +groups. Most of the rows of the head-to-tail dipoles are parallel to +each other and the system adopts a bulk anti-ferroelectric state. We +believe that this is the first time the organization of the head +groups in ripple phases has been addressed. + +Although the size-mismatch between the heads and molecular bodies +appears to be the primary driving force for surface convexity, the +persistence of the bilayer through the use of rippled structures is a +function of the strong, attractive interactions between the heads. +One important prediction we can make using the results from this +simple model is that if the dipole-dipole interaction is the leading +contributor to the head group attractions, the wave vectors for the +ripples should always be found {\it perpendicular} to the dipole +director axis. This echoes the prediction we made earlier for simple +elastic dipolar membranes, and may suggest experimental designs which +will test whether this is really the case in the phosphatidylcholine +$P_{\beta'}$ phases. The dipole director axis should also be easily +computable for the all-atom and coarse-grained simulations that have +been published in the literature.\cite{deVries05} + +Experimental verification of our predictions of dipolar orientation +correlating with the ripple direction would require knowing both the +local orientation of a rippled region of the membrane (available via +AFM studies of supported bilayers) as well as the local ordering of +the membrane dipoles. Obtaining information about the local +orientations of the membrane dipoles may be available from +fluorescence detected linear dichroism (LD). Benninger {\it et al.} +have recently used axially-specific chromophores +2-(4,4-difluoro-5,7-dimethyl-4-bora-3a,4a-diaza-s-indacene-3-pentanoyl)-1-hexadecanoyl-sn-glycero-3-phospocholine +($\beta$-BODIPY FL C5-HPC or BODIPY-PC) and 3,3' +dioctadecyloxacarbocyanine perchlorate (DiO) in their +fluorescence-detected linear dichroism (LD) studies of plasma +membranes of living cells.\cite{Benninger:2005qy} The DiO dye aligns +its transition moment perpendicular to the membrane normal, while the +BODIPY-PC transition dipole is parallel with the membrane normal. +Without a doubt, using fluorescence detection of linear dichroism in +concert with AFM surface scanning would be difficult experiments to +carry out. However, there is some hope of performing experiments to +either verify or falsify the predictions of our simulations. + Although our model is simple, it exhibits some rich and unexpected -behaviors. It would clearly be a closer approximation to the reality -if we allowed greater translational freedom to the dipoles and -replaced the somewhat artificial lattice packing and the harmonic -elastic tension with more realistic molecular modeling potentials. -What we have done is to present a simple model which exhibits bulk -non-thermal corrugation, and our explanation of this rippling +behaviors. It would clearly be a closer approximation to reality if +we allowed bending motions between the dipoles and the molecular +bodies, and if we replaced the rigid ellipsoids with ball-and-chain +tails. However, the advantages of this simple model (large system +sizes, 50 fs time steps) allow us to rapidly explore the phase diagram +for a wide range of parameters. Our explanation of this rippling phenomenon will help us design more accurate molecular models for -corrugated membranes and experiments to test whether rippling is -dipole-driven or not. - +corrugated membranes and experiments to test whether or not +dipole-dipole interactions exert an influence on membrane rippling. \newpage \bibliography{mdripple} \end{document}