Pattern formation in reaction-diffusion system on membrane with mechanochemical feedback.

Shapes of biological membranes are dynamically regulated in living cells. Although membrane shape deformation by proteins at thermal equilibrium has been extensively studied, nonequilibrium dynamics have been much less explored. Recently, chemical reaction propagation has been experimentally observed in plasma membranes. Thus, it is important to understand how the reaction-diffusion dynamics are modified on deformable curved membranes. Here, we investigated nonequilibrium pattern formation on vesicles induced by mechanochemical feedback between membrane deformation and chemical reactions, using dynamically triangulated membrane simulations combined with the Brusselator model. We found that membrane deformation changes stable patterns relative to those that occur on a non-deformable curved surface, as determined by linear stability analysis. We further found that budding and multi-spindle shapes are induced by Turing patterns, and we also observed the transition from oscillation patterns to stable spot patterns. Our results demonstrate the importance of mechanochemical feedback in pattern formation on deforming membranes.

Introduction

Membrane deformation is a fundamental biological process involved in many cellular functions such as vesicular transport[1], cell division[2], and cell motility[3]. To understand these phenomena, the mechanism of membrane deformation by intracellular proteins has been investigated in detail[4-8]. Recently, it has been shown that the deformation of biological membranes is not just a passive phenomenon but also plays physiological roles[8-15]. For example, membrane curvature induces localization of membrane proteins in highly curved domains[9] and phase separation of lipid membranes[10-12]. This clustering can lead to the emergence of lipid rafts, which are believed to play important roles in cell signaling and membrane trafficking[12,13,16]. Membrane binding by curvature-inducing proteins that are involved in vesicular transport is also regulated by membrane curvature and by various proteins[7,17,18]. For example, recruitment of curvature-inducing protein FBP17, involved in endocytosis, onto the membrane is regulated by the local membrane curvature, membrane tension, and endocytic proteins[17-19]. This mechanism is suggested to play important roles in cell polarization[19], endocytosis[20], and cell division[21].

To understand pattern formation on curved surfaces, several types of studies have been conducted[22-32]. One typical approach is to analyze pattern formation at thermal equilibrium based on phase separation[22-26]. This type of study has shown that membrane shapes and domain patterns of equilibrium states are affected by the line tension of domain boundaries, bending rigidity, and local curvatures[22-26]. Such studies have successfully described the experimentally observed patterns of multi-component lipid vesicles. However, studies pertaining to kinetics are limited to the dynamics of relaxation toward an equilibrium state[22-26,31,33].

Most of previously conducted theoretical and numerical studies have examined only the effects of protein binding; however, in living cells, it is known that many proteins typically work in concert to regulate biological functions. Propagation waves in membranes are often observed during cell migration, spreading, growth, or division[34-41]. Such waves and chemical patterns can be reproduced through activator-inhibitor systems of reaction–diffusion models[42]. The reaction–diffusion system was first proposed by Turing to describe the symmetry breaking of morphogenesis[43], and has been applied to curved surfaces such as animal skins and tissues[44-46]. These studies have shown that geometry affects pattern formation and domain localization[29]; however, the conclusions of such studies are limited by the fact that the surface shape is fixed, although the effects of size increase have been investigated[27,28]. Recently, the propagating waves of F-BAR protein and actin growth have been explained by the reaction–diffusion systems of five chemical reactants on a quasi-flat membrane[18]. As large membrane deformations caused by the coupling of curvature and reaction–diffusion systems have not yet been studied[41], the effects of membrane deformation on reaction–diffusion systems have not been elucidated.

In this study, we investigated the coupling effects between membrane deformation and reaction–diffusion systems by simulating vesicle deformation through curvature-inducing proteins and also chemical reactions using a reaction–diffusion model. Our model accounts for the mechanochemical feedback between membrane curvature and protein concentration. We employed a dynamically triangulated surface model to represent the membrane and calculated the curvature energy to solve the membrane deformation dynamics[47-49]. We employed the Brusselator model[50], one of the simplest reaction–diffusion systems, modifying it to include the mechanochemical feedback from membrane curvature. As the dynamics of a non-deformable surface are well understood, we were able to analyze the evident membrane-deformation effects. We describe how this coupling changes the vesicle shape and pattern formation.

Results

Reaction–diffusion model and stability analysis

A two-dimensional reaction–diffusion system with two reactants is written as where is a time constant, and are diffusion coefficients of reactants and , and is a two-dimensional Laplace–Beltrami operator. In this study, we consider the Brusselator model, which is described by the reaction scheme below:

The reaction equations are given by and where and are positive parameters[50].

In the coupling of the reaction–diffusion system with the change in membrane curvature, represents the local area fraction covered by curvature-inducing binding proteins on the membrane and is the concentration of a protein to regulate the protein binding. The free energy in relation to curvature is expressed as withwhere and represent the bending rigidity without or with the bound proteins, respectively; is the spontaneous curvature; is the surface area; and is the mean curvature, where and are two principal curvatures. The corresponding curvature term is added to the reaction equation thus the reaction–diffusion equations are written aswhere is the mechanochemical feedback magnitude of the reaction (), and is a normalization factor expressed as used to obtain Turing and oscillation phases at To maintain is restricted between the lower and upper bounds: it is set to or when the time evolution of Eq. (2) crosses those bounds. The first reaction becomes which can be considered to represent the binding of protein from the solution surrounding the membrane. Thus, the binding of is enhanced at a membrane curvature where so that On the other hand, the time evolution of is not directly dependent on the local membrane curvature. Note that the mixing-entropy term of the protein concentration is not accounted to reproduce the normal Brusselator dynamics when the membrane shape is fixed. In this study, we use and for all simulations.

Based on the linear stability analysis around the fixed point, [51], the conditions for Hopf and Turing bifurcations with a membrane curvature effect A′ on a fixed spherical surface are and respectively; and temporal oscillations and spatial patterns appear above these. The membrane curvatures for Hopf and Turing bifurcations are given below, respectively:where . At , a homogeneous phase is formed, because . The phase stability diagram is shown in Fig. 1. This diagram shows that bifurcations occur as the magnitude of the spontaneous curvature and mechanochemical coupling magnitude increase at .

Figure 1

The phase diagram for the Brusselator, modified to include a membrane curvature effect, on a surface of a constant mean curvature at , , and . The purple and green lines are the Turing bifurcation curve and the Hopf bifurcation, respectively. These curves separate the regions in which the homogeneous stable patterns (H), stationary Turing patterns (T), or temporal oscillation patterns (O) occur. The red line indicates .

Pattern formation on membrane

The membrane motion is solved by the Langevin dynamics of dynamically triangulated surface model, which formed a triangular network of spherical topology with vertices, as described previously[47]. In this study, the presence of curvature-inducing proteins is considered in addition to the model as given in Eq. (1). We use and where is the thermal energy (see “Methods” for more details). The results are displayed with the length unit , energy unit , and time unit .

First, we analyzed the pattern formation on the fixed surface of a spherical vesicle at the reduced volume, , where is the vesicle volume (Fig. 2a,b,g). The results are consistent with those of the linear stability analysis (Fig. 2g). The effects of thermal fluctuations are discussed in the “Supplementary Material”. Figure 2a,b show typical snapshots. One large circular Turing domain appears at and (Fig. 2b).

Figure 2

(a–f) Snapshots of the vesicles and (g,h) phase diagrams for , , and . (a,b,g) (fixed shape). (c–f,h) . (a,c) and . (b) and . (d) and . (e) and . (f) and . The color in snapshots indicates the concentration of the curvature-inducing protein, . Purple and green lines on the phase diagrams represent the Turing bifurcation curve and Hopf bifurcation, respectively, and the symbols represent the simulation results. The red line indicates . Two or three overlapped symbols indicate the coexistence of multiple patterns.

In contrast, membrane deformation changes the chemical patterns in deformable vesicles at (Fig. 2c–f,h). The oscillation phase is suppressed, and Turing pattern is observed in a wider parameter region (Fig. 2h). At high spontaneous curvature , budding and spicule shapes are formed, accompanied by Turing patterns (Fig. 2d,e). These spicule shapes only appear under conditions of Turing pattern formation, while budding can occur in homogeneous membranes. Moreover, budded spheres typically have a high value of that is homogeneously distributed and form a Turing domain boundary separating two phases with higher or lower value of at the narrow connective neck, as shown in Fig. 2f, because of the reduction in diffusion through the neck. Thus, the Turing pattern is modified by the membrane shape deformation. Bud formation is obtained for at (Fig. 2h). This is reasonable, as the curvature energy of a spherically shaped bud with a radius which is fully covered by the curvature-inducing protein () is minimal. The condition of bud formation is given by , since the volume of the rest of a vesicle of a spherical shape is maximal. In the case of , the threshold is .

For high values of , different shapes can be formed depending on the initial shapes, such as the prolate and budded shapes shown in Fig. 2h. Figure 3 shows another example. Vesicles of three or four spicules are formed from prolate and oblate vesicles, respectively, with (Fig. 3a,b). As pattern formation progresses, the vesicle shape changes according to the chemical pattern (Supplementary Movie S1). In order to evaluate the non-uniformity of and the smoothed local curvature , we calculated separation metrics, and , where and are the between-class variance and within-class variance, respectively[52] (The curvature smoothing method is described in the “Supplementary Material”). Each variance is calculated as below:where is the probability of each class, is the class mean value, and is the class variance. The threshold value to divide into two classes is determined to maximize the metric value. Therefore, the metrics becomes large when Turing patterns are formed clearly, whereas becomes small when the two phases are gently separated or not separated (i.e., homogeneous patterns). Figure 3c,d show that increases as the Turing pattern develops, followed by an increase in ; this sequence is consistent with that depicted by the sequential snapshots and indicates that non-uniformity can be distinguished by calculating the separation metrics. We also calculated the time development of asphericity, , to evaluate vesicle deformation (Fig. 3e). Asphericity is the degree of deviation from a spherical shape, calculated as below:where is the eigenvalue of the gyration tensor of the vesicle[47,53,54]. For a sphere, , and for the thin-rod limit, ( and . As the vesicle forms three or four spindles, decreases (Fig. 3e,f).

Figure 3

Examples of pattern formation and membrane deformation. (a,b) Sequential snapshots of the vesicles for , , , , , and starting from (a) prolate and (b) oblate shapes. The color indicates the concentration of curvature-inducing protein, . (c–e) Time evolution of (c) the separation metric of the protein concentration, , (d) that of the local curvature, , (e) asphericity, , and (f) the number of domains, . The purple and green lines indicate the simulations starting from prolate and oblate shapes, respectively. Results are presented as the mean ± standard error (n = 8).

To further investigate the effect of coupling between the Brusselator and vesicle deformation, we conducted the simulation with different and values at (Fig. 4 and Supplementary Fig. S3). As decreases, the number of domains, and decrease, whereas increases (Fig. 4e–g). In addition, the domain size increases as increases. Therefore, higher and values and a lower are obtained at a lower (Fig. 4h–j). When , convex regions are formed in various directions and the vesicle becomes nearly spherical, but when , the vesicle becomes prolate in shape, and increases (Fig. 4a–d). Thus, chemical pattern formation affects vesicle deformation and the relation between and the preferred curvature of the domains is important in determining the stable shapes. The results do not significantly differ between simulations that start from prolate or oblate shapes, except under the condition at and (Fig. 4a,f). Under that condition, with starting from a prolate-shaped vesicle, two domains arise at the pole of prolate, and the vesicle shape remains in the prolate shape. However, when the simulation starts from the oblate-shaped vesicle, multiple domains arise at the edge of oblate, and vesicle shape morphs into a multi-spindle shape (Fig. 4a,f,g). As well as the effect of chemical pattern formation to the vesicle deformation, vesicle shape can also affect chemical pattern formation.

Figure 4

(a–d) Snapshots of vesicles for , , , , and for two values of and starting from prolate and oblate shapes. (a) and . (b) and . (c) and . (d) and . The color indicates the concentration of curvature-inducing protein, . (e–j) Time evolution of (e,h) the separation metric of the local curvature, , (f,i) asphericity, , and (g,j) the number of domains, . The data for and at are shown in (e–g), and the data for and at are shown in (h–j). The orange and red lines indicate simulations starting from prolate shapes, and the purple and green lines indicate simulations starting from oblate shapes. Results are presented as mean ± standard error (n = 8).

A comparison of Fig. 2g with Fig. 2h shows that the region encompassing Turing patterns is enlarged in the phase diagram at from , as increases. To investigate this change, we performed simulations at and with different and (Fig. 5). For or , Turing patterns occur instead of oscillations, whereas for , an oscillation occurs at , and the oscillating patterns transition to the Turing pattern at and (Fig. 5d,g,j, and Supplementary Movie S2). As shown in Fig. 5e,h,k, the maximum values of the local curvature at and eventually increase over time; this does not occur at . As the local curvature increases, the position on the phase diagram shifts toward the upper left, as shown in Supplementary Fig. S5. Therefore, the transitions from an oscillation pattern to a Turing pattern is induced by a local increase in .

Figure 5

(a–c) Snapshots of the vesicles for , , , , and for three values of and . (a) . (b) . (c) . The color indicates the concentration of curvature-inducing protein, . (d–l) Time development of (d,g,j) the separation metric of the protein concentration, , (e,h,k) the maximum value of the local curvature, , and (f,i,l) the mean area ratio of one domain . The data for , , and are shown in (d–f), (g–i), and (j–l), respectively. The purple, green, and orange lines indicate the simulation data for , and , respectively. Results of one typical simulation run are shown. The results averaged from eight independent simulations are shown in Supplementary Fig. S4.

At , a small domain is generated and stabilized by the local deformation of the vesicle at . In contrast, a large domain is temporarily generated at but is not stabilized, since the stable domain size is much larger than the sphere of preferred curvature thus the vesicle cannot sufficiently deform (Fig. 6 and Supplementary Movies S3 and S4). In addition, the oscillation period for is significantly longer for than for or for (Fig. 5). The oscillation period is calculated from the peak of the Fourier spectrum of for the eight independent runs: , , and at , , and , respectively. This fact and the time evolution of indicate that membrane deformation is suppressed by the volume restriction for (Fig. 5e). In contrast, substantial membrane deformation occurs at the reduced volumes of and , which enables frequent generation of domains. Thus, membrane deformation can change both oscillation period and stability.

Figure 6

Sequential snapshots of the vesicles for , , , , , and for (a) and (b) . The videos are shown in Supplementary Movies S3 and S4.

Discussion

In this study, we have examined the coupling effects between a reaction–diffusion system and membrane deformation by simulating membrane deformation using a dynamically triangulated surface model. We adapted the Brusselator model to include mechanochemical feedback between local membrane curvature and the concentration of curvature-inducing proteins. Based on the linear stability analysis of the reaction–diffusion system with a membrane curvature effect on a fixed spherical surface, we have clarified that bifurcation curves depend on the mechanochemical coupling magnitude and the value of spontaneous curvature of curvature-inducing proteins with respect to the local membrane curvature (Fig. 1). Thus, the stability of both Turing and oscillation dynamics depend on the membrane shape. We have shown that various shapes, such as buds and multi-spindles, depend on , , and the diffusion constant (Figs. 2, 3, 4). In addition, since the domain formation of curvature-inducing proteins is promoted at regions with high local curvature, the initial shape of the vesicles affects the dynamics of pattern formation (Fig. 4a). Therefore, the dynamics of protein pattern formation change the shape of vesicles, while membrane deformation simultaneously affects pattern formation. This feedback loop can drastically alter the chemical reaction patterns from those on non-deformable surfaces (Fig. 2g,h). A dynamic transition from an oscillating pattern to a Turing pattern is induced by membrane deformation (Fig. 5g–i, and Supplementary Movie S2). Such transitions have not been reported in previous studies.

In the context of living cells, many kinds of proteins and other molecules function interdependently on membranes, where the function of one protein is often activated or inhibited by those of others. Membrane deformation brought about by competing forces of protein-induced curvature changes and surface tension changes impelled by actin growth has been studied[4,8,18,19]. By choosing not to consider the dynamics of actin in this study, we demonstrated that various membrane deformations, accompanied by Turing patterns and oscillations, can be produced by one curvature-inducing protein and one or a small number of regulatory proteins without actin interactions.

Here, we analyzed the coupling of a reaction–diffusion system with membrane deformation utilizing the fixed parameters , , and , focusing primarily on Turing patterns, and oscillatory conditions to a lesser extent. The experimental results indicate that observed patterns, which include a feedback loop between curvature-inducing proteins and membrane deformation, are not only stable spot patterns, such as those observed during cell polarization[19], but are also propagating waves[18]. Similarly, the reconstituted Min system in liposomes, which regulates bacterial cell division, has been shown to exhibit propagating wave patterns[38-40]. These patterns, which induce oscillating membrane deformation, are also described by reaction–diffusion systems. The system developed in this paper can also be applied to these patterns observed in living systems, by adjusting the parameters. Other chemical reaction models, such as the Oregonator[55], which was developed to model the Belousov–Zhabotinsky reaction, and the F-BAR-actin model[18], are also easily applied. Thus, the present model system is a powerful tool that can be used to study a wide range of chemical reaction systems that are coupled with membrane deformation.

Methods

Membrane model

Membrane contains vertices connected by bonds of an average length , with volumes and masses, , excluded. The local curvature energy in Eq. (1) is discretized using dual lattices. The surface area and volume of a vesicle are kept constant at about 0.01% accuracy by harmonic constraint potentials. Details of the potentials are described in Ref.[47]. For the coefficients of area and volume constraint potentials, four times greater values are employed than those in Ref.[47]. To produce membrane fluidity, bonds are flipped to the diagonal of two adjacent triangles using the Monte Carlo method. The membrane motion is solved by molecular dynamics (MD) with the Langevin thermostat:where is the friction coefficient, and is Gaussian white noise, which obeys the fluctuation–dissipation theorem. The hydrodynamic interactions are not considered. The time unit in MD is based on diffusion, and is used. To allow membrane deformation followed by concentration changes in , is employed. Equation (9) is numerically integrated by the leapfrog algorithm with time steps .

Discretization of reaction–diffusion equations

We developed a finite volume method to discretize Eqs. (2) and (3). Since the Kelvin–Stokes theorem holds for curved surfaces, it is straightforwardly applicable, as employed on a flat surface. A vertex-centered finite volume approach is applied to the dual lattices used for the calculation of membrane curvature[47]. The time evolution of of the ith vertex is discretized using the following forward difference method:where is the vertex area, is the side length between neighboring vertex cells, and is the distance between the neighboring vertices. The effect of curvature on diffusion is included as the variation of side lengths. Similarly, Eq. (3) is discretized. In this study, is used. The initial concentrations for the simulations are set around the fixed point , with small random perturbations. When or or 1 are taken, respectively, as the initial concentration instead.

Supplementary Information 1.

Supplementary Video S1.

Supplementary Video S2.

Supplementary Video S3.

Supplementary Video S4.