CurD: A Tool for Diffusion Analyses on Curved Membranes

Curved membranes are abundant and functionally relevant in living matter, yet they have eluded computational studies due to methodological limitations. With multiple approaches available for setting up simulations on such curved membranes, there is a growing need for efficient and versatile tools to analyze their outcomes. Here, we present CurD, a freely available tool for analyzing the diffusion of membrane constituents along curved surfaces. The tool efficiently uses the Vertex-oriented Triangle Propagation algorithm to compute geodesic distances significantly faster than conventional implementations of path-finding algorithms, while providing a friendly command-line interface. With CurD, we resolve the effects of curvature and lipid packing densities on the diffusion of lipids on two curved membranes—one with only mean curvature and another with both Mean and Gaussian curvature. We find that Gaussian curvature plays a surprisingly small role, whereas mean curvature or packing of lipid headgroups largely dictates their mobility.

Molecular simulations are starting to achieve sufficient length-and time-scales that allow us to efficiently probe the compositional and topological complexity of real cellular membranes. [1][2][3] Notably, many key functions related to these membranes encapsulating either the entire cell or its organelles involve significant local membrane curvature. 4 The plasma membrane has specific invaginated signaling platforms, 5 the complex and dynamic topography of the mitochondrial inner membrane is involved in numerous processes, 6 and the cells store energy in lipid droplets that bud from endoplasmic reticulum. 7 In these processes, curvature is generated by both lipids 8 and proteins, 9 yet they both also sense curvature, leading to their spatial sorting with functional implications. 10,11 Adequate modelling 12 and analysis 13 of such non-planar lipid bilayers requires the development of algorithms that take into account the geometry of the membrane. A particularly fascinating property impacted by membrane curvature is the lateral diffusion of membrane components, which has a direct influence on the interpretation of experimental results obtained from techniques such as Fluorescence Correlation Spectroscopy (FCS) and Single Particle Tracking (SPT) that often assume planar geometries. [14][15][16] Notably, the effects of the underlying curvature can be falsely interpreted by slowed-down or even anomalous diffusion, [16][17][18] calling for further clarification by alternative approaches such as computer simulations.
Biological membranes undergo thermal fluctuations affecting their curvature. These fluctuations are generally considered to slow diffusion down by increasing the thickness or extending the geometric path. 19 Taking them into account would require the treatment of changing surfaces. Instead, here, we restrain ourselves to work in the static membrane limit, 20 which is an excellent approximation for structures such as the cristae in the endoplasmic reticulum (ER) or auditory outer hair cells. 21,22 Another aspect of the classification of membrane surfaces is based on their curvature. The curvature of a surface can be completely described by two scalar fields, the Gaussian curvature K(r) and the mean curvature H(r). Caveolae and budding vesicles with nonzero Gaussian curvature have radii of few dozen to a hundred nanometers. 23 Tubular structures such as those in the mitochondrion are examples of developable surfaces, that only have non-vanishing mean curvature. 24 The ER is also of great interest due to its complex membrane topology of folded membranes and highly curved regions. 25,26 Several studies employing continuum methods have investigated the effect of both mean and Gaussian curvature on biological membranes, 20,27-29 resulting in two major conclusions. First, at the continuum level, the surface diffusion of particles can only depend on Gaussian curvature K, and not on mean curvature H. This is related to the fact that developable surface-surfaces with vanishing K-are isometric to a plane. 27 Obviously, this is only true for diffusion along the surface, whereas any curvature, K or H, will affect the diffusion coefficients measured in experiments that assume a planar membrane. Secondly, on surface regions with K > 0 (elliptic paraboloids) diffusion is slowed down, while on region with K < 0 (hyperbolic paraboloids) it is sped up. 27,28 In addition to these observations, the ratio of the real and projected long-time diffusion coefficients are shown to closely follow a so called area-scaling law 20 under a broad range of conditions. Details about the area-scaling law can be found in the Supporting Information.
A major shortcoming of all these continuum methods is that they cannot account for the effects of lipid packing or the inclusion of proteins or other membrane heterogeneities.
These effects carry the possibility of mean curvature influencing surface diffusion through heterogeneous membrane structures. 29 While such effects are explicitly inherent in molecular simulations, their analyses presents another kind of drawback. Namely, the majority of the current analysis tools fall short of dealing with the changing membrane normal in curved membranes. 30 Here, we present an algorithm and a software that implements it, coined CurD, for the calculation of Mean Square Displacement (MSD) in molecular simulations of curved membranes. We apply the method to coarse-grained Martini 31,32 simulations of phospholipid bilayers forming a vesicular bud-like membrane protrusion ("Budded") and an undulating wave-like surface ("Wave"). Although slightly over-emphasized, the mean and Gaussian curvatures of our simulated systems are not far from those present in some biological systems.
All relevant details of the simulation can be found in the Supporting Information. The simulated systems and their curvatures are shown in Fig. 1. The membrane in the "Wave" system has only mean curvature, H, hence it is isomorphic to a plane. This is not true for the lipid bilayer in the "Budded" system, as it also possesses Gaussian curvature, K. An important attribute of Gaussian curvature is its insensitivity to the orientation of the constituents of the bilayer. A bowl-shaped (saddle-shaped) region always has positive (negative) Gaussian curvature, irrespective of whether the bud bulges into the cytoplasmic or the extracellular space. The curvatures (bottom panels of Fig. 1), while properly capture the topology of the system, exhibit some irregularities due to the use of mesh surfaces (see the Supporting Information). Even though their relative sign depends on the chosen leaflet, a strong correlation between the magnitudes of H and K is clearly apparent.
The actual displacements of particles constrained to move on a surface are best described using geodesic distances. 27,28 Geodesic distances can be efficiently computed on mesh surfaces using specialized methods, even without explicitly storing or constructing the corresponding shortest path. 33 One such algorithm is the Vertex-oriented Triangle Propagation (VTP) of Qin et al. 34 that simultaneously computes all distances from a source vertex to all other vertices on the mesh. Therefore, as the first step of our approach, we create a triangular mesh surface from each of the leaflets of the simulated membranes, using the center of mass of the lipid. Then, we use the position of the closest vertex as an approximation of the surface position of the lipid centers of mass. For sufficiently fine meshes, the error introduced by this discretization is negligible. To make full use of VTP algorithm and avoid multiple evaluation of the individual source vertices, the displacements must be grouped in a particular way. Triplets of integers (start vertex, end vertex, lagtime) representing the path of particles are written to a binary file and suitably sorted. By performing the grouping, all paths starting or ending at a given vertex form a contiguous block, irrespective of the Radial distance (nm) upper lower After preparing a triangular mesh surface, all particle displacements, that is, initial and final particle positions (white circles) separated by lagtime are mapped onto the surface ensuring that both the starting and end point are on the mesh, as close to the origin as possible (gray circles). Particle displacements at large lagtimes require taking into account more periodic images. The resulting v i and v j pairs of mesh vertices are written to a binary file along with the corresponding lagtime, ∆. To enable the efficient use of the VTP algorithm, all pairs are ordered so that the lower vertex index appears first (see for example the numbers highlighted in blue), and subsequently ordered along the first column. This way, all particle displacements involving a given v i appear as a single contiguous block. Finally, the VTP algorithm is needed to be called on each source vertex v i only once to evaluate all distances in the block. These evaluations are independent, therefore are readily parallelized, as done in the current implementation. The computed displacements are assigned to both vertices while also taking account of the associated lagtime.
To investigate the differences between the various diffusion measurement methods, we evaluated the Mean Square Displacements both by projecting the motion of the particles onto the macroscopic plane of the membrane-as typical in numerous experimental approachesand by using our approach presented above. The projected and surface MSD values as a function of the position along the y axis or the radial distance in case of the "Wave" and "Budded" systems are shown at a few selected lagtimes in Fig. 3 and 4, respectively.  Figure 3: Two dimensional ("Flat", discarded z coordinate) and geodesic ("Upper" (for upper leaflet)) Mean Square Displacement values at selected lagtimes (∆ = 2, 53, 453 ns from bottom upward) as a function of the position along the y axis in the "Wave" system. The values in the "Flat" system are averaged across both leaflets. The "Lower" leaflet is a shifted version of the "Upper" leaflet, and as such, it is omitted. The images on the right illustrate the distributions of MSD values at the corresponding lagtimes.
In both systems, using the projected displacements has a major effect on the apparent motion of the molecules due to ignoring their motion along the axis perpendicular to the plane of the membrane. The magnitude of this effect is directly related to local orientation of the membrane segments, thus the projected displacements in the almost vertical regions of the "Wave" and "Budded" membranes underestimate the actual displacements the most. This phenomenon is purely geometric, and can have a significant impact on results obtained from experiments. The 2D MSD profiles (black lines) get progressively smoother with increasing lagtime due to the mixing of molecules originating from regions of different curvatures. gradually mix as the lagtime increases, producing uniform MSD profiles at large enough lagtimes. To better quantify these effects, we computed the diffusion coefficient distributions using the conventional 2D calculation method and the geodesic-based approach presented here. Furthermore, we also decomposed the latter based on local curvature (see Supporting Information for details). The diffusion coefficients at the individual mesh points were calculated at the lag time of 53 ns, instead of separately applying linear regression to them.
It must be noted, that computing a diffusion coefficient distribution on surface mesh points is not strictly equivalent to the per-particle distribution. However, the almost uniform density of lipid centers of mass (see while the Flat region of the "Budded" system corresponds to the "apparently fastest" domain in the Projected distribution (both are around 6 · 10 −7 cm 2 /s), the Flat part of the "Wave" is conclusively faster than the projection. This is in agreement with the orientation of the planar membrane regions in the two systems (cf. Fig. 1).
Contrary to continuum theories predicting the lack of influence of mean curvature H on surface diffusion, 27-29 the "Wave" system seems to exhibit clear correlations between where the headgroups are more compressed. The observed discrepancy between continuum predictions and molecular simulations is a direct consequence of lipid packing effects, which are usually not included in the continuum models. 29 Even though the differences are minor in the systems studied here, they can be more significant with e.g. a more complex lipid mixture where the different lipid species are sorted by curvature.
Because there is a fundamental asymmetry in the curvature of its leaflets, in the analysis of the "Budded" system the individual leaflets must be treated separately. In case of the upper leaflet (blue curves in Fig. 4), positive Gaussian curvature, K > 0, seems to correlate with faster diffusion, while regions with K < 0 exhibit slower diffusion. This behavior is verified by the analysis demonstrated in Fig. 5, and again goes against the conclusion drawn from continuum theories, where K > 0 results in slower diffusion and K < 0 causes faster diffusion. 28 Similarly to the mean curvature, the speedup in the upper layer can be ascribed to a less dense packing of lipids, at least at the level of headgroups.
Importantly, Gaussian curvature is insensitive to the direction of the membrane normal: bowls have positive and saddles have negative Gaussian curvature. Therefore, if the diffusion only depended on the Gaussian curvature, one would expect similar tendencies in both upper and lower leaflets. This is clearly not the case, as the lower leaflet seems to follow the theoretical prediction presented above. ? Consequently, the theoretically predicted role of Gaussian curvature is not the only factor determining surface diffusion in molecular systems, mean curvature also must play an important role. What is more, in our "Budded" system the absolute magnitude of Gaussian and mean curvatures positively correlate, while the mean curvature also encodes the orientation of the lipids. Thus, the mean curvature seems sufficient to explain the behavior observed in our simulations. Joined with the headgroup densities, the arising picture nicely follows that of the "Wave" system possessing only mean curvature. This line of reasoning is further supported by Fig.?? containing the correlations of the variables of interest.
To conclude, we have developed an algorithm and implemented it to analyze the diffusion dynamics along curved membrane surfaces. We have applied our method to two simulated membranes with different topologies-one with only mean curvature and another with ad-