Mesoscale simulation of biomembranes with FreeDTS

We present FreeDTS software for performing computational research on biomembranes at the mesoscale. In this software, a membrane is represented by a dynamically triangulated surface equipped with vertex-based inclusions to integrate the effects of integral and peripheral membrane proteins. Several algorithms are included in the software to simulate complex membranes at different conditions such as framed membranes with constant tension, vesicles and high-genus membranes with various fixed volumes or constant pressure differences and applying external forces to membrane regions. Furthermore, the software allows the user to turn off the shape evolution of the membrane and focus solely on the organization of proteins. As a result, we can take realistic membrane shapes obtained from, for example, cryo-electron tomography and backmap them into a finer simulation model. In addition to many biomembrane applications, this software brings us a step closer to simulating realistic biomembranes with molecular resolution. Here we provide several interesting showcases of the power of the software but leave a wide range of potential applications for interested users.


Introduction
Fluid artificial and biological membranes can adapt to a diverse range of morphologies from multispherical structures to the astonishing forms seen in the subcellular organelles such as endoplasmic reticulum, Golgi apparatus, and mitochondria [1,2].These shapes are often dynamic and constantly undergo significant changes that are crucial for cell function e.g., such as endocytosis, cellular respiration.Also, biological membrane shapes contain information about overall health of an organism, physiological state of the cell and abnormal membrane architectures are implicated in many diseases such as Parkinson's Disease [3][4][5][6][7].The shape classes of simple lipid membranes are well understood and they can be described by a few macroscopic parameters such as spontaneous membrane curvature and pressure difference [1].The membranes of cells, however, are much more complex and distinct mechanisms govern their overall arrangement at various scales, ranging from the molecular to the macroscale [8][9][10].Thus, the study of membrane organization continues to be an active and essential field of science.Computer simulations are an effective tool for studying biomembrane architecture and the mechanisms involved in its remodeling [11].For analyzing membrane shape at length scales of up to 100nm, molecular simulations, such as molecular dynamics and dissipative particle dynamics simulations, have been very effective.Examples of exciting developments include the local curvature of membranes caused by a single protein or by the assembly of proteins [12][13][14], membrane shape induced lipid sorting [15], wetting induced membrane deformation [16], and even curvature sensing of proteins [17].Due to limitations in accessibility of large time and length scales, molecular-based simulations alone cannot provide a comprehensive picture of membrane remodeling.On the other end, macroscopic modeling that incorporates protein effects as a mean field density function enables the description of large-scale and generic membrane remodeling behaviors.However, it overlooks many phenomena which are the result of membrane fluctuations and rotational and translational entropy of single and (few) proteins [8,18,19].In between these two scales, mesoscopic modes in which large biomolecules such as proteins are explicitly considered, while lipids are modeled in a mean field manner (but allows for undulations and shape fluctuations) are needed to fully map out the organization of complex membranes [20,21].
Several mesoscopic models have been used to explore diverse range of membrane associated processes such as protein clustering via membrane-mediated interaction [22], membrane shape remolding by crowding of intrinsically disordered proteins [23], membrane neck constriction by assembly of proteins [24] and even activity-driven membrane remodeling [25] (for more see [20] and the reference within).In spite of this, these studies are often conducted with in-house software, or/and the software is limited to those specific applications that are difficult to apply to new research questions, which has hampered progress in mesoscopic membrane modeling.Recently, a software package called TriMem has been released that facilitates the simulation of triangulated surfaces with a strong focus on performance, in order to integrate the evolution of a system efficiently [26].However, this package currently lacks the protein model and only applicable to simple lipid bilayers.
Here, we present FreeDTS, a software package for the computational investigation of biomembranes, at the mesoscopic length scale which can also be used for macroscale modeling of membranes.The software is designed to cover a diverse range of biological processes and is easy to expand to cover more.In FreeDTS a membrane is represented by a dynamically triangulated surface equipped with vertex-based inclusions to integrate the effects of integral and peripheral membrane proteins.Several algorithms are implemented in FreeDTS to perform complex membrane simulations in different conditions, e.g., framed membranes with constant tension.Membranes also can be confined into a fixed region of the space to explore the effect of the environment on the membrane shape and fluctuations.In addition, FreeDTS allows one to turn off the shape evolution of the membrane and only explore organization of membrane proteins.This allows us to take realistic membrane shapes obtained, for instance, from cryo-electron tomography and obtain heterogeneous organization of biomolecules which can be backmapped to finer simulations models.This feature with helps from backmapping schemes e.g., TS2CG, brings us a step closer to simulating realistic biomembranes with molecular resolution [27,28].In the following, we report several interesting examples to show the power of the software and the detail information of how to use the software is included in the manual.

The model
Membranes: FreeDTS represents a membrane as a triangulated surface containing   vertices,   edges and   triangles (Figure SI-1).During the system evaluation, the vertices position gets updated and the mutual link between two adjacent triangulates can be flipped (See below and [29]).These two movies allow us to sample through all possible triangulations for a given   ,   and   and therefore such representation is often referred to as dynamically triangulated surfaces (DTS).DTS is a very popular and successful model to describe shape configurations of interfaces, surfaces, and lipid membranes [30,31].For the purpose of a mesoscopic model, each vertex represents a membrane patch containing hundreds of lipids [29].Discrete geometric operations are used to determine the geometric properties of the surface at each vertex.Several methods are available, each with its own advantages and disadvantages [32,33].In the current version we are using a method based on Shape Operator described in [34].The reason for this choice is that this scheme allows us to obtain, principal curvatures ( 1 ,  2 ), and principal directions ( ̂1,  ̂2), in addition to an associated area, surface normal ( ̂) on each vertex (see SI section 1 and Figure SI-1).This is particularly important when modeling anisotropic proteins and protein-protein interactions (see below).Due to the fluid nature of membranes (the type of membrane that we are considering which encompass most of biological membranes), bending energy will be a function of mean and Gaussian curvatures.In the current version, we will use a discretized version of the Helfrich Hamiltonian that is a function of surface curvature up to the second order (equation 1).
Where  =  1 + 2 2 and  =  1  2 and ,   and  ̅ are bending rigidity, Gaussian modulus and spontaneous membrane curvature respectively (the model parameters of the membrane).As a note, different bending energy functions can be easily adopted in FreeDTS without any significant performance reductions.
To ensure self-avoidance, the edge length is   ≤   ≤   .It has been tested that   =   and   = √3  with a mild constraint on the dihedral angle between two neighboring triangles is enough to ensure self-avoidance of the surface (  is the minimum distance between any two vertices).Nevertheless,   ,   and the dihedral angle can be set by the user in FreeDTS.Due to the flexibility of the edge size, the area of each triangle (in principle) can vary between √3  2 /4 ≤   ≤ √3  2 /4, but on average the total area of the surface remains constant (See Figure SI-2 and table SI 1).However, to control the surface area, the area constraint algorithm can be loaded which couple the system energy to additional term as equation 2.
Where  0 is the targeted area and   is compressibility modules.For membrane patches in a periodic box, the possibility of the box change (dynamic box) is crucial to allow membrane shape remodeling.FreeDTS is equipped with a tension-controlling algorithm that allows for performing simulation at constant (or zero) tension.This algorithm couples the system energy to   = −  where  is frame tension and   is the projected area.We have described the details of this algorithm in [35] and the SI section 2.
Due to the osmolarity effect, closed surfaces, such as vesicles, will be bound by the volume of the media they encompass.It is possible to apply an energy potential as equation 3, which supports both first-and second-order couplings of the volume in FreeDTS.
Where ∆ is the pressure difference between in and outside,  is the vesicle volume,  0 = 1 6 1/2  3/2 (i.e., the volume of a perfect sphere with the area of ) and  is the targeted reduced volume (0 <  ≤ 1).
FreeDTS also provide a more sophisticated volume coupling through the Jacobus van 't Hoff equation to better model osmotic pressure.In this way, the energy of the system is coupled to Where   is the initial volume and ̅  (̅  ) is effective concertation of the inside (outside) solute (more details can be found in SI section 3).
While equation 1 allows for controlling the spontaneous membrane curvature, spontaneous membrane curvature could be global, originating from for example the number of lipid molecules in each monolayer.
FreeDTS allows for controlling global membrane curvature energy to a potential as where  = ∑ 2

𝑁 𝜐 𝜐
,  0 is average membrane global curvature, and   is the coupling constant.Since the area difference between inner and outer membrane (∆) is proportional to  (∆ = ℎ), equation 5 can also be used to control area difference.
Proteins: Proteins are modelling as in-plane inclusions with an in-plane orientation.They interact with the membrane and locally modify the membrane physical and mechanical properties e.g., spontaneous and rigidity.They also interact with one another.In FreeDTS two types of membrane-inclusion interactions exist.First are for proteins with symmetric (or almost symmetric) interactions (Figure SI-3A), such as the one seen for Shiga and cholera toxins [36,37].The interaction energy between a vertex and an inclusion of this type is given by equation 6.
Where   ,   and   are mean, gaussian and area associated with the vertex  and Δ,  0 Δ  are the inclusion model parameters.These parameters has physical meaning, Δ and Δ  can be seen as an increase the in membrane bending rigidity and gaussian modulus and  0 as a local curvature (for a discussion on this, see [8]).
Second protein type are the one which break in-plane symmetry of the membrane (Figure SI-3B), such as dynamin and BAR protein family.Membrane-protein interaction between this type of inclusions is given by equation 7.
|| and  ⊥ are membrane curvature in a direction parallel and perpendicular to the inclusion orientation.Since the used discretized geometric operations allows us to obtain  ̂1,  ̂2 and  1 ,  2 .We can obtain  || and  ⊥ using Eulers curvature formula as where  is the angle between the orientation of the inclusion and principal direction curvature ( ̂1) (Figure 1-A). 1 ,  2 , ||0 ,  ⊥0 are protein model parameters. 1 and  2 are the directional bending rigidities imposed by the inclusion on the membrane. ||0 ( ⊥0 ) is inclusion preferred membrane curvature in the direction of (perpendicular to) its in-plane orientation (for a discussion on this, see [8,29]).
Protein-Protein interactions will be a function of the angle between the two inclusions along the geodesic (Θ = Θ 2 − Θ 1 and see Figure 1-B) and the normal angle between the vertices that the proteins occupy (cos γ =  ̂1 . ̂2 ).However, this angle is defined positive if tip to tip distance is larger than the edge length and negative otherwise.As the results, the simplest interaction between two neighboring inclusions  and  can be written as The first term ( , ) models the isotropic part of the interaction, while the second term models anisotropic interactions (Figure 1-C) and the third is a new term allowing for additional membrane curvature imprint by two proteins related to the value of γ 0 (See figure 1-B). is the least common multiple of the degree of the ,  proteins symmetry in the plane of the membrane. , ,  , and  , are strength coefficients of each energy term.Θ 0 is phase shift.The third term in equation 9 models proteins that induce membrane curvature through dimerization or oligomerization (see Figure 1-D) [38,39].Also, this interaction energy can also be used to model lipid domain formations in multicomponent membranes.More complicated interaction function can be easily adopted in FreeDTS [29].Currently, we do not see any specific need for any other function equation 9 suffice within the resolution of the mesoscopic model.
System evolution: System evolution and the equilibrium properties of the membranes are evaluated by Monte Carlo sampling of Boltzmann's probability distribution.Every Monte Carlo move consists of,   Alexander moves, a trial to flip the mutual link between two neighboring triangles,   vertex positional updates and   inclusion moves (where   is the number of the inclusions in the system) [29].If a system is coupled to tension controlling scheme, then, with the a given probability, there will be also one trial move to change the box size (see SI section 2) [35].There is also a possibility to activate parallel tempering algorithm.This will run multiple replicas at different temperature using OpenMP parallelization to better sampling of the system (for more detail, see SI section 4).

Results
In this section, we present a number of results to demonstrate the power of FreeDTS in capturing membrane mesoscale organization.These include capturing membrane undulation spectrum, proteininduced membrane deformation, tether pulling, curvature-induced protein sorting, and finally protein sorting by real mitochondrial membrane shape (the proteins are hypothetical and do not represent any realistic mitochondrial proteins).
Undulation of framed membranes: Flat membranes are a very common and pragmatic model of membrane segments as the curvature of typical cellular or model membranes is small and can be locally considered flat.FreeDTS allows simulations of flat membranes within a periodic box with constant frame tension.Theoretical analysis, confirmed by experiment and molecular simulations shows that a free membrane undulation spectrum (with the exception of membranes with small bending rigidities [40]) follows equation 10 [41][42][43][44].
Where   and  are effective bending rigidity (renormalized bending rigidity) and frame tension respectively.Accordingly, undulation spectrum for membranes under tension; follows  −4 for large  (small wavelength) and  −2 for small , while for tensionless membranes it follows  −4 .(see Figure 2-A).It is, however, more difficult to obtain the spectrum for membranes in different conditions, such as complex membranes or membranes in confinement.Using FreeDTS, we are able to obtain membrane undulation spectrums under a variety of different conditions.Figure 2-B shows undulation spectrums for three different flat membranes (1) framed membrane with zero tension and 20% coverage by inclusions (4,5) framed membrane with zero tension confined between two sandwiching walls.

Membrane shape deformation by proteins:
Proteins are one of the main drivers of membrane deformations [45,46].FreeDTS allows for exploring the membrane shaping and sorting capacity of membrane proteins.Figure 3 shows several examples of membrane deformation by different membrane proteins at the mesoscale.For small protein-protein interactions (above a certain  threshold [35]) buds form and the membrane is divided into two domains of protein rich and protein poor.However, for larger protein-protein interactions, proteins cluster and form buds, leaving insignificant number of proteins on the flat surface.
Pulling a membrane tether: One of the common procedures to deform membranes is pulling a tether (nanotube) using, for example, optical tweezers [12,47,48].Since FreeDTS allows for membrane simulations with constant tension, it can also be used to pull a membrane tether by applying a harmonic force to a vertex (see video 1 and Figure 4-SI).Such a system can be used, for example, to study protein sorting [12] or to obtain membrane bending rigidity [49].Figure 4 illustrates how FreeDTS can be used to pull a membrane tether as well as to study curvature-driven protein sorting that is dependent not only on protein curvature but also on protein-protein interactions.
Other shapes: Vesicles, high genus membranes and tubes: The type of membranes (surfaces) that can be simulated with FreeDTS is not limited to flat membranes.In principle, any triangulated surface that is closed (even if it is through a period box) can be handled.A prime example is a spherical surface such as vesicles.Figure SI-5 shows the transition of a spherical vesicle from prolate-to-oblate and oblate-to-stomatocyte by volume reduction [50].Figure 5 A-B also shows how both types of inclusions could induce tubular membrane invaginations in a vesicle.An important characteristic of membranes is the topology of their surfaces.While membranes are flexible and easy to bend, their surface topology tends to remain constant, as topological changes require membrane fission and fusion processes that are restricted by a high energy barrier.Topology is characterized by the topological genus g, which counts the number of handles attached to a sphere.For instance, g=0 for a sphere and g=1 for a coffee mug.Although organelle membranes, such as the ones found in mitochondria and Golgi apparatus, exhibit high-genus topology, there has been only a limited number of studies conducted on membranes with non-spherical topologies, which have considered only simple, low-genus membranes [51][52][53].FreeDTS allows for the exploration of these surfaces.Just as an example, Figure 6-C-D shows how a closed surface with topological genus transformed by some inclusions into a stomatocyte structure.This process might be relevant for nuclear membrane formation and assembly of nuclear pore complex driven by membrane curvature [54].Another type of membrane structure is tube, a surface periodic in one direction (Figure 6-A).Tubes can be used to study processes on the segments of, for example, smooth endoplasmic reticulum.Figure shows how FreeDTS can be used to explore tube deformations by different means.

Confined membranes:
In FreeDTS, membranes can be contained within spaces of various shapes and sizes.In the current version, a membrane can be confined between two sandwiching walls, an ellipsoidal shape, an ellipsoidal shell, or a block (Figure 2-B and Figure SI-6).
Protein sorting on realistic membranes: Advances in experimental techniques such as cryo-electron tomography now allows for resolving membrane shape of a full organelle or even a cell [2,55].When the density and mesoscale parameters of each biomolecule are provided, FreeDTS can use these structures as input files to determine biomolecular organization.Figure 7 shows a DTS simulation of a real inner mitochondrial membrane with two types of proteins.It is important to note that these proteins are not representative of realistic mitochondrial proteins and have only been used as a proof of concept.Proteins are sorted in different regions of the membrane based on the type of protein and their interactions with one another.In this simulation, only the proteins have the possibility of organizing and the triangulated mesh does not change (shape change).

About the source code
FreeDTS is implemented in C++.It is self-contained and does not depend on any additional library.It is very well-link to TS2CG through input and output meshes and can be expanded to be linked with TS2CG on the fly.Also, it is easy to compile and run (for more see the tutorial section in the manual or the tutorial section on the GitHub).The energy unit is    and the unit length is   which is the minimum distance between any pair of vertices.

Discussion
Mesoscopic simulations are a necessary tool to explore many important features of membrane-involved cellular processes [20].In particular DTS approaches have shown to be very robust in capturing spatial and lateral organization of biomembranes.FreeDTS as shown in the current manuscript and our previous works [12,24,35,40] offer a great technology for exploring these processes.Also, as shown in Figure 7, FreeDTS can be used to obtain the organization of biomolecules on experimentally resolved membrane shapes.In principle, protein density and protein-protein interactions of main membraneshaping proteins can be obtained from experiments such as cross-linking mass spectrometry [56] or from higher-resolution simulations [29].Other mesoscopic model parameters (such as induced local curvature) can be obtained from molecular simulations [37] or from a top-down approach through fitting to experimental data in controlled biophysical experiments [12,57].Once these parameters are available, FreeDTS can be used to obtain the organization of these biomolecules.Furthermore, the equilibrated output of the FreeDTS can be processed directly by TS2CG to create the system structure with a molecular resolution [27,28] for full cell or organelle simulations [11,58].
To obtain biologically relevant information, models such as DTS may appear to be highly dependent on calibrating their parameters to start with.It should be noted, however, that this is not entirely accurate.
Even without any knowledge of membrane shaping protein structure, DTS simulation can provide some knowledge about their structure by tuning the model parameters against macroscopic biophysical experiments [24,57].
FreeDTS can be expanded (and is in our agenda) in many different ways to become a versatile tool to explore membrane remolding processes.For example, during a simulation in FreeDTS, the number of vertices, triangles, and edges of the triangulated mesh remains constant.In addition, the membrane surface is considered closed, i.e., no possibility of an open edge, such as a hole, exists.Due to these limitations, FreeDTS does not allow for several membrane-related cellular processes, including membrane fission, fusion and membranes with large holes [59,60].Fission and fusion can be captured by developing a dynamic topology algorithm such as ref.[61][62][63].It is, however, more challenging for a membrane with open edges, especially since the energy of the edge also depends on the geodesic curvature of the edge [64].This indeed requires theoretical developments on how to evaluate these geometrical variables, e.g., edge geodesic curvature on triangulated surfaces.
The Helfrich Hamiltonian (the membrane bending energy used by FreeDTS) is a function of principal curvatures up to the second order of principal curvatures assuming the energy should be invariant under in-plane rotation.It has been demonstrated that this bending energy is very excellent for describing large-scale membrane shapes [1].Nevertheless, at the mesoscale, the effects of higher order may become significant as cellular membranes often exhibit curvature radii that are comparable to the thickness of the membrane [45,65].However, the challenge of resolving this problem is more theoretical in nature, i.e., finding a suitable bending energy function.Once obtained, its implementation within FreeDTS is rather straightforward without significant increase in computational cost.
FreeDTS currently runs most movies sequentially and only uses one CPU core (with the exception of parallel tempering where each replica is run on a separate core).The moves can in fact be parallelized to some extent [26] in order to explore membranes with more vertices or obtain results more quickly for smaller systems.However, this may not lead to better membrane exploration.It is, for example, not amendable to properly capture configurations of membranes larger than a few folds just by allocating more resources, since the required steps to properly sample the configurations increase rapidly (nonlinearly) with system size.Additionally, membranes exhibit shapes with very different configurational structures, separated by many large energy barriers, but energetically close (degenerate).Therefore, often it is more practical to perform many replicas to capture these possible configurations, or use enhanced sampling methods e.g., parallel tempering and Hamiltonian replica exchange methods that are much easier to parallelize efficiently for DTS-like methods.
In biomembranes, different phases with specific molecular compositions coexist.These phases can recruit membrane proteins which have been postulated to be essential for many biological processes such as signaling and endocytosis [66,67].Currently, FreeDTS is capable of capturing such phenomena, but the total surface area of each phase remains constant.In a biological system, however, these phases are dynamic and can appear, disappear, shrink, or expand depending on the physiological conditions.In the future developments we will address this limitation.
Finally, the current software can be enhanced by equipping it with models of, for example, cell walls, cortical actin cytoskeleton and active processes (breaking detailed balance) or by coupling the membrane mechanics to kinetic models of the cellular processes [68] to provide a kinetic-mechanic model of the cell membranes.

Conclusion
We have presented FreeDTS software to perform computational research on biomembranes at the mesoscale.The model parameters of proteins have physical meaning and in principle can be calibrated using finer scale simulation techniques e.g., all-atom and coarse grain molecular dynamics, or through a top-down approach through experimental data.Several algorithms are included in the software that allows for simulations of framed membranes with constant tension, vesicles with various fixed volume or constant pressure difference, confined membranes into the fixed region of the space, constant fixed global curvature, and application for external forces on regions of the membrane.In addition, the software allows one to turn off the shape evolution of the membrane and only explore inclusions organization.This allows us to take realistic membrane shapes obtained from cryo-electron tomography and obtain heterogeneous organization of biomolecules which can be backmapped to finer simulations models.In addition to its use for exploring many biomembrane application, this software brings us a step closer to simulating realistic biomembranes with molecular resolution.

Figure 2 :
Figure 2: undulations spectrum for different membranes obtained using FreeDTS.(A) undulation spectrums for membranes with different tensions; for large  (small wavelength) it follows  −4 , while for small , it follows  −2 .(B) undulations spectrum for membranes in different conditions; (blue) framed membrane with zero tension and 20% inclusion coverage (red) membranes between two sandwiching walls with the thickness of 2ddts, (green) membranes between two sandwiching walls with the thickness of 4ddts.

Figure 4 :
Figure 4: Tether pulling from membranes.(A) A tube was pulled from a free membrane coupled with a constant tension of  = 2[    −2 ]; (B) protein sorting due to curvature.Protein coverage 10%, type

Figure 7 :
Figure 7: Experimentally resolved real mitochondrial membranes can be load in FreeDTS to obtain protein organizations.