Simulating realistic membrane shapes

Biological membranes exhibit diversity in their shapes and complexity in chemical compositions that are linked to many cellular functions. These two central features of biomembranes have been the subject of numerous simulation studies, using a diverse range of computational techniques. Currently, the field is able to capture this complexity at increasing levels of realism and connect the microscopic view on protein – lipid interactions to cellular morphologies at the level of entire organelles. Here we highlight recent advances in this topic, identify current bottlenecks, and sketch possible ways ahead.


Introduction
Biomembranes are essential cellular components. Together with membrane-adhered structures, such as the cytoskeleton, cell membranes constitute incredibly heterogeneous and crowded environments, containing hundreds of different lipid types and being densely packed with a large variety of membrane proteins. They provide identity not only to the cell as a whole, through the enveloping plasma membrane (PM), but also to many internal organelles. Membrane spatial organization ranges from a locally flat, lamellar geometry to highly curved ones as can be found, for instance, in the endoplasmic reticulum, Golgi apparatus, thylakoids, or mitochondria. Moreover, these structures are dynamic, involved in many remodeling processes often with highly curved intermediates [1].
Cellular and artificial membranes thus undergo constant shape transformations at a length scale much larger than their thickness. These large-scale membrane remodeling events are driven by three major phenomena: (i) alteration in macroscopic variables e.g., osmotic pressure or pH changes, (ii) external forces e.g., actin polymerization, or (iii) cooperative action of lipids and membrane proteins. Here, we focus on the ability of lipids and proteins to shape membranes. Understanding the coupling between membrane composition and curvature generation and sensing is of fundamental importance and remains an active field of research.
A qualitative understanding is obtained from theoretical concepts, such as the lipid shape model, in which the intrinsic curvature of a lipid is linked to the effective ratio between its head and tail volumes. Similarly, the propensity of proteins to induce membrane curvature can be explained by a number of basic mechanisms [2], including insertion of amphipathic or hydrophobic domains (wedging), entropic repulsion between soluble domains, asymmetric membrane protein shape, or scaffolding. Although powerful, these concepts treat membranes in a mean-field approximation and are limited to describe generic behavior. To obtain more quantitative insights in multicomponent and heterogeneous systems, computer simulations are an indispensable tool [3e5]. To unravel the complex interplay between proteins and lipids giving rise to membrane shape, a multitude of simulation techniques is required, ranging from detailed all-atom (aa) and coarse-grain (cg) molecular dynamics (MD) to mesoscopic simulation techniques.
Over recent years, important progress in this field has been made; in our view exciting new developments are the ability to capture the complexity of realistic membranes with detailed aaMD and cgMD and to simulate complex membrane shapes and shape transformations with mesoscale models. The field is currently heading toward connecting these scales, to eventually enable modeling of realistic membranes at the level of full cells and cell organelles.
The scope of this opiniated minireview is as follows. We first provide a brief overview of the state-of-the-art of the multiscale modeling methods used in this field, followed by selected recent examples of the major progress along the lines sketched above, culminating with current attempts to bridge to the cell (organelle) level. We finish with a short concluding section with our ideas on possible future roads ahead.

Multiscale modeling framework
To capture the broad range of time and length scales involved in membrane shape regulation, researchers in this field rely on a hierarchical modeling framework (Figure 1). At the high-resolution end, aaMD is a popular technique, with typical simulations containing a single protein solvated by a few thousand lipids, probing dynamics into the microsecond timescale. In parallel to technological advances in computational hard and software and development of better force fields, the potential of aaMD in the field of membrane simulations has increased rapidly, with simulations of entire small organelles already reported [6]. There is a growing range of complex membranes that have been simulated by aaMD, mimicking cell membranes of a variety of cell types and providing a better connection to experimental in vivo work [7,8]. An important development is also the ability of aaMD to quantify the effect of proteins on membrane curvature, with recent exciting examples revealing the curvature imprint of an Annexin trimer [9], membrane bending due to ESCRT proteins [10], and curvature generation due to protein crowding [11]. Such data are not only important to unravel the molecular driving forces of membrane reshaping by diverse classes of proteins, but are also essential to calibrate cg and mesoscale models in the hierarchical modeling framework. However, large-scale membrane shape changes are still beyond the aaMD capacities.
Coarse-grain models, in which groups of atoms are represented by effective interaction sites, can overcome this limitation. Owing to the inherent speedup, cgMD has become a powerful computational technique for biomembrane investigations; the feasible time and length scale allows for capturing lateral heterogeneities, as well as large-scale membrane shape deformations while keeping track of molecular details. Typical systems may contain up to hundreds of individual proteins and probe the dynamics in the microsecond to millisecond range. State-of-the-art examples are simulations of crowded PM models [12], tissue-specific complex Schematic view of a hierarchical modeling framework. At a length scale around 10 nm, the detailed membrane interaction of single proteins can be investigated using aaMD on the microsecond timescale. Protein-protein interactions can be probed using cgMD, at spatio-temporal scales of 100 nm and milliseconds. The results from aaMD and cgMD can be used to calibrate mesoscale models, which allow simulating the membrane remodeling induced by many proteins at micrometer length and macroscopic timescale (seconds). Backmapping to coarse-grain resolution can provide a near-atomistic view of these remodeled membrane shapes, with the possibility to further backmap local subregions to atomistic detail. The example shown illustrates the collective effect of Shiga toxin B subunit proteins causing the formation of a micrometer-long tubular membrane invagination, as detailed in Pezeshkian et al. [42]. membranes [13], and multicomponent bacterial membranes [14,15], to name but a few. The importance of these works is that they start mimicking the compositional complexity and nanoscale structuring observed in realistic cellular membranes. In addition, cgMD enables simulations of membrane remodeling by the collective action of multiple components. Exemplary recent work in this area includes the stabilization of mitochondrial folds by ATP synthase [16] Biomembranes show certain behaviors that take place at length and time scales that are beyond the capacity of both aaMD and cgMD. To investigate these systems, computational scientists have adopted even more aggressive coarse-graining approaches, in which proteins and lipids are represented by a few interaction sites only [24,25], or, field-based modeling techniques that ignore the molecular details altogether and treat the membrane, typically, as a continuous surface decorated by vector fields representing proteins or heterogeneity in the membrane [26e28]. The accuracy of these latter models is dependent on the validity of the macroscopic Hamiltonians that are used to describe the membrane shape. A popular and successful energy function describing large-scale membrane shape is the Helfrich Hamiltonian that has been adopted into many computational modeling schemes. Fiorin et al. [29], using sophisticated collective variables to access shape deformations in MD simulations, showed that the Helfrich Hamiltonian works nicely for length scales larger than 10 nm. The timescale of mesoscale methods, however, is generally less well defined but can be improved by coupling to continuum-based hydrodynamics [30]. Among mesoscopic models, dynamic triangulated surface (DTS) simulations have been very powerful tools to investigate generic physical principles in biological membranes [31,32]. Pezeshkian and Ipsen [32] recently provided an algorithm for DTS simulation of membrane patches at constant frame tension allowing themto study the sorting affinity of Annexin proteins into nanotube membrane tethers [9]. A grand canonical ensemble of DTS simulation also has been developed to obtain membranes with constant tension [33].
In recent years, attempts have been made to exploit mesoscale models for specific biological processes such as filament induced membrane tubulation [ Overall, mesoscale simulation techniques show a great power at capturing realistic membrane shapes and shape remodeling. A limitation is that, unlike aaMD and cgMD, there is no general-purpose open-source software and force field (generic mapping procedure) available to make it accessible to the larger computational community. Another challenge for mesoscale models is to realistically include dynamic aspects, e.g., how to deal with inhomogeneous friction due to curvature gradients. The dynamics of biomolecules embedded in a membrane are strongly affected by each monolayer viscosity and intermonolayer friction between the two leaflets, both of which can be curvature dependent as lipids residing in a curved membrane pack differently. As a result, diffusion of lipids and proteins will be affected by membrane shape [40].
An important, ongoing, effort is to systematically bridge the various levels of resolution into hybrid or multiscale setups. An ideal way of combining scales would be in a concurrent scheme, in which areas of high resolution are interfacing with a lower resolution surrounding, either in a static or adaptive approach. Although there have been several attempts to create concurrent multiscale schemes for membrane systems, as in the recent work of Yiu et al. [41], the application of such methods appears still limited by artefacts arising at the resolution interface or lack of computational efficiency. Arguably a more fruitful approach is that of consecutive multiscaling. In consecutive schemes, typically, modeling starts with aaMD to provide mapping constraints for cgMD. In turn, aaMD and cgMD trajectories are used to calibrate mesoscale models, which then can be used to obtain large complex equilibrated membrane structures. These structures can be converted back to the cg and atomistic levels in a process denoted backmapping, providing a connection between all the scales ( Still, there is room for better integration of aa/cgMD and mesoscale models. The major limitation of mesoscale approaches is that their resolution is above the membrane thickness, whereas the size of many interacting proteins is within this scale. Although the global effect of membraneeprotein interactions on the large-scale membrane shape can be captured by adding a few additional terms into the Helfrich Hamiltonian, the local effects of these proteins, through lipid-mediated and direct proteineprotein interactions, are missing. Complicating factors are that such interactions are not only strongly dependent on the details of the environment but also time-dependent. A promising solution to this limitation is a dynamic coupling between the levels of resolution, i.e. performing aaMD or cgMD in parallel to mesoscale simulations during which the proteine protein interactions handled by the higher resolution methods are feeding back into the mesoscale simulations [44,45].

Lipid sorting in curvature gradients
The ability of aaMD and cgMD to model multicomponent membranes, as well as large enough patches to consider different morphologies has opened the way to quantify curvature-induced sorting of lipids in complex environments. This is important to aid in the interpretation of experiments, which are becoming more and more sophisticated in probing the local membrane organization. Furthermore, quantification of lipid (and protein) sorting preferences can also accelerate the simulation of highly curved and large membrane patches, by providing close to equilibrium lateral distributions of the constituent biomolecules as starting configurations.
Simulating large PM patches composed of seven different lipids types, Koldsø et al. [46] revealed that shape fluctuations lead to curvature induced sorting of lipids. They observed strongest sorting of phosphatidylethanolamine (PE) and cholesterol, preferring regions of negative curvature as expected from the shape concept, but also of phosphatidylinositol 4,5bisphosphate lipids and gangliosides. Similar sorting of gangliosides and other PM lipids was observed by Jefferies and Khalid [47] and by Sharma and Lindau [48], but work from Lipowsky et al. [22] indicates that gangliosides in fact cause positive curvature, seemingly at odds with the other results mentioned and pointing to an important role for the size of the ganglioside headgroup which differed between these studies, or to concentration effects. Simulations of tethers pulled from asymmetric PM membranes containing >60 lipid types show a very rich and complex sorting behavior, with often subtle enrichments/depletions of lipid types between the tether and basal membrane [49]. Notably, the sorting was found to strongly depend on the pulling direction, and to be affected more by the lipid tails than the headgroup type. An interesting role of lipids with polyunsaturated tails was found, being enriched in the strongly curved tether in both leaflets and for both pulling directions (Figure 2-ii). This promiscuous character results from the highly flexible polyunsaturated tails which can sustain both positive and negative curvature gradients, as is also apparent from the work of Gautier and Antonny [50]. Nonintuitive sorting can also arise from competitive effects, as is observed for instance when simulating mixtures of PE and cardiolipin (CL), which both prefer negative curvature. At high enough concentrations, CL outcompetes PE and forces the latter to sort to regions with zero or positive curvature (Figure 2-i). Interesting sorting behavior can also arise in lipid mixtures which phase separate into liquidordered (Lo) and liquid-disordered domains, with the stiffer Lo domains even able to occupy highly curved regions when put under tension (Figure 2-iii) [51e53].
The emerging picture from these, and related, simulations is that, even though driving forces are small, these can lead to significant compositional gradients between differently curved membrane regions. Lipid sorting can be rationalized to a large extent by the shape concept, even in complex environments. However, exceptions to this simplistic view arise due to d at least d three distinct mechanisms ( Figure 2): (i) nonadditive effects, such as competitive sorting where lipids with the strongest curvature preference may push other lipids to less favorable regions, (ii) the presence of promiscuous lipids, such as polyunsaturated lipids which can adapt to different curvature gradients, and (iii) tensionecurvature interplay, which can cause a change in membrane properties such as softening of Lo domains.
To further unravel these nonadditive effects on curvature-induced lipid sorting and lipid-induced membrane curvature generation (two sides of the same coin), more systematic studies are required as recently undertaken by several groups [54e56]. Efficient approaches to create and analyze arbitrary membrane morphologies, such as TS2CG [43], BUMPy [57], and MemSurfer [58], as well as methods to stabilize these, e.g., using layers of virtual particles [59] or collective variables [60,61], will be very helpful in this respect. These tools would also allow simulations of curvature-induced protein sorting at (near) atomic detail [62], which is so far mostly in the realm of mesoscale models.
An emerging problem when simulating complex multicomponent membranes, either with aaMD or cgMD, is dealing with trans-bilayer asymmetry. Simply matching the area of each monolayer does not guarantee a tensionless membrane [63e65]. A surplus of lipids on one side by itself can generate curvature stress. It has been found that membranes containing fast flipflopping molecules, such as cholesterol, generates a tensionless membrane while retaining trans-bilayer asymmetry, at least when the chemical potential of cholesterol is equal in both leaflets [64]. However, rapid flip-flop of small molecules can only be achieved in cg models. The problem of proper choice of leaflet composition is compounded when considering nonlamellar membrane shapes. This is still largely unexplored and requires careful consideration both from the theoretical and simulation side.

Toward full cell membrane shapes
Growing experimental evidence from 3-dimensional reconstructions of biomolecular structures, e.g., electron tomography and confocal microscopy, indicates that biological membranes take up a diverse range of shape conformations. For instance, although the PM typically is planar on the nanometer scale, at the greater length scale and during many cellular processes, large morphological membrane changes occur. In addition, eukaryotic cells contain various membrane-bound organelles, e.g., mitochondria and the endoplasmic reticulum, that show highly complex shapes that are dynamic in both morphology and surface topology. For instance, the inner membrane of a small mitochondrion exhibits a strongly curved surface morphology with a high topological genus of 46 (the number of handles on the surface, for example, for a torus the genus is one and for a sphere it is zero) [43,66].
Although aaMD and cgMD simulations of viruses and small organelles are now achievable [43, 67,68], the required simulation time for membrane remodeling processes at realistic organelle's scale is still beyond the capacity of these techniques. On the other hand, mesoscopic models have proven to be powerful for capturing large-scale membrane shape and lateral organization remodeling. Biologically relevant processes at length scales observable by optical microscopy, e.g., the formation of membrane tubular protrusions by polymerizing filaments [34,33], membrane constriction by chiral assembly of proteins [69] and tubular membrane invaginations by membrane proteins [28] are now readily achievable (Figure 3a-i,ii). These methods are close to capturing the level of entire organelles. Tachikawa et al., for instance, developed a DTS scheme to create a minimal model of the self-organization of the Golgi apparatus through assembly of small vesicles (Figure 3aeiii) [36].
Nevertheless, mesoscopic models can provide very little information about the molecular organizations as the size of many biomolecules is below the resolution of these models. To provide a realistic (i.e. realistic in both size and composition) model of organelle membranes, it is clear that an integrative multiscale simulation framework is required. A promising protocol toward this aim is as follows: From experimental data e.g., cryoelectron tomography (cryoET), a 3-dimensional structure of a realistic membrane shape is obtained. The experimental data are then converted to mashed surfaces [70] which are typically the input for mesoscopic simulation techniques that can be exploited to obtain local organization of biomolecules (proteins and lipid density) within their resolutions. The output of the mesoscale model will then be backmapped, e.g., using TS2CG, to a full organelle structure at the cg level. For such systems, although very large, microsecond cgMD simulations are still achievable allowing equilibration of the local organization of the constituents. Any part of the full system that requires more detailed information, can be extracted and backmapped to atomistic resolution. As proof of concept, we recently [43] used this protocol to simulate the membranes of an entire mitochondrion with realistic size and lipid compositions (Figure 3b).
To further progress on modeling realistic membrane morphologies at the cellular scale, an important driver will be the availability of detailed experimental data. In particular, data on local structure and dynamics, as most published measurements (e.g. proteomics and lipidomics analysis) are static ensemble averages representing a population mean. Obtaining this information at high spatial and dynamic resolutions is extremely difficult, particularly in living cells. Fortunately, exciting progress in experimental assays is being made, which implies an immense step forward in elucidating the architecture and stoichiometry of cellular components at unprecedented resolutions, setting the stage for spatially detailed simulations of whole cells [71,72].

Conclusion
During the past years, significant progress has been made in our ability to simulate the complexity of realistic membrane shapes, based on a hierarchical modeling framework. On the higher resolution side, we are now able to study curvature-induced lipid sorting mechanisms in multicomponent membrane systems with aaMD and cgMD, and soon will be able to extend this systematically toward protein sorting. At the lower resolution end, we are now capable of modeling realistic morphologies at the organelle scale with a variety of mesoscale modeling techniques. Exciting current efforts are directed to combine these models of different resolution, to fully capture the complexity and plasticity of cell membranes.
To add to this challenge, cell membranes are not isolated entities: they are in continuous interaction with their surroundings. Signaling pathways connect the inside of cells or cell organelles to the outside and may trigger local structural changes, such as the formation of raft platforms, all the way to large-scale membrane remodeling events during fusion or fission. In turn, dynamic processes happening at the membrane trigger numerous downstream processes in other parts of the cell. The increasing recognition that the cytoplasm also constitutes a very crowded and heterogeneous environment makes it clear that we are only at the beginning of our quest for modeling membranes in a native context. In the end, whole-cell models at molecular resolution are key.
Moreover, an essential characteristic of living cell membranes is that they are nonequilibrium structures, subject to and driven by active processes, specifically constant exchange of membrane-bound materials and the action of active proteins that shuttle between different states, driving changes in the local structure of the membrane upon energy consumption. In addition, they are in contact with the intracellular milieu that is the source of many out of equilibrium processes. Therefore, membrane shape and lateral organizations are affected by far from equilibrium processes. This feature of biomembranes has rarely been considered in simulations to date.
New methods are therefore needed to unravel the molecular driving forces that underlie membrane-mediated cellular processes. Major remaining challenges, at all levels of resolution, are as follows: (i) to further increase complexity in terms of membrane composition, making use of advanced experimental data d in particular with respect to local compositions in living cells; (ii) to add the cellular context, for instance, the cytoskeleton framework [73] and realistic models for the cytoplasm [11]; (iii) and to include out of equilibrium aspects, for instance, using open boundary methods [74] and incorporating reactive schemes [75] into classical and mesoscale methods. These challenges are compounded by the need to capture these processes over a multitude of spatio-temporal scales, requiring further progress in bridging the various computational methods into a consistent modeling framework. Advances in incorporating machine learning and artificial intelligence may prove to be a key asset in this respect [45].

Conflict of interest statement
Nothing declared.