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.

Pattern formation in reaction-diffusion system on membrane with mechanochemical feedback Naoki Tamemoto & Hiroshi Noguchi * 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.
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][5][6][7][8] . Recently, it has been shown that the deformation of biological membranes is not just a passive phenomenon but also plays physiological roles [8][9][10][11][12][13][14][15] . For example, membrane curvature induces localization of membrane proteins in highly curved domains 9 and phase separation of lipid membranes [10][11][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][18][19] . This mechanism is suggested to play important roles in cell polarization 19 , endocytosis 20 , and cell division 21 .
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][35][36][37][38][39][40][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][45][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 Scientific Reports | (2020) 10:19582 | https://doi.org/10.1038/s41598-020-76695-x www.nature.com/scientificreports/ 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][48][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 τ ∂u ∂t = D u �u + f (u, v) and τ ∂v ∂t = D v �v + g(u, v), where τ is a time constant, D u and D v are diffusion coefficients of reactants u and v , 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 f (u, v) = A − (B + 1)u + u 2 v and g(u, v) = Bu − u 2 v, where A and B are positive parameters 50 .
In the coupling of the reaction-diffusion system with the change in membrane curvature, u represents the local area fraction covered by curvature-inducing binding proteins on the membrane (u ∈ [0, 1]), and v is the concentration of a protein to regulate the protein binding. The free energy in relation to curvature is expressed as F cv = f cv dS, with where κ 0 and κ 1 represent the bending rigidity without or with the bound proteins, respectively; C 0 is the spontaneous curvature; S is the surface area; and H is the mean curvature, H = (C 1 + C 2 )/2, where C 1 and C 2 are two principal curvatures. The corresponding curvature term A ′ is added to the reaction equation (u, v); thus the reaction-diffusion equations are written as where G is the mechanochemical feedback magnitude of the reaction ( G ≥ 0 ), and k u is a normalization factor expressed as k u u, used to obtain Turing and oscillation phases at u ∈ [0, 1]. To maintain 0 ≤ u ≤ 1, u is restricted between the lower and upper bounds: it is set to u = 0 or u = 1 when the time evolution of Eq. (2) crosses those bounds. The first reaction becomes A + A ′ → u, which can be considered to represent the binding of protein u from the solution surrounding the membrane. Thus, the binding of u is enhanced at a membrane curvature On the other hand, the time evolution of v 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 A = 4.5, , the conditions for Hopf and Turing bifurcations with a membrane curvature effect A′ on a fixed spherical surface are B > 1 + A + A ′ 2 and B > 1 + A + A ′ η 2 , respectively; and temporal oscillations and spatial patterns appear above these. The membrane curvatures for Hopf and Turing bifurcations are given below, respectively: i.e., 2A < GH 2 E cv , a homogeneous phase is formed, because u s < 0 . The phase stability diagram is shown in Fig. 1. This diagram shows that bifurcations occur as the magnitude of the spontaneous curvature C 0 and mechanochemical coupling magnitude G increase at A + A ′ ≥ 0.
(1) www.nature.com/scientificreports/ 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 N 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 κ 0 /k B T = 20 and κ 1 /k B T = 40, where k B T is the thermal energy (see "Methods" for more details). The results are displayed with the length unit R = √ S/4π , energy unit κ 0 , and time unit τ. First, we analyzed the pattern formation on the fixed surface of a spherical vesicle at the reduced 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 Gκ 0 /R 2 = 0.061 and C 0 R = 8 (Fig. 2b).
In contrast, membrane deformation changes the chemical patterns in deformable vesicles at V * = 0.8 ( 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 C 0 , 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 u that is homogeneously distributed and form a Turing domain boundary separating two phases with higher or lower value of u 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 C 0 R ≥ 3 at V * = 0.8 (Fig. 2h). This is reasonable, as the curvature energy of a spherically shaped bud with a radius r b = 2/C 0 , which is fully covered by the curvature-inducing protein ( u = 1 ) is minimal. The condition of bud formation is given by V * ≤ (r b /R) 3 + 1 − (r b /R) 2 3/2 , since the volume of the rest of a vesicle of a spherical shape is maximal. In the case of V * = 0.8 , the threshold is R/r b ≥ 2.2.
For high values of C 0 , 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 (u, v) ≃ (u s , v s ) (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 u and the smoothed local curvature where ρ i is the probability of each class, µ i is the class mean value, and σ i 2 is the class variance. The threshold value to divide into two classes is determined to maximize the metric value. Therefore, the metrics s u becomes large when Turing patterns are formed clearly, whereas s u becomes small when the two phases are gently separated or not separated (i.e., homogeneous patterns). Figure 3c,d show that s u increases as the Turing pattern develops, followed by an increase in s H ; 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: , and for the thin-rod limit, α = 1 ( 1 = 1 and 2 = 3 = 0) . As the vesicle forms three or four spindles, α decreases (Fig. 3e,f). To further investigate the effect of coupling between the Brusselator and vesicle deformation, we conducted the simulation with different C 0 and D u values at Gκ 0 /R 2 = 0.046 ( Fig. 4 and Supplementary Fig. S3). As C 0 decreases, the number of domains, N d and s H decrease, whereas α increases (Fig. 4e-g). In addition, the domain size increases as D u increases. Therefore, higher N d and s H values and a lower α are obtained at a lower D u (Fig. 4h-j). When N d > 2 , convex regions are formed in various directions and the vesicle becomes nearly spherical, but when N d = 2 , the vesicle becomes prolate in shape, and α increases (Fig. 4a-d). Thus, chemical pattern formation affects vesicle deformation and the relation between N d 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 C 0 R = 7 and D u = 20 (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.
A comparison of Fig. 2g with Fig. 2h shows that the region encompassing Turing patterns is enlarged in the phase diagram at V * = 0.8 from V * = 1 , as G increases. To investigate this change, we performed simulations at Gκ 0 /R 2 = 0.077 and C 0 R = 10 with different V * and D u (Fig. 5). For D u = 10 or 20 , Turing patterns occur instead of oscillations, whereas for D u = 50 , an oscillation occurs at V * = 0.95 , and the oscillating patterns transition Figure 2.  Supplementary Fig. S5. Therefore, the transitions from an oscillation pattern to a Turing pattern is induced by a local increase in H.
At V * = 0.95 , a small domain is generated and stabilized by the local deformation of the vesicle at D u = 10 . In contrast, a large domain is temporarily generated at D u = 50, but is not stabilized, since the stable domain size is much larger than the sphere of preferred curvature C 0 /2; thus the vesicle cannot sufficiently deform ( Fig. 6 and Supplementary Movies S3 and S4). In addition, the oscillation period for D u = 50 is significantly longer for V * = 0.95 than for V * = 0.8 or for V * = 0.65 (Fig. 5). The oscillation period τ os is calculated from the peak of the Fourier spectrum of s u for the eight independent runs: τ os /τ = 100 , 11 , and 8 at V * = 0.95 , 0.8 , and 0.65 , respectively. This fact and the time evolution of ∼ Hmax indicate that membrane deformation is suppressed by the volume restriction for V * = 0.95 (Fig. 5e). In contrast, substantial membrane deformation occurs at the reduced volumes of V * = 0.8 and 0.65 , which enables frequent generation of domains. Thus, membrane deformation can change both oscillation period and stability.

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  www.nature.com/scientificreports/ on the mechanochemical coupling magnitude G and the value of spontaneous curvature of curvature-inducing proteins C 0 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 multispindles, depend on G , C 0 , and the diffusion constant D u (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-  Here, we analyzed the coupling of a reaction-diffusion system with membrane deformation utilizing the fixed parameters A , B , 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 curvatureinducing 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][39][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 N = 4000 vertices connected by bonds of an average length a , with volumes and masses, m , excluded. The local curvature energy f cv in Eq. (1) is discretized using dual lattices. The surface area S = 0.41a 2 (2N − 4) ≃ 3280a 2 and volume V 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 g i (t) is Gaussian white noise, which obeys the fluctuation-dissipation theorem. The hydrodynamic interactions are not considered. The time unit in MD is τ md = ζ a 2 /k B T based on diffusion, and m = ζ τ md is used. To allow membrane deformation followed by concentration changes in u , τ md = 0.1τ is employed. Equation (9) is numerically integrated by the leapfrog algorithm with time steps �t md = 0.001τ md .

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 u of the ith vertex is discretized using the following forward difference method: where S i is the vertex area, l ij is the side length between neighboring vertex cells, and r ij 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, t rd = 0.1 t md is used. The initial concentrations for the simulations are set around the fixed point (u s , v s ) , with small random perturbations. When u s < 0 or u s > 1, u = 0 or u = 1 are taken, respectively, as the initial concentration instead.