A tutorial for mesoscale computer simulations of lipid membranes: tether pulling, tubulation and fluctuations

Abstract

Lipid membranes and membrane deformations are a long-standing area of research in soft matter and biophysics. Computer simulations have complemented analytical and experimental approaches as one of the pillars in the field. However, setting up and using membrane simulations can come with barriers due to the multidisciplinary effort involved and the vast choice of existing simulations models. In this review, we introduce the non-expert reader to coarse-grained membrane simulations at the mesoscale. Firstly, we give a concise overview of the modelling approaches to study fluid membranes, together with guidance to more specialized references. Secondly, we provide a conceptual guide on how to develop mesoscale membrane simulations. Lastly, we construct a hands-on tutorial on how to apply mesoscale membrane simulations, by providing a pedagogical examination of membrane tether pulling, shape and mechanics of membrane tubes, and membrane fluctuations with three different membrane models, and discussing them in terms of their scope and how resource-intensive they are. To ease the reader's venture into the field, we provide a repository with ready-to-run tutorials.

Fig. 1. Classification and applications of fluid membrane models. Continuum models have been used to study membrane fluctuation spectra, particle uptake, uptake dynamics, vesicle shapes, shape transitions, tether pulling, adhering vesicles, vesicle domains, and domain induced budding among others. Examples of available software packages to solve continuum models are Mathematica, MATLAB, Julia or FEniCS. Applications of mesh models are particle uptake, membrane–filament interactions, membranes and active particles and tether pulling. Available mesh model software are the Surface Evolver, Mem3DG, PyMembrane, TriMem or Flippy. Particle based models representing 1 lipid as a collection of 10–100 particles have been applied to study lipid phase behavior, studies on lipid types, membrane channels and ion transport and transmembrane proteins. These models can be simulated using software like CHARMM-GUI, GROMACS or NAMD. Simulations with lipids represented as 1–10 particles (which we refer to as several-beads-per-lipid models in Section 2.2.3), and focused on the 10 nm–10 μm scale have been applied to studying the self-assembly of lipids and other biomolecules into structured complexes, mechanical properties, membrane domain formation and membrane fusion and pore formation. Simulation models where a patch of membrane is represented as a single particle (one-bead-per-lipid-patch models in Section 2.2.4) have been exploited to probe phenomena at larger scales (20 nm–10 μm) such as particle wrapping and uptake, vesicle shapes, whole cell membranes and in-plane mechanics, among others. Available software to simulate these types of models are LAMMPS, HOOMD or ESPResSO. Images in the figure reproduced with permission from ref. 36, 39, 44, 45, 47, 49, 51, 53, 60, 61, 65, 67, 71, 72, 74, 75, 76, 80, 83 and 89.

Fig. 2. Tether pulling Monte Carlo simulation using a mesh model. (A) In a dynamically triangulated mesh, membrane fluidity is ensured by regularly swapping bonds in the mesh. To extrude a membrane tube, we set-up a flat membrane patch consisting of edge (fixed) and bulk vertices and tether the central membrane vertex to a bead (blue particle) via a harmonic spring. (B) To measure the force required to extrude a tether, we develop a two-stage protocol: first, we move the blue particle to generate simulation checkpoints; second, we reinitialize the system at each checkpoint to let the membrane relax for a fixed position of the blue bead. (C) Example of relaxation curves for the force pulling on the membrane (shades of black) and membrane elongation (shades of blue) obtained during the second protocol step for a membrane with κ = 20k_BT, γ = 1k_BT/σ² and average initial bond length 〈l〉/σ = 1.5. The blue bead is fixed at z ≈ 80σ from the initial membrane plane. As the simulation progresses, the pulling force on the membrane equilibrates to a non-zero value. (D) Force–elongation profile. The curve is obtained by averaging over the final equilibrated force in three replica simulations. The curve is normalised by the plateauing force f₀ measured in simulations. (E) Representative simulation snapshots for different membrane elongations shown as coloured bars in panel (D).

Fig. 3. Membrane tube equilibration using the YLZ potential. (A) The YLZ model models the membrane as a single layer of interacting beads. A tube of an arbitrary radius is constructed out of a regular trigonal lattice in a simulation box with periodic boundary conditions with ((↔)) showing the barostat-coupled z box dimension L_z, which can vary. (B) The size of the box along the cylinder axis is allowed to vary using a modified NPH Nose–Hoover barostat. Regardless of the initial state, under the action of the barostat, the membrane always equilibrates to the same radius if we set the tension and membrane parameters to be the same. (C) As predicted by the theory, the same membrane tube equilibrates to a different radius under different tensions with radius decreasing with increasing tension. (The shaded part shows a fixed-box equilibration.) (D) The parameter μ of the YLZ model changes the bending rigidity of the membrane, κ. The value of κ was estimated for three different values of μ using the tube radii at different tensions (red). This is compared to the estimates from the membrane fluctuation spectrum (blue). (E) The Helfrich theory predicts that the radius of a thermodynamically stable tube is proportional to . The simulations adhere to a power law nearly perfectly, with the exponent being close to that given by the theory. The equilibrium value of the radius for each set of parameters was obtained by averaging over time and three different random seeds. The three seeds were also used for an estimation of the error, which, however, is negligible (see the barely visible shaded area). All quantities in the figure are expressed in simulation units.

Fig. 4. Membrane fluctuation simulation using the Cooke and Deserno model. (A) (top) Diagram of Cooke bilayer lipid and interaction matrix showing pair potential sketches. (bottom) Simulation setup for a small membrane patch (small size chosen for clarity), with ((↔)) showing the barostat-coupled x, y box dimensions L_x = L_y = L, which can vary. (B) Analysis pipeline: skip over initial equilibration determined by the evolution of the box length L in time, then for the rest of the trajectory consider L to be 〈L〉; process equally spaced frames to obtain for each the instantaneous height fluctuation spectrum h_ij; finally for each mode of the spectrum compute its average power 〈h_q²〉 in Fourier space and corresponding uncertainty as the adjusted standard deviation . (C) Ensemble average of the height fluctuation power spectrum, in log–log space, with snapshot of simulation. The spectrum is linear in these coordinates (fit as dashed black line), up to the cut-off shown as a vertical line at the 12σ wavelength.

Fig. 5. Model comparison by computational time. (A) CPU time (hours) spent in simulating different membrane deformations using various membrane models. Mesh simulations: cargo budding simulations with a mesh model were performed with the parallelised software TriLMP, using two vesicles with N₁ = 2562 (with diameter σ₁/σ ≈ 30) and N₂ = 10 242 (σ₂/σ ≈ 60) vertices in the mesh; two cargo particles were tested with diameters σ_c,1/σ₁ = 6 and σ_c,2/σ₂ = 6, and the strength of the interaction was chosen to ensure sufficiently strong adhesion (ε/k_BT = 10). Dashed lines around symbols indicate that the model cannot accommodate the topological changes required to complete the deformation. The tube pulling results for the mesh MC model correspond to the time required to equilibrate the membrane and obtain the curves presented in Section 3.2. YLZ simulations: the CPU time for the tube equilibration was that for the simulations presented in Section 3.3. The time for the cell division was for a vesicle of 50 000 YLZ particles following a reliable division protocol requiring ∼800 000 time steps. Cooke simulations: the CPU time ranges were obtained from simulation datasets appropriate for each measurement of shape fluctuations, tube equilibration and cargo budding; for this latter case we included the CPU time required for the equilibration of the final state. For tube pulling and cell division, we used estimates based on scaling the budding simulation by a factor of respectively 4 and 140, the latter corresponding to a cell of area 1 μm². (B) Comparison of the performance of three openly available, parallelised membrane models. The drastic fall in efficiency for the mesh model TriLMP is due to the partially serial nature of the code: bond update moves, although computed in parallel, are currently serially implemented in TriLMP. Parallelisation in LAMMPS depends on how the simulation box is subdivided. The YLZ tube simulation was parallelised only by dividing the box along one axis, hence the monotonic curve. The Cooke fluctuation simulation was divided along two axes, which explains the dip at 8 CPUs.