Active membrane recycling induced morphology changes in vesicles

Membranes of organelles in the intracellular trafficking pathway continuously undergo recycling through fission and fusion processes. The effect of these recycling processes on the large-scale morphology of organelles is not well understood. Using a dynamically triangulated surface model, we developed a membrane morphology simulator that allows for membrane trafficking, and analyzed the steady state shape of vesicles subjected to such active remodeling. We study a two-component vesicle composed of 1) active species which can have nonzero spontaneous curvature and participate in the recycling and 2) inactive species which do not participate in the recycling. We obtain a plethora of steady state morphologies as a function of the activity rate, spontaneous curvature, and the strength of interaction between species. We observe that morphology changes, as a function of rate of activity, are diametrically opposite for the two signs of the spontaneous curvature, but only have a weak effect on its magnitude. The interplay between the in-plane diffusion, the activity rate, and the spontaneous curvature are shown to determine the vesicle morphology at the steady state. It is shown that the spontaneous curvature and activity inhibits the formation of clusters of active species on the surface. We carry out linear stability analysis of a continuum model and show that the spherical shape of a vesicle is indeed unstable when subjected to active membrane recycling above a certain activity rate.

1 Introduction

Organelle membranes exhibit a wide variety of morphologies that are believed to be crucial for their function [1–3]. While peroxysomes and lysosomes are spherical, the endoplasmic reticulum (ER) has sheet-like and tube-like regions [4]. Adjacent sheets of the rough ER are connected through parking-garage like structures called “Terasaki spiral ramps” [5]. The Golgi apparatus consists of flattened membrane compartments called the cisternae and the mitochondria is made of an outer membrane and a highly convoluted inner membrane. There are a number of structural motifs that are commonly observed in these organelle membranes such as sheets, tubules, or interconnected networks [6, 7]. The dimensions of such features have managed to remain conserved throughout evolution. For instance, the diameter of the tubules in the cristae junction of the inner mitochondrial membrane are strictly maintained at 28 nms [7]. To know the functional organization within a cell, it is essential to understand the mechanisms that are responsible for generation and maintenance of organelle morphologies. Together with the fact that organelles dissolve and reassemble to take nearly the same structure during each cell cycle suggests that the generation and stabilization of organelle morphology is the product of a concerted effort. Although there are several molecular players and complex feedback loops involved in the processes of organelle morphogenesis, it is possible that these process can be understood in terms of a few self-organizing principles [8–10].

Different models and mechanisms, that can spontaneously generate such complex shapes have been proposed [6, 7, 11, 12]. One such mechanism, which is specially relevant for the organelle membranes in the membrane trafficking pathway, is the regulated rates of fission and fusion of vesicles. This mechanism is reported to be a crucial factor for maintaining the morphology and composition of these organelles [13]. This process of fission and fusion are active in the sense that they require external energy input through hydrolysis of ATP/GTP and are assisted by special proteins such as SNAREs [14]. It is accepted that such active recycling of essential proteins and lipids are imperative to maintain the composition and function of the membranes. However, an understanding of the effect of such active processes on the morphology and organisation of organelle membranes is still lacking.

In a planar membrane subjected to active deposition and evaporation of membrane materials, it was found that the activity can result in a dynamical instability resulting in long tubular shapes of the membrane [15]. Ramakrishnan et al. [16] studied the spontaneous morphological changes of a vesicle under adsorption and desorption of active curvature modifying agents. It was seen that spontaneous curvature fluctuations can lead to shape transitions through an activity-dependent tension. On a similar line, Goutaland et al. explored the stability of a free standing membrane subjected to active and thermal spontaneous curvature fluctuations induced by binding/unbinding of curvature active components [17]. Both these models however do not consider actual area changes of the membrane due to deposition and evaporation or binding and unbinding of curvature active agents. Another recent work investigates the role of nonequilibrium recycling involving actual material exchange on the conformations of 1D active polymer models [18].

In a work that is closely related to what is presented here Tachikawa and Mochizuki carried out coarse-grained simulations that capture the reassembly of Golgi apparatus as a result of vesicle aggregation, fusion, and membrane shape relaxation [19]. While their work is focused on reproducing the shape of Golgi, our emphasis is on understanding the steady states shapes that results from energy consuming transport of membrane components and curvature modifying agents.

In-plane organisation of membrane components, that are curvature active, can have significant effect on the morphology of the membrane. Coarsening dynamics and domain size distributions of membrane domains undergoing active recycling is investigated in Ref. [20]. It was shown that active stresses in a thin sheet, that depend on the concentration of a stress-regulating molecular species which are dynamically changing in response to flows on and deformations of the surface, can result in spontaneous generation of nontrivial surface shapes and shape oscillations [21]. Another factor that can contribute to the morphology of a vesicle is the rate of growth of volume relative to that of area, which was addressed recently [22].

In contrast, to the above mentioned works, the work presented in this paper addresses the question of the morphology of vesicles, with continuous flux of materials into and out of its surface, at the steady state. We explore the effects of fluctuations in local density and spontaneous curvature using a new method based on Dynamical Triangulation Monte Carlo (DTMC) and local density based potentials. DTMC is a useful technique to explore the non-axisymmetric shapes of membranes which are analytically intractable [23]. Here, we introduce a modified DTMC method, that now allows for addition and removal of membrane material, and explore the steady-state morphology of the vesicle as a function of material addition and removal rates. Here, the recycling of membrane material is considered as an active (energy consuming) process. We observe that such nonequilibrium processes lead to shape instabilities when material exchange is coupled to local membrane curvature.

2 Simulation method

The vesicle is represented by an infinitely thin triangulated surface where each triangle represents a patch of lipid bilayer. The length unit of the simulation is set by the diameter of the vertex beads. The configuration of the system is specified by the number of vertices, their position vectors $(\vec{\mathbf{R}})$, the connectivity of the mesh $(\mathcal{T})$, and the nature of each species $(\vec{\mathit{\rho }})$. From the position vectors for the vertices, one can directly calculate the area and normal for each triangle on the surface. Using this, an area associated with the vertex A_i and a surface normal n_i can be evaluated by appropriate weighted averaging. One can then define mean curvature at each vertex using [23],

$H_i=\frac{1}{{\sigma }_i}{\mathbf{n}}_i\cdot \sum _{j\left(i\right)}\frac{{\sigma }_{ij}}{l_{ij}}\left({\mathbf{R}}_i-{\mathbf{R}}_j\right).(1)$

Here, l_ij is the distance between the two vertices located at R_i and R_j, σ_ij is the length of the bond in the dual lattice, σ_i is the area of the dual cell of vertex i. In Eq. 1, j(i) denotes all vertices j that are neighbors of vertex i.

The equilibrium dynamics of the membrane is simulated using the vertex displacement and bond flip MC moves, which respectively captures the local radial fluctuations and the in-plane lipid diffusion. More details about the DTMC technique is available in Refs. [23, 24]. The process of fusion of a smaller vesicle is considered to add more membrane material (lipids and proteins). For the large vesicle, the process of fusion essentially increases its area. Similarly, the process of fission removes a small amount of area from the vesicle. In simulations, the addition and removal of area from the vesicle is captured with the help of MC moves for addition and removal of a vertex respectively. The conventional DTMC method has a restriction on the length of the tether connecting neighboring vertices. In order to introduce vertex addition and removal MC moves, we need to relax this restriction on the maximum length of the tether. We use a local density based attractive interaction potential between the vertices to keep the triangulated surface intact [25]. The local density at each vertex is defined using,

${\rho }_i=\sum _{j\left(i\right)}{f}_{\text{cut}}\left(l_{ij}\right)(2)$

where the f_cut is a cutoff function defined as

$f_{\text{cut}}\left(r\right)=\left\{\begin{array}{cc}\exp \left\{A\left(1+\frac{1}{{\left(r/r_{\text{cut}}\right)}^n-1}\right)\right\}\hfill & \left(r<r_{\text{cut}}\right)\hfill \\ 0\hfill & \left(r\ge r_{\text{cut}}\right)\hfill \end{array}\right.(3)$

with n = 12, r_cut = 2.1, r_half = 1.8, $A=\log (2)\left[{\left(r_{\text{cut}}/r_{\text{half}}\right)}^n-1\right]\approx 3.715$. The total energy of the surface includes three contributions—the bending energy (E_bend), the energy from local density based attractive potential (E_att), and the interaction energy (E_int),

$\mathcal{H}=E_{\text{bend}}+E_{\text{att}}+E_{\text{int}}(4)$

$\begin{array}{ll}\hfill & =\frac{\kappa }{2}\sum _{i=1}^{N_v}{\left(2H_i-C_p{\varphi }_i\right)}^2{A}_i\hfill \\ \hfill & +\frac{ϵ}{4}\sum _{i=1}^{N_v}\log \left(1+e^{\left[-4\left({\rho }_i-{\rho }^{*}\right)\right]}\right)-J\sum _{i=1}^{N_v}\sum _{j\in i}{\varphi }_i{\varphi }_j\hfill \end{array}(5)$

where ρ* = 8 is taken as the preferred local density and ϵ = 6 is used for the strength of the potential. The two parameters of the attractive potential, ϵ and ρ* was chosen such that simulations reproduce the edge length and coordination number distributions similar to those obtained using the tether potential (Figures 1A,B). The last term in Eq. 5 allows for interaction between the species on the surface. Here, ϕ_i denotes the nature of the species occupying the vertex; ϕ_i = 1 indicates that an active species occupies vertex i and ϕ_i = 0 indicates that a passive species is occupying vertex i. The interaction is attractive if J > 0 and repulsive otherwise.

FIGURE 1

FIGURE 1. The histogram for edge length (A) and coordination number (B) distribution for different parameter choices for the local density based potential in our modified DTMC technique. The histograms for the conventional DTMC technique is shown as “Tether” for comparison. (C) Schematic representation of fission/fusion process in our model. The fission/fusion process takes place with rates r⁻/r⁺ respectively. (D) Schematic showing the local configurational changes during vertex addition and removal MC moves in the simulation.

We assume that the fission and fusion processes are actively assisted by and happens only with the help of protein machineries that are present on the smaller vesicle. In simulations, the presence of the protein machinery is taken into account by tagging a vertex as active. Thus, a two component membrane is considered by tagging each vertex as either active or passive to indicate two types of species. Since the protein machinery residing on the small vesicle can also induce membrane curvature, the active species is assigned a preferred curvature C_p. On the other hand, the passive species is assumed to have no preferred curvature. The equilibrium dynamics of the membrane is simulated using the vertex displacement, edge flip, and protein diffusion MC moves. To simulate active processes of fusion and fission, we introduce two new MC moves—vertex addition and vertex removal respectively (Figure 1D). Addition of a vertex to the triangulated surface is considered here as equivalent to fusion of a small vesicle to a large membrane surface. Since the typical size of transport vesicles in the intracellular environment is around 50 nms, this is considered as the size of a vertex. Similarly, removal of a vertex is considered equivalent to fission of a small vesicle from the membrane surface. Note that these 2 MC moves break the detailed balance condition and therefore drives the simulations away from equilibrium. The details of the implementation of the two active MC moves are described below.

Each vertex addition MC move increases the number of vertices in the triangulated surface by one and thereby increases the total area. We perform the vertex addition move in two ways—three neighbor addition and six neighbor addition. In the former, we choose a random face on the surface and attempt to add a vertex at a distance d ∈ (0, d_add) above (outside the surface) it along the face normal. The initially chosen face is removed and three new faces are added along with the new vertex. The coordination number of the newly added vertex in this case is three.

In six neighbor addition, a random face is chosen and a new vertex placement is attempted above it at a distance d ∈ (0, d_add) along the face normal. If the move is accepted, the chosen face and the three edges of the chosen face are removed. This results in six vertices surrounding the newly added vertex. Six new faces and edges are added by connecting the newly added vertex with the six vertices. Thus the coordination number of the newly added vertex is six in this algorithm. The vertex addition attempt is rejected if it violates the self- avoidance criteria. We observed that when d_add < = 0.7, almost all the addition attempts were rejected due to self-intersection. For d_add > 1.1, most of the newly added vertices results in highly obtuse triangle faces wherein the curvature computations are not accurate. Within the range 0.7 < d_add < 1.1, we observed that the results of our simulation do not vary significantly. For the results presented here, we used parameter value d_add = 0.9. Note that vertex addition adds “bumps” on the surface (Figure 2) that models the fusion of a small vesicle to the membrane surface. The sharp edges of these bumps are short-lived compared to the timescale for morphology change of the vesicle.

FIGURE 2

FIGURE 2. (A) Activity makes the spherical shape unstable when C_p ≠ 0.0. When C_p > 0, the reduced volume of the vesicle undergoes a sharp decrease as the activity rate is increased, indicating a sharp change in the vesicle morphology from spherical to flattened shapes. When C_p < 0, the reduced volume undergoes a diametrically opposite transition with activity rate, corresponding to a morphology change from tubular to spherical shapes. The blue filled star, which is overlapping with the red symbols, corresponds to simulations where C_p = −0.1 with no in-plane diffusion of active species. The interaction strength J is set to zero in all these cases. (B) A qualitative morphology phase diagram indicating the parameter ranges where the vesicle assumes a flattened/tubular shape (dark regions) and spherical shape (light regions). The dashed vertical lines show three cross sections corresponding to three different regimes of morphology changes observed in our model—tubular to spherical when C_p < 0 (blue), always spherical when C_p = 0 (black), and spherical to tubular/stomatocyte when C_p > 0 (red). (C) Snapshots of morphologies observed at the steady state for—C_p = −0.1 (blue outline), C_p = 0 (black outline), and C_p = 0.1 (red outline)—at different activity rates. The parameter values used in the simulation are $\bar{N}=1200$, ζ = 0.1, and κ = 10k_BT.

As vesicle fusion is an active process driven by energy supply from the hydrolysis of ATP molecules, we assume that the statistics of vertex addition is not controlled by the Hamiltonian. Rather, we use the following probability for accepting the addition move

${\mathcal{P}}_{\text{add}}={\alpha }_{+}\frac{1}{\eta +e^{\zeta \left(N_v-\bar{N}\right)}}(6)$

where N_v is the number of vertices in the vesicle at any instant, $\bar{N}$ is a parameter determining the mean size of the vesicle, and ζ controls the width of size fluctuations. Such an assumption is reasonable because it is generally considered that there are regulatory mechanisms which are responsible for controlling the size of organelles. The parameter α₊ N_v is the number of attempts to add a vertex in 1 MC Sweep (MCS). Fusion process is always assumed to add an active species. Since a vertex is always added from outside, this also implies a positive curvature at the added vertex independent of the value of preferred curvature. The fission process is captured by the vertex removal MC move, wherein a randomly chosen vertex with an active species is removed from the surface. The resultant polygon shaped opening on the surface is then retriangulated while maintaining self-avoidance. The probability for vertex removal is taken as

${\mathcal{P}}_{\text{rem}}={\alpha }_{-}\frac{1}{1+e^{-\zeta \left(N_v-\bar{N}\right)}},(7)$

where α₋ N_v is the number of attempts to remove a vertex in one MCS. The parameters α₊ and α₋ respectively affects the rate of fusion and fission processes with respect to the rate of bending relaxation. In order to satisfy the condition that the addition and removal probabilities are equal when $N_v=\bar{N}$, we choose the parameter η = 2α₊/α₋ − 1 in Eq. 6. For simplicity, we assume that α₊ = α₋ = α throughout this work.

Note that, we have neglected the Gaussian curvature term in the free energy of the vesicle as we are only allowing vesicles of spherical topology. In order to maintain self-avoidance of the surface during addition and removal moves, we performed collision checks between vertices, edges, and faces during each active MC attempt. Collision checks are optimized by using an Octree data structure to represent and perform spatial queries using the geometric objects.

3 Linear stability analysis in the presence of membrane recycling

In this section, we derive the dynamical equations for a nearly spherical vesicle subjected to active membrane recycling. A two component vesicle, which has an active component that takes part in the active processes of fission/fusion and an inactive component that does not, is considered. The active component may represent either curvature inducing proteins or small lipid vesicles having a preferred curvature C_p. We derive the coupled dynamical equations for the vesicle radius field and the density of the active component, that includes the equilibrium bending relaxation, the diffusion of the active components and the effects of recycling. The stability of the spherical shape of such a vesicle is analyzed with respect to the activity strength.

The vesicle is modeled as a nearly spherical 2D surface. The configuration of the vesicle is represented using spherical parametrization where the distance from the centre of mass, to each point on the surface, is given by,

$\vec{r}\left(\theta ,\phi \right)=r_0\left[1+u\left(\theta ,\phi \right)\right]\widehat{e_r}.(8)$

Here r₀ is the radius of a reference sphere and u is the scaled deviation in radius at point (θ, φ). Assuming that the deviation (u) from the reference spherical shape is small, the relevant local geometric quantities—square root of determinant of the metric tensor and the mean curvature—upto 2^nd order in u are,

$\sqrt{g}\approx \left(1+2u+u^2+\frac{{\left({\nabla }_{\perp }u\right)}^2}{2}\right),(9)$

$H\approx -\frac{1}{2r_0}\left(-2+2u+{\nabla }_{\perp }^{2}u-2u^2-2u{\nabla }_{\perp }^{2}u\right),(10)$

where ∇_⊥ is the surface gradient operator defined as ${\nabla }_{\perp }={\widehat{e}}_{\theta }{\partial }_{\theta }+(1/\sin \theta ){\widehat{e}}_{\phi }{\partial }_{\phi }$. With the overall minus sign in H, we follow the convention that a sphere has positive curvature. The local density of the active component is denoted by $\tilde{\psi }(\theta ,\phi )$ and the projected local density is obtained as $\psi (\theta ,\phi )=\sqrt{g}\tilde{\psi }(\theta ,\phi )$. We consider small fluctuations (ϕ ≪ 1) of this projected density about its steady state value ψ_st,

$\psi ={\psi }_{\text{st}}+\varphi .(11)$

An effective energy functional can be written for the vesicle based on the Helfrich Hamiltonian as [26],

$\mathcal{H}=\int \mathrm{d}A\sqrt{g}\left\{\sigma +\frac{\kappa }{2}{\left(2H-C_p\tilde{\psi }\right)}^2+\frac{\chi }{2}{\tilde{\psi }}^2\right\}(12)$

where κ is the bending rigidity of the membrane and C_p is the preferred curvature of the active component. The integral is defined over the entire area of the reference sphere, where $dA=r_0^2\sin \theta d\theta d\phi$. The equilibrium dynamical equation of this system is given by,

${\partial }_{t}u=-\mathrm{\Gamma }\frac{\delta \mathcal{H}}{\delta u}+{\xi }_u\left(t\right),(13)$

${\partial }_t\varphi =M{\nabla }_{\perp }^2\frac{\delta \mathcal{H}}{\delta \varphi }+{\xi }_{\varphi }\left(t\right),(14)$

where Γ is the membrane mobility and M is the lateral mobility of the active species. The first term on the right hand side of Eq. 13 and Eq. 14 represent the relaxation of the local radius and the density field respectively. The second term represents the equilibrium random noise. In this paper we only look at the regime where random forces from the surroundings on the system are small compared to the active noise (to be introduced below) and set ξ_u(t) = 0 and ξ_ϕ(t) = 0 for the rest of the analysis.

The active processes shown schematically in Figure 1C are modeled as follows. Fusion can happen anywhere on the surface regardless of whether active species are present in the region or not. In a very general case, the rates of fusion may also depend on the local curvature and composition of the active species. Thus, we assume the following form for the rate of fusion

$r^{+}={\lambda }_1+{\lambda }_{2}H+{\lambda }_3\tilde{\psi }+{\lambda }_4\tilde{\psi }H(15)$

where λ₁ represents a uniform rate of fusion and λ₂, λ₃, and λ₄ are appropriate coupling constants. In our model, fission is allowed only from regions where active species are present on the surface, therefore it must always be coupled to $\tilde{\psi }$. The rate of fission is taken as,

$r^{-}={\mu }_1\tilde{\psi }+{\mu }_2\tilde{\psi }H(16)$

where μ₁ and μ₂ are coupling constants. Note that, in the special case of spatially uniform fission and fusion, all the coupling parameters except λ₁ and μ₁ are zero. For a perfectly spherical (u = 0) and homogeneous (ϕ = 0) vesicle, the steady state is attained when the rates are balanced, ${\left(r^{+}-r^{-}\right)}_{u=0,\varphi =0}=0$. From this condition, the steady state concentration can be written in terms of the coupling parameters as,

$\begin{array}{ll}\hfill {\psi }_{\text{st}}& =\frac{-\left({\lambda }_1+\frac{{\lambda }_2}{r_0}\right)}{\left({\lambda }_3+\frac{{\lambda }_4}{r_0}\right)-\left({\mu }_1+\frac{{\mu }_2}{r_0}\right)}\hfill \\ \hfill & =-\frac{{\lambda }_1+\frac{{\lambda }_2}{r_0}}{{\delta }_1+\frac{{\delta }_2}{r_0}},\hfill \end{array}(17)$

where δ₁ = λ₃ − μ₁ and δ₂ = λ₄ − μ₂. We now have four parameters related to the activity—1) λ₁: the uniform rate of fusion, 2) λ₂: the local curvature dependent rate of fusion, 3) δ₁: the activity rate that depends on the local density of active species, and 4) δ₂: the activity rate that depends both on the local density and the local curvature. They are not all independently controllable as the steady-state concentration imposes a constraint. Thus, we can write the parameter δ₂ in terms of the other three as,

${\delta }_2=-\left[r_0{\delta }_1+\frac{1}{{\psi }_{\text{st}}}\left(r_0{\lambda }_1+{\lambda }_2\right)\right].(18)$

For a configuration that deviates from the reference state by u and ϕ respectively in the radius and local density fields, the net change due to activity is,

${\left(r^{+}-r^{-}\right)}_{u,\varphi }=-\frac{1}{{\psi }_{\text{st}}}\left({\lambda }_1+\frac{{\lambda }_2}{r_0}\right)\varphi +2\left({\lambda }_1+\frac{{\lambda }_2}{r_0}\right)u+\frac{1}{2}\left({\lambda }_1+{\psi }_{\text{st}}{\delta }_1\right)\left(2u+{\nabla }_{\perp }^{2}u\right).(19)$

The dynamical equations for the local radius and the local density fields of the active vesicle are,

${\partial }_{t}u=-\mathrm{\Gamma }\frac{\delta \mathcal{H}}{\delta u}+\frac{l}{\tau }\left(r^{+}-r^{-}\right)(20)$

${\partial }_t\varphi =M{\nabla }_{\perp }^2\frac{\delta \mathcal{H}}{\delta \varphi }+\frac{a}{\tau }\left(r^{+}-r^{-}\right),(21)$

where the quantities l and a represent the temporal variation of local radius and local density field, respectively, due to activity and τ represents the activity timescale. In both Eqs 20, 21, the first term of the right hand side captures the equilibrium relaxation and the second term corresponds to the nonequilibrium dynamics. So as to carry out a linear stability analysis, the variables u and ϕ are expanded in the spherical harmonic basis as

$u\left(\theta ,\phi ,t\right)=\sum _{n=2}^{\infty }\sum _{m=-n}^n{u}_{nm}\left(t\right)Y_{nm}\left(\theta ,\phi \right)(22)$

$\varphi \left(\theta ,\phi ,t\right)=\sum _{n=1}^{\infty }\sum _{m=-n}^n{\varphi }_{nm}\left(t\right)Y_{nm}\left(\theta ,\phi \right),(23)$

where Y_nm (θ, φ) represents a spherical harmonic mode. The modes corresponding to n = 0 and n = 1 with respect to variable u capture the radial expansion and the lateral translation modes respectively. As these two modes do not change the energy of the system, they are neglected in the spherical harmonic expansion. From Eq. 20 and Eq. 21, we get the linearized coupled dynamical equation for the spherical harmonic coefficients u_nm and ϕ_nm as,

$\begin{array}{ll}\hfill {\partial }_t\left(\begin{array}{c}\hfill u_{nm}\hfill \\ \hfill {\varphi }_{nm}\hfill \end{array}\right)& =\frac{1}{{\tau }_{\kappa }}\left(\begin{array}{cc}\hfill -\bar{\sigma }\left[n\left(n+1\right)+2\right]-n\left(n+1\right)\left(n-1\right)\left(n+2\right)-4\overline{C_p}{\psi }_{\text{st}}\left[n\left(n+1\right)-1\right]\hfill & \hfill \overline{C_p}\left[n\left(n+1\right)-2\right]+2{\psi }_{\text{st}}\left({\overline{C_p}}^2+\bar{\chi }\right)\hfill \\ \hfill +\frac{1}{2}{\psi }_{\text{st}}^2\left({\overline{C_p}}^2+\bar{\chi }\right)\left[n\left(n+1\right)-6\right]-\frac{\bar{l}}{2}\left({\lambda }_1+{\psi }_{\text{st}}{\delta }_1\right)\left[n\left(n+1\right)-2\right]+2\bar{l}\left({\lambda }_1+{\bar{\lambda }}_2\right)\hfill & \hfill -\frac{\bar{l}}{{\psi }_{\text{st}}}\left({\lambda }_1+{\bar{\lambda }}_2\right)\hfill \\ \hfill \bar{M}\overline{C_p}n\left(n+1\right)\left(n-1\right)\left(n+2\right)+2\bar{M}{\psi }_{\text{st}}\left({\overline{C_p}}^2+\bar{\chi }\right)n\left(n+1\right)\hfill & \hfill -\bar{M}\left({\overline{C_p}}^2+\bar{\chi }\right)n\left(n+1\right)\hfill \\ \hfill -\frac{\bar{a}}{2}\left({\lambda }_1+{\psi }_{\text{st}}{\delta }_1\right)\left[n\left(n+1\right)-2\right]+2\bar{a}\left({\lambda }_1+{\bar{\lambda }}_2\right)\hfill & \hfill -\frac{\bar{a}}{{\psi }_{\text{st}}}\left({\lambda }_1+{\bar{\lambda }}_2\right)\hfill \end{array}\right)\hfill \\ \hfill & \left(\begin{array}{c}\hfill u_{nm}\hfill \\ \hfill {\varphi }_{nm}\hfill \end{array}\right),\hfill \end{array}(24)$

where, 1/τ_κ = Γκ. The dimensionless quantities in Eq. 24 are defined as given in Table 1.

TABLE 1

TABLE 1. Parameters considered in the model, their dimensionless form, and the (range of) scaled values used for the stability analysis.

In the absence of activity [τ_κ/τ = 0], the dynamical equations represent equilibrium relaxation of the membrane. In addition to this, when C_p = 0, the dynamical matrix becomes diagonal and the eigenvalues are ∼ − n⁴ and $-\bar{M}\bar{\chi }{n}^2$. These are respectively the bending relaxation and diffusion modes of a membrane in equilibrium [27]. A non-zero preferred curvature (and a non-zero protein density) will result in coupling of the two relaxation modes.

4 Results

4.1 Activity-induced steady state morphologies

We use the modified DTMC method (described in 2) to explore and characterize the steady state morphologies of a two-component vesicle with a minority component that is curvature inducing. Only the curvature inducing components on the membrane undergo active recycling, when vertices labelled with these components are added on the outer side or removed from the vesicle. The simulation ensures that the mean and the width of fluctuation in the number of these active species is set to specific values. A nearly spherical vesicle is used as the initial configuration in the simulation runs. The reduced volume, defined as the ratio of the vesicle volume to the volume of a sphere having the same area, is used to broadly characterize the shapes. A value of reduced volume close to unity indicates that the vesicle shape is nearly spherical, whereas values smaller than unity indicates that the vesicle is deformed. The organization of active sites on the surface is quantified with the help of cluster size probability distributions. In all the simulations, the spherical topology of the membrane is maintained, the vesicle volume is allowed to change, and the net pressure difference between the inside and outside is set to be zero. We choose the bending modulus and spontaneous curvature values such that in the absence of activity (no vertex addition or vertex removal MC moves) the vesicle remains spherical.

For the results presented in this subsection, we assume that there is no direct interaction between the active components (i.e., J = 0 in Eq. 5). The dependence of the average reduced volume of steady state vesicle configurations on activity rate is shown in Figure 2A. When the preferred curvature of the active species is zero (C_p = 0.0), the reduced volume does not deviate significantly from unity even at higher activity rates suggesting that the vesicle remains spherical (see black curve in Figure 2A). For non-zero spontaneous curvatures, Figure 2A shows a sharp transition in the average reduced volume with increasing activity rate. Such a change in the average reduced volume indicates that there is a drastic change in the vesicle morphology. When the preferred curvature of the active species is positive, we observe a transition from spherical shape to flattened sacs and stomatocytes. On the other hand, when the preferred curvature is negative, we see tubular shapes at lower activity which transitions sharply to spherical shape at higher activity. The steady state shapes obtained at various values for the activity rate and spontaneous curvatures are shown in Figure 2C. A morphology phase diagram as a function of activity rate and preferred curvature of the active components is shown in Figure 2B.

The morphology changes can be understood vis a vis three competing timescales involved in the process; 1) curvature relaxation, 2) membrane recycling or activity and 3) in-plane diffusion. In our simulations, we observed that it requires $\sim 500$ MCS for the curvature, added through vertex addition, to relax. Membrane recycling/activity timescale τ is defined as the average time between successive active events (fusion or fission). The membrane recycling timescale is set by the activity rate in our simulations. We find that the two time scales are comparable when the activity rate is 0.2/MCS. The activity rate is referred to as low/high in comparison with this rate. Since the addition of material is always from outside the added curvature is always positive irrespective of the value of C_p. For the case of C_p > 0.0, since the preferred curvature is closer to the curvature of the newly added vertex, the curvature active component takes longer to diffuse out from the region, i.e., the in-plane diffusion is slow in this case. This implies that the probability of material removal from the same region where it was added is higher. At low activity rates, the membrane gets enough time to relax the added curvature to the C_p value. Whereas at high activity rates, i.e., when the curvature relaxation timescale is much larger than the activity timescale, the curvature fluctuation due to addition/removal is high and leads to spherical shape of the vesicle being destabilized as the result of the dynamic tension [16]. Since the membrane-recycling time dominates this regime, the added vertex do not get sufficient time to sense the preferred curvature. This explains the low-sensitivity to the magnitude of C_p (see red curves in Figure 2A). In the case of negative spontaneous curvature C_p < 0.0, due to a mismatch between the added curvature, which is always positive due to the addition of vertices only on the outer side of the vesicle, and the preferred curvature, the active species quickly diffuses out of the addition site (i.e., fast in-plane diffusion) and eventually gets removed from a different site. At low activity rates, the membrane has time to relax to the preferred curvature of the active species. Hence, what gets removed is the negatively curved region, while what is added is a positively curved region. This biased curvature fluctuation drives an instability leading to tubular shapes as seen in Figure 2C. At high activity rates the membrane recycling dominates over the curvature relaxation. What makes C_p < 0 different from that of C_p > 0 at this high activity rate is that in the former the vertex removed will have the average curvature value of the vesicle while in the latter it is most likely that the added vertex with a higher curvature itself gets removed. This disparity arises due to two factors 1) the curvature added is always more than the average curvature of the vesicle and 2) the active species moves faster out of the regions with unfavorable curvature. Thus, in-plane diffusion of the active species is a deciding factor for the stability of the spherical shape. In order to cross check the role played by in-plane diffusion, in deciding the steady state shapes, we ran simulations without in-plane diffusion for the case of C_p = −0.1. Since there is no in-plane diffusion of active species, the site of vertex addition and removal remains the same and we see a transition from spherical to flattened shapes similar to the case of C_p > 0.0 (see blue filled stars in Figure 2A). Thus, we see that an interplay between the in-plane diffusion, the activity rate, and the added curvature offers a mechanism to spontaneously generate and stabilize various morphologies of the vesicle.

The transition from low activity to high activity is independent of the value of C_p and is set by the typical shape relaxation time. The shape relaxation in our simulations is $\sim 500$ MCS and at steady state, we have $\sim 120$ active species on the membrane surface. In order for each one of the active species to have a chance for removal within the curvature relaxation time, the rate of activity should be greater than 0.2/MCS. Since the diffusion time scale is larger, this correlates with the region of transition shown in Figure 2 for C_p > 0. On the other hand, when C_p < 0 diffusion dominates and the competition is between diffusion and activity rate which we verify explicitly in our simulation as described above.

4.2 Active recycling inhibits formation of large clusters

When the interaction strength J, between the active species, is greater than zero, it is energetically favourable for the active species to form clusters. In this subsection, we focus on the effect of activity rate and the spontaneous curvature on the size of such clusters. Figure 3A shows, a log-log plot between the probability of observing a cluster of a particular size versus the cluster size. The dependence on activity rate and the preferred curvature of the active species C_p are shown separately. The shoulder like regions in both the plots show that there is a preference for clusters having around 10 active species for low activity and low C_p values.

FIGURE 3

FIGURE 3. (A) The dependence of the cluster size distribution on activity rate. Increasing the activity rate suppresses the formation of large clusters. The spontaneous curvature for all the curves is C_p = 0.6. (B) The dependence of the cluster size distribution on the spontaneous curvature. Increasing the spontaneous curvature suppresses the formation of large clusters. The activity rate for all the curves is 0.06 per MCS. The parameter values used in the simulation are $\bar{N}=1200$, ζ = 0.1, κ = 10k_BT, and J = 1.0.

Figure 3A compares the cluster size distributions at different activity rates when the spontaneous curvature of the active species is fixed at 0.6. We see that as activity rate is increased, the probability of larger clusters diminishes. Similarly, Figure 3B shows the cluster size distributions at a low activity rate (0.06 per MCS) for different spontaneous curvatures. In this case too, we see that the probability of large cluster sizes diminishes with increase in spontaneous curvature. Thus, we conclude that both activity rate and its coupling to the curvature of the membrane inhibit the formation of clusters on the vesicle surface.

4.3 Stability of the spherical shape

We have carried out a linear shape stability analysis of a two-component vesicle subjected to active membrane recycling. The derivation of the linearized dynamical equation about the homogeneous spherical reference state is given in Section 3. The spherical shape of the vesicle remains stable at all C_p values in the absence of activity when the steady state density of the active component ψ_st is fixed at 0.1. The stability phase diagram of the vesicle in the (C_pr₀, ψ_st) plane, in the absence of activity, is shown in Figure 4A. Here r₀ is the radius of a reference sphere. The dark regions in Figure 4A corresponds to the case where curvature-composition coupling is strong enough to destabilize the spherical shape [28]. Values for other dimensionless parameters are chosen as shown in Table 1.

FIGURE 4

FIGURE 4. (A) Stability phase diagram for the vesicle without activity. Spherical morphology is unstable in the high C_p and high ψ_st regime (dark blue). The dashed line marks the case of ψ_st = 0.1, where the spherical shape of the vesicle is stable for the range of spontaneous curvature values considered here. (B,C) The stability diagram in the presence of activity in ${\bar{C}}_p$–γ plane for two values of ${\bar{\lambda }}_2$. (D–F) The stability diagram in the presence of activity in ${\bar{\lambda }}_2$–γ plane for two values of ${\bar{C}}_p$. The values of all other parameters are as shown in Table 1.

In the presence of activity, the stability is also governed by the coupling parameters, namely the local curvature dependent rate of fusion $({\bar{\lambda }}_2)$ and the local projected density dependent activity rate (δ₁). The activity rate is controlled by changing the ratio of timescale for equilibrium membrane relaxation and typical activity timescale, defined as γ = τ_κ/τ. For the analysis, we fix the coefficient representing the bare rate of fusion λ₁ as unity. When both the coupling parameters $({\bar{\lambda }}_2,{\delta }_1)$ are negative (Figure 4B), we observe that the spherical shape is stable for all negative spontaneous curvatures and unstable for positive spontaneous curvatures above a certain activity rate. On the other hand, when the curvature coupling parameter ${\bar{\lambda }}_2$ is positive, i.e., more fusion in positively curved regions, we see that the spherical shape is unstable in the negative half-plane of ${\bar{C}}_p$ and stable up to a certain activity rate for positive ${\bar{C}}_p$ (Figure 4C). Figure 4 (bottom panel) shows the ${\bar{\lambda }}_2$–γ plane of stability phase diagram for three values of ${\bar{C}}_p$ while the δ₁ is fixed. The stability phase diagrams for δ₁ = 5 (not reported in Figure 4) were qualitatively similar to those for δ₁ = −5.

The linear stability analysis predicts the conditions under which the spherical shape is unstable in a certain parameter space. As a general trend, we observe that, when activity is large and non-uniform $({\bar{\lambda }}_2\ne 0)$, the spherical shape is unstable. On the other hand, when the activity is independent of local curvature (i.e., ${\bar{\lambda }}_2=0$) the spherical shape is seen to be stable at all activity rates when ${\bar{C}}_p\ge 0$. A more direct mapping of the coupling parameters to those in simulations is complicated and is not pursued here. Moreover, it is difficult to predict using such a stability analysis, whether the unstable vesicle achieves a non-spherical steady state shape and if yes, discern what those shapes are.

5 Discussion

Active transport of membrane materials is considered essential for maintaining the composition of intracellular membranes, particularly for those in the endocytic pathway. Such transport is also expected to influence the large scale morphology of organelles. We have developed a nonequilibrium simulation method based on DTMC to study the effect of fluctuations in area and spontaneous curvature on the steady state morphology of a vesicle. The spherical shape of the vesicle was found to be stable at all activity rates as long as the spontaneous curvature of the active component was zero. At positive spontaneous curvatures, the vesicle is spherical at low activity rates ($<10^{-1}$/MCS). At very high activity rates ($>1$/MCS), the vesicle assumes a stomatocyte shape (Figure 2A). For intermediate activity rates between 0.1–1.0/MCS, the vesicle assumes a tubular or a flattened disc morphology. On the other hand, for C_p < 0.0 we see almost the opposite trend, with tubular shapes for low activity and spherical vesicle for high activity. Thus the steady state vesicle morphology was found to depend on the spontaneous curvature of the active species, the interaction between active species, and the activity rate. A qualitative morphology phase diagram is presented in Figure 2B, which shows the parameter range in which the three classes of steady state shapes are observed. The results presented here is qualitatively different from that shown earlier for membranes with only active curvature fluctuation [16], where the transitions were found to be insensitive to the sign of the spontaneous curvature of the active species. The primary difference between the work reported here and that in Ref. [16] is that we allow for local membrane area fluctuations along with active curvature fluctuations. The fact that the fusion considered here is in/out asymmetric could be the cause for this sign dependence of shape transition as active species with positive spontaneous curvature are more likely to stay in the region of fusion, while active species with negative C_p is more likely to diffuse out of the region of fusion.

The distribution on the membrane, of active proteins that can cause membrane curvature, will have serious implications on the morphology of the vesicle. We expect a non-zero attractive interaction between the active species to result in formation of large clusters on the surface. While at low activity rates we do find the formation of such large clusters, higher activity rates and high spontaneous curvatures were seen to inhibit the formation of large clusters. Our results show that regulation of the membrane recycling rate could offer a robust mechanism to maintain the vesicle in spherical, tubular, or a stomatocyte morphology, without having large clusters of curvature active components.

In order to gain a little more understanding of the factors that drive the membrane shape instability, we derived the dynamical equations for the local shape and composition fields of a nearly spherical vesicle that is subjected to active material transport. The coupled dynamical equations took into account the bending relaxation of the membrane, the diffusion of the components, and changes in the local radius and local composition due to the active processes. The vesicle was composed of passive and active components, with the passive species having no preferred curvature while the active species possessed a spontaneous curvature and took part in membrane fission and fusion. The conditions under which the active recycling of membrane material can destabilize a spherical vesicle was identified using a linear stability analysis.

Though the biological membranes, with its chemical diversity, is far more complex than the simplified models used in this work, we are able to see that membrane area and curvature fluctuations, resulting from active transport of materials, can lead to non-trivial vesicle shapes. Future work will explore the effect of multiple active agents that undergo reaction and diffusion on the vesicle surface.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

PBSK and TVSK conceived the presented idea. PBSK developed the numerical code and performed the computations. PBSK supervised the findings of this work. All authors discussed the results and contributed to the final manuscript.

Acknowledgments

TS thanks IIT Palakkad for the hospitality and for providing high performance computing resources. The authors thank Madan Rao and John Ipsen for fruitful discussions.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P. Molecular biology of the cell. 2 edn. New York: Garland Science (2008).

Google Scholar

2. Phillips R, Kondev J, Theriot J, Garcia H. Physical biology of the cell. New York, NY: Garland Science (2012).

Google Scholar

3. Heald R, Cohen-Fix O. Morphology and function of membrane-bound organelles. Curr Opin Cell BiolCell architecture (2014) 26:79–86. doi:10.1016/j.ceb.2013.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Shibata Y, Voeltz GK, Rapoport TA. Rough sheets and smooth tubules. Cell (2006) 126:435–9. doi:10.1016/j.cell.2006.07.019

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Guven J, Huber G, Valencia DM. Terasaki spiral ramps in the rough endoplasmic reticulum. Phys Rev Lett (2014) 113:188101. doi:10.1103/PhysRevLett.113.188101

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Voeltz GK, Prinz WA. Sheets, ribbons and tubules—How organelles get their shape. Nat Rev Mol Cell Biol (2007) 8:258–64. doi:10.1038/nrm2119

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Shibata Y, Hu J, Kozlov MM, Rapoport TA. Mechanisms shaping the membranes of cellular organelles. Annu Rev Cell Dev Biol (2009) 25:329–54. doi:10.1146/annurev.cellbio.042308.113324

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Misteli T. The concept of self-organization in cellular architecture. J Cell Biol (2001) 155:181–6. doi:10.1083/jcb.200108110

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Karsenti E. Self-organization in cell biology: A brief history. Nat Rev Mol Cell Biol (2008) 9:255–62. doi:10.1038/nrm2357

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Rafelski SM, Marshall WF. Building the cell: Design principles of cellular architecture. Nat Rev Mol Cell Biol (2008) 9:593–602. doi:10.1038/nrm2460

PubMed Abstract | CrossRef Full Text | Google Scholar

11. McMahon HT, Gallop JL. Membrane curvature and mechanisms of dynamic cell membrane remodelling. Nature (2005) 438:590–6. doi:10.1038/nature04396

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zimmerberg J, Kozlov MM. How proteins produce cellular membrane curvature. Nat Rev Mol Cell Biol (2006) 7:9–19. doi:10.1038/nrm1784

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Karbowski M, Youle R. Dynamics of mitochondrial morphology in healthy cells and during apoptosis. Cell Death Differ (2003) 10:870–80. doi:10.1038/sj.cdd.4401260

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Chen YA, Scheller RH. Snare-mediated membrane fusion. Nat Rev Mol Cell Biol (2001) 2:98–106. doi:10.1038/35052017

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Rao MRCS, R. C. S. Active fusion and fission processes on a fluid membrane. Phys Rev Lett (2001) 87:128101. doi:10.1103/PhysRevLett.87.128101

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Ramakrishnan N, Ipsen JH, Rao M, Kumar PBS. Organelle morphogenesis by active membrane remodeling. Soft Matter (2015) 11:2387–93. doi:10.1039/C4SM02311K

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Goutaland Q, van Wijland F, Fournier JB, Noguchi H. Binding of thermalized and active membrane curvature-inducing proteins. Soft Matter (2021) 17:5560–73. doi:10.1039/D1SM00027F

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Natesan R, Gowrishankar K, Kuttippurathu L, Kumar PBS, Rao M. Active remodeling of chromatin and implications for in vivo folding. J Phys Chem B (2022) 126:100–9. doi:10.1021/acs.jpcb.1c08655

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Tachikawa M, Mochizuki A. Golgi apparatus self-organizes into the characteristic shape via postmitotic reassembly dynamics. Proc Natl Acad Sci U S A (2017) 114:5177–82. doi:10.1073/pnas.1619264114

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Rautu SA, Rowlands G, Turner MS. Size-dependent recycling of membrane clusters. EPL (Europhysics Letters) (2018) 121:58004. doi:10.1209/0295-5075/121/58004

CrossRef Full Text | Google Scholar

21. Mietke A, Jülicher F, Sbalzarini IF. Self-organized shape dynamics of active surfaces. Proc Natl Acad Sci U S A (2019) 116:29–34. doi:10.1073/pnas.1810896115

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Ruiz-Herrero T, Fai TG, Mahadevan L. Dynamics of growth and form in prebiotic vesicles. Phys Rev Lett (2019) 123:038102. doi:10.1103/PhysRevLett.123.038102

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Nelson D, Piran T, Weinberg S. Statistical mechanics of membranes and surfaces. 2 edn. Singapore: World Scientific Singapore (2004).

Google Scholar

24. Ramakrishnan N, Kumar PS, Ipsen JH. Membrane-mediated aggregation of curvature-inducing nematogens and membrane tubulation. Biophysical J (2013) 104:1018–28. doi:10.1016/j.bpj.2012.12.045

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Noguchi H, Gompper G. Meshless membrane model based on the moving least-squares method. Phys Rev E (2006) 73:021903. doi:10.1103/PhysRevE.73.021903

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Helfrich W. Elastic properties of lipid bilayers: Theory and possible experiments. Z für Naturforschung C (1973) 28:693–703. doi:10.1515/znc-1973-11-1209

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Seifert U, Langer SA. Viscous modes of fluid bilayer membranes. Europhys Lett (1993) 23:71–6. doi:10.1209/0295-5075/23/1/012

CrossRef Full Text | Google Scholar

28. Leibler S. Curvature instability in membranes. J Phys France (1986) 47:507–16. doi:10.1051/jphys:01986004703050700

CrossRef Full Text | Google Scholar

29. Gauthier NC, Masters TA, Sheetz MP. Mechanical feedback between membrane tension and dynamics. Trends Cell Biology (2012) 22:527–35. doi:10.1016/j.tcb.2012.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar