Membrane lateral structure: The influence of immobilized particles on domain size

Membrane lateral structure: The influence of immobilized particles on domain size

Timo Fischer Georg-August-Universität Göttingen, Institute of Theoretical Physics, Friedrich-Hund-Platz 1, D-37077 Göttingen, Germany    H. Jelger Risselada Max Planck Institute for Biophysical Chemistry, Department of Theoretical and Computational Biophysics, Am Faßberg 11, D-37077 Göttingen, Germany.    Richard L. C. Vink Georg-August-Universität Göttingen, Institute of Theoretical Physics, Friedrich-Hund-Platz 1, D-37077 Göttingen, Germany

In experiments on model membranes, a formation of large domains of different lipid composition is readily observed. However, no such phase separation is observed in the membranes of intact cells. Instead, a structure of small transient inhomogeneities called lipid rafts are expected in these systems. One of the numerous attempts to explain small domains refers to the coupling of the membrane to its surroundings, which leads to the immobilization of some of the membrane molecules. These immobilized molecules then act as static obstacles for the remaining mobile ones. We present detailed Molecular Dynamics simulations demonstrating that this can indeed account for small domains. This confirms previous Monte Carlo studies based on simplified models. Furthermore, by directly comparing domain structures obtained using Molecular Dynamics to Monte Carlo simulations of the Ising model, we demonstrate that domain formation in the presence of obstacles is remarkably insensitive to the details of the molecular interactions.

65.20.De, 87.16.dt, 64.60.De, 64.70.Ja, 61.20.Ja

I Introduction

Membranes are two-dimensional (2D) fluid environments, consisting of many different types of lipids and proteins citeulike:1595800 (). The membrane constituents are not arranged randomly, but are expected to be spatially organized into domains of different size, composition, and dynamics citeulike:5470515 (); citeulike:4308274 (); citeulike:4312704 (); citeulike:3042594 (); citeulike:6562069 (). In biological membranes this is important because it links to key processes in cells, such as signaling, endocytosis, and adhesion citeulike:6505558 (); citeulike:6499105 (). In artificial membranes, domain formation is relevant for applications, ranging from photolithographic patterning, spatial addressing, microcontact printing, and microfluidic patterning citeulike:6528926 (); citeulike:6549481 (). To identify the factors that control domain formation, and to understand their underlying physical mechanisms, is therefore of key practical importance.

One challenge is to account for a heterogeneous domain structure in thermal equilibrium citeulike:5470515 (). Due to line tension, one would naively expect a multi-domain structure to be unstable as the free energy would be reduced by domain coalescence. Indeed, in model membranes, this is precisely what is observed citeulike:3850881 (); Veatch20033074 (). At high temperature, these systems are in a mixed state, but they phase separate into (macroscopically large) liquid-ordered and liquid-disordered domains at low temperature. At the temperature where the phase separation begins to occur (the so-called critical temperature), the transition passes through a critical point citeulike:3850881 (); citeulike:3367620 (); citeulike:3850871 (). At the critical point, the membrane is highly heterogeneous, featuring domains of all sizes. Hence, at least in model membranes, a non-trivial heterogeneous domain structure arises only near the critical point. In line with this, it has been proposed that critical behavior might also explain a heterogeneous domain structure in biological membranes citeulike:3852117 (). We note that critical behavior is not the only feasible mechanism. Also the coupling between the elastic properties of the membrane and its composition has been identified to be a key factor affecting domain shape and size citeulike:8893054 (); schick:2011 (); citeulike:5045047 (); citeulike:9473614 (); citeulike:2802537 (); citeulike:4025142 (). In addition, the presence of hybrid lipids that collect at the interface between liquid-ordered and liquid-disordered domains could also induce a heterogeneous domain structure citeulike:10418449 ().

i.1 Particle immobilization in membranes

Figure 1: Schematic sketches (not to scale) showing possible origins of quenched disorder in membranes. (a) In biological membranes, quenched disorder may arise from the anchoring of proteins and other molecules to the cytoskeleton network. Such anchored proteins are effectively immobilized. The typical distance between immobilized particles can become very small  nm citeulike:6859704 (). (b) An analogous in vitro example is provided by supported membranes. In this case, particle immobilization may arise from surface roughness or, in the case of larger molecules, from surface friction. (c) In both examples (a) and (b), a top-view of the membrane bilayer qualitatively resembles a 2D array of immobilized obstacles (Q), through which the mobile particles (M) diffuse.

The above mechanisms share in common that they do not require the membrane to be coupled to its environment in any specific way. However, an essential difference between model membranes and real cells may well be the fact that the latter are coupled to their surroundings. For example, the cytoskeleton of a cell may induce a tension on the membrane affecting its domain structure uline:2012 (). Another possible effect – the one that we focus on in this work – is that the coupling of the membrane to its environment causes some of its components to become immobilized. A possible mechanism leading to particle immobilization in a real cell may be an actin network underlying the membrane. The actin network is attached to the membrane via anchoring proteins and, since the actin network is disordered, their positions will be random. The anchoring proteins are essentially immobile compared to the freely diffusing membrane components, and thus act as obstacles for the mobile components [fig. 1(a)]. In addition to the anchoring proteins themselves, also other molecules can attach to the actin strands, which thereby also become obstacles. In fact, the distance between obstacles anchored to the actin strands can become very small  nm citeulike:6859704 ().

The immobilization of molecules can also arise in vitro in supported membranesciteulike:6609138 (); citeulike:6528926 (); citeulike:6116695 (); citeulike:5806549 (); citeulike:6127454 (); citeulike:6598182 (); citeulike:4071068 (). Fluorescence measurements have revealed that surface roughness profoundly limits lipid diffusion, even leading to complete immobilization. For example, the immobilized lipid fraction is around 15% on alumina substrates, and 5% on silica citeulike:6116695 (). Furthermore, due to surface friction, embedded proteins may also become immobilized citeulike:6127454 () [fig. 1(b)]. Hence, both in biological membranes as well as in artificially supported lipid bilayers, lateral diffusion may be envisioned as taking place inside a random network of static obstacles [fig. 1(c)]. In physical terms, the presence of a random network of static obstacles constitutes a form of quenched disorder.

i.2 Quenched disorder in membranes:
Summary of known results

The above examples show that the presence of quenched disorder in membranes is inevitable in many practical situations. This naturally raises the question to what extent this can account for the properties of a membrane. Regarding single particle diffusion, it is known that lateral diffusion constants are drastically reduced in the presence of obstacles citeulike:918597 (); citeulike:6505316 (). The diffusion constant of band 3 in mouse erythrocytes without cytoskeleton, i.e. in the absence of quenched disorder, is 50 times larger compared to that in healthy cells citeulike:6505316 (). Additionally, the diffusion becomes anomalous (non-Brownian) on intermediate time-scales. These results can be explained by assuming that the quenched obstacles divide the membrane into compartments citeulike:6859704 (); citeulike:918597 (); citeulike:9836137 (); citeulike:5393528 (); citeulike:6505558 (); citeulike:6195415 (); citeulike:6116283 (); citeulike:4304224 (); citeulike:139732 (). Within a single compartment, the diffusion is fast and Brownian, while diffusion between compartments is governed by a much slower “hopping” dynamics citeulike:918597 ().

In addition to single particle diffusion, quenched obstacles also affect the collective behavior of the mobile constituents. By collective behavior we mean in this context the size and shape of the liquid-ordered and liquid-disordered domains that can form in the membrane. Precisely this question was addressed in a series of recent Monte Carlo (MC) simulations citeulike:7115548 (); citeulike:6599228 (); citeulike:8864903 (); Machta20111668 (); Ehrig201180 (). A remarkable finding is that, provided (1) the quenched obstacles are randomly distributed, and (2) show a preferred attraction to one of the mobile components, macroscopic domain formation is entirely prevented citeulike:8864903 (). This behavior directly connects to the fundamentals of fluid phase behavior in quenched porous media, for which numerous theoretical results are available madden.glandt:1988 (); citeulike:4892193 (); citeulike:10004339 (). In particular, when conditions (1) and (2) above are met, the obstacles induce a change in the universality class, from Ising toward random-field Ising gennes:1984 (). As is well known, the latter does not support macroscopic domain formation in two dimensions (); citeulike:2841682 (). Hence, based on the fundamentals of statistical physics, particle immobilization provides a robust mechanism to account for a heterogeneous membrane domain structure.

i.3 Aim of paper

In this paper, we put the recent MC findings citeulike:8864903 (); Machta20111668 (); Ehrig201180 () pertaining to particle immobilization and its effect on membrane domain formation to a stringent test. It should be noted that these previous studies are based on rather simple models: all exclude membrane height fluctuations, most are defined on a lattice, lipid tail degrees of freedom are either ignored or only very crudely modeled, while solvent and cholesterol are absent. In addition, since all deployed the MC method, the dynamics could only be approximated or was omitted altogether. Hence, it is important to verify if the observed trends in the MC studies also survive in more realistic membrane models. The aim of this paper is to provide this verification. We present a high-resolution Molecular Dynamics (MD) study of a phase-separating DPPC/DLiPC/cholesterol lipid bilayer in the presence of an immobilized component, which captures important details left out in earlier works.

The outline of our paper is as follows: In Section II, we present our MD simulation results, and show that the presence of quenched obstacles in a membrane indeed results in a heterogeneous structure of small domains. In Section III, we then compare our MD results to the predictions of a simple MC model (the 2D Ising model), and demonstrate that the latter is already well-suited to describe domain structures. We end with a summary and conclusion in Section IV.

Ii Molecular Dynamics simulations

ii.1 MD simulation setup

For the MD simulations presented here, a setup from a previous workciteulike:3483967 () was used, adapted to our needs by changing the size and shape of the membrane patch and adding quenched disorder to the system. In the setup without disorder, which is simulated first for comparison, the membrane consists of cholesterol molecules and two species of phospholipids, namely molecules of dipalmitoyl-phosphatidylcholine with fully saturated carbon tails (DPPC), and molecules of the polyunsaturated phospholipid dilinoleyl-phosphatidylcholine (DLiPC). This composition lies deep within the coexistence region of the phase diagram VeatchPhasediagram (); citeulike:5629089 (); citeulike:9049927 (). Hence, our lipids have a strong tendency to demix.

Run Type of obstacles
I None  s
IIa Same on both sides, mobile in -direction  s
IIb Same on both sides, mobile in -direction  s
III On one side only, pinned in -direction  s
Table 1: Overview of the four MD simulation runs performed in this work, listing the immobilized obstacle type, and the total simulation time . Runs IIa and IIb use the same obstacle configuration, but different random numbers for the initialization of the mobile lipids were used. Note that diffusion rates in the MARTINI force field are times larger compared to atomistic simulations and experiments doi:10.1021/jp071097f (). Values for time given in this work refer to the simulation time coordinate, not to the rescaled (i.e. larger) “physical time”.

In the initial state, the lipids form a flat bilayer patch in the plane, using a periodic simulation box of size  nm with the same number of DPPC and DLiPC lipids in both layers. Within each layer, the phospholipids are initially randomly mixed. The molecules are modeled with the MARTINI force fielddoi:10.1021/jp071097f (), which provides a near-atomic resolution of - non-hydrogen atoms per simulation bead. In order to capture the effect of the solvent, solvent beads were included. The simulation has been performed with GROMACS JCC:JCC20291 (); doi:10.1021/ct700301q () Version  at temperature  K and a pressure of  bar (for further details we refer the reader to the original setup citeulike:3483967 ()). In total, four different simulation runs have been performed, summarized in Table 1, and to be discussed further in the following sections.

ii.2 The free membrane

We first consider the “free” membrane (run I in Table 1), by which we mean a membrane without any quenched disorder. In the absence of quenched disorder, model ternary membranes of phospholipids and cholesterol readily phase separate into macroscopic (m) domains Veatch20033074 (). The domains observed in experiments are typically round, as this shape minimizes the length of the line interface. The key difference with computer simulations using periodic boundary conditions is that the equilibrium domain shape is not necessarily round. This is a crucial point in our analysis and to avoid possible confusion with related studies Ehrig201180 (); citeulike:9234255 () we explain it explicitly. In a 2D system with periodic boundary conditions, the two candidate domain shapes are (1) a circle, or (2) a stripe that spans the system along the shorter edge of the simulation box citeulike:7134501 (); citeulike:7237424 (). The shape that prevails is the one with the shortest line interface, which depends on the area fraction of the two coexisting phases and the aspect ratio of the simulation box. For the simulations presented here, the hallmark of macroscopic domain formation in a finite simulation box with periodic boundaries is the appearance of the latter stripe geometry. Consequently, in our MD simulations, the appearance of such a stripe domain will be regarded as evidence for macroscopic domain formation.

Figure 2: Membrane without quenched disorder (run I) at  s, showing the upper leaflet. Due to inter-leaflet coupling, the domain structure on the lower leaflet is almost identical. From the initially mixed state, two domains of DPPC (dark blue) and cholesterol (yellow) have formed in a background of DLiPC (bright brown). The DPPC domain in the foreground has assumed the stripe shape, which is the simulation analogue of large round domains seen in experiments.

We started the MD simulation of the free membrane with a mixed configuration. During the first few s of simulating the membrane, we observe the characteristics of binary demixing, namely the formation, coarsening, and merging of domains. These domains reflect the liquid-ordered phase (rich in DPPC and cholesterol) and the liquid-disordered phase (rich in DLiPC). Due to inter-leaflet couplingciteulike:3483967 (), domains in both leaflets arrange in similar shapes and at the same locations. After  s, one of the domains has assumed the stripe geometry, which remains stable for the rest of the simulation [fig. 2]. Following the explanation above, this indicates macroscopic demixing into large DPPC-rich and DLiPC-rich domains, such as seen experimentally in free-standing giant unilamellar vesicles Bagatolli2000290 (); Veatch20033074 (). Ultimately, the small DPPC domain on the right of fig. 2 will merge with the stripe, but this process is very slow and beyond our computational resources.

ii.3 Membranes with quenched disorder

Figure 3: “Early time” domain structures after  s (a),  s (b), and  s (c) as observed in simulation run IIa (cf. Table 1). The obstacles are shown as white; the remaining color-coding is the same as in fig. 2.
Figure 4: “Late time” domain structures of runs IIa (a) and IIb (b) at  s. FPPC molecules are shown as white and are at the same locations in both systems (the slightly different positions of the white beads between the snapshots result from only the phosphate atom of the FPPC being quenched). Otherwise, the color-coding is the same as in fig. 2. In both simulations, the stripe domain characterizing macroscopic domain formation is absent. Instead, structures consisting of small domains are seen.

To demonstrate how quenched disorder affects structure formation in membranes, obstacles have been added to the free membrane of Section II.2. Contrary to previous studies of membranes with quenched disorderciteulike:7115548 (); citeulike:6599228 (); citeulike:8864903 (); Machta20111668 (); Ehrig201180 () that only took into account a flat monolayer our model takes into account the full bilayer and allows for membrane undulations. Hence, there are two additional degrees of freedom we can investigate. First of all, the obstacles may be allowed to follow the membrane undulations in the -direction, or be fixed in this direction. Secondly, the obstacles may be placed at the same lateral locations in both leaflets, which may be conceived to describe a trans-membrane protein, or they could be chosen to reside on only one leaflet.

ii.3.1 Obstacles that follow membrane fluctuations

We first consider the case whereby the obstacles have been placed at the same locations in both leaflets and are allowed to move freely in the -direction. In this scenario, the obstacles correspond to trans-membrane proteins that are free to move in the -direction themselves, or so long that the membrane can fluctuate around them. To model the obstacles, we chose to take DPPC molecules that are not allowed to diffuse laterally. Such obstacles, which we refer to as FPPC molecules in what follows, naturally attract the DPPC-rich phase. The motivation to use DPPC obstacles stems from experimental considerations: cytoskeleton anchor proteins typically favor the liquid-ordered phase citeulike:10759324 (); Ehrig201180 (). In principle, we could have chosen DLiPC-affine obstacles or a mix of DPPC-affine and DLiPC-affine obstacles. Similar results as those presented here would be expected in these cases also.

To generate an initial state of a membrane with obstacles, we take the initial state from Section II.2 and chose random lateral locations in the membrane. Then, in each leaflet the nearest phosphate atoms belonging to a DPPC molecule are identified and moved to the target positions by a series of small translations and subsequent energy minimizations (so that the DPPC molecules remain intact and that there is no overlap between lipids). Finally, the so-moved DPPC molecules are labeled FPPC and the -coordinates of their phosphate atoms are bound to their target locations via a hard harmonic potential. During the course of the simulation, the lateral positions to which the phosphates of the FPPC are bound scale with the volume changes of the simulation box. Otherwise they are static, effectively making the FPPC quenched.

Two simulation runs with the same FPPC obstacle configuration have been performed (IIa and IIb of Table 1), both starting from a state of “mixed” mobile lipids. In both simulations, a formation and coarsening of small DPPC-rich domains is seen during the first few microseconds (as shown in fig. 3 for run IIa). After approximately  s, however, domain growth and merging of domains into larger ones is arrested. Beyond this time, the domains exhibit shape fluctuations, but do not significantly change in size or location for the rest of the simulation. Most importantly, as shown in fig. 4, neither of the simulation runs reveals the stripe geometry seen in fig. 2 that would indicate domains larger than the size of the simulation box. We conclude that the obstacles indeed prevent macroscopic demixing, and replace it with the formation of small domains that appear stable over times at least in the order of  s.

In recent MC studies, it was shown that quenched obstacles induce regions in the membrane where either a liquid-ordered or a liquid-disordered domain is preferredciteulike:8864903 (); Machta20111668 (). At zero temperature this would essentially fix the domain structure into the groundstate configuration, while at finite temperature thermal fluctuations still allow some freedom as to where the domains form. This explains why, despite identical FPPC locations, the domain structures in fig. 4(a) and (b) are not the same. We return to this point in detail in Section III.

ii.3.2 Obstacles on a single leaflet

Figure 5: Top-view of a membrane with FPPC obstacles in only one bilayer leaflet (run III) at  s. For clarity, this figure does not show the DLiPC or cholesterol molecules. DPPC belonging to the leaflet with the FPPC are drawn in blue, those of the “free” leaflet in yellow. The domains form “in registry” in both leaflets, but they are larger compared to the case where FPPC resides in both leaflets. The structure is dominated by a single large domain in both leaflets, but it is distinctively different from the stripe geometry that would indicate macroscopic domains.

We now consider the case where the FPPC obstacles reside in only one of the bilayer leaflets, and do not follow the membrane height fluctuations (run III). Such obstacles thus locally “pin” the membrane height, which could mimic a cytoskeleton anchoring site or, in the case of an adhered membrane, a (stiff) ligand-receptor bond citeulike:10571682 (). This scenario is physically interesting because it probes to what extent inter-leaflet coupling transmits the influence of quenched disorder from one leaflet to the other. To address this scenario in our MD simulations, we perform the same setup steps as in Section II.3.1 with the same target locations for the obstacles, but we place the obstacles on only one of the leaflets and also bind the -coordinate of the FPPC phosphate atoms to the hard harmonic potential.

Our simulations indicate that inter-leaflet coupling is indeed an important mechanism in domain formation, which is consistent with previous studies putzel:2011 (); schick:2011 (). During the first few microseconds, the formation, coarsening, and merging of domains is observed in both leaflets. The domains form “in registry” in both leaflets, occurring at the same (lateral) positions. After  s, most of the DPPC molecules are contained in a single domain as shown in fig. 5. We thus obtain larger domains compared to the case of FPPC residing in both leaflets. At the same time, we do not obtain the stripe geometry, suggesting that the domains are not macroscopic. Fig. 5 thus presents a “hybrid” structure, between the stripe domain of fig. 2 and the small domains of fig. 4. This result shows that the effect of quenched disorder is transmitted from the leaflet with FPPC to the one without, but also indicates that the tendency of the leaflet without FPPC to form large domains “back-transmits”. The net result is a structure of finite-sized domains (i.e. not macroscopic) but of a typical size that exceeds the case where FPPC resides in both leaflets.

ii.4 Interface properties of domains

Figure 6: Ratio of the interface length and domain area of the largest domain for the free membrane (run I), the membrane with FPPC in both leaflets (runs IIa and IIb), and for a membrane with FPPC in only one leaflet (run III). For each run, results are shown for each leaflet separately (hence 8 curves in total). Already at early times, for the free membrane differs from the membranes with quenched disorder, and saturates at a significantly lower value. All curves shown are interpolations of the raw simulation data. For clarity, the raw data are only shown explicitly for one curve. The scatter in the raw data of the other curves is comparable.

We now consider a more quantitative measure to distinguish membranes with quenched disorder from those without. To this end, we consider the ratio of the domain interface length to domain area . In the free membrane, the formation of macroscopic domains is driven by the desire of the system to minimize the interface length, yielding a low value of . In the membrane with quenched disorder, where stable small domains are expected, will be larger. To measure , a grid of rectangles is superimposed on the membrane. The positions of the phosphate atoms of the DPPC molecules are then projected onto the grid. If a grid cell contains the phosphate atom of a DPPC molecule, then that cell is considered to be part of a DPPC-rich domain. Cells with less than four neighbors belonging to a DPPC-rich domain are defined as interface.

In fig. 6, the ratio of the number of interface squares to the number of total squares for the largest DPPC-rich domains are shown for both leaflets and for all our simulations. In all cases, the curves first show a rapid drop resulting from the initial formation of domains, followed by a saturation as the system equilibrates. The saturation values for the free membrane, however, are distinctly below those of the membranes with quenched disorder. This is consistent with macroscopic domain formation in the former, and stable finite domains in the latter. Note also that, for the membrane with FPPC obstacles in only one leaflet, is largest in the FPPC containing leaflet.

ii.5 Choice of order parameter

In our analysis, domains were identified by their phospholipid content. An alternative choice is to identify domains via the tail ordering of the phospholipids, i.e. the identification of liquid-ordered and liquid-disordered domains citeulike:9234255 (). To this end, we have also quantified the elongation of the phospholipid tails via a nematic order parameter. DPPC and DLiPC lipids typically have high and low nematic order parameters, respectively. Simulation snapshots in which the phospholipids are color-coded according to their nematic order parameter reveal the same domain structures as those presented here. The identification of domains according to lipid content or tail conformation is therefore equivalent, and will not be explicitly presented here.

Iii Comparison to the 2D Ising model

In Section II, we have shown using a detailed membrane model that quenched obstacles prevent the formation of large (macroscopic) domains. This confirms the validity of the trends observed in earlier MC studies on simple membrane models citeulike:6599228 (); citeulike:8864903 (); Machta20111668 (); Ehrig201180 (). One issue that remains open is whether the apparent elimination of macroscopic domains really is the correct equilibrium behavior, or merely caused by the obstacles slowing down the dynamics of domain formation such that the relevant timescales become inaccessible in the simulation. Because of their high computational demand111Each of the MD runs of Table 1 required several months of simulation with 32 processors in parallel. our MD simulations alone cannot resolve this issue. However, in MC simulations of simple models, where the proper equilibrium behavior (of these simple models) can be accessed222In the case of MC simulations of the 2D Ising model, the equivalent of one MD run takes only a few seconds on a single processor., it is predicted that macroscopic domain formation is not merely delayed, but in fact eliminated (the theoretical motivation is the Imry-Ma argument (); citeulike:8864903 ()). In this section, we will show that these MC models make further predictions that can be tested against our MD results. Provided one uses the same positions for the obstacles, the resulting domain structure is remarkably insensitive to the model details. In fact, the minimal model requirement is that the unlike lipid species repel each other, which is already contained in the 2D Ising model. In what follows, we will show that the domain structures obtained in fig. 4 with MD, can indeed be reproduced in MC simulations of the 2D Ising model.

In the 2D Ising model, each lattice site  of a 2D square lattice is occupied by a spin variable . In the context of membranes, means that site  contains a DPPC (DLiPC) lipid, while cholesterol is ignored on this level Machta20111668 (). The energy of an Ising spin configuration is given by , with coupling constant (the sum runs over nearest neighbors). We use the Ising model to mimic the domain structures of fig. 4 by simulating a lattice. This approximates the aspect ratio of the MD simulation box, as well as the total number of phospholipids in a single leaflet. To represent the FPPC obstacles, the lattice is superimposed on the initial state of the MD simulations with obstacles. Those sites containing an FPPC phosphate atom subsequently have their spin value frozen to . The remaining sites are randomly initialized with values , under the constraint that the resulting (DPPC:DLiPC) composition matches that of the MD simulation as well as possible.

We then use MC simulations with Kawasaki dynamics newman.barkema:1999 () to determine equilibrium domain structures. That is, in each MC move, two non-FPPC lattice sites are chosen at random. Then, the corresponding spin values are exchanged with the Metropolis probability, , with the energy difference caused by the swap, the Boltzmann constant, and the temperature. In the absence of quenched disorder, the Ising model exhibits macroscopic phase separation when , as was shown by Onsager PhysRev.65.117 (). Note that is the only free parameter in the context of this work.

Figure 7: Typical domain structures of the “Ising” membrane MC simulation at . (a) Domain structure in the absence of quenched obstacles, and with DPPC (blue) and DLiPC lipids (orange). The stripe geometry is recovered, in agreement with the discussion of Section II.2. The snapshots (b,c,d) show typical structures in the presence of FPPC obstacles (red). In each of these snapshots, the FPPC configuration is the same, and chosen to reflect the MD configuration of fig. 4 (runs IIa and IIb). The differences in the snapshots (b,c,d) reflect the thermal fluctuations.

At , the parameter we have chosen to use for the comparison with the MD results, the system therefore assumes the stripe geometry [fig. 7(a)]. In the presence of quenched disorder, macroscopic phase separation is eliminated [fig. 7(b,c,d)]. Instead, finite domains that qualitatively resemble those of fig. 4 are seen. This similarity is quite remarkable, given that the simulation models used are entirely different (only the positions of the obstacles are the same). It shows that the generic features of the domain structure are essentially encoded in the obstacle configuration, and rather insensitive to the details of the interactions. To put this statement on a more solid basis, we now consider a quantitative measure of the correspondence between the domain structures seen in the two models.

iii.1 Correlation between the models

Figure 8: The ensemble averaged snapshot (i.e. the ’s) of the “Ising” membrane simulation corresponding to the obstacle configuration of fig. 7(b,c,d) (which, in turn, was derived from the MD obstacle configuration of fig. 4). Preferential locations for the appearance of domains are seen, which coincide well with the domain structures obtained in the MD simulations of fig. 4.

Because of the Ising model’s computational simplicity, many different equilibrium states can be generated quickly by means of MC simulation. In fact, for the lattice size used here, a representative set of the full equilibrium ensemble can be generated. This allows us to create an ensemble averaged snapshot by generating a lot of snapshots such as shown in fig. 7(b,c,d) and evaluating the average spin value at each site . As shown in fig. 8, this average snapshot reveals structure, which is due to the fact that the obstacles induce regions that locally prefer one of the lipid species citeulike:8864903 (); Machta20111668 (). With this average structure, we demonstrate the correspondence between the MD simulation results and those of the Ising model in two steps: First, we define a compatibility measure that quantifies how much a given snapshot of either model resembles the average spin values . Then, we evaluate the range of compatibilities within the Ising model itself, and compatibilities obtained from a structure of domains at random locations (i.e. we want to rule out that an observed correspondence is due to pure chance). It will be shown that the compatibilities of our MD results fall into the region of the former, and distinctively disagree with the predictions for the latter.

We define the compatibility of an Ising model snapshot with the as


where the summation runs over the lattice sites that do not contain an obstacle. This measure quantifies the degree of correlation between a single snapshot (characterized by spin values ) and the ensemble-averaged snapshot of fig. 8 (characterized by averages ). To calculate the compatibility of an MD snapshot with the averaged Ising structure of fig. 8, we formally use eq. (1), but with proper MD equivalents of the appearing terms. We consider each leaflet of the membrane separately, so the sum in eq. (1) runs over all DPPC and DLiPC molecules in the leaflet, where is their total number. If lipid  is a DPPC (DLiPC) lipid, then the equivalent spin value (). The corresponding value of is determined by projecting the leaflet onto the plane, and to then superimpose the averaged Ising structure of fig. 8. For each lipid , the corresponding is taken from the grid cell in the averaged structure that contains the lipid’s center of mass.

Figure 9: Compatibility of our MD simulations with the average domain structure of the Ising model in fig. 8. The compatibility is normalized by , which is the average compatibility of the Ising model with its own average structure. The upper shaded region corresponds to one standard deviation around this average. The inner and outer lower shaded areas represent one standard deviation regions that are expected for randomly located domains (corresponding to  nm and  nm, respectively; see details in text). The MD simulations with FPPC obstacles in both leaflets (runs IIa and IIb) equilibrate to states consistent with the Ising model (i.e. around the upper shaded region) and are clearly incompatible with random domains. In the case of FPPC in only one leaflet (run III), the -values are still distinctly different from random domains (i.e. outside the lower shaded region), but neither the FPPC leaflet (solid curve) nor the “free leaflet” (dashed curve) appear to be fully compatible with the Ising model.

Since the will in general have values between and , and since typical snapshots [cf. fig. 7(b,c,d)] do not exactly have the domain structure of the average snapshot [fig. 8], we do not expect the compatibilities of our MD simulation to reach . Instead, the proper values to compare against are those of the Ising model snapshots themselves, which are straightforward to obtain: We simply let the MC simulation that was used to produce fig. 8 run a second time and collect a histogram of values. The histogram is Gaussian, with mean and width ; the upper shaded region in fig. 9 marks the interval . In case the domains obtained in the MD simulations are compatible with the averaged Ising structure of fig. 8, we expect the corresponding compatibilities to be inside or at least in the proximity of the upper shaded area. Indeed, as fig. 9 shows, the MD compatibilities profoundly increase with time, indicating that the corresponding domain structures increasingly correlate with the equilibrium domains of the Ising simulation. As may be expected, the correlation is most pronounced in the case where the FPPC obstacles reside in both leaflets, where an excellent agreement with the Ising model is found. However, also when FPPC resides in only one leaflet, the correlation remains. In this case, the correlation is most pronounced in the leaflet containing the obstacles, as might be expected.

To demonstrate that the agreement between the domain structures of the Ising model and the MD simulations does not merely occur by chance, we also measure the typical compatibilities that one obtains in case the domains appear at random locations. To generate a random domain structure, we take the initial state of the MD simulation, but with all phospholipids (except FPPC) being DLiPC. Then, we pick a random location in the membrane and change the identity of all DLiPC lipids within a radius of this location to DPPC. This step is repeated until the (DPPC:DLiPC) ratio equals that of the simulation (note that the resulting domains are usually not circular, since the circular regions picked at each step can overlap). Since this process is fast, we can generate many such random domain structures and collect a histogram of compatibilities. This histogram is also Gaussian, but centered around a lower average value, and with a different width that depends on . The range of compatibilities expected for these domain structures is marked in fig. 9 by the lower shaded areas. As can be seen in the figure, they are clearly distinct from the MD curves (which holds for all values of ). We therefore conclude that the correspondence between domains obtained in our MD simulations and the Ising model is not coincidental. Despite the very different interaction potentials of the two models, both yield equivalent domain structures. Consequently, we expect that the domains seen in our MD simulations are indeed equilibrated states, as opposed to kinetically trapped states that would eventually macroscopically phase separate.

We emphasize that our findings are not unique to considered here. We have also performed the analysis for , and always found the same scenario: A reasonable agreement between the MD simulations with obstacles on both sides and the Ising model, while the “pure chance” regime is always excluded. The reason for this agreement is that over a wide range of temperature, the average domain structure [fig. 8] qualitatively remains the same, and merely becomes sharper as the temperature is decreased citeulike:8864903 (); Machta20111668 ().

Iv Summary

In this work, we have presented results from MD simulations of a DPPC/DLiPC/cholesterol membrane at near-atomic resolution in the absence and in the presence of quenched disorder. Our results show that static components in the membrane destroy the formation of large domains. Instead, a heterogeneous structure of small domains is seen [fig. 4]. For this to happen, it suffices that the obstacles are present in only one of the leaflets, since quenched disorder effects can be transmitted to the other leaflet via inter-leaflet coupling [fig. 5]. As the results have been obtained for a system with a strong demixing tendency, the elimination of macroscopic domain formation should occur for other lipid compositions also. A coupling of a cell to its surroundings or its cytoskeleton that causes an immobilization of membrane components may therefore explain the appearance of raft-like domains in thermal equilibrium.

By means of a suitably constructed compatibility measure we have also shown that when the obstacles reside in both leaflets, the domain structure can be predicted with very simple models, such as the Ising model (cf. fig. 9). We do not expect this to be a feature unique to the Ising model, but instead consider this result a signal that structure formation in the presence of quenched disorder is universal and does not strongly depend on the microscopic details of the system. Other similarly simple membrane models citeulike:9234255 (); citeulike:8864903 (); pinkmodel () are expected to be equally suited to predict domain structure. The exception is in cases where the obstacle configurations differ between the leaflets. In this situation, a simple monolayer model may be inadequate due to the neglect of inter-leaflet coupling. The construction and use of a related bilayer model is conceivable but has not been tested.

If one accepts the evidence that a MC simulation of a simple model can adequately describe the equilibrium domain structure of membranes in the presence of quenched disorder, their comparably low computational demands allow considering questions that are inaccessible in detailed MD simulations. The possibility to probe larger length scales permits the investigation of obstacle configurations with a structure on their own Ehrig201180 (); Machta20111668 (), and connect to the length scales of experimental visualization methods Hell:STED (); Eggeling:STED (), while the option to simulate many different setups allows to compare different positionings and different types of static obstaclesciteulike:8864903 (); citeulike:7115548 ().


This work was supported by the Deutsche Forschungsgemeinschaft via the Emmy Noether program (VI 483) and the SFB 803 (project B2).


  • (1) Singer SJ and Nicolson GL, The fluid mosaic model of the structure of cell membranes., Science (New York, N.Y.) 175, 720 (1972)
  • (2) Lenne PF and Nicolas A, Physics puzzles on membrane domains posed by cell biology, Soft Matter 5, 2841 (2009)
  • (3) Veatch S and Keller S, Seeing spots: Complex phase behavior in simple membranes, Biochim. Biophys. Acta 1746, 172 (2005)
  • (4) Pike LJ, The challenge of lipid rafts, Journal of Lipid Research 50, S323 (2009)
  • (5) Simons K and Ikonen E, Functional rafts in cell membranes, Nature 387, 569 (1997)
  • (6) Ramachandran S, Laradji M, and Kumar PBS, Lateral Organization of Lipids in Multi-component Liposomes, J. Phys. Soc. Jpn. 78, 041006 (2009)
  • (7) Saxton M, Chapter 8: Lateral Diffusion of Lipids and Proteins, vol. 48 of Current Topics in Membranes and Transport, chap. 8, 229–282 (Elsevier, 1999)
  • (8) Lenaz G, Lipid fluidity and membrane protein dynamics, Biosci. Rep. 7, 823 (1987)
  • (9) Castellana E and Cremer P, Solid supported lipid bilayers: From biophysical studies to sensor design, Surf. Sci. Rep. 61, 429 (2006)
  • (10) Yoon TY, Jeong C, Lee SW, Kim JH, Choi MC, Kim SJ, Kim MW, and Lee SD, Topographic control of lipid-raft reconstitution in model membranes, Nat. Mater. 5, 281 (2006)
  • (11) Veatch SL, Cicuta P, Sengupta P, Honerkamp-Smith A, Holowka D, and Baird B, Critical Fluctuations in Plasma Membrane Vesicles, ACS Chem. Biol. 3, 287 (2008)
  • (12) Veatch SL and Keller SL, Separation of Liquid Phases in Giant Vesicles of Ternary Mixtures of Phospholipids and Cholesterol, Biophys. J. 85, 3074 (2003)
  • (13) Honerkamp-Smith A, Veatch S, and Keller S, An introduction to critical points for biophysicists; observations of compositional heterogeneity in lipid membranes, Biochim. Biophys. Acta 1788, 53 (2009)
  • (14) Veatch SL, Soubias O, Keller SL, and Gawrisch K, Critical fluctuations in domain-forming lipid mixtures., PNAS 104, 17650 (2007)
  • (15) Veatch S, From small fluctuations to large-scale phase separation: lateral organization in model membranes containing cholesterol., Semin. Cell Dev. Biol. 18, 573 (2007)
  • (16) Semrau S, Idema T, Holtzer L, Schmidt T, and Storm C, Accurate Determination of Elastic Parameters for Multicomponent Membranes, Phys. Rev. Lett. 100, 088101 (2008)
  • (17) Schick M, Membrane heterogeneity: Manifestation of a curvature-induced microemulsion, Phys. Rev. E 85, 031902 (2012)
  • (18) Semrau S and Schmidt T, Membrane heterogeneity - from lipid domains to curvature effects, Soft Matter 5, 3174 (2009)
  • (19) Hu J, Weikl T, and Lipowsky R, Vesicles with multiple membrane domains, Soft Matter 7, 6092 (2011)
  • (20) Baumgart T, Hess ST, and Webb WW, Imaging coexisting fluid domains in biomembrane models coupling curvature and line tension, Nature 425, 821 (2003)
  • (21) Parthasarathy R and Groves JT, Curvature and spatial organization in biological membranes, Soft Matter 3, 24 (2007)
  • (22) Yamamoto T and Safran SA, Line tension between domains in multicomponent membranes is sensitive to degree of unsaturation of hybrid lipids, Soft Matter 7, 7021 (2011)
  • (23) Kusumi A, Shirai YM, Koyama-Honda I, Suzuki KGN, and Fujiwara TK, Hierarchical organization of the plasma membrane: Investigations by single-molecule tracking vs. fluorescence correlation spectroscopy, FEBS Lett. 584, 1814 (2010)
  • (24) Uline MJ, Schick M, and Szleifer I, Phase Behavior of Lipid Bilayers under Tension, Biophys J 102, 517 (2012)
  • (25) Tamm L and Groves JT, Supported Membranes, in: J. Struct. Biol., vol. 168, 1–222 (Elsevier, 2009)
  • (26) Mager MD, Almquist B, and Melosh NA, Formation and Characterization of Fluid Lipid Bilayers on Alumina, Langmuir 24, 12734 (2008)
  • (27) Goksu EI, Nellis BA, Lin WC, Satcher JH, Groves JT, Risbud SH, and Longo ML, Effect of Support Corrugation on Silica Xerogel-Supported Phase-Separated Lipid Bilayers, Langmuir 25, 3713 (2009)
  • (28) Goennenwein S, Tanaka M, Hu B, Moroder L, and Sackmann E, Functional Incorporation of Integrins into Solid Supported Membranes on Ultrathin Films of Cellulose: Impact on Adhesion, Biophys. J. 85, 646 (2003)
  • (29) Steinem C and Janshoff A, Modellmembranen auf Oberflächen. Verankert und doch mobil, Chem. unserer Zeit 42, 116 (2008)
  • (30) Nováková E, Mitrea G, Peth C, Thieme J, Mann K, and Salditt T, Solid supported multicomponent lipid membranes studied by x-ray spectromicroscopy, Biointerphases 3, FB44 (2008)
  • (31) Kusumi A, Nakada C, Ritchie K, Murase K, Suzuki K, Murakoshi H, Kasai RS, Kondo J, and Fujiwara T, PARADIGM SHIFT OF THE PLASMA MEMBRANE CONCEPT FROM THE TWO-DIMENSIONAL CONTINUUM FLUID TO THE PARTITIONED FLUID: High-Speed Single-Molecule Tracking of Membrane Molecules, Annu. Rev. Biophys. Biomol. Struct. 34, 351 (2005)
  • (32) Sheetz MP, Schindler M, and Koppel DE, Lateral mobility of integral membrane proteins is increased in spherocytic erythrocytes, Nature 285, 510 (1980)
  • (33) Kusumi A, Suzuki KGN, Kasai RS, Ritchie K, and Fujiwara TK, Hierarchical mesoscale domain organization of the plasma membrane, Trends Biochem Sci 36, 604 (2011)
  • (34) Sung BJ and Yethiraj A, Computer Simulations of Protein Diffusion in Compartmentalized Cell Membranes, Biophys. J. 97, 472 (2009)
  • (35) Auth T and Gov NS, Diffusion in a Fluid Membrane with a Flexible Cortical Cytoskeleton, Biophys. J. 96, 818 (2009)
  • (36) Saxton MJ, The spectrin network as a barrier to lateral diffusion in erythrocytes. A percolation analysis, Biophys. J. 55, 21 (1989)
  • (37) Saxton M, Anomalous diffusion due to obstacles: a Monte Carlo study, Biophys. J. 66, 394 (1994)
  • (38) Bouchaud JP and Georges A, Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications, Phys. Rep. 195, 127 (1990)
  • (39) Gómez J, Sagués F, and Reigada R, Effect of integral proteins in the phase stability of a lipid bilayer: Application to raft formation in cell membranes, J. Chem. Phys. 132, 135104 (2010)
  • (40) Yethiraj A and Weisshaar JC, Why Are Lipid Rafts Not Observed In Vivo?, Biophys. J. 93, 3113 (2007)
  • (41) Fischer T and Vink RLC, Domain formation in membranes with quenched protein obstacles: Lateral heterogeneity and the connection to universality classes, J. Chem. Phys. 134, 055106 (2011)
  • (42) Machta BB, Papanikolaou S, Sethna JP, and Veatch SL, Minimal Model of Plasma Membrane Heterogeneity Requires Coupling Cortical Actin to Criticality, Biophys. J. 100, 1668 (2011)
  • (43) Ehrig J, Petrov EP, and Schwille P, Near-Critical Fluctuations and Cytoskeleton-Assisted Phase Separation Lead to Subdiffusion in Cell Membranes, Biophys. J. 100, 80 (2011)
  • (44) Madden WG and Glandt ED, Distribution functions for fluids in random media, J. Stat. Phys. 51, 537 (1988)
  • (45) Gordon PA and Glandt ED, Liquid–liquid equilibrium for fluids confined within random porous materials, J. Chem. Phys. 105, 4257 (1996)
  • (46) Vink RLC, Fischer T, and Binder K, Finite-size scaling in Ising-like systems with quenched random fields: Evidence of hyperscaling violation, Phys. Rev. E 82, 051134 (2010)
  • (47) De Gennes PG, Liquid-liquid demixing inside a rigid network: qualitative features, J. Phys. Chem. 88, 6469 (1984)
  • (48) Imry Y and Ma SK, Random-Field Instability of the Ordered State of Continuous Symmetry, Phys. Rev. Lett. 35, 1399 (1975)
  • (49) Binder K, Imry Y, and Pytte E, Interface roughening and random-field instabilities in Ising systems in three or less dimensions, Phys. Rev. B 24, 6736 (1981)
  • (50) Risselada HJ and Marrink SJ, The molecular face of lipid rafts in model membranes, PNAS 105, 17367 (2008)
  • (51) Veatch S, Polozov I, Gawrisch K, and Keller S, Liquid Domains in Vesicles Investigated by NMR and Fluorescence Microscopy, Biophys. J. 86, 2910 (2004)
  • (52) Marsh D, Cholesterol-induced fluid membrane domains: A compendium of lipid-raft ternary phase diagrams, Biochim. Biophys. Acta 1788, 2114 (2009)
  • (53) Wolff J, Marques CM, and Thalmann F, Thermodynamic Approach to Phase Coexistence in Ternary Phospholipid-Cholesterol Mixtures, Phys. Rev. Lett. 106, 128104 (2011)
  • (54) Marrink SJ, Risselada HJ, Yefimov S, Tieleman DP, and de Vries AH, The MARTINI Force Field:  Coarse Grained Model for Biomolecular Simulations, J. Phys. Chem. B 111, 7812 (2007)
  • (55) Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, and Berendsen HJC, GROMACS: Fast, flexible, and free, Journal of Computational Chemistry 26, 1701 (2005)
  • (56) Hess B, Kutzner C, van der Spoel D, and Lindahl E, GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation, Journal of Chemical Theory and Computation 4, 435 (2008)
  • (57) Ehrig J, Petrov EP, and Schwille P, Phase separation and near-critical fluctuations in two-component lipid membranes: Monte Carlo simulations on experimentally relevant scales, New J. Phys. 13, 045019 (2011)
  • (58) Fischer T and Vink RLC, The Widom-Rowlinson mixture on a sphere: elimination of exponential slowing down at first-order phase transitions, J. Phys.: Condens. Matter 22, 104123 (2010)
  • (59) Grossmann B and Laursen ML, The confined-deconfined interface tension in quenched QCD using the histogram method, Nucl. Phys. B 408, 637 (1993)
  • (60) Bagatolli LA and Gratton E, Two Photon Fluorescence Microscopy of Coexisting Lipid Domains in Giant Unilamellar Vesicles of Binary Phospholipid Mixtures, Biophys. J. 78, 290 (2000)
  • (61) Kinnunen PKJ, On the principles of functional ordering in biological membranes, Chem. Phys. Lipids 57, 375 (1991)
  • (62) Speck T and Vink RLC, Random pinning limits the size of membrane adhesion domains (2012)
  • (63) Putzel GG, Uline MJ, Szleifer I, and Schick M, Interleaflet Coupling and Domain Registry in Phase-Separated Lipid Bilayers, Biophys. J. 100, 996 (2011)
  • (64) Newman MEJ and Barkema GT, Monte Carlo Methods in Statistical Physics (Clarendon Press, Oxford, 1999)
  • (65) Onsager L, Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition, Phys. Rev. 65, 117 (1944)
  • (66) Mouritsen OG, Boothroyd A, Harris R, Jan N, Lookman T, MacDonald L, Pink DA, and Zuckermann MJ, Computer simulation of the main gel-fluid phase transition of lipid bilayers, J. Chem. Phys. 79, 2027 (1983)
  • (67) Hell SW and Wichmann J, Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy, Optics Lett. 19, 780 (1994)
  • (68) Eggeling C, Ringemann C, Medda R, Schwarzmann G, Sandhoff K, Polyakova S, Belov VN, Hein B, von Middendorff C, Schönle A, and Hell SW, Direct observation of the nanoscale dynamics of membrane lipids in a living cell, Nature 457, 1159 (2009)
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description