25 |
|
\bibliographystyle{achemso} |
26 |
|
|
27 |
|
\title{Dipolar ordering in the ripple phases of molecular-scale models |
28 |
< |
of lipid membranes} |
28 |
> |
of lipid membranes} |
29 |
|
\author{Xiuquan Sun and J. Daniel Gezelter \\ |
30 |
|
Department of Chemistry and Biochemistry,\\ |
31 |
|
University of Notre Dame, \\ |
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 |
78 |
> |
within the gel phase,~\cite{Cevc87} although Tenchov {\it et al.} have |
79 |
> |
recently observed near-hexagonal packing in some phosphatidylcholine |
80 |
> |
(PC) gel phases.\cite{Tenchov2001} The X-ray diffraction work by |
81 |
|
Katsaras {\it et al.} showed that a rich phase diagram exhibiting both |
82 |
|
{\it asymmetric} and {\it symmetric} ripples is possible for lecithin |
83 |
|
bilayers.\cite{Katsaras00} |
85 |
|
A number of theoretical models have been presented to explain the |
86 |
|
formation of the ripple phase. Marder {\it et al.} used a |
87 |
|
curvature-dependent Landau-de~Gennes free-energy functional to predict |
88 |
< |
a rippled phase.~\cite{Marder84} This model and other related continuum |
89 |
< |
models predict higher fluidity in convex regions and that concave |
90 |
< |
portions of the membrane correspond to more solid-like regions. |
91 |
< |
Carlson and Sethna used a packing-competition model (in which head |
92 |
< |
groups and chains have competing packing energetics) to predict the |
93 |
< |
formation of a ripple-like phase. Their model predicted that the |
94 |
< |
high-curvature portions have lower-chain packing and correspond to |
95 |
< |
more fluid-like regions. Goldstein and Leibler used a mean-field |
96 |
< |
approach with a planar model for {\em inter-lamellar} interactions to |
97 |
< |
predict rippling in multilamellar phases.~\cite{Goldstein88} McCullough |
98 |
< |
and Scott proposed that the {\em anisotropy of the nearest-neighbor |
99 |
< |
interactions} coupled to hydrophobic constraining forces which |
100 |
< |
restrict height differences between nearest neighbors is the origin of |
101 |
< |
the ripple phase.~\cite{McCullough90} Lubensky and MacKintosh |
102 |
< |
introduced a Landau theory for tilt order and curvature of a single |
103 |
< |
membrane and concluded that {\em coupling of molecular tilt to membrane |
104 |
< |
curvature} is responsible for the production of |
105 |
< |
ripples.~\cite{Lubensky93} Misbah, Duplat and Houchmandzadeh proposed |
106 |
< |
that {\em inter-layer dipolar interactions} can lead to ripple |
107 |
< |
instabilities.~\cite{Misbah98} Heimburg presented a {\em coexistence |
108 |
< |
model} for ripple formation in which he postulates that fluid-phase |
109 |
< |
line defects cause sharp curvature between relatively flat gel-phase |
110 |
< |
regions.~\cite{Heimburg00} Kubica has suggested that a lattice model of |
111 |
< |
polar head groups could be valuable in trying to understand bilayer |
112 |
< |
phase formation.~\cite{Kubica02} Bannerjee used Monte Carlo simulations |
113 |
< |
of lamellar stacks of hexagonal lattices to show that large head groups |
88 |
> |
a rippled phase.~\cite{Marder84} This model and other related |
89 |
> |
continuum models predict higher fluidity in convex regions and that |
90 |
> |
concave portions of the membrane correspond to more solid-like |
91 |
> |
regions. Carlson and Sethna used a packing-competition model (in |
92 |
> |
which head groups and chains have competing packing energetics) to |
93 |
> |
predict the formation of a ripple-like phase. Their model predicted |
94 |
> |
that the high-curvature portions have lower-chain packing and |
95 |
> |
correspond to more fluid-like regions. Goldstein and Leibler used a |
96 |
> |
mean-field approach with a planar model for {\em inter-lamellar} |
97 |
> |
interactions to predict rippling in multilamellar |
98 |
> |
phases.~\cite{Goldstein88} McCullough and Scott proposed that the {\em |
99 |
> |
anisotropy of the nearest-neighbor interactions} coupled to |
100 |
> |
hydrophobic constraining forces which restrict height differences |
101 |
> |
between nearest neighbors is the origin of the ripple |
102 |
> |
phase.~\cite{McCullough90} Lubensky and MacKintosh introduced a Landau |
103 |
> |
theory for tilt order and curvature of a single membrane and concluded |
104 |
> |
that {\em coupling of molecular tilt to membrane curvature} is |
105 |
> |
responsible for the production of ripples.~\cite{Lubensky93} Misbah, |
106 |
> |
Duplat and Houchmandzadeh proposed that {\em inter-layer dipolar |
107 |
> |
interactions} can lead to ripple instabilities.~\cite{Misbah98} |
108 |
> |
Heimburg presented a {\em coexistence model} for ripple formation in |
109 |
> |
which he postulates that fluid-phase line defects cause sharp |
110 |
> |
curvature between relatively flat gel-phase regions.~\cite{Heimburg00} |
111 |
> |
Kubica has suggested that a lattice model of polar head groups could |
112 |
> |
be valuable in trying to understand bilayer phase |
113 |
> |
formation.~\cite{Kubica02} Bannerjee used Monte Carlo simulations of |
114 |
> |
lamellar stacks of hexagonal lattices to show that large head groups |
115 |
|
and molecular tilt with respect to the membrane normal vector can |
116 |
< |
cause bulk rippling.~\cite{Bannerjee02} |
116 |
> |
cause bulk rippling.~\cite{Bannerjee02} Recently, Kranenburg and Smit |
117 |
> |
described the formation of symmetric ripple-like structures using a |
118 |
> |
coarse grained solvent-head-tail bead model.\cite{Kranenburg2005} |
119 |
> |
Their lipids consisted of a short chain of head beads tied to the two |
120 |
> |
longer ``chains''. |
121 |
|
|
122 |
|
In contrast, few large-scale molecular modeling studies have been |
123 |
|
done due to the large size of the resulting structures and the time |
349 |
|
along dipole $i$ and $\mathbf{\hat r}_{ij}$ is the unit vector |
350 |
|
pointing along the inter-dipole vector $\mathbf{r}_{ij}$. |
351 |
|
|
352 |
+ |
Since the charge separation distance is so large in zwitterionic head |
353 |
+ |
groups (like the PC head groups), it would also be possible to use |
354 |
+ |
either point charges or a ``split dipole'' approximation, |
355 |
+ |
\begin{equation} |
356 |
+ |
V_{ij}({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf |
357 |
+ |
\hat{r}}_{ij})) = \frac{1}{{4\pi \epsilon_0 }}\left[ {\frac{{\mu _i \cdot \mu _j }}{{R_{ij}^3 }} - |
358 |
+ |
\frac{{3\left( {\mu _i \cdot r_{ij} } \right)\left( {\mu _i \cdot |
359 |
+ |
r_{ij} } \right)}}{{R_{ij}^5 }}} \right] |
360 |
+ |
\end{equation} |
361 |
+ |
where $\mu _i$ and $\mu _j$ are the dipole moments of sites $i$ and |
362 |
+ |
$j$, $r_{ij}$ is vector between the two sites, and $R_{ij}$ is given |
363 |
+ |
by, |
364 |
+ |
\begin{equation} |
365 |
+ |
R_{ij} = \sqrt {r_{ij}^2 + \frac{{d_i^2 }}{4} + \frac{{d_j^2 |
366 |
+ |
}}{4}}. |
367 |
+ |
\end{equation} |
368 |
+ |
Here, $d_i$ and $d_j$ are effect charge separation distances |
369 |
+ |
associated with each of the two dipolar sites. This approximation to |
370 |
+ |
the multipole expansion maintains the fast fall-off of the multipole |
371 |
+ |
potentials but lacks the normal divergences when two polar groups get |
372 |
+ |
close to one another. |
373 |
+ |
|
374 |
|
For the interaction between nonequivalent uniaxial ellipsoids (in this |
375 |
|
case, between spheres and ellipsoids), the spheres are treated as |
376 |
|
ellipsoids with an aspect ratio of 1 ($d = l$) and with an well depth |
379 |
|
et al.} and is appropriate for dissimilar uniaxial |
380 |
|
ellipsoids.\cite{Cleaver96} |
381 |
|
|
382 |
< |
The solvent model in our simulations is identical to one used by |
383 |
< |
Marrink {\it et al.} in their dissipative particle dynamics (DPD) |
384 |
< |
simulation of lipid bilayers.\cite{Marrink04} This solvent bead is a |
385 |
< |
single site that represents four water molecules (m = 72 amu) and has |
386 |
< |
comparable density and diffusive behavior to liquid water. However, |
387 |
< |
since there are no electrostatic sites on these beads, this solvent |
388 |
< |
model cannot replicate the dielectric properties of water. |
382 |
> |
The solvent model in our simulations is similar to the one used by |
383 |
> |
Marrink {\it et al.} in their coarse grained simulations of lipid |
384 |
> |
bilayers.\cite{Marrink04} The solvent bead is a single site that |
385 |
> |
represents four water molecules (m = 72 amu) and has comparable |
386 |
> |
density and diffusive behavior to liquid water. However, since there |
387 |
> |
are no electrostatic sites on these beads, this solvent model cannot |
388 |
> |
replicate the dielectric properties of water. Note that although we |
389 |
> |
are using larger cutoff and switching radii than Marrink {\it et al.}, |
390 |
> |
our solvent density at 300 K remains at 0.944 g cm$^{-3}$, and the |
391 |
> |
solvent diffuses at 0.43 $\AA^2 ps^{-1}$ (roughly twice as fast as |
392 |
> |
liquid water). |
393 |
|
|
394 |
|
\begin{table*} |
395 |
|
\begin{minipage}{\linewidth} |
447 |
|
solvent beads (24 water molecules) per lipid. These configurations |
448 |
|
were then equilibrated for another $30$ ns. All simulations utilizing |
449 |
|
the solvent were carried out at constant pressure ($P=1$ atm) with |
450 |
< |
$3$D anisotropic coupling, and constant surface tension |
450 |
> |
$3$D anisotropic coupling, and small constant surface tension |
451 |
|
($\gamma=0.015$ N/m). Given the absence of fast degrees of freedom in |
452 |
|
this model, a time step of $50$ fs was utilized with excellent energy |
453 |
|
conservation. Data collection for structural properties of the |
454 |
|
bilayers was carried out during a final 5 ns run following the solvent |
455 |
< |
equilibration. All simulations were performed using the OOPSE |
456 |
< |
molecular modeling program.\cite{Meineke05} |
455 |
> |
equilibration. Orientational correlation functions and diffusion |
456 |
> |
constants were computed from 30 ns simulations in the microcanonical |
457 |
> |
(NVE) ensemble using the average volume from the end of the constant |
458 |
> |
pressure and surface tension runs. The timestep on these final |
459 |
> |
molecular dynamics runs was 25 fs. No appreciable changes in phase |
460 |
> |
structure were noticed upon switching to a microcanonical ensemble. |
461 |
> |
All simulations were performed using the {\sc oopse} molecular |
462 |
> |
modeling program.\cite{Meineke05} |
463 |
|
|
464 |
|
A switching function was applied to all potentials to smoothly turn |
465 |
< |
off the interactions between a range of $22$ and $25$ \AA. |
465 |
> |
off the interactions between a range of $22$ and $25$ \AA. The |
466 |
> |
switching function was the standard (cubic) function, |
467 |
> |
\begin{equation} |
468 |
> |
s(r) = |
469 |
> |
\begin{cases} |
470 |
> |
1 & \text{if $r \le r_{\text{sw}}$},\\ |
471 |
> |
\frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2} |
472 |
> |
{(r_{\text{cut}} - r_{\text{sw}})^3} |
473 |
> |
& \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\ |
474 |
> |
0 & \text{if $r > r_{\text{cut}}$.} |
475 |
> |
\end{cases} |
476 |
> |
\label{eq:dipoleSwitching} |
477 |
> |
\end{equation} |
478 |
|
|
479 |
|
\section{Results} |
480 |
|
\label{sec:results} |
785 |
|
($\mu$).\label{fig:phaseDiagram}} |
786 |
|
\end{figure} |
787 |
|
|
788 |
< |
We have computed translational diffusion coefficients for lipid |
789 |
< |
molecules from the mean square displacement, |
790 |
< |
\begin{eqnarray} |
791 |
< |
\langle {|\left({\bf r}_{i}(t) - {\bt r}_{i}(0) \right)|}^2 \rangle \\ \\ |
792 |
< |
& = & 6Dt |
793 |
< |
\end{eqnarray} |
794 |
< |
of the lipid bodies. The values of the translational diffusion |
795 |
< |
coefficient for different head-to-tail size ratio are shown in table |
796 |
< |
\ref{tab:relaxation}. |
797 |
< |
|
798 |
< |
We have also computed orientational diffusion constants for the head |
799 |
< |
groups from the relaxation of the second-order Legendre polynomial |
800 |
< |
correlation function, |
801 |
< |
\begin{eqnarray} |
802 |
< |
C_{\ell}(t) & = & \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf |
803 |
< |
\mu}_{i}(0) \right) \rangle \\ \\ |
804 |
< |
& \approx & e^{-\ell(\ell + 1) \theta t}, |
805 |
< |
\end{eqnarray} |
806 |
< |
of the head group dipoles. In this last line, we have used a simple |
807 |
< |
``Debye''-like model for the relaxation of the correlation function, |
808 |
< |
specifically in the case when $\ell = 2$. The computed orientational |
809 |
< |
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... |
788 |
> |
We have computed translational diffusion constants for lipid molecules |
789 |
> |
from the mean-square displacement, |
790 |
> |
\begin{equation} |
791 |
> |
D = \lim_{t\rightarrow \infty} \frac{1}{6 t} \langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \rangle, |
792 |
> |
\end{equation} |
793 |
> |
of the lipid bodies. Translational diffusion constants for the |
794 |
> |
different head-to-tail size ratios (all at 300 K) are shown in table |
795 |
> |
\ref{tab:relaxation}. We have also computed orientational correlation |
796 |
> |
times for the head groups from fits of the second-order Legendre |
797 |
> |
polynomial correlation function, |
798 |
> |
\begin{equation} |
799 |
> |
C_{\ell}(t) = \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf |
800 |
> |
\mu}_{i}(0) \right) |
801 |
> |
\end{equation} |
802 |
> |
of the head group dipoles. The orientational correlation functions |
803 |
> |
appear to have multiple components in their decay: a fast ($12 \pm 2$ |
804 |
> |
ps) decay due to librational motion of the head groups, as well as |
805 |
> |
moderate ($\tau^h_{\rm mid}$) and slow ($\tau^h_{\rm slow}$) |
806 |
> |
components. The fit values for the moderate and slow correlation |
807 |
> |
times are listed in table \ref{tab:relaxation}. Standard deviations |
808 |
> |
in the fit time constants are quite large (on the order of the values |
809 |
> |
themselves). |
810 |
|
|
811 |
+ |
Sparrman and Westlund used a multi-relaxation model for NMR lineshapes |
812 |
+ |
observed in gel, fluid, and ripple phases of DPPC and obtained |
813 |
+ |
estimates of a correlation time for water translational diffusion |
814 |
+ |
($\tau_c$) of 20 ns.\cite{Sparrman2003} This correlation time |
815 |
+ |
corresponds to water bound to small regions of the lipid membrane. |
816 |
+ |
They further assume that the lipids can explore only a single period |
817 |
+ |
of the ripple (essentially moving in a nearly one-dimensional path to |
818 |
+ |
do so), and the correlation time can therefore be used to estimate a |
819 |
+ |
value for the translational diffusion constant of $2.25 \times |
820 |
+ |
10^{-11} m^2 s^{-1}$. The translational diffusion constants we obtain |
821 |
+ |
are in reasonable agreement with this experimentally determined |
822 |
+ |
value. However, the $T_2$ relaxation times obtained by Sparrman and |
823 |
+ |
Westlund are consistent with P-N vector reorientation timescales of |
824 |
+ |
2-5 ms. This is substantially slower than even the slowest component |
825 |
+ |
we observe in the decay of our orientational correlation functions. |
826 |
+ |
Other than the dipole-dipole interactions, our head groups have no |
827 |
+ |
shape anisotropy which would force them to move as a unit with |
828 |
+ |
neighboring molecules. This would naturally lead to P-N reorientation |
829 |
+ |
times that are too fast when compared with experimental measurements. |
830 |
|
|
831 |
|
\begin{table*} |
832 |
|
\begin{minipage}{\linewidth} |
833 |
|
\begin{center} |
834 |
< |
\caption{Rotational diffusion constants for the head groups |
835 |
< |
($\theta_h$) and molecular bodies ($\theta_b$) as well as the |
836 |
< |
translational diffusion coefficients for the molecule as a function of |
837 |
< |
the head-to-body width ratio. The orientational mobility of the head |
838 |
< |
groups experiences an {\it order of magnitude decrease} upon entering |
839 |
< |
the rippled phase, which suggests that the rippling is tied to a |
840 |
< |
freezing out of head group orientational freedom. Uncertainties in |
841 |
< |
the last digit are indicated by the values in parentheses.} |
842 |
< |
\begin{tabular}{lccc} |
834 |
> |
\caption{Fit values for the rotational correlation times for the head |
835 |
> |
groups ($\tau^h$) and molecular bodies ($\tau^b$) as well as the |
836 |
> |
translational diffusion constants for the molecule as a function of |
837 |
> |
the head-to-body width ratio (all at 300 K). In all of the phases, |
838 |
> |
the head group correlation functions decay with an fast librational |
839 |
> |
contribution ($12 \pm 1$ ps). There are additional moderate |
840 |
> |
($\tau^h_{\rm mid}$) and slow $\tau^h_{\rm slow}$ contributions to |
841 |
> |
orientational decay that depend strongly on the phase exhibited by the |
842 |
> |
lipids. The symmetric ripple phase ($\sigma_h / d = 1.35$) appears to |
843 |
> |
exhibit the slowest molecular reorientation.} |
844 |
> |
\begin{tabular}{lcccc} |
845 |
|
\hline |
846 |
< |
$\sigma_h / d$ & $\theta_h (\mu s^{-1})$ & $\theta_b (1/fs)$ & $D ( |
847 |
< |
\times 10^{-11} m^2 s^{-1} \\ |
846 |
> |
$\sigma_h / d$ & $\tau^h_{\rm mid} (ns)$ & $\tau^h_{\rm |
847 |
> |
slow} (\mu s)$ & $\tau_b (\mu s)$ & $D (\times 10^{-11} m^2 s^{-1})$ \\ |
848 |
|
\hline |
849 |
< |
1.20 & $0.206(1) $ & $0.0175(5) $ & $0.43(1)$ \\ |
850 |
< |
1.28 & $0.179(2) $ & $0.055(2) $ & $5.91(3)$ \\ |
851 |
< |
1.35 & $0.025(1) $ & $0.195(3) $ & $3.42(1)$ \\ |
852 |
< |
1.41 & $0.023(1) $ & $0.024(3) $ & $7.16(1)$ \\ |
849 |
> |
1.20 & $0.4$ & $9.6$ & $9.5$ & $0.43(1)$ \\ |
850 |
> |
1.28 & $2.0$ & $13.5$ & $3.0$ & $5.91(3)$ \\ |
851 |
> |
1.35 & $3.2$ & $4.0$ & $0.9$ & $3.42(1)$ \\ |
852 |
> |
1.41 & $0.3$ & $23.8$ & $6.9$ & $7.16(1)$ \\ |
853 |
|
\end{tabular} |
854 |
|
\label{tab:relaxation} |
855 |
|
\end{center} |
879 |
|
is convex. These structures are held together by the extremely strong |
880 |
|
and directional interactions between the head groups. |
881 |
|
|
882 |
< |
Dipolar head groups are key for the maintaining the bilayer structures |
883 |
< |
exhibited by this model. The dipoles are likely to form head-to-tail |
884 |
< |
configurations even in flat configurations, but the temperatures are |
885 |
< |
high enough that vortex defects become prevalent in the flat phase. |
886 |
< |
The flat phase we observed therefore appears to be substantially above |
887 |
< |
the Kosterlitz-Thouless transition temperature for a planar system of |
888 |
< |
dipoles with this set of parameters. For this reason, it would be |
889 |
< |
interesting to observe the thermal behavior of the flat phase at |
890 |
< |
substantially lower temperatures. |
882 |
> |
The attractive forces holding the bilayer together could either be |
883 |
> |
non-directional (as in the work of Kranenburg and |
884 |
> |
Smit),\cite{Kranenburg2005} or directional (as we have utilized in |
885 |
> |
these simulations). The dipolar head groups are key for the |
886 |
> |
maintaining the bilayer structures exhibited by this particular model; |
887 |
> |
reducing the strength of the dipole has the tendency to make the |
888 |
> |
rippled phase disappear. The dipoles are likely to form attractive |
889 |
> |
head-to-tail configurations even in flat configurations, but the |
890 |
> |
temperatures are high enough that vortex defects become prevalent in |
891 |
> |
the flat phase. The flat phase we observed therefore appears to be |
892 |
> |
substantially above the Kosterlitz-Thouless transition temperature for |
893 |
> |
a planar system of dipoles with this set of parameters. For this |
894 |
> |
reason, it would be interesting to observe the thermal behavior of the |
895 |
> |
flat phase at substantially lower temperatures. |
896 |
|
|
897 |
|
One feature of this model is that an energetically favorable |
898 |
|
orientational ordering of the dipoles can be achieved by forming |