*Perspective* **Why a Large-Scale Mode Can Be Essential for Understanding Intracellular Actin Waves**

#### **Carsten Beta 1, Nir S. Gov <sup>2</sup> and Arik Yochelis 3,4,\***


Received: 16 April 2020; Accepted: 18 June 2020; Published: 23 June 2020

**Abstract:** During the last decade, intracellular actin waves have attracted much attention due to their essential role in various cellular functions, ranging from motility to cytokinesis. Experimental methods have advanced significantly and can capture the dynamics of actin waves over a large range of spatio-temporal scales. However, the corresponding coarse-grained theory mostly avoids the full complexity of this multi-scale phenomenon. In this perspective, we focus on a minimal continuum model of activator–inhibitor type and highlight the qualitative role of mass conservation, which is typically overlooked. Specifically, our interest is to connect between the mathematical mechanisms of pattern formation in the presence of a large-scale mode, due to mass conservation, and distinct behaviors of actin waves.

**Keywords:** nonlinear waves; actin polymerization; bifurcation theory; mass conservation; spatial localization; pattern formation; activator–inhibitor models

#### **1. Introduction**

Biological pattern formation refers to the emergence of complex spatiotemporal variations in living systems that are typically far from thermodynamic equilibrium [1–3]. Even though these systems can differ in composition and scales, they share many similarities and generic phenomena that are observed in a wide variety of natural settings, such as stationary periodic patterns of pigments on animal skins, spiral waves in biological cells and cardiac arrhythmia, or swarming phenomena in bacterial colonies and in flocks of birds or fish. The theoretical study of biological pattern formation can be roughly divided into two time periods: (i) The second half of the twentieth century, following the seminal works by Turing on morphogenesis [4] and by Hodgkin and Huxley (HH) on action potentials in the giant squid axon [5], and (ii) the beginning of the twenty-first century, where an ever-increasing amount of quantitative biological data provided the basis for more detailed mechanistic models of biological systems.

During the first period, theoretical studies were largely limited to a few prototypical reaction–diffusion (RD) or activator–inhibitor (AI) models [6], such as the FitzHugh–Nagumo (FHN) [7,8], Gierer–Meinhardt [9], and Keller–Segel equations [10]. Based on the relative simplicity of these models, e.g., the FHN model as compared to the HH equations, and their relation with models of inanimate matter, e.g., the Swift-Hohenberg model of thermal fluid convection [11] and the Gray–Scott model of chemical reactions [12–14], several pattern formation methodologies have been advanced [1], such as weakly nonlinear and singular-perturbation methods. These models provided deep insights into universal aspects of pattern formation phenomena and generic relations to applications were

substantiated, such as frequency locking and spiral waves in the cardiac system. The second time period has manifested a gradual shift of research interests towards specific detailed biological and medical systems [15,16], including, for example, micron-scale intracellular waves, the development of tissues and organs, sound discrimination in the auditory system, and pathologies such as cancer metastasis. In particular, systems in these contexts are generally described by elaborate, system-specific models that are less amenable to mathematical analysis than earlier toy models of pattern formation.

Consequently, despite the common pattern formation thread that connects these different biological and medical applications, it became difficult to navigate through the vast number of distinct models and approaches, particularly in cases where technical jargon makes it difficult to adopt cross-disciplinary integration between different communities including biophysics, computational biology, mathematical biology, biochemistry, dynamical systems, and numerical analysis. Even though the formulation of complete models for biological systems is currently unrealistic, uncovering partial mechanisms that drive pattern formation phenomena remains of utmost importance for understanding functional aspects of living systems and for developing technological and medical applications, such as drugs or implants. Moreover, mechanistic studies of pattern forming systems are also fertile sources of new mathematical questions that advance the development of analytical and numerical methods [2,17–25], which, in turn, contribute new insights into the original applications.

In this perspective, we will focus on intracellular actin waves (IAW), a topic that recently gained much interest not only in the biological context but also as an inspiring showcase of active matter. More specifically, we are interested in IAW that are affected by a large-scale mode—a situation that arises due to the conservation of actin monomers (over the time-scale of the IAW phenomenon). We note that phenomena such as Ca2<sup>+</sup> waves are, in general, beyond the scope of this perspective as they involve the transport of ions between the cell interior and the extracellular space (which acts as an infinite reservoir) [26–28], unless conservation can be accounted for [29]. Moreover, we emphasize that we aim to provide a perspective and not a comprehensive review, as such reviews are already available, see, e.g., in [30–36].

The perspective is organized as follows. In Section 2, we introduce the rich phenomenology of IAW and the modeling aspects that are associated with mass conservation. Then, in Section 3 we present the theoretical aspects of a large scale mode in the context of physicochemical settings and also indicate its significance to IAW as a consequence of the conservation of actin monomers on the time scales at which many actin-driven processes operate. Finally, in Section 4 we discuss why incorporation of mass conservation is a plausible qualitative step in unfolding the robustness of IAW mechanisms, and in Section 5 we conclude by emphasizing the theoretical strategies for modeling and control of wave persistence as a potential roadmap toward applications in synthetic biology.

#### **2. Intracellular Actin Waves**

The functions of many cells are tied to their ability to dynamically change their shape, mostly via the spatiotemporal organization of their actin cytoskeleton. Examples of this include diverse cell types, such as human neutrophils, fish keratocytes, or the social amoeba *Dictyostelium discoideum*. Among the most prominent dynamical patterns in the actin cytoskeleton are IAW, which have attracted much attention over the past decade [31]. These waves are assumed to play a role in several essential cellular functions, among them cell locomotion, cytokinesis, and phagocytic uptake of extracellular matter. Many competing models at different levels of complexity have been developed to describe cortical actin waves, mostly relying on coupled nonlinear AI equations. Even though intracellular actin waves involve a large number of interacting molecular species as well as multiple local and global interactions, prototypical AI models have been shown to capture many features of the overall dynamics. However, important effects due to mass conservation constraints have been hitherto largely neglected.

#### *2.1. Phenomenology from Experiments*

Actin waves are characterized by propagating cytoskeletal regions that are enriched in filamentous actin and actin-related proteins. Depending on the cell type, IAW may differ in their biochemical composition and dynamics, including different wave morphologies and propagation speeds. One of the earliest examples of IAW was reported from cultured neurons that show the propagation of fin-like actin-filled membrane protrusions along their axon [37]. They were found to depend on actin polymerization and have been associated with neural polarization [38,39]. Similar fin-like actin waves also emerge in non-neural cell types when cultured on thin fibers [40]. Moreover, adherent cells that are attached to flat substrates may display traveling wave-like protrusions of their cell shape. They are particularly prominent when moving laterally along the cell border, such as in mouse embryonic fibroblasts [41] or at the leading edge of fish keratocytes [42].

Traveling actin waves have also been observed at the dorsal and ventral sides of adherent cells. In neutrophils, small dynamic wave fragments emerge that organize cell polarity and leading edge formation [43]. Larger ring-shaped waves were found to travel across the substrate-attached bottom membrane of *D. discoideum* cells [44]. They enclose a region that is structurally distinct from the cortical area outside the actin ring [45,46], and their dynamics often shows rotating spiral cores and mutual annihilation upon collision [47,48], but they could not be initiated by external receptor stimuli [49]. While understanding the rich dynamics of IAW is challenging on its own right, there are prominent applications and functional properties that stimulate further studies of IAW in different contexts:


**Figure 1.** Examples of intracellular actin waves. (**A**) Actin wave nucleation and propagation in a migrating immature dendritic cell. Red arrows indicate the origin of the wave, scale bar 20 μm. (**B**) Overlay of contours of representative actin waves shown in panel (**A**) during propagation. (**C**) Wave-mediated binary cytofission in a *Dictyostelium discoideum* cell, scale bar 10 μm. An actin wave in a cell with two nuclei becomes unstable and splits into two independent segments that move in opposite directions and induce a cytofission event. Panels (**A**) and (**B**) are reproduced from [50] and panel (**C**) is reproduced from [53]. Copyright 2020 National Academy of Sciences.

Despite intense studies over the past years, the molecular details of IAW mechanisms remain largely unclear and most likely vary between different cell types.

#### *2.2. Modeling Approaches of Actin Waves*

Following the numerous experimental observations of IAW in different cell types and during different cellular functions, many model equations have been proposed to describe this phenomenon. Here, we will briefly describe the main types and features of theoretical models that have been employed while referring the reader to [31,34,35,63,64] for more details.

The growth of the cortical actin network within IAW is a complex dynamical process that involves many components that perform a coordinated set of functions, giving rise to the formation of a three-dimensional network of actin filaments that propagates along the cell membrane. This process involves the activation of actin-associated proteins some of them membrane bound, that initiate the nucleation of actin polymerization, branching of actin filaments, cross-linking, and bundling, as well as severing and depolymerization. There are very few theoretical models that attempt to give a molecular-scale description of the IAW phenomenon, where all of these processes are described. One example for such a model that attempts to describe the waves at the scale of the individual actin filaments is given in [65]. Although providing detailed pictures of the actin network, it is difficult and time-consuming to use such modeling to extract understanding regarding the large-scale dynamics of the IAW. Such modeling efforts could in the future include more molecular components [66,67], on larger length and time scales, and provide a platform for theoretical advances in this field, that works in conjunction with filament-scale experimental data [46].

As the IAW have characteristic spatial scales in the range of hundreds of nanometers, propagate over tens of microns, and persist over hours, it is natural to describe them using coarse-grained models that avoid prescribing the molecular-scale details of the actin network. As will be shown, many of these models agree with some qualitative or even quantitative features of the observed IAW in cells. It is therefore difficult at present to reach a consensus regarding the validity of these models. Comparisons between such models is also complicated as they often include different components, and it is unclear if and which one of those components play a fundamental role in the emergence of IAW or can be neglected, otherwise.

Among the coarse-grained models, we can find a small class of models that contain biophysical elements, such as forces and/or the membrane shape, which play a key role in the mechanism that drives the propagation of the IAW. One example is well demonstrated by Gholami et al. [68], who show that the dynamics of the actin polymerization/depolymerization drive the oscillatory propagation of waves. When actin filaments polymerize against the cell membrane, they exert a protrusive pressure on the membrane, which pushes the membrane forward and the actin network backwards. The interplay between the rate of actin polymerization and the rate at which the actin filaments are cross-linked into a stable gel-like network, determines if the cortical actin is stable or whether it exhibits an unstable oscillatory regime.

Another group of biophysics-based models contains curved membrane proteins that nucleate the cortical actin polymerization [69–73]. In these models, the curved proteins flow/adsorb to the membrane regions that have a curvature similar to their intrinsic shape, and their concentration is therefore affected by the membrane deformations that are induced by the forces exerted by the actin cytoskeleton. These forces include the protrusive force of actin polymerization, as well as contractile forces due to myosin-II mediated contractility. Recently, also models combining an RD kinetics coupled to mechanical properties through the impact of curved actin nucleators and/or membrane shape and tension were introduced [42,74]. Other models combine the RD dynamics with a physical effect, such that the directed or random lateral actin polymerization can physically drive the treadmilling of the IAW components along the membrane [75]. The advantage of the biophysical class of models is that they can naturally account for the observed effects of physical parameters on the IAW, such as membrane tension [68,74] or the contractile forces of myosin-II motors [76].

In many cases, however, RD equations that include both positive and negative feedback loops, are sufficient to demonstrate the formation of propagating waves, fronts, or pulses. These models exhibit different levels of complexity and different numbers of components. In the simplest cases, generic activator–inhibitor models of FHN-type were proposed. In particular, they were used together with a local-excitation, global-inhibition (LEGI) mechanism to account for the response of the receptor-mediated signaling pathway and the downstream actin cytoskeleton to external cues [77–79]. Other basic RD models describe the actin dynamics, including the monomeric and filamentous species, and one form of an actin activator, using the filamentous actin itself as a source of negative [65] or positive [80] feedback. More complex models include different numbers of activators of actin polymerization, inhibitors, and their complex network of interactions [55,81,82]. Yet, in general, RD equations are not subjected to the conservation of mass constraint although often some of the components are conserved, for example, when they represent two different forms of the same protein [83]. In other cases, the actin is conserved as it is converted from monomeric to filamentous forms and back, see for example [58,84]. In what follows, we address the qualitative role of conservation, which is reflected by the existence of a large scale mode, on the dynamics of IAW, using as much as possible generic principles, i.e., extracting conclusions that are qualitatively independent of the specific molecular details that are included in the model.

#### **3. Actin Dynamics as a Constrained Continuous Medium: Implications and Applications**

The phenomenology of dissipative waves can be conveniently demonstrated through a dynamical systems approach via prototypical models, such as FHN. As summarized above, many variants of such activator–inhibitor models have been used to describe different aspects of cytoskeletal dynamics and in particular the formation of actin waves. Although these are heuristic models, they are analytically tractable and thus allow for fundamental insights into spatiotemporal behavior, which cannot be obtained through the analysis of more realistic multivariable equation sets. Propagating waves are traditionally classified into three universality classes [1,2,6,85]:


While the mathematical mechanisms are distinct, the emerging patterns can show similar characteristics, for example, all classes may display the formation of spiral waves [2,85]. Consequently, comparisons to experimental observations can often only be qualitative, making insights uncertain. Moreover, it is not always clear whether the simplified models comprise the minimal set of qualitative ingredients, e.g., interactions (local vs. non-local), spatial coupling, essential degrees of freedom and feedback loops, finite domain effects, or existence of conserved observable(s). In a broader context, IAW can be classified as AI type media [31,86], although unlike the typical RD media, the number of actin monomers is conserved over the time scales of wave dynamics. As such, mass conservation is an inherent constraint of the modeling framework [35,58,86,87], which is generically reflected by coupling to a large scale mode in the dispersion relation, as illustrated in Figure 2.

**Figure 2.** Schematic representation of a dispersion relation obtained from infinitesimal periodic perturbations, proportional to exp (*σt* + *ikx*), about a uniform steady state; *Re*[*σ*] is the perturbation growth rate and *k* its wavenumber. The right-hand part of the dispersion relation represents the onset of an instability of a finite wavenumber type (often also referred to as Turing instability), while the left-hand part reflects a conserved quantity and stays always neutral; both parts are model-independent. The curves may connect as typically occurs in systems such as (1) or belong to different curves, such as for (5). The imaginary part of *σ* corresponds to stationary nonuniform patterns if zero, and otherwise describes time-dependent solutions.

#### *3.1. Conservation in Physicochemical Systems*

It is convenient to first consider total conservation of an observable, described by the continuity equation

$$\frac{\partial u}{\partial t} = \nabla \cdot \left[ M(u) \, \nabla \frac{\delta F(u)}{\delta u} \right],\tag{1}$$

where *u* is a scalar field, *M* is a mobility function, and *F* is a free energy. If the free energy contains an intrinsic length scale, like in the phase field crystal model or in wetting, stationary periodic and localized patterns may emerge [88–95]. The mutual aspect is coupling between the large-scale mode (*k* = 0), which is model independent and remains always neutral due to conservation (also known as the Goldstone mode), and the pattern forming instability of finite wavenumber (Turing) type [96], as shown in Figure 2. The impact of the conserved quantity has been analyzed mostly via a weakly nonlinear reduction to a set of two amplitude equations: One is the complex Ginzburg–Landau equation for the finite wave-length mode, while the other is for the neutral large scale mode [97–106]. Both super- and sub-critical bifurcation types have been studied, and showed that indeed inclusion of the large-scale mode may qualitatively change the nature of the solution in terms of organization and

stability [107,108] (and references therein). For example, in the absence of the large-scale mode, spatially localized solutions form in coexistence with a periodic (Turing-type) solution and are organized in a vertical snakes-and-ladders homoclinic structure. In the presence of a large-scale mode, these solutions can also form outside the existence region of periodic solutions and only partially overlap, i.e., the homoclinic snaking structure becomes slanted [101,109].

In fact, similar asymptotic intuition and analysis methods apply also if the observable has a velocity-like behavior (often Galilean invariance) [110], obeying the symmetry *x* → −*x* and *u* → −*u*. Such behavior arises in systems that are being driven out of equilibrium, such as convection [111–114], propagation of flames [115,116], surface waves [117–119], and electro-diffusion in ion channels [120,121]. In such cases, leading order approximations show that the dynamics can still be enslaved to an oscillatory (Hopf) finite wavenumber mode and a large scale mode [122–126]. While many fundamental advances have been made in understanding the coupling between the complex Ginzburg–Landau equation and the large scale mode, e.g., in terms of stability of periodic and solitary waves in one space dimension and dynamics of spiral waves in two-dimensional systems, several pattern formation issues remain open [127]. Consequently, as over the time scales on which IAW occur the system is far from equilibrium, it is natural to assume that a large scale mode due to mass conservation alters the pattern formation mechanism, even without explicit flux conservation.

#### *3.2. Activator–Inhibitor Patterns with Conservation*

In general, AI systems are modeled in a similar fashion as chemical reactions [6,15,16,36,128–131], which are not limited by supply of new substrates into the reactor:

$$\frac{\partial \mu}{\partial t} = -f(u, v) + D\_{\mu} \nabla^{2} u\_{\prime} \tag{2a}$$

$$\frac{\partial v}{\partial t} = -g(u, v) + D\_v \nabla^2 v,\tag{2b}$$

where *u* is the activator that typically contains an autocatalytic or enzymatic kinetic term and a diffusion constant *Du*, and *v* is an inhibitor that rapidly diffuses with a diffusion constant *Dv*, where typically *Dv Du*. Note that we do not discuss transport by cross-diffusion here. As intracellular processes often take place on very different time scales, effective mass conservation may arise, for example, in cases where protein synthesis and/or degradation occurs much slower than a particular biochemical reaction of interest. Conservation in AI models is associated with a local conservation of mass

$$\int\_{\Omega} \left[ u(\mathbf{x}, t) + v(\mathbf{x}, t) \right] d\mathbf{x} = \text{constant},\tag{3}$$

where Ω is the physical domain, or by defining in (2)

$$
\mathcal{g}(u,v) = -f(u,v). \tag{4}
$$

Linear stability analysis about uniform solutions leads to dispersion relations that contain the persistent neutral (large scale) mode, as shown in Figure 2. As in the case of Equation (1), also Equation (4) supports multiplicity of uniform solutions since *u* depends on an arbitrarily chosen value of *v* (or vise versa), and this degenerate degree of freedom appears as the *k* = 0 mode. This constraint plays effectively the role of a chemical potential. However, in the pattern forming case, where an additional bifurcation is present (e.g., a Turing bifurcation), the dispersion relations may contain both the neutral mode at *k* = 0 and another at a finite wavenumber. In this formulation, models for cell polarity [132] and molecular motors [133] had inspired several mathematical works in the context of existence and emergence of stationary [134–138] and time-dependent [139–141] patterns.

However, as has been described in Section 2.2, IAW are multicomponent processes and essentially comprise involve a large number activators and inhibitors. Moreover, in such an AI network, not all

the components obey conservation [86,142], namely, to Equations (2) and (4) can be added at least one additional non-conserved observable *w*,

$$\frac{\partial u}{\partial t}^{\prime} = -f(u, v, w) + D\_{\text{ul}} \nabla^2 u\_{\prime} \tag{5a}$$

$$\frac{\partial \upsilon}{\partial t} = -f(u, v, w) + D\_{\upsilon} \nabla^2 \upsilon,\tag{5b}$$

$$\frac{\partial w}{\partial t} = -h(u, v, w) + D\_w \nabla^2 w,\tag{5c}$$

where *h* can be either a linear or a nonlinear functional and essentially does not have to include transport of *w* via diffusion; these details are naturally determined by the characteristics of the biological system. Equation (5) thus reflects only a partial conservation and has been employed to study the emergence of IAW in the context of CDR [58], where a variety of complex behaviors have been observed experimentally, ranging from distinct types of propagating fronts to apparently chaotic spiral waves.

#### **4. Discussion and Example**

The complex pattern formation exhibited by CDR raises the question about the modeling strategy, specifically, with respect to the minimal set of equations and the necessity of a conserved quantity. As has already been indicated in Section 2.2, there are many ways to model IAW but all of them are prone to subjective interpretations.

In the absence of a clear physical intuition, as IAW are far from equilibrium phenomena, dynamical systems offer an efficient platform for creating an appropriate qualitative framework. More specifically, the study of bifurcations may provide the minimal qualitative set of constraints, exactly as phase-transitions allow us to classify many types of physical phenomena. On the other hand, bifurcation analysis can also be a tedious task as there may be many local and global bifurcations that coexist in a given parameter range (as an example we refer the reader to a systematic extension of excitable media by Champneys et al. [143]). Nevertheless, utilizing recent advances in nonlinear perturbations [83,144] and numerical path continuation methods [145–149] it might be possible to navigate between coexisting bifurcations and a multiplicity of emerging stable and unstable solutions [144,150]. Next, we turn to conservation and ask whether it may prescribe a fundamental and robust qualitative change, as compared to typical local RD modeling in the absence of conserved quantities. To exemplify this case, we exploit a reduced CDR model (of Equations (5) type), which has been used to examine solitary wave collisions in the context of IAW [151]. In the reduced CDR model, the conserved AI system of Equations (5) is replaced by the conservation of the actin monomers, as they are converted from the monomeric to the filamentous form (and back) while the IAW propagates.

Observation of solitary waves dates back to John S. Russell (1834), yet only after the work in Zabusky and Kruskal [152] were solitary waves distinguished by their collision properties [153,154]: *solitons* if after a collision of two pulses, two pulses emerge (particle-like identity), and *dissipative solitons* or *excitable pulses* if they are annihilated. Solitons are often being discussed in the context of conservative media, which mathematically means exploiting the integrable nature of the governing model equations [107, 155], while excitable pulses often arise in RD type systems. Although collisions of solitons may involve high spatiotemporal complexity, the outcome of two colliding solitons remains unchanged (i.e., elastic particle-like dynamics) [156,157]. On the other hand, the annihilation of excitable pulses after the collision is recognized as paramount for electrophysiological function, i.e., it would be impossible to maintain directionality, and thus rhythmic behavior, under the reflection of action potentials [158].

Importantly, collision of pulses implies merging of the pulses in space, i.e, through the formation of a *collision zone*. This behaviour is distinct from the *interaction* between excitable pulses that is due to repulsion and can exhibit dynamics that may resemble a solitonic behavior [159,160]. Moreover, more complex scattering scenarios have been observed in generic RD models such as, for example, the Gray–Scott model [161,162]. Note that there exists a vast literature on the latter topic that we do not intend to review in total here. Taken together, the distinction between solitons and excitable pulses is important for numerous applications.

Yochelis et al. [151] showed that the minimal IAW model in one space dimension, in the setting of Equations (5), may indeed support rich and robust spatiotemporal dynamics following pulse collisions, in contrast to IAW models which do not contain explicit mass conservation [28,55,81,163,164]: annihilation, reflection, and "birth" of new pulses after reflection, as shown in Figure 3. In a broader RD context, where similar aspects have been also observed, these dynamics do not require special properties, such as non-locality [165–168], cross-diffusion [169], and heterogeneity [170–172]. Moreover, the phenomenon is robust and occurs over a wide range of parameter values, whereas for a typical RD model without mass conservation, such as FHN, somewhat similar dynamics of propagating pulses are observed only in a narrow range near the onset of an oscillatory Hopf bifurcation about a uniform steady state [173,174]. The distinction between the FHN model and a system of Equations (5) type can be elaborated by geometrical intuition, as pulses are of large amplitude and thus cannot be unfolded using weakly nonlinear analysis such as in Section 3. Argentina et al. [173] showed that in the FHN model a manifold construction about the collision state of two pulses ("collision droplet", Figure 4A) can explain why a Hopf bifurcation may impact the collision zone and thus generate crossover of pulses (soliton-like behavior). A similar geometric picture shows that mass conservation in Equations (5) changes the nature of the collision zone by addition of a generic two-dimensional neutral manifold (Figure 4B), relating the pulse crossover behavior to a localized unstable mode and does not require any Hopf bifurcation of the uniform state [173,174]. In other words, for the colliding pulses to avoid annihilation, there has to be a mechanism for recovery—a spontaneous regrowth of the fields after the collision. In the FHN model [173], the proximity to the oscillatory onset can re-initiate the pulses. In the case of actin conservation, the colliding pulses first disintegrate the polymerized actin, thereby releasing a large local pool of monomers. If these monomers do not diffuse too fast, they are available to re-initiate the pulses by polymerization. For more details, we refer the reader to the work in [151].

**Figure 3.** Space–time plots showing (from left to right) annihilation, reflection/crossover, and "birth" of new pulses following collision (a behavior that resembles backfiring), respectively, as obtained from direct numerical integration of the minimal CDR model equations [151] that have the same structure as Equations (5). No-flux boundary conditions were used. From left to right the amount of actin monomers increases (see details in [151]). The dark shaded color indicates higher values of filamentous actin in the IAW. Reprinted figure with permission from the work in [151] Copyright 2020 by the American Physical Society.

**Figure 4.** Excitable solitons, geometric schematics of the dynamics during collision of two pulses. (**A**) FitzHugh–Nagumo model and (**B**) an reaction–diffusion model with mass conservation, of Equation (5) type. (**A**) Reprinted from Publication [173], with permission from Elsevier and (**B**) from [151], Copyright 2020 by the American Physical Society.

Note that the nucleation of new pulses after a collision should not be confused with the well-known scenario of backfiring, an instability that appears when a localized propagating pulse becomes unstable and splits into two new counterpropagating pulses that, upon collision, annihilate [175]. Backfiring has been observed in a wide range of model systems [166,176,177], and also in recent experiments of CO electrooxidation on Pt [178]. However, in contrast to backfiring, the nucleation of new pulses that we addressed here and that is shown in Figure 3 always requires a preceding collision event and thus has to be distinguished from the classical backfiring scenario.

#### **5. Conclusions**

The case of the reduced CDR model discussed above provides a glimpse of the profound impact of mass conservation on the dynamics. In conventional FHN-type AI models without mass conservation, colliding pulses typically annihilate upon collision. Here, a soliton-like crossover occurs only under special conditions, e.g., near a Hopf point, and thus requires fine-tuning of the parameters. In contrast, if mass conservation is taken into account, propagating pulses robustly exhibit rich collision scenarios over a wide range of parameters, including crossover and formation of new pulses following a collision. Even though this has only been demonstrated for a simple toy model, the universal nature of the underlying bifurcations suggests that similar behavior will be observed also in more detailed, high-dimensional models of IAW, provided that mass conservation is included, e.g., for mechanochemical waves under the conservation of calcium [29].

The impact of mass conservation on pattern formation in biological systems has recently attracted increasing attention, in particular in the context of well-controlled, confined systems such as the bacterial Min protein oscillator [87]. However, many biological systems involve multiple components, not all of which are conserved, so that the consequences of strict mass conservation as implied by Equations (2) and (4) are often relaxed and require a more general view. This is provided, in the simplest case, by adding a third dynamical variable to the system that is coupled to the conserved quantities but does not obey mass conservation itself, see Equations (5). It demonstrates that a large scale mode is a key feature that mass conservation introduces to the system and that triggers specific dynamical properties, such as soliton-like crossover of pulses and the collision-induced birth of new pulses in a wide range of parameters. Equations (5), and its resulting dynamics, can serve as motivation for future studies of the synthesis between classical AI models and models with complete mass conservation (such as those used in the context of the Min and Par systems [179]).

Similar to neural systems, where annihilation of colliding pulses is essential to maintain directionality of information transport, we conjecture that also in the case of IAW, the crossover of colliding pulses, which is favored due to the mass conservation constraint, plays an important functional role. This may be particularly true when sustained wave activity is a key requirement for proper cell functions, as, for example, in cases where cell locomotion or nutrient uptake depend on IAW (see Section 2.1). For traditional excitable pulses that annihilate upon collision, wave activity is likely to get extinguished regularly, thus hampering cellular activities that rely on persistent IAW. In contrast, soliton-like crossover and collision-induced nucleation of new pulses that are robust properties of a mass-conserved system may ensure prolonged wave activity even in the absence of actively triggered pulse nucleation or local heterogeneities that may serve as pacemakers. Moreover, cells may also actively exploit shifts between parameter regimes of pulse annihilation and soliton-like behavior to control their level of IAW activity, as shown in Figure 3.

Finally, the study of simplified models to elucidate generic properties of IAW patterns may also prove useful for the future design of synthetic cellular systems. A current focus of bottom-up approaches in synthetic biology is to introduce artificial cytoskeletal structures into membrane vesicles, thus assembling the essential building blocks of a primitive cell [180,181]. The next logical step along this line of research will be to endow the artificial cytoskeletal components with the simple pattern forming properties that may ultimately serve as a basis for essential cellular functions, such as motility and cytokinesis. This requires a thorough understanding of the key properties that are necessary to reconstitute the desired wave patterns in a minimal model system. We thus expect that an understanding of the essential bifurcations and instabilities that govern the dynamics of IAW will provide a useful guideline for the future design of artificial cell cortices.

**Author Contributions:** Conceptualization and writing—original draft preparation, A.Y.; writing—review and editing, C.B., N.S.G., A.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** N.S.G. is the incumbent of the Lee and William Abramowitz Professorial Chair of Biophysics and thanks the support by the Israel Science Foundation (Grant No. 1459/17).

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
