# Excitable reaction-diffusion waves of curvature-inducing proteins on
deformable membrane tubes

Naoki Tamemoto  Institute for Solid State Physics, University of Tokyo,
Kashiwa, Chiba 277-8581, Japan  Hiroshi Noguchi<EMAIL_ADDRESS>Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba
277-8581, Japan

###### Abstract

Living cells employ excitable reaction-diffusion waves for internal cellular
functions, in which curvature-inducing proteins are often involved. However,
the role of their mechanochemical coupling is not well understood. Here, we
report the membrane deformation induced by the excitable reaction-diffusion
waves of curvature-inducing proteins and the alternation in the waves due to
the deformation, using a coarse-grained simulation of tubular membranes with a
modified FitzHugh–Nagumo model. Protein-propagating waves deform tubular
membranes, and large deformations induce budding and erase waves. The wave
speed and shape are determined by a combination of membrane deformation and
spatial distribution of the curvature-inducing protein. Waves are also
undulated in the azimuthal direction depending on the condition. Rotationally
symmetric waves locally deform the tubes into a symmetric shape but maintain a
straight shape on average. Our simulation method can be applied to other
chemical reaction models and used to investigate various biomembrane
phenomena.

##  I Introduction

Membrane deformation is a common phenomenon in living cells Zimmerberg and
Kozlov (2006); McMahon and Gallop (2005); Suetsugu  _et al._ (2014); Voeltz
and Prinz (2007). The shapes of biomembranes are regulated by various proteins
including curvature-inducing proteins Zimmerberg and Kozlov (2006); McMahon
and Gallop (2005); Suetsugu  _et al._ (2014); Baumgart  _et al._ (2011). For
example, Bin/Amphiphysin/Rvs (BAR) superfamily proteins, which have a banana-
shaped dimer structure, bind to membranes and bend them according to intrinsic
protein shapes Suetsugu  _et al._ (2014); Baumgart  _et al._ (2011); Peter
_et al._ (2004). Clathrin, and coat protein complexes (COPI and COPII) also
bend membranes and induce membrane budding. These proteins cooperatively work.
For example, during clathrin-mediated endocytosis, the BAR proteins bind to
membranes, and subsequently, clathrin-assembly forms a spherical bud; the
binding of dynamin to the bud neck induces membrane fission McMahon and
Boucrot (2011). While proteins bend membranes as described above, the membrane
curvature affects protein binding (curvature sensing) McMahon and Gallop
(2005); Baumgart  _et al._ (2011); Peter  _et al._ (2004); Zhao  _et al._
(2017); Sorre  _et al._ (2012); Has and Das (2021). Various curvature-inducing
proteins have different preferences for high membrane curvatures depending on
protein shape. Thus, membrane curvature also affects their dynamics Zhao  _et
al._ (2017).

Recently, there have been experimental and theoretical studies on the feedback
regulation between protein dynamics and membrane mechanics. Reaction-diffusion
waves of proteins coupled with membrane deformation have been observed and
theoretically analyzed under various conditions Wu  _et al._ (2018); Wu and
Liu (2021); Peleg  _et al._ (2011); Gov (2018). Further, the membrane shape
affects intracellular proteins’ wave propagation, and they subsequently
regulate cell motion Weiner  _et al._ (2007); Taniguchi  _et al._ (2013); Miao
_et al._ (2019). The dynamics of the curvature-inducing proteins are
associated with cellular functions, such as cell migration Tsujita  _et al._
(2015), and cell division Xiao  _et al._ (2017). In addition, giant
unilamellar vesicles (GUVs) containing Min proteins, which generate reaction-
diffusion waves in Escherichia coli and regulate cell division, exhibit cyclic
deformation accompanying the Min protein waves Godino  _et al._ (2019); Christ
_et al._ (2021); Litschel  _et al._ (2018). In signal transduction, receptors
on cell membranes receive extracellular signals, transfer information
downstream in the intracellular reaction network, and regulate cellular
functions according to the surrounding environment, such as chemotaxis Xiong
_et al._ (2010). Some of these cell signaling are modeled using excitable
reaction-diffusion models, and various studies have shown that the signaling
network regulates cell migration Huang  _et al._ (2013); Zhan  _et al._
(2020). Researchers have analyzed excitable reaction-diffusion waves on fixed
curved surfaces and found that the surface geometry affects the stability,
speed, and shape of the waves Horibe  _et al._ (2019); Kneer  _et al._ (2014).
However, studies that have considered membrane deformation and reaction-
diffusion dynamics are limited, and their biological implications are not
understood yet.

In previous studies, we examined the coupling between the large deformation of
a fluid vesicle and reaction-diffusion dynamics with the binding of curvature-
inducing proteins using a dynamically triangulated membrane model and a non-
excitable reaction-diffusion model, Brusselator model, to investigate the
mechanochemical feedback Tamemoto and Noguchi (2020, 2021). We understood that
membrane deformation stabilizes Turing patterns, simultaneously Turing
patterns induce membrane deformation, resulting in the formation of budded and
multispindle-shaped vesicles Tamemoto and Noguchi (2020). Moreover, we
successfully reproduced the self-oscillation of vesicle shapes as seen in the
reconstructed Min system in GUVs, and showed that the starting point of
traveling waves is determined by the local membrane curvature Tamemoto and
Noguchi (2021).

In this study, we simulated a coupled system between membrane deformation and
excitable reaction-diffusion waves of curvature-inducing proteins. Excitable
reaction-diffusion waves are fundamental biological phenomena, and models,
such as the FitzHugh–Nagumo model for nerve signaling, have been used for
studying them FitzHugh (1961); Nagumo  _et al._ (1962). Moreover, various
protein waves interacting with membrane mechanics are demonstrated using
excitable reaction-diffusion models Wu  _et al._ (2018); Wu and Liu (2021). We
employed the modified FitzHugh–Nagumo model to be coupled with membrane
deformation. Here, we consider that the curvature-inducing proteins have a
laterally isotropic shape. We employed tubular membranes and waves propagating
along the tube axis to reduce the complexity of waves with higher dimensions.
Tubular membrane structures are observed in cell membranes (e.g., axons of
neurons) as well as intracellular organelles (e.g., endoplasmic reticulum (ER)
and mitochondria) McMahon and Gallop (2005); Voeltz and Prinz (2007). Tubular
membranes can be produced in vitro by pulling a vesicle using optical tweezers
and a micropipette, and the tube radius is controlled by the amplitude of the
imposed force Dimova (2014); Roux (2013). Experimentally, tubular membranes
are often used to investigate the curvature dependence of protein binding to
membranes Baumgart  _et al._ (2011); Sorre  _et al._ (2012); Has and Das
(2021). We investigated the dynamics of protein-propagating waves using a
combination of the FitzHugh–Nagumo model and tubular membranes.

##  II Methods

###  II.1 Membrane model

We applied the membrane model used in our previous studies Tamemoto and
Noguchi (2020, 2021) with a few modifications. We triangulated a membrane
surface to discretize the membrane tube. In this triangular network, $N$
vertices are connected by bonds of average length $a$, and the vertices have
excluded volumes and masses $m$. The surface area $S$ is constrained by the
harmonic potentials Noguchi and Gompper (2005). The tube lies along the
$x$-axis and is connected using a periodic boundary condition. Unlike our
previous studies, the volume enclosed by the membrane is left unconstrained in
this study. This tube is considered to be extended from a vesicle by optical
tweezers or protein filaments, so that the tube is connected to a reservoir of
a large volume. The triangulated bond network is re-meshed using the bond-flip
Monte Carlo method to allow lateral diffusion of vertices Gommper and Kroll
(2004).

The deformation of the tubular membrane is expressed by free energy of the
curvature, which is modified to include the effect of curvature-inducing
proteins. The curvature energy is expressed as $F_{\mathrm{cv}}=\int
f_{\mathrm{cv}}dS$ with

$f_{\mathrm{cv}}=(1-u)\frac{\kappa_{0}}{2}(2H)^{2}+u\frac{\kappa_{1}}{2}(2H-C_{0})^{2},$
(1)

where $u$ is the concentration of curvature-inducing proteins on the membranes
($u\in[0,1]$); $H$ is the mean curvature; $\kappa_{1}$ and $\kappa_{0}$ are
the bending rigidity with and without the proteins on the membranes,
respectively; and $C_{0}$ is the spontaneous curvature of the curvature-
inducing proteins Tamemoto and Noguchi (2020, 2021); Noguchi (2021a, 2022).
Although the difference in the saddle-splay modulus $\bar{\kappa}$ (also
called the Gaussian modulus) can modify the protein binding Noguchi (2021a,
2022), it is not considered here for simplicity. For the budding of a vesicle,
the $\bar{\kappa}$ difference can be accounted for by the rescaling of
$\kappa_{1}$ and $C_{0}$ Noguchi (2021a). Similarly, the $\bar{\kappa}$
difference does not likely cause a qualitative change in the present case.
When $\bar{\kappa}$ is constant for an entire membrane, $\bar{\kappa}$ does
not contribute to membrane shape transformation for a fixed topology because
of the Gauss–Bonnet theorem. According to this curvature energy expression,
the larger the protein concentration $u$, the more strongly the membrane bends
when $C_{0}$ is sufficiently large.

Membrane deformation is solved through molecular dynamics simulation using a
Langevin thermostat:

$m\frac{\partial^{2}\bm{r}_{i}}{\partial t^{2}}=-\frac{\partial
U}{\partial\bm{r}_{i}}-\zeta\frac{\partial\bm{r}_{i}}{\partial
t}+\bm{g}_{i}(t),$ (2)

where $\zeta$ is the friction coefficient and $\bm{g}_{i}$ is Gaussian white
noise, which obeys the fluctuation-dissipation theorem. We use the potential
$U=U_{\mathrm{S}}+U_{\mathrm{b}}+U_{\mathrm{r}}+U_{\mathrm{cv}}$, where
$U_{\mathrm{S}}$ is the constraint potential for the surface area $S$;
$U_{\mathrm{b}}$ and $U_{\mathrm{r}}$ are the bond and repulsive potentials,
respectively; and $U_{\mathrm{cv}}$ is the discretized potential for the
bending energy, $F_{\mathrm{cv}}$, using the dual lattice.

$\displaystyle U_{\mathrm{S}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}k_{\mathrm{S}}(S-S_{0})^{2},$ (3) $\displaystyle
U_{\mathrm{b}}$ $\displaystyle=$
$\displaystyle\sum_{\mathrm{bond}}\frac{b\exp(1/(l_{\mathrm{c0}}-r_{i,j}))}{l_{\mathrm{max}}-r_{i,j}}\Theta(r_{i,j}-l_{\mathrm{c0}}),$
(4) $\displaystyle U_{\mathrm{r}}$ $\displaystyle=$
$\displaystyle\sum_{\mathrm{all\
pairs}}\frac{b\exp(1/(r_{i,j}-l_{\mathrm{c1}}))}{r_{i,j}-l_{\mathrm{min}}}\Theta(l_{\mathrm{c1}}-r_{i,j}),$
(5) $\displaystyle U_{\mathrm{cv}}$ $\displaystyle=$
$\displaystyle(1-u)\frac{\kappa_{0}}{2}\sum_{i}\frac{1}{\sigma_{i}}\bigg{(}\sum_{j(i)}\frac{\sigma_{i,j}{\bf
r}_{i,j}}{r_{i,j}}\bigg{)}^{2}$
$\displaystyle+u\frac{\kappa_{1}}{2}\sum_{i}\sigma_{i}\bigg{(}\frac{1}{\sigma_{i}}\sum_{j(i)}\frac{\sigma_{i,j}{\bf
r}_{i,j}}{r_{i,j}}-C_{0}{\bf n}_{i}\bigg{)}^{2},$

where $\Theta(x)$ is the unit step function, $r_{i,j}$ is the distance between
two vertices $i$ and $j$, and ${\bf n}_{i}$ is the normal vector at the $i$-th
vertex. The sum over $j(i)$ is over the neighbors of the $i$-th vertex, which
are connected by tethers. The length of a bond in the dual lattice is
$\sigma_{i,j}=r_{i,j}[\cot(\theta_{1})+\cot(\theta_{2})]/2$, where the angles
$\theta_{1}$ and $\theta_{2}$ are opposite to bond $ij$ in the two triangles
sharing this bond. $\sigma_{i}=0.25\sum_{j(i)}\sigma_{i,j}r_{i,j}$ is the area
of the dual cell of vertex $i$. Here, we use
$k_{\mathrm{S}}=4k_{\mathrm{B}}T$, $b=80k_{\mathrm{B}}T$,
$l_{\mathrm{max}}=1.33a$, $l_{\mathrm{c0}}=1.15a$, $l_{\mathrm{c1}}=0.85a$,
and $l_{\mathrm{min}}=0.67a$. The details of these potentials are provided in
Ref. Noguchi and Gompper (2005).

###  II.2 Reaction-diffusion model

We employed the FitzHugh–Nagumo model FitzHugh (1961); Nagumo  _et al._ (1962)
as an excitable reaction-diffusion model to investigate stable propagating
waves with excitable dynamics. Reaction-diffusion equations for concentrations
of curvature-inducing proteins $u$ and regulatory proteins $v$ are given by

$\displaystyle\begin{split}\tau\frac{\partial u}{\partial
t}&=-(u-0.2)(u-A)(u-0.9)-v-G\frac{\partial f_{\mathrm{cv}}}{\partial u}\\\
&\quad+D_{u}\nabla^{2}u+Q_{u},\end{split}$ (7)
$\displaystyle\tau\frac{\partial v}{\partial t}$
$\displaystyle=B(C(u-0.2)-v)+D_{v}\nabla^{2}v+Q_{v},$ (8)

where $\tau$ is a time constant; $D_{u}$ and $D_{v}$ are the diffusion
coefficients of the two types of proteins; $A$, $B$, and $C$ are positive
reaction parameters; $Q_{u}$ and $Q_{v}$ are the stimuli for wave initiation
at a thin slice of the tube; and $G$ is the magnitude of mechanochemical
coupling. To introduce the effect of membrane curvature on the reaction-
diffusion model, we assumed that the binding rate of curvature-inducing
proteins onto membranes is proportional to the partial derivative of the local
curvature energy with respect to the concentration of curvature-inducing
proteins, $\partial f_{\mathrm{cv}}/\partial u$, as in our previous study
Tamemoto and Noguchi (2020). According to this coupling, the binding of
curvature-inducing proteins is enhanced around the preferred curvature
$H\simeq C_{0}/2$, where $\partial f_{\mathrm{cv}}/\partial u<0$. We constrain
the range of $u$ during the time evolution to maintain $u\in[0,1]$; when Eq.
(7) gives $u$ out of this range, it is set to $u=0$ or $u=1$. When a small
stimulus is given in the equilibrium state, the state quickly recovers to the
fixed point. However, the system shows excitability when a larger stimulus is
given (see Fig. 1).

Figure 1:  Phase plane of the FitzHugh–Nagumo model with the modified protein
unbinding process for $A=0.4$, $B=0.001$, $C=0.1$, and $G=0$. Purple and green
lines are nullclines for $\partial u/\partial t=0$ and $\partial v/\partial
t=0$, respectively. Black lines indicate trajectories after the equilibrium
state is stimulated.

###  II.3 Parameters

In this study, we used membrane tubes with a length of $L_{x}=200a$, where the
number of vertices is $N=15~{}000$ or $30~{}000$. For these numbers of
vertices, the tubular surface areas are $S\simeq 12~{}300a^{2}$ and
$24~{}600a^{2}$, and the corresponding radii of the tubes are
$R_{\mathrm{t}}=S/2\pi L_{x}\simeq 9.79a$ and $19.6a$, respectively. Unless
otherwise mentioned, we applied the conditions
$\kappa_{0}/k_{\mathrm{B}}T=20$, $A=0.4$, $B=0.001$, $C=0.1$,
$D_{u}/a^{2}=0.1$, $D_{v}=0$, $Gk_{\mathrm{B}}T/a^{2}=0.002$, and
$\tau_{\mathrm{md}}=\tau$, where $\tau_{\mathrm{md}}=\zeta
a^{2}/k_{\mathrm{B}}T$. For this ratio of the time units of the reaction and
membrane motion, the membrane motion can follow the reaction dynamics. The
starting protein concentration values were set around the fixed point
calculated using a fixed value of the membrane curvature
$(H=1/2R_{\mathrm{t}})$, with small, random perturbations.

For the initiation of propagating waves, we added stimuli at the thin slice of
the tubular membrane, given by

$\displaystyle Q_{u}$ $\displaystyle=\exp\left(-0.0625(x-10)^{2}\right)\text{
and}$ (9) $\displaystyle Q_{v}$
$\displaystyle=0.02\exp\left(-0.25(x-5)^{2}\right)$ (10)

during the initial short period ($0\leq t/\tau\leq 1$). The profiles of
$Q_{u}$ and $Q_{v}$ are shown in Fig. 2(a), and they induce a one-dimensional
wave that propagates on a flat membrane with $G=0$ as shown in Fig. 2(b).

Figure 2:  Initiation of propagating waves. (a) The profiles of the stimuli
$Q_{u}$ and $Q_{v}$ to initiate the propagating wave. (b) Propagating waves
induced by the stimuli shown in (a). The color variation indicates different
concentration of the curvature-inducing protein, $u$.

Equations (7) and (8) were integrated numerically using the forward difference
method for time and the finite-volume scheme for space Tamemoto and Noguchi
(2020). Equation (2) was integrated using the leapfrog method. Error bars were
calculated based on three independent runs.

In the estimation of local curvature, we calculated a smoothed local curvature
$\tilde{H}$ by averaging the local curvature $H$ up to the second-order
adjacent vertices. The details are described in the Supplementary Material of
Ref. Tamemoto and Noguchi (2020).

###  II.4 Stability of Nonexcited Membrane

In thermal equilibrium, a cylindrical membrane tube generates an force
$f_{\mathrm{ex}}$ along the tube axis Noguchi (2021b):

$f_{\mathrm{ex}}=2\pi\kappa_{0}C_{0}\bigg{\\{}\frac{1}{R_{\mathrm{t}}C_{0}}\Big{[}\Big{(}\frac{\kappa_{1}}{\kappa_{0}}-1\Big{)}u+1\Big{]}-\frac{\kappa_{1}}{\kappa_{0}}u\bigg{\\}}.$
(11)

For $f_{\mathrm{ex}}>0$, the membrane tube is stable and this force is
balanced with an external force (here imposed by the periodic boundary
condition). In contrast, for $f_{\mathrm{ex}}\leq 0$, the cylindrical tubes
are unstable, so that the unstable condition is given by

$C_{0}R_{\mathrm{t}}\geq 1+\frac{\kappa_{0}(1-u)}{\kappa_{1}u}.$ (12)

This condition is identical to that in Ref. Ou-Yang and Helfrich (1989) (the
membrane spontaneous curvature $\geq$ the tube curvature). In particular,
$f_{\mathrm{ex}}=0$, the membrane deforms spherical buds via unduloid
deformation, maintaining the mean curvature of the membrane everywhere
Kenmotsu (2003); Naito  _et al._ (1995). Previously, budding of a tubular
membrane with a positive spontaneous curvature has been reported Tsafrir  _et
al._ (2001); Campelo and Hernández-Machado (2007); Tozzi  _et al._ (2019);
Noguchi (2021b). Note that the budding also occurs slightly below this
threshold under thermal fluctuations Noguchi (2021b). Buckling of the tube is
induced by a large excess membrane area when the tube is compressed (negative
surface tension) Liu  _et al._ (2016); den Otter  _et al._ (2003). In our
simulation condition, an unexcited membrane has $u\simeq 0.2$ (see Fig. 1).
For example, the cylindrical tube is unstable at $C_{0}R_{\mathrm{t}}\gtrsim
3$ for $\kappa_{1}/\kappa_{0}=2$.

##  III Results and Discussions

###  III.1 Propagating waves with mechanochemical coupling

We investigated the propagating waves on tubular membranes for
$R_{\mathrm{t}}=9.79a$ and $19.6a$ (narrow tubes and wide tubes,
respectively). Before examining the deformable tubes, we evaluated the effects
of constant positive or negative inputs on waves on a non-deformable narrow
membrane tube. Instead of $-G\partial f_{\mathrm{cv}}/\partial u$, we added a
constant input $I$ to Eq. (7) over the entire membrane-tube region. Figures
3(a)–(c) show typical snapshots of the propagating waves for $I=0.002$, $0$,
and $-0.002$, respectively. Among them, the width and speed of the waves were
largest for $I=0.002$ and smallest for $I=-0.002$. The wave is generated once
but shrinks and disappears without propagating for $I=-0.01$ (see Fig. 3(d)).
Figure 4 shows the changes in wave width, $w_{\mathrm{wave}}$ and speed,
$s_{\mathrm{wave}}$. Both of wave speed and width decrease with decreasing
$I$, and the waves disappear when $I$ is smaller than the threshold.

Figure 3:  Typical sequential snapshots of propagating waves on non-deformable
tubular membranes for (a) $I=0.002$, (b) $I=0$, (c) $I=-0.002$, and (d)
$I=-0.01$. The color variation indicates different concentration of the
curvature-inducing protein, $u$.  Figure 4:  Effect of constant inputs $I$ on
the wave properties with non-deformable membranes. In the gray shaded area,
propagating waves disappear.

First, we investigated the effects of mechanochemical coupling on propagating
waves in deformable tubular membranes with a small radius
($R_{\mathrm{t}}=9.79a$). A straight wave propagates at a small spontaneous
curvature $C_{0}$, similar to non-deformable tubes, as the coupling effect is
weak (see Fig. 5(a) and Supplementary Movie S1). However, undulations of the
membrane shape and wavefront are observed at an intermediate $C_{0}$ (see Fig.
5(b) and Supplementary Movie S2). At $C_{0}R_{\mathrm{t}}=3$ and
$\kappa_{1}/\kappa_{0}=3$, the tube is deformed in to a meandering form, as
shown in the first snapshot of Fig. 5(b), owing to the instability of
cylindrical shapes. When a wave is stimulated in this meandering tube, the
wave propagates on the tube. Owing to the tube deformation, the shape of the
protein wave is deformed, which could also affect the membrane shape. In
addition, the waves are unstable and eventually disappear at large $C_{0}$
(see Fig. 6 and Supplementary Movie S3 for $C_{0}R_{\mathrm{t}}=4$ and
$\kappa_{1}/\kappa_{0}=4$). These results indicate that the interaction
between membrane deformation and propagating waves changes the stability of
the membrane shapes and protein waves.

Figure 5:  Typical sequential snapshots of propagating waves on deformable
tubular membranes with a small radius ($R_{\mathrm{t}}=9.79a$) and small
spontaneous curvatures $C_{0}$. (a) $C_{0}R_{\mathrm{t}}=1$ and
$\kappa_{1}/\kappa_{0}=2$. (b) $C_{0}R_{\mathrm{t}}=3$ and
$\kappa_{1}/\kappa_{0}=3$. The color variation indicates different
concentration of the curvature-inducing protein, $u$. The corresponding movies
are provided in Supplementary Movies S1 and S2.  Figure 6:  Typical sequential
snapshots of propagating waves on deformable tubular membranes with the small
radius ($R_{\mathrm{t}}=9.79a$) with large spontaneous curvatures $C_{0}$.
$C_{0}R_{\mathrm{t}}=4$ and $\kappa_{1}/\kappa_{0}=4$. The corresponding movie
is provided in Supplementary Movie S3.

Figure 7 shows the movement of the center of mass of the waves
($g_{\mathrm{wave},x}$) along the tube ($x$) axis. The wave moves almost
linearly for small spontaneous-curvature conditions (purple and green lines in
Fig. 7); however, the wave is weakly wobbled by the membrane deformation. The
inputs from the mechanochemical effect on the reaction of $u$ for
$H=1/2R_{\mathrm{t}}$ are $-G\partial f_{\mathrm{cv}}/\partial u=0.0002$ and
$-0.0023$ for the purple and green lines in Fig. 7, respectively. These
results indicate that the waves move at an approximate constant speed even on
deformable tubes. The wave speed decreases with decreasing mechanochemical
inputs $-G\partial f_{\mathrm{cv}}/\partial u$, as in the case of non-
deformable tubes (Fig. 4).

At large $C_{0}$ and $\kappa_{1}/\kappa_{0}$ values, $g_{\mathrm{wave},x}$
moves slightly but the wave disappears (Fig. 6). However, at
$C_{0}R_{\mathrm{t}}=4$ and $\kappa_{1}/\kappa_{0}=4$, the mechanochemical
input $-G\partial f_{\mathrm{cv}}/\partial u|_{H=1/2R_{\mathrm{t}}}=-0.0073$,
which is still in the range of the stable wave for the non-deformable tubular
membranes (see Fig. 4). Therefore, the disappearance of the wave is partly
caused by the reduction in diffusion through the narrow neck.

Figure 7:  Time evolution of $x$-coordinate of the center of mass of the wave
($g_{\mathrm{wave},x}$), on deformable tubular membranes with the small radius
($R_{\mathrm{t}}=9.79a$). The simulation condition is the same as that of
Figs. 5 and 6.

The wave width also decreases with decreasing mechanochemical input
$-G\partial f_{\mathrm{cv}}/\partial u$ (Fig. 8). Figure 8(a) shows the time
evolution of the wave width at $C_{0}R_{\mathrm{t}}=3$ for
$\kappa_{1}/\kappa_{0}=1$, $2$, and $6$ ($-G\partial f_{\mathrm{cv}}/\partial
u=-0.0006$, $-0.0015$, and $-0.0048$ at $H=1/2R_{\mathrm{t}}$, respectively).
Straight waves propagate at $\kappa_{1}/\kappa_{0}=1$ and $2$. The wave widths
are also approximately the same everywhere and maintain a constant value in
time. Similar changes were observed at $\kappa_{1}/\kappa_{0}=6$ for
$C_{0}R_{\mathrm{t}}=1$, $2$, and $3$ ($-G\partial f_{\mathrm{cv}}/\partial
u=0.0002$, $-0.001$, and $-0.0048$ at $H=1/2R_{\mathrm{t}}$, respectively)
(Fig. 8(b)).

Figure 8:  Time evolution of the wave width on deformable tubular membranes
with a small radius ($R_{\mathrm{t}}=9.79a$). (a) $\kappa_{1}/\kappa_{0}$
varied at $C_{0}R_{\mathrm{t}}=3$. (b) $C_{0}R_{\mathrm{t}}$ varied at
$\kappa_{1}/\kappa_{0}=6$.

Mechanochemical effects were also observed for the wide tubes (Figs. 9 and
10). At $C_{0}R_{\mathrm{t}}=4$ and $\kappa_{1}/\kappa_{0}=4$, the wavefront
is undulated with a small azimuthal deformation of the membrane. In contrast,
the other tube regions maintain a straight shape (see Fig. 9 and Supplementary
Movie S4). The wavefront undulation also occurs even for small $C_{0}$ and
large $\kappa_{1}/\kappa_{0}$ (see Supplementary Movie S5)). The detection and
evaluation methods of the wavefront undulation are described in Appendix A. A
large membrane deformation is induced when the value of $C_{0}$ is large, and
the wave locally disappears and breaks due to the winding membrane shape (Fig.
10, $t/\tau=1500\text{--}2000$, see also Supplementary Movie S6). Once the
wave breaks, the shapes of the wave and membrane become largely disordered
(Fig. 10, $t/\tau=3000\text{--}8000$).

Figure 9:  Typical sequential snapshots of propagating waves on deformable
wide tubes ($R_{\mathrm{t}}=19.6a$) at $C_{0}R_{\mathrm{t}}=4$ and
$\kappa_{1}/\kappa_{0}=4$. The color variation indicates different
concentration of the curvature-inducing protein, $u$. The corresponding movie
is provided in Supplementary Movie S4.  Figure 10:  Typical sequential
snapshots of propagating waves on deformable wide tubes
($R_{\mathrm{t}}=19.6a$) at $C_{0}R_{\mathrm{t}}=8$ and
$\kappa_{1}/\kappa_{0}=8$. The color variation indicates different
concentration of the curvature-inducing protein, $u$. The corresponding movie
is provided in Supplementary Movie S6.

Figure 11 shows the movements of the $x$-coordinate of the vertices in the
wave ($r_{\mathrm{wave},x}$) on a deformable wide tube
($R_{\mathrm{t}}=19.6a$) at $C_{0}R_{\mathrm{t}}=8$ and
$\kappa_{1}/\kappa_{0}=8$. The wave breaks and splits into two or three pieces
(Fig. 11, $t/\tau=3500\text{--}7000$). Even though the shapes of the waves and
membranes are disordered, the waves move almost linearly along the membrane
tube after division.

Figure 11:  Time evolution of the $x$-coordinate of vertices in the wave
($r_{\mathrm{wave},x}$) on a deformable wide tube ($R_{\mathrm{t}}=19.6a$) at
$C_{0}R_{\mathrm{t}}=8$ and $\kappa_{1}/\kappa_{0}=8$. The simulation
condition is the same as that of Fig. 10.

Based on these results, we classified the waves into three phases, straight,
undulated, and unstable, as shown in Fig. 12. The shapes of the tubes and
waves remain straight in the straight wave phase, similar to the case of non-
deformable tubes (see Figs. 3 and 5(a)). This phase is observed when the
effect of the mechanochemical coupling is weak. In the unstable wave phase at
large $C_{0}$, the waves disappear and break. At large
$\kappa_{1}/\kappa_{0}$, the wave phase appears at smaller values of $C_{0}$,
since the mechanochemical coupling is strong. Large membrane deformations,
such as the formation of narrow necks, are also induced due to large values of
$C_{0}$ (see Fig. 6). The wave’s shape is disordered when the wave breaks.
When the disordered wave does not disappear, the tube’s shape is largely
deformed due to the disordered wave (see Fig. 10). In the intermediate
condition, tube shape and/or the wavefront is undulated but the wave
propagates. We call it undulated wave phase. The membrane shape and wavefront
are undulated and unstable, depending on the spontaneous curvature, bending
rigidity, and tube radius. At $C_{0}=0$, a cylindrical tube is stable at any
value of $\kappa_{1}/\kappa_{0}$ in thermal equilibrium. However, the wave
deforms the tube in azimuthal direction at large $\kappa_{1}/\kappa_{0}$, as
explained in Sec. III.2.

Figure 12:  Pattern diagrams with deformable tubular membranes for
$R_{\mathrm{t}}=9.79a$ and $19.6a$ (narrow tube and wide tube, respectively).
The symbols represent the simulation results. Overlapped symbols indicate the
coexistence of multiple patterns. The dashed lines are given by Eq. (12) with
$u=0.2$.

###  III.2 Analysis of undulated waves

For further investigation, we focused on the undulated wave phase (green
triangles in Fig. 12). Figure 13 shows typical snapshots of the undulated
waves and their cross-section views. For a wide tube ($R_{\mathrm{t}}=19.6a$),
at $C_{0}R_{\mathrm{t}}=4$ and $\kappa_{1}/\kappa_{0}=4$, the global shape of
the tube is not undulated, but the wavefront is winding (Fig. 13(a)). The
front and rear edges of the wave elongate perpendicularly (major axes of their
elliptical shapes at the cross sections are perpendicular). Curvature-inducing
proteins exist in two separate regions in the cross-sectional view at the
edges of the wave owing to this wavefront shape, which induce azimuthal
deformation in the membrane shape. In the azimuthal spatial distribution, the
dominant azimuthal wavenumber is two. At large $\kappa_{1}/\kappa_{0}$ and
$C_{0}=0$, a winding wavefront appears (e.g. at $C_{0}R_{\mathrm{t}}=0$ and
$\kappa_{1}/\kappa_{0}=32$, see Fig. 13(b) and Supplementary Movie S5). In
contrast, for a narrow tube ($R_{\mathrm{t}}=9.79a$), such winding wavefront
with the locally elliptic tube shape are not observed. At
$C_{0}R_{\mathrm{t}}=0$ and $\kappa_{1}/\kappa_{0}=12$, the wavefront is
winding (Fig. 13(c)). In this condition, the dominant azimuthal wavenumber of
protein distribution is one.

Figure 13:  Typical snapshots of undulated waves and cross-section views of
the waves at front and rear edges. (a) $C_{0}R=4$, $\kappa_{1}/\kappa_{0}=4$,
and $R_{\mathrm{t}}=19.6a$ (wide tube). (b) $C_{0}R_{\mathrm{t}}=0$ and
$\kappa_{1}/\kappa_{0}=32$, and $R_{\mathrm{t}}=19.6a$ (wide tube). (c)
$C_{0}R_{\mathrm{t}}=0$ and $\kappa_{1}/\kappa_{0}=12$, and
$R_{\mathrm{t}}=9.8a$ (narrow tube) The color variation indicates different
concentration of the curvature-inducing proteins, $u$.

We focused on the wavefront of the propagating wave and performed a discrete
Fourier transformation (DFT) on the cross-sections at the wavefront to
quantify the azimuthal deformation ($R_{\mathrm{t}}=19.6a$,
$C_{0}R_{\mathrm{t}}=4$, and $\kappa_{1}/\kappa_{0}=4$). $A_{u_{q}}$ and
$A_{H_{q}}$ indicate the $q$th amplitudes of DFT for the concentration of the
curvature-inducing protein $u$ and smoothed local membrane curvature
$\tilde{H}$ in an azimuthal direction, respectively. The magnitude of
undulation of the wavefront is quantified by the value $\langle
d_{\mathrm{RMSD}}\rangle$, which is the average of root mean squared
deviations between the wavefront and local geodesic lines, described in
Appendix A. Figure 14 shows that the undulation of the wavefront increases,
and the $q=2$ modes of $u$ and $\tilde{H}$ dominantly grow together.

Figure 14:  Time evolution of the magnitude of undulation of wavefronts,
$\langle d_{\mathrm{RMSD}}\rangle$ for $C_{0}R_{\mathrm{t}}=4$,
$\kappa_{1}/\kappa_{0}=4$, and $R_{\mathrm{t}}=19.6a$ (wide tube), with $q$th
amplitudes of DFT for (a) $u$ and (b) $\tilde{H}$ in an azimuthal direction,
$A_{u_{q}}$ and $A_{H_{q}}$, respectively ($q=1$ and $2$). The simulation
condition is the same as that of Figs. 9 and 13(a).

It is considered that the azimuthal wavenumber is determined by a combination
of the spatial structures of the reaction-diffusion waves and tubular membrane
deformation through mechanochemical coupling. For the wide tube, the
circumference is long enough to form the winding wavefront with azimuthal
wavenumber $q=2$. The tubes deform in to the locally elliptic shape. For the
narrow tube, the winding wavefront with azimuthal wavenumber $q=1$ is formed.

##  IV Conclusions

In this study, we have studied the coupled dynamics between the deformation of
tubular membranes and excitable reaction-diffusion waves of curvature-inducing
proteins. Three dynamic modes appear depending on spontaneous curvature,
bending rigidity, and tube radius: the straight wave phase with small membrane
deformation, undulated wave phase, and unstable wave phase. Membrane
deformation changes the protein-binding rate and diffusion along the tube
axis, and they modify the wave shape. The waves disappear at narrow membrane
necks in the unstable phase.

The waves with azimuthal wavenumber $q=2$ (twofold rotationally symmetric)
grow with the undulation of the wavefront for the wide tubes, in which the
membrane locally deforms in an elliptic shape, and the front and rear edges of
the wave elongate perpendicularly. Three (or more) -fold symmetric waves may
propagate for much wider tubes. Although axisymmetric assumption has often
been used for membrane simulations, full three-dimensional simulations are
needed to investigate non-axisymmetric waves, as demonstrated in this study.

The shape of organelles is dynamically regulated in living cells. For example,
the stack structure of the Golgi apparatus is formed by the cis side and
deconstructed from the opposite side Matsuura-Tokita  _et al._ (2006);
Papanikou and Glick (2014). In this dynamic equilibrium state, curved membrane
structures such as the cisternae and tubular networks are regulated by various
proteins at different maturation stages Kulkarni-Gosavi  _et al._ (2019).
Pearled tubular membranes are formed as transport intermediates for cargo
delivery from ER to Golgi. The formation of these structures is associated
with COPI and COPII Weigel  _et al._ (2021); Zucker and Kozlov (2022). In
addition, the COPI travels during the cargo delivery Weigel  _et al._ (2021).
Reaction-diffusion waves may play a role in these structural changes. Further,
endocytosis in hotspots (emerging from the same position) is enhanced in the
presence of propagating waves of curvature-inducing proteins Yang  _et al._
(2017). The propagation of membrane undulations of a certain wavenumber with
protein propagation waves might be related to the mechanism of these spatial
pattern generation. We employed the FitzHugh–Nagumo model to generate an
excitable chemical wave. However, our simulation method can be applied to
other chemical reaction models and used to investigate various biomembrane
phenomena.

This study considered the proteins that induce laterally isotropic spontaneous
curvatures. BAR superfamily proteins and other proteins have laterally
anisotropic shapes (often modeled as an elliptic shape) and bend the bound
membrane anisotropically Baumgart  _et al._ (2011); Suetsugu  _et al._ (2014);
they exhibit different behaviors from isotropic proteins. The BAR proteins
induce membrane tubulation Baumgart  _et al._ (2011); Suetsugu  _et al._
(2014); Mim and Unger (2012); Kabaso  _et al._ (2011); Noguchi (2019), whereas
the isotropic proteins induce the formation of spherical buds Hurley  _et al._
(2010); Foret (2014); Goutaland  _et al._ (2021); Noguchi (2021a). The maximum
binding (sensing) curvature of anisotropic proteins depends on the binding
chemical potential, whereas that of isotropic proteins is constant Noguchi
_et al._ (2022); Noguchi (2022). Chemical waves involving BAR proteins have
been observed Wu  _et al._ (2018); Yang  _et al._ (2017); Xiao  _et al._
(2017), but the effects of anisotropy in the waves are not understood so far.
It is essential to further studies how such differences in curvature-inducing
proteins change coupling with reaction-diffusion dynamics.

Figure 15:  Calculation of local geodesic lines and $\langle
d_{\mathrm{RMSD}}\rangle$ for a narrow tube ($R_{\mathrm{t}}=9.79a$,
$C_{0}R_{\mathrm{t}}=2$, and $\kappa_{1}/\kappa_{0}=6$). (a) Lines and points
obtained in the calculation process of the local geodesic line. (b) Snapshots
of waves and local geodesic lines at $t/\tau=2000$ and $t/\tau=3000$. The
color in snapshots indicates the concentration of curvature-inducing proteins,
$u$. (c) Time development of $\langle d_{\mathrm{RMSD}}\rangle$.

###### Acknowledgements.

We would like to thank S. Ishihara and M. Yanagisawa for valuable discussions.
This work was supported by JSPS KAKENHI Grant Number JP21K03481.

##  Appendix A Algorithm to evaluate the wavefront undulation

First, we determine that a spatial wave exists if $u_{\mathrm{max}}\geq 0.55$
and $u_{\mathrm{max}}-u_{\mathrm{min}}>0.4$, where $u_{\mathrm{max}}$ and
$u_{\mathrm{min}}$ are the maximum and minimum values of $u$, respectively.
Vertices for $u\geq 0.55$ are classified to the wave region when the wave
exists, and the other vertices to the non-wave region. When vertices in the
wave and non-wave regions are adjacent to each other, the internally dividing
point, in which $u=0.55$ on the bond between the vertices, is considered a
wave edge. The wavefront is determined by the condition on $v$ at the edge,
$v<v_{\mathrm{min}}+0.4(v_{\mathrm{max}}-v_{\mathrm{min}})$, where
$v_{\mathrm{max}}$ and $v_{\mathrm{min}}$ are the maximum and minimum values
of $v$, respectively, because when the reaction-diffusion wave propagates, the
oscillation phase of $v$ should shift from that of $u$ toward the reverse
direction of wave propagation.

Let $V_{\mathrm{wf}}$ be a vertex on the wavefront (black point in Fig.
15(a)). Since $V_{\mathrm{wf}}$ is on a bond of the original triangular
network graph, it faces two triangles and is adjacent to the two other
vertices of the wave edge. Let $L_{\mathrm{wf}}$ be a line connecting those
two vertices. We roughly estimate the direction of the line tangent to the
wavefront at $V_{\mathrm{wf}}$ as the direction of $L_{\mathrm{wf}}$. Next, we
calculate the cutting line of the tube along the $x$-axis through
$V_{\mathrm{wf}}$ using a greedy algorithm to advance in the $x$-axis
direction with a lookahead search with a depth of ten on a discrete graph
(black dotted line ($L_{\mathrm{cl0}}$) in Fig. 15(a)). Then, we determine the
shortest path from $V_{\mathrm{wf}}$ to $V_{\mathrm{wf}}^{\prime}$ (the
corresponding vertex to $V_{\mathrm{wf}}$ duplicated by cutting) on the graph
cut by $L_{\mathrm{cl0}}$. Thus, we obtain a path that starts from
$V_{\mathrm{wf}}$, circles the tube, and returns to itself (green circle line
($C_{\mathrm{wf}}$) in Fig. 15(a)). As well as $L_{\mathrm{cl0}}$, we obtain
two cutting lines $L_{\mathrm{cl1}}$ and $L_{\mathrm{cl2}}$, which are through
vertices where $C_{\mathrm{wf}}$ is divided by one quarter and three quarters,
respectively (yellow and red dotted lines in Fig. 15(a)). Then, we select
several vertices on $L_{\mathrm{cl1}}$ and $L_{\mathrm{cl2}}$
($V_{\mathrm{cl1}}$ and $V_{\mathrm{cl2}}$) and calculate the approximate
shortest path on the graph from $V_{\mathrm{cl1}}$ and $V_{\mathrm{cl2}}$ to
$V_{\mathrm{wf}}$ ($L_{\mathrm{geo1}}$ and $L_{\mathrm{geo2}}$, depicted by
cyan and orange lines in Fig. 15(a)). The algorithm for obtaining geodesic
lines is an iterative refinement of the discrete graph and Dijkstra’s
algorithm, the details of which can be found in Ref. Kanai and Suzuki (2000).
We find the best two geodesic lines, one from $L_{\mathrm{geo1}}$ and one from
$L_{\mathrm{geo2}}$, which are (i) smoothly connected at $V_{\mathrm{wf}}$,
and (ii) nearly parallel to $L_{\mathrm{wf}}$. These local geodesic lines are
calculated for approximately $20$ points on the wavefront, and the root mean
squared deviation between the geodesic line and the wavefront
($d_{\mathrm{RMSD}}$) is calculated as below:

$d_{\mathrm{RMSD}}=\sqrt{\frac{1}{N_{\mathrm{g}}}\sum_{i}\left|\bm{r}_{\mathrm{g},i}-\bm{r}_{\mathrm{wf,c}}\right|^{2}},$
(13)

where $N_{\mathrm{g}}$ is the number of vertices of the local geodesic line;
$\bm{r}_{\mathrm{g},i}$ is a position vector of the $i$th vertex of the local
geodesic line; and $\bm{r}_{\mathrm{wf,c}}$ is the position vector of the
vertex of the wavefront closest to the $i$th vertex of the local geodesic
line. The average of $d_{\mathrm{RMSD}}$ ($\langle d_{\mathrm{RMSD}}\rangle$)
is used to evaluate the wavefront undulation at each time step.

Examples of the calculation results are shown in Fig. 15(b) and (c). At
$t/\tau=2000$, the wavefront is not undulated and the local geodesic line is
aligned along the wavefront, whereas at $t/\tau=3000$, the wavefront is
undulated and a local geodesic line deviates slightly from the wavefront. (see
Fig. 15(b)). Thus, $\langle d_{\mathrm{RMSD}}\rangle$ has a small value at
$t/\tau\leq 2000$ but increases around $t/\tau=3000$ (see Fig. 15(c)). These
results demonstrate that the undulation of the wavefront can be evaluated
using $\langle d_{\mathrm{RMSD}}\rangle$.

## References

  * Zimmerberg and Kozlov (2006) J. Zimmerberg and M. M. Kozlov, Nat. Rev. Mol. Cell Biol. 7, 9 (2006).
  * McMahon and Gallop (2005) H. T. McMahon and J. L. Gallop, Nature 438, 590 (2005).
  * Suetsugu  _et al._ (2014) S. Suetsugu, S. Kurisu, and T. Takenawa, Physiol. Rev. 94, 1219 (2014).
  * Voeltz and Prinz (2007) G. K. Voeltz and W. A. Prinz, Nat. Rev. Mol. Cell Biol. 8, 258 (2007).
  * Baumgart  _et al._ (2011) T. Baumgart, B. R. Capraro, C. Zhu, and S. L. Das, Annu. Rev. Phys. Chem. 62, 483 (2011).
  * Peter  _et al._ (2004) B. J. Peter, H. M. Kent, I. G. Mills, Y. Vallis, P. J. G. Butler, P. R. Evans, and H. T. McMahon, Science 303, 495 (2004).
  * McMahon and Boucrot (2011) H. T. McMahon and E. Boucrot, Nat. Rev. Mol. Cell Biol. 12, 517 (2011).
  * Zhao  _et al._ (2017) W. Zhao, L. Hanson, H.-Y. Lou, M. Akamatsu, P. D. Chowdary, F. Santoro, J. R. Marks, A. Grassart, D. G. Drubin, Y. Cui, and B. Cui, Nat. Nanotech. 12, 750 (2017).
  * Sorre  _et al._ (2012) B. Sorre, A. Callan-Jones, J. Manzi, B. Goud, J. Prost, P. Bassereau, and A. Roux, Proc. Natl. Acad. Sci. USA 109, 173 (2012).
  * Has and Das (2021) C. Has and S. L. Das, Biochim. Biophys. Acta 1865, 129971 (2021).
  * Wu  _et al._ (2018) Z. Wu, M. Su, C. Tong, M. Wu, and J. Liu, Nat. Commun. 9, 136 (2018).
  * Wu and Liu (2021) M. Wu and J. Liu, Curr. Opin. Cell Biol. 68, 45 (2021).
  * Peleg  _et al._ (2011) B. Peleg, A. Disanza, G. Scita, and N. S. Gov, PLoS ONE 6, e18635 (2011).
  * Gov (2018) N. S. Gov, Phil. Trans. R. Soc. B 373, 20170115 (2018).
  * Weiner  _et al._ (2007) O. D. Weiner, W. A. Marganski, L. F. Wu, S. J. Altschuler, and M. W. Kirschner, PLoS Biol. 5, e221 (2007).
  * Taniguchi  _et al._ (2013) D. Taniguchi, S. Ishihara, T. Oonuki, M. Honda-Kitahara, K. Kaneko, and S. Sawai, Proc. Natl. Acad. Sci. USA 110, 5016 (2013).
  * Miao  _et al._ (2019) Y. Miao, S. Bhattacharya, T. Banerjee, B. Abubaker-Sharif, Y. Long, T. Inoue, P. A. Iglesias, and P. N. Devreotes, Mol. Syst. Biol. 15, e8585 (2019).
  * Tsujita  _et al._ (2015) K. Tsujita, T. Takenawa, and T. Itoh, Nat. Cell Biol. 17, 749 (2015).
  * Xiao  _et al._ (2017) S. Xiao, C. Tong, Y. Yang, and M. Wu, Dev. Cell 43, 493 (2017).
  * Godino  _et al._ (2019) E. Godino, J. N. López, D. Foschepoth, C. Cleij, A. Doerr, C. F. Castellá, and C. Danelon, Nat. Commun. 10, 4969 (2019).
  * Christ  _et al._ (2021) S. Christ, T. Litschel, P. Schwille, and R. Lipowsky, Soft Matter 17, 319 (2021).
  * Litschel  _et al._ (2018) T. Litschel, B. Ramm, R. Maas, M. Heymann, and P. Schwille, Angew. Chem. 57, 16286 (2018).
  * Xiong  _et al._ (2010) Y. Xiong, C.-H. Huang, P. A. Iglesias, and P. N. Devreotes, Proc. Natl. Acad. Sci. USA 107, 17079 (2010).
  * Huang  _et al._ (2013) C.-H. Huang, M. Tang, C. Shi, P. A. Iglesias, and P. N. Devreotes, Nat. Cell Biol. 15, 1307 (2013).
  * Zhan  _et al._ (2020) H. Zhan, S. Bhattacharya, H. Cai, P. A. Iglesias, C.-H. Huang, and P. N. Devreotes, Dev. Cell 54, 608 (2020).
  * Horibe  _et al._ (2019) K. Horibe, K.-i. Hironaka, K. Matsushita, and K. Fujimoto, Chaos 29, 093120 (2019).
  * Kneer  _et al._ (2014) F. Kneer, E. Schöll, and M. A. Dahlem, New J. Phys. 16, 053010 (2014).
  * Tamemoto and Noguchi (2020) N. Tamemoto and H. Noguchi, Sci. Rep. 10, 19582 (2020).
  * Tamemoto and Noguchi (2021) N. Tamemoto and H. Noguchi, Soft Matter 17, 6589 (2021).
  * FitzHugh (1961) R. FitzHugh, Biophys. J. 1, 445 (1961).
  * Nagumo  _et al._ (1962) J. Nagumo, S. Arimoto, and S. Yoshizawa, Proc. IRE 50, 2061 (1962).
  * Dimova (2014) R. Dimova, Adv. Colloid Interface Sci. 208, 225 (2014).
  * Roux (2013) A. Roux, Soft Matter 9, 6726 (2013).
  * Noguchi and Gompper (2005) H. Noguchi and G. Gompper, Phys. Rev. E 72, 011901 (2005).
  * Gommper and Kroll (2004) G. Gommper and D. M. Kroll, in  _Statistical Mechanics of Membranes and Surfaces_ , edited by D. Nelson, T. Piran, and S. Weinberg (World Scientific, Singapore, 2004) 2nd ed., pp. 359–426.
  * Noguchi (2021a) H. Noguchi, Phys. Rev. E 104, 014410 (2021a).
  * Noguchi (2022) H. Noguchi, Int. J. Mod. Phys. B36, 2230002 (2022).
  * Noguchi (2021b) H. Noguchi, Soft Matter 17, 10469 (2021b).
  * Ou-Yang and Helfrich (1989) Z. C. Ou-Yang and W. Helfrich, Phys. Rev. A 39, 5280 (1989).
  * Kenmotsu (2003) K. Kenmotsu,  _Surfaces with constant mean curvature_ (American Mathematical Society, Providence, R.I., 2003).
  * Naito  _et al._ (1995) H. Naito, M. Okuda, and O.-Y. Zhong-can, Phys. Rev. Lett. 74, 4345 (1995).
  * Tsafrir  _et al._ (2001) I. Tsafrir, D. Sagi, T. Arzi, M.-A. Guedeau-Boudeville, V. Frette, D. Kandel, and J. Stavans, Phys. Rev. Lett. 86, 1138 (2001).
  * Campelo and Hernández-Machado (2007) F. Campelo and A. Hernández-Machado, Phys. Rev. Lett. 99, 088101 (2007).
  * Tozzi  _et al._ (2019) C. Tozzi, N. Walani, and M. Arroyo, New J. Phys. 21, 093004 (2019).
  * Liu  _et al._ (2016) X. Liu, F. Tian, T. Yue, X. Zhang, and C. Zhong, Soft Matter 12, 9077 (2016).
  * den Otter  _et al._ (2003) W. K. den Otter, S. A. Shkulipa, and W. J. Briels, J. Chem. Phys. 119, 2363 (2003).
  * Matsuura-Tokita  _et al._ (2006) K. Matsuura-Tokita, M. Takeuchi, A. Ichihara, K. Mikuriya, and A. Nakano, Nature 441, 1007 (2006).
  * Papanikou and Glick (2014) E. Papanikou and B. S. Glick, Curr. Opin. Cell. Biol. 29, 74 (2014).
  * Kulkarni-Gosavi  _et al._ (2019) P. Kulkarni-Gosavi, C. Makhoul, and P. A. Gleeson, FEBS Letters 593, 2289 (2019).
  * Weigel  _et al._ (2021) A. V. Weigel, C.-L. Chang, G. Shtengel, C. S. Xu, D. P. Hoffman, M. Freeman, N. Iyer, J. Aaron, S. Khuon, J. Bogovic, W. Qiu, H. F. Hess, and J. Lippincott-Schwartz, Cell 184, 2412 (2021).
  * Zucker and Kozlov (2022) B. Zucker and M. M. Kozlov, Proc. Natl. Acad. Sci. USA 119, e2116142119 (2022).
  * Yang  _et al._ (2017) Y. Yang, D. Xiong, A. Pipathsouk, O. D. Weiner, and M. Wu, Dev. Cell 43, 507 (2017).
  * Mim and Unger (2012) C. Mim and V. M. Unger, Trends Biochem. Sci. 37, 526 (2012).
  * Kabaso  _et al._ (2011) D. Kabaso, E. Gongadze, P. Elter, U. van Rienen, J. Gimsa, V. Kralj-Iglič, and A. Iglič, Mini Rev. Med. Chem. 11, 272 (2011).
  * Noguchi (2019) H. Noguchi, Sci. Rep. 9, 11721 (2019).
  * Hurley  _et al._ (2010) J. H. Hurley, E. Boura, L.-A. Carlson, and B. Różycki, Cell 143, 875 (2010).
  * Foret (2014) L. Foret, Eur. Phys. J. E 37, 42 (2014).
  * Goutaland  _et al._ (2021) Q. Goutaland, F. van Wijland, J.-B. Fournier, and H. Noguchi, Soft Matter 17, 5560 (2021).
  * Noguchi  _et al._ (2022) H. Noguchi, C. Tozzi, and M. Arroyo, Soft Matter 18, 3384 (2022).
  * Kanai and Suzuki (2000) T. Kanai and H. Suzuki, in  _Proceedings Geometric Modeling and Processing 2000._ (2000) pp. 241–250.