Disappearance, division, and route change of excitable reaction-diffusion waves in deformable membranes

Shapes of biomembrane in living cells are regulated by curvature-inducing proteins. However, the effects of membrane deformation on signal transductions such as chemical waves have not been researched adequately. Here, we report that membrane deformation can alter the propagation of excitable reaction-diffusion waves using state-of-the-art simulations. Reaction waves can induce large shape transformations, such as membrane budding and necking, that erase or divide the wave, depending on the curvature generated by the waves, feedback to the wave propagation, and the ratio of the reaction and deformation times. In genus-2 vesicles, wave division occurs at branching points and collided waves disappear together. We demonstrate that the occasional disappearance of the waves can alter the pathway of wave propagation. Our findings suggest that membrane deformation and reaction waves can together regulate signal transductions on biomembranes.


Results
Membrane modeling. A fluid membrane is represented by a triangulated surface as employed in our previous study [21][22][23] . The binding of two types of proteins is considered: curvature-inducing proteins and regulatory proteins (Fig. 1a). The bending energy of the membrane is represented by Canham-Helfrich model 24,25 modified for two-component systems 26 and is expressed as F cv = f cv dS with where u is the concentration of the curvature-inducing proteins ( u = 1 at maximum coverage) and H is the mean curvature of the membrane. The first and second terms in Eq. (1) give the bending energies of the unbound and bound membranes, respectively. The unbound membrane has the bending rigidity κ 0 = 20k B T and zero spontaneous curvature, where k B T is the thermal energy. The bound membrane has the bending rigidity κ 1 and a finite spontaneous curvature C 0 . This value of κ 0 is typical for lipid membranes 27 and κ 1 > κ 0 is required for curvature generation 28 . The membrane properties are independent of the concentration v of regulatory proteins. We incorporated a modified FitzHugh-Nagumo model as RD equations for u and v on the membrane as follows: where τ rd is the RD time unit and G is the magnitude of mechanochemical coupling. That is, the curvatureinducing proteins binds more in the case the binding reduces the bending energy. Unless specified, G = 0.002 is used. The diffusion coefficient D u = 0.1σ 2 , and ∇ 2 is the two-dimensional Laplace-Beltrami operator where σ is the average edge length of the triangular mesh. The FitzHugh-Nagumo model 29,30 is modified to maintain 0 < u < 1 as follows: . The proteins of u and v are the activator and inhibitor, respectively, and an excitation occurs along the nullcline for u as shown in Fig. 1b. In the absence of the excitation, the proteins maintain the constant concentrations (the fixed point depicted in Fig. 1b). Once excited, the binding of the curvature-inducing proteins (u) induces the binding of the regulatory proteins (v) and subsequently, a high concentration of v leads to the unbinding of the curvature-inducing proteins (see the arrows in Fig. 1b).
To initiate a wave, the proteins are bound on a narrow region of the membrane at initial states, as indicated by the snapshot in the top row of Fig. 2b. The regulatory proteins are set at the left side of the curvature-inducing proteins to generate a wave in the right direction. The membrane motion is solved by molecular dynamics (MD) with the Langevin thermostat, that is, thermal conduction is assumed to be sufficiently fast, so that a constant temperature is maintained. The time τ md , in which a free particle diffuses the distance σ , is the MD time unit. A detailed description of the methods is provided in the "Methods" section and Supplementary Information. Waves in straight tubes. First, we describe the wave propagation in a straight membrane tube (Figs. 2, 3, 4, 5, 6). We investigate the conditions, in which the cylindrical shapes of unbound ( u = 0.07 , unexcited) and bound ( u = 0.9 , excited) membranes are stable and unstable, respectively; i.e., the region between two dashed lines in Fig. 3a. These two thresholds are determined by the condition that the homogeneous membrane generates a shrinking force along the tube axis given by 23,31 : www.nature.com/scientificreports/   Figure 3. Phase diagrams at G = 0.002 . The symbols + and × indicate the wave propagation and disappearance, respectively. (a) C 0 R tube vs. κ 1 /κ 0 . The data for τ rd /τ md = 1 , 2.5, 5, and 20 are overlayed at each point. The upper and lower dashed lines represent the upper thresholds, indicating that the cylindrical membrane shape is stable for the unexcited ( u = 0.07 ) and excited ( u = 0.9 ) membranes, respectively.
Both the symbols are overlayed when both dynamics are observed in six simulation runs.  www.nature.com/scientificreports/ Under the unstable condition (caused by high spontaneous curvature C 0 and/or more rigid proteins), the membrane exhibits buckling and buddings to a pearl-necklace shape. In our simulations, the tube radius R tube = 9.79σ.
Since the membrane ends are connected by a periodic boundary, a wave repeatedly propagates on the membrane, when the membrane deformation is small (see Figs. 2a and 4a and Supplementary Movie 1). However, the protein-binding induced deformation can generate a narrow neck in front of the wave and erases the wave (see Fig. 2b). The condition of this wave disappearance strongly depends on the time ratio τ rd /τ md of RD to membrane deformation. For very fast wave propagation at τ rd /τ md = 1 , the membrane is only slightly deformed even at a high bending rigidity of the bound membrane ( κ 1 /κ 0 = 5 ) for C 0 R tube = 3 (see Fig. 3b). Interestingly, a reentrant transition is observed; as τ rd /τ md increases, waves disappear first and then propagate again (see Fig. 3b and c). A slow wave propagation allows the bound membrane to form a narrow tube with a roughly constant length because the membrane can relax into a thermal-equilibrium shape (see Figs. 2c and 4b and Supplementary Movie 2). A more rigid bound membrane forms a shorter tube (see Fig. 4b).
In undeformed tubes, the wave velocity v wave is proportional to 1/τ rd . However, the membrane deformation slows down the waves so that v wave decreases with increasing κ 1 and C 0 (see Fig. 4c). Note that the wave propagation has hysteresis around the phase boundaries; that is, when the wave of a narrow tube, as shown in the bottom three snapshots in Fig. 2c, is used for an initial state, the wave propagates for a slightly faster RD of τ rd /τ md = 10 in most cases (five of six runs), but the wave disappears for τ rd /τ md = 5 at κ 1 /κ 0 = 4 and C 0 R tube = 3 (see Supplementary Movie 3).
At high C 0 and long τ rd , several spherical buds are formed before the disappearance (see Fig. 5a and Supplementary Movie 4). The excited membranes divide into several small domains and each domain forms a bud. www.nature.com/scientificreports/ Since the excited domains are separated from the membrane tube by narrow necks, the waves are stopped and membrane returns back into a cylindrical shape. Although the membrane topology is fixed in our simulation, vesicle formation from the buds are expected if topological changes are allowed as in Refs. 32 and 33 . Next, we discuss the dependence on the feedback constant G. Since the nullcline for u shifts downwards with increasing G, waves become slower (see Supplementary Fig. 1) and eventually disappear at high G even for τ rd /τ md = 1 (see Fig. 6). Interestingly, a wave division is observed at high C 0 (indicated by • in Fig. 6c): a wave changes its shape from a ring to a strip and subsequently, a spot-shaped excited domain is pinched off; the spot grows to a long strip and rotates around the tube axis (see Fig. 5b and Supplementary Movie 5). Thus, wave propagation can change into an azimuthal direction through membrane deformation.
To examine the robustness of these dynamics, we simulate a tubular vesicle with the reduced volume V * = V /(4πR ves 3 /3) = 0.4 where R ves = √ A/4π = 25.5σ , and V and A are the volume and surface area of the vesicle, respectively. Similar wave dynamics are observed in a tubular vesicle as shown in Supplementary Figs. 2 and 3 and Supplementary Movies 6 and 7. A difference is that the vesicle length is finite so that the wave is terminated at the vesicle end. Since a free-standing vesicle can bend significantly in fluctuation, the cooccurrence region in the phase diagram is larger than the membrane tube (compare Supplementary Fig. 3 and Fig. 6). The wave rotation around the tube axis before the disappearance is also obtained at high G as in the membrane tube (see Supplementary Figs.

2b and Supplementary Movie 6).
Waves in tubular networks. In living cells, organelles, such as ER and Golgi apparatus, have branched tubular networks [34][35][36] . They can be considered as a circuit for RD waves. Theoretically, multiple states can coexist in such networks 37,38 . As the simplest tubular network, we consider a genus-2 vesicle, which has two branching points, at V * = 0.4 , as shown in Fig. 7. Vesicles with two and higher genus can be produced experimentally [39][40][41][42] , and their shapes can be reproduced using the simulation method used in this study 42 .
For small-deformation conditions, a wave is divided into two at the branching points. Both the waves propagate separately (see Fig. 7a); subsequently, they collide in the left tube branch and disappear together (see the right snapshots in Fig. 7a); the remaining wave returns to the original position. These division and collision dynamics occur repeatedly (see Fig. 8a and Supplementary Movie 8). Similar to straight tubes, the waves in tubular networks disappear before reaching the first branching point at high κ 1 and/or high C 0 (see the line for C 0 R ves = 10 in Fig. 8a and the phase diagram in Fig. 8b and c). In intermediate conditions, waves occasionally disappear; however, some of the waves stochastically continue to propagate (see Figs. 7b and 8). In this case, waves often propagate in the opposite direction: for example, Fig 7b shows that the wave propagates upwards in the middle branch and clockwise in the right branch. Thus, the wave pathways are stochastically changed by membrane deformation. www.nature.com/scientificreports/

Discussion
In this study, we have shown that the membrane deformation generated by the wave propagation and feedback to the wave exhibit characteristic behaviors, which are not seen in undeformable surfaces. When the wave propagates faster than the membrane deformation, the deformation is minimal. Conversely, for slower wave propagation, the membrane forms a thermal-equilibrium shape for the excited domain. When the shape of equilibrium is a narrow cylindrical tube locally, the excited domain moves translationally while retaining its shape. However, when spherical buds are formed, the wave disappears. In the case of the compatible speeds, the formation of membrane neck erases the wave. When the feedback of deformation to the protein binding is strong, the wave can be destabilized owing to the modification of reaction equation. In this case, a strip shape of waves can be formed, and subsequently, wave division and azimuthal rotation occur. Thus, reaction waves exhibit various dynamics due to the coupling of membrane deformation. In this study, unbound membranes are considered as laterally homogeneous. However, biomembranes are heterogeneous and lipid rafts are formed as a platform of the functions 43 . Hence, different wave modes can be generated depending on the position of membrane in living cells. For example, the wave propagation induces membrane budding and subsequent vesicle formation at specific membrane locations. It was reported that endocytosis is enhanced in the presence of propagating waves of curvature-inducing proteins 44 . In living cells, vesicles are generated in ER, Golgi apparatus, and plasma membranes in intracellular membrane traffic 15,16,35 . In particular, ER exhibits necking, and its contacts with mitochondria may trigger chemical waves 35,36 . Thus, the wave-induced membrane deformation can play a role in ER and other organelles.
The network structures of membranes are observed in vitro and in vivo. They can work as the circuits for RD waves. We have demonstrated that membrane deformation can change the propagation pathways in genus-2 vesicles. This suggests that the deformation of membrane tubes or other soft substrates, such as gels, can be used as a tool to switch the wave circuits. Moreover, such switches may play a role in control signal transductions in living cells.

Methods
The membrane tubes and vesicles are represented by triangular meshes of N = 15000 and 10000 vertices, respectively. The membrane potentials are discretized as in Ref. 45 (see Supplementary Information). To produce a lateral fluidity, the meshes are stochastically reconstructed by a bond-flip Monte Carlo method 45,46 . Straight membrane tubes lie along the x axis with a length of L x = 200σ . An initial state is set up as follows: first, the membrane tube is equilibrated in the absence of excitation. Then, the protein concentration u is locally raised by 0.8 for 20 < (x − x end )/σ < 40 , and v by 0.00125(x − x end )/σ for 0 < (x − x end )/σ ≤ 20 and by 0.00125(40 − x + x end )/σ for 20 < (x − x end )/σ < 40 , where x end is the left end position of the membrane tube. For tubular and genus-2 vesicles, initial states are similarly set up (see Supplementary Information).
We consider the membrane vertices of u > 0.5 as the excited wave region. The center position x wave and width w wave of the wave are calculated by x wave = u i >0.5 x i /N wave and w 2 wave = u i >0.5 (x i − x wave ) 2 /N wave , where N wave is the number of vertices belonging to the wave region. The area ratio of the wave is obtained as φ wave = N wave /N . Error bars are calculated from three or more independent runs. To calculate the phase diagrams, six and ten simulation runs are conducted for regions close to the phase boundaries of wave disappearance and division, respectively. Additional details about the methods are described in Supplementary Information.

Data availability
All data generated or analysed during this study are included in this published article [and its supplementary information files]. The data for C 0 R ves = 6 and 8 correspond to those shown in Fig. 7a and b, respectively. (b) Phase diagram for κ 1 /κ 0 = 4 . (c) Phase diagram for C 0 R ves = 8 . The symbol + indicates the wave propagation along both branches. The symbol × indicates the immediate wave disappearance before reaching the first branching point. The symbol indicates the partial wave propagation, in which one of the waves disappears at least before the wave reaches the original position.