Efficient tunable generic model for fluid bilayer membranes

We present a model for the efficient simulation of generic bilayer membranes. Individual lipids are represented by one head- and two tail-beads. By means of simple pair potentials these robustly self-assemble to a fluid bilayer state over a wide range of parameters, without the need for an explicit solvent. The model shows the expected elastic behavior on large length scales, and its physical properties (eg fluidity or bending stiffness) can be widely tuned via a single parameter. In particular, bending rigidities in the experimentally relevant range are obtained, at least within $3-30 k_{\text{B}}T$. The model is naturally suited to study many physical topics, including self-assembly, fusion, bilayer melting, lipid mixtures, rafts, and protein-bilayer interactions.

Lipid molecules in aqueous solution spontaneously assemble into bilayer membranes. In biological systems, such membranes are involved in tasks over an extraordinary range of length scales, from transport of water and ions at the scale of nm, up to phagocytosis, amoebal motion and cell budding at the scale of µm [1]. Computer simulations designed to understand some aspects of this structural and functional range must accordingly be tailored to the specific length and timescales involved. Techniques that probe both the smallest [2] and largest [3] of these length and time scales are now comparatively well established; however, accessing intermediate regimes has proven far more difficult, and it is only recently that significant progress has been made in this regard. The need for a comprehensive suite of techniques to study lipid bilayers at the mesoscale is highlighted by the sheer number of relevant problems in this regime, which include viral budding, raft formation, fusion, phase separation of multicomponent systems and protein sorting during vesiculation.
Most existing approaches to mesoscale simulation employ coarse grained lipids and require explicit solvent particles to stabilize the bilayer. This strategy is convenient and natural, yet it comes at a high price: Already for small flat systems the solvent accounts for most of the computational effort, but the problem gets significantly worse when three dimensional objects such as vesicles are to be simulated. Membrane and solvent play the role of surface and bulk, respectively, hence the solvent degrees of freedom vastly outnumber the lipids even for rather modest sized vesicles. An obvious solution to this problem is to replace the solvent by effective lipid interactions. Given the great success of this approach in polymer physics it is perhaps surprising that solvent-free bilayer simulations have so far failed to find widespread acceptance. The problem appears to be that naive choices for interparticle potentials (e. g. Lennard-Jones) do not lead to a fluid bilayer phase but only to ordered "solid" bilayers at low temperature and low density phases at high temperature. Many attempts to obtain a broadly stable fluid phase have been made with varying degrees of success. So far, however, none of these has resulted in a model that is sufficiently robust, simple, or versatile for general use. For example, the original solvent-free model of Drouffe et al. [4] and later modifications by Noguchi [5] rely on density dependent interaction potentials to stabilize the fluid phase. But the multibody nature of interactions poses serious problems for interpretation and measurement of thermodynamic quantities. Other models do not exhibit the crucial property of unassisted self-assembly [6,7] and require the use of angular dependent potentials [7] or a large set of highly tuned interaction parameters [6]. Thus, there remains a clear need for an efficient, robust and tuneable solvent-free bilayer model.
In this letter we present a very simple method for simulating solvent-free fluid bilayer membranes that is based on pair potentials, is highly robust and tuneable, and reproduces macroscopic properties of real bilayers. The key ingredient is an attractive potential between lipid tails for which the range of attraction, w c , can be varied. First we map out the properties of a lipid system as a function of w c and temperature and find that a sufficiently large w c is the key to obtaining a fluid phase. Moreover, w c can be used to tune bilayer properties and to construct multi-component systems with interfacial tension. This is demonstrated in the final part of our Letter where we examine the kinetics of domain formation in a two-component vesicle.
Let us now describe a model that demonstrates the principle mentioned above. Each lipid molecule is represented by one "head" bead followed by two "tail" beads. Their size is fixed via a Weeks-Chandler-Andersen potential with r c = 2 1/6 b. We use ǫ as our unit of energy. In order to ensure an effective cylindrical lipid shape we choose b head,head = b head,tail = 0.95 σ and b tail,tail = σ, where σ is the unit of length. The three beads are linked by two FENE bonds with stiffness k bond = 30 ǫ/σ 2 and divergence length r ∞ = 1.5 σ. Lipids are straightened by a harmonic spring with rest length 4σ between head-bead and second tail-bead which corresponds in lowest order to a harmonic bending potential 1 2 k bend σ 2 ϑ 2 for the angle π − ϑ between the three beads. We fixed the bending stiffness at k bend σ 2 = 10 ǫ. Finally, all tail beads attract each other according to This describes an attractive potential with a depth of ǫ which for r > r c smoothly tapers to zero. Its decay range w c is the key tuning parameter in our model. We performed Molecular Dynamics (MD) simulations using the ESPResSo package [8]. A fixed number of lipids was simulated in a cubic box of side length L subject to periodic boundary conditions. The canonical state was reached by means of a Langevin thermostat [9] (with a time step δt = 0.01 τ and a friction constant Γ = τ −1 in Lennard-Jones units). If needed, constant tension conditions were also implemented via a modified Andersen barostat [10] (with a box friction Γ box = 2 × 10 −4 τ −1 and box mass within the range Q = 10 −5 . . . × 10 −4 ).
One of our principal objectives is to map out the conditions under which the fluid bilayer is stable. We identified the fluid phase in two different ways. First, a box-spanning bilayer was pre-assembled from 1000 lipids, and its equilibration under zero lateral tension was attempted (requiring -if successful -box lengths of L ≈ 25 σ). Three qualitatively different outcomes were observed: (i) At sufficiently low temperature the bilayer adopted a "gel" phase; (ii) within a more elevated temperature range a fluid phase can be reached, provided w c 0.8 σ; (iii) at sufficiently high temperature a bilayer under zero tension always fell apart. Fluid and gel phases are already clearly distinct under visual inspection (long range order in the gel phase and none in the fluid). Across the gel-fluid boundary [11] we observed a sharp increase in in-plane diffusion constant D, an abrupt decrease in orientational order S = 1 2 3(a i ·n) 2 −1 i (where a i is the vector along the axis of the i th lipid and n is the average bilayer normal) and the emergence of a nonzero flip-flop rate r f (the probability per unit time that a single lipid changes its monolayer). Typical values for these parameters along the k B T = 1.1ǫ isotherm are D = 0.06 − 0.03σ 2 /τ , S = 0.5 − 0.8 and r f = 2 − 90 × 10 −5 τ −1 . For w c = 1.6σ and k B T = 1.1ǫ we also measured the bilayer rupture tension Σ r ≈ 4mN/m and the compressibility modulus at zero tension, K ≈ 50mN/m (mapping to real length scales by assuming a bilayer thickness of 5nm). These quantities are within their respective measured ranges for synthetic and biological membranes, Σ = 0.1 − 12mN/m [12] and K = 50 − 1700mN/m (depending on cholesterol content) [13]. Second, lipids were simulated at constant volume, but starting from a random "gas" configuration. Under all conditions which previously gave stable tensionless membranes, a bilayer patch quickly self-assembled which at the right L could zip up with its open ends to span the box. If the box was too big, the patch either remained free, or (sometimes) closed upon itself to form a vesicle. Box spanning bilayers could also occur somewhat above the evaporation boundary of Fig. 1 [14], indicating that this line should be viewed as the location where the rupture tension for a bilayer approaches zero. We finally remark that upon approaching this line from below, the bilayer becomes increasingly disordered, and flipflop rates and diffusion constants increase strongly.
Looking at figure 1 we see that the temperature range over which a fluid bilayer is stable increases as w c is increased and that for relatively narrow potentials w c < 0.7σ the fluid region disappears completely [15]. It is noteworthy that these general features are not restricted to the present functional form of Eqn. (4). Indeed, we have also obtained a qualitatively similar phase diagram for a Lennard-Jones like potential with variable width. It should therefore be emphasized that it is not the precise functional form of our tail attractions that is important for the stabilization of fluid bilayers but rather the length scale over which these attractions are effective.
Next we illustrate that the fluid bilayers approach the correct elastic continuum limit. If one expands the bilayer in modes h(r) = q h q e iq·r with q = 2π L (n x , n y ), (linearized) Helfrich theory predicts the mode spectrum [16] where κ is the bending modulus and σ the lateral tension. Below the crossover wave vector q c = σ/κ one has |h 2 q | ∼ q −2 (tension regime), while |h 2 q | ∼ q −4 holds above (bending regime). Once 1/q approaches length scales comparable to the bilayer thickness, continuum theory breaks down and further effects (e. g. protrusion modes [17]) set in. In order to see the characteristic q −4 scaling of the bending regime over a sufficiently wide range requires q c to be as small as possible, hence we simulated at zero tension. Furthermore, to get away from microscopic lengths (such as the bilayer thickness) we took systems four times as big as the ones we used for mapping the phase diagram (4000 lipids, L ≃ 50 σ). Note that reaching the continuum limit in MD simulations is not trivial, since the relaxation time of bending modes scales as q −4 .
After setting up the bilayer, we first waited until tension, box length, and energy had equilibrated (which took typically 10 5 τ for fluid systems). Then on the order of 100 configurations separated by 2000τ were used to measure the mode spectrum. The bilayer mid-plane was identified by tracking the tail-beads and interpolating their vertical position onto a 16 × 16 grid. Possible stray lipids had to be excluded from this procedure. An FFT then yields the power spectrum |h 2 q | . Fig. 2 illustrates (for the system with w c = 1.6σ and k B T = 1.1ǫ) the q −4 scaling. Notice that length scales exceeding L ≈ 20 σ (i. e., about four times the bilayer thickness) are required to reach the asymptotic regime, outside of which a fit to Eqn. (5) and the subsequent extraction of a bending constant is meaningless. The inset shows the corresponding values of κ for the six simulated stable bilayer systems at this particular value of the temperature (see Fig. 1). We emphasize that its value is easily tuneable over the experimentally interesting range, at least from κ ≈ 3 k B T up to κ ≈ 30 k B T .
So far we have used the parameter w c only to tune the bilayer stiffness. Let us now illustrate another application, namely, the possibility of creating mixed lipid systems. We consider a simple example in which two lipid types A and B are present and where w AA c = w BB c . A line tension can now be generated by choosing the cross interaction w AB c smaller than the homogeneous ones. Provided a sufficiently small cross term has been chosen, A and B domains will form. Figure  3 illustrates qualitative features of this process for a critical quench (A:B = 1:1). First we note that matching domains always appear on inner and outer leaflets. As is typical for a critical system, domains are initially elongated in shape and then coarsen to form circular patches or stripes. In this case the domain size and A-B line tension are sufficient to induce budding [18]. This process has recently been studied in detail via DPD simulations [19]. For matching lipid proportions in outer and inner leaflets, budding was induced by cleavage of the domain boundary. We also observed this "flaking mechanism" for very high line tensions (w AB c ≪ w AA,BB c ), while for less pronounced line tensions a more continuous type budding occured (see Fig. 3), provided the vesicle was big enough.
Next we investigate the kinetics of domain formation on a larger vesicle in the off critical regime. Although the kinetics of domain formation has been extensively studied for supported membranes, only a small number of recent simulations [19,20,21] and experiments [22] have tackled the problem for free vesicles. The importance of these studies is underscored by the fact that many features of a vesicular system are absent in a flat supported membrane. These include bending fluctuations, volume constraints, and kinetic effects of curvature-composition coupling [23].
To study the kinetics of phase separation we performed simulations using pre-equilibrated vesicles with an A:B ratio of 3:7 and measured the number n(t) of clusters of A lipids as a function of time. This quantity displays two distinct kinetic regimes (see Fig. 4). In the first regime n(t) decays exponentially, corresponding to a conversion of an initially exponential cluster size distribution to an arrangement where meso-sized clusters begin to dominate. At later times such clusters display an n(t) ∼ t θ power law decay with θ ≈ −0.4, which agrees with the value expected for coarsening via patch-coalescence in the non-hydrodynamical regime [24]. This scenario is confirmed visually, as well as by characteristic jumps in the (size weighted) average cluster size. Interestingly, we do not observe evidence for a Lifshitz-Slyozov type ripening, for which one would have θ = −2/3 [25].
Having demonstrated the key features of the present model, including its tunability and its application to multi-component systems, we now turn to performance aspects. Although one of the clear advantages of our approach is speed, it is difficult to make meaningful comparisons across different computer architectures and implementations. One basic quantity that provides a rough implementation-independent comparison is the particle number. Simulating a vesicle with 16000 coarse grained lipids, as in our example, would require roughly 25 times as many solvent particles [19,20]. (Actually this factor scales with the vesicle radius). Obviously, one must also account for the time-step which is often chosen somewhat larger in DPD simulations. Still, this leads us to an at least 5-fold speed-up for the present method over DPD, and a crude comparison based on direct CPU time usage agrees with this estimate [26].
In summary, we have presented a method that is fundamentally different from existing coarse grained lipid membrane simulations that use explicit solvent. Our method should be seen as complementary to these techniques since it presents a compelling speed advantage, especially for systems in three dimensions (vesicles, bicontinuous phases). It does not naturally include volume constraints or hydrodynamics. Indeed, the model is inspired by the great success of simple "beadspring" models used to study polymer systems in the Rouse regime. We have presented Langevin Dynamics simulations here, but the simplicity of the model and concept easily permit implementation of other integrators or even Monte Carlo. Volume constraints should be relatively simple and efficient to include via the concept of a "phantom solvent" [27]. The particle based nature of the model and its tunable interactions allow for great flexibility such that components with different pre-determined stiffness and lipid shape present no difficulty, and topological changes occur naturally. Despite its astounding simplicity our approach allows one to capture all these essential aspects of lipid bilayer physics within a single unified framework.