Traveling waves of an FKPP-type model for self-organized growth.

We consider a reaction-diffusion system of densities of two types of particles, introduced by Hannezo et al. (Cell 171(1):242-255.e27, 2017). It is a simple model for a growth process: active, branching particles form the growing boundary layer of an otherwise static tissue, represented by inactive particles. The active particles diffuse, branch and become irreversibly inactive upon collision with a particle of arbitrary type. In absence of active particles, this system is in a steady state, without any a priori restriction on the amount of remaining inactive particles. Thus, while related to the well-studied FKPP-equation, this system features a game-changing continuum of steady state solutions, where each corresponds to a possible outcome of the growth process. However, simulations indicate that this system self-organizes: traveling fronts with fixed shape arise under a wide range of initial data. In the present work, we describe all positive and bounded traveling wave solutions, and obtain necessary and sufficient conditions for their existence. We find a surprisingly simple symmetry in the pairs of steady states which are joined via heteroclinic wave orbits. Our approach is constructive: we first prove the existence of almost constant solutions and then extend our results via a continuity argument along the continuum of limiting points.

Motivation and result

The mechanics of tissue-growth have drawn the attention of the scientific community. A central question is, how the cells are organized, how they react to and communicate with their environment on the microscopic level, and how their behavior during the growth phase gives rise to distinct macroscopic structures. Mathematical models can help to understand these processes and works regarding organoids, wound healing or tumor growth are abundant (Montes-Olivas et al. 2019; Mammoto and Ingber 2010; Falco et al. 2021; d’Alessandro et al. 2021; Jain et al. 2021). However, for most of these models, our numerical skills far predominate the possibility to analyze them rigorously. Hence, for understanding the basic mechanics of the underlying biological processes, the need for simplified models arises.

Especially when studying spatiotemporal effects and macroscopic pattern formation, reaction–diffusion systems and their traveling waves have proven insightful. One of the oldest and most studied models is the FKPP-equation (Fisher 1937; Kolmogorov et al. 1937), describing the advance of an advantageous population. The arise of more complex spatial patterns due to the instability of a homogeneous state was first described in Turings groundbreaking paper The Chemical basis of Morphogenesis (Turing 1952). More recently, systems of Keller-Segel type have been studied extensively, where growth, movement and self-organization of a population are driven by chemotactic guidance (Keller and Segel 1971; Perthame 2004; Painter 2019). Broad introductions to mathematical modeling of pattern formation in developmental biology have been written by Painter (2019) and Othmer et al. (2009), among others.

The group of Hannezo et al. proposed A Unifying Theory of Branching Morphogenesis in epithelial tissues (Hannezo et al. 2017). They introduced a stochastic model, related to branching and annihilating random walks (Cardy and Täuber 1996). In this model, a branched structure is represented by a network. This network undergoes stochastic growth dynamics, where each branch of the network grows independently from the others and follows a set of simple, local rules. At its tip, each branch elongates or splits up at certain rates and these tips are called active. When an active tip comes too close to a different branch, it irreversibly ceases any activity and becomes inactive. The numerical results of Hannezo et al. reveal that this stochastic growth process self-organizes: the active tips concentrate at the boundary of the network and form a rather sharp layer of growth. The center of the network is static and - rather surprisingly - exhibits a homogenous geometry, in particular a constant density of branches. Remarkably, as mentioned by the authors, this model self-organizes without any signaling gradients. Even a directional bias of the branches can be achieved, as the result of an appropriate spatial boundary. Moreover, the authors observed that their simulations were in good agreement with biological data from mammary glands, kidneys and the human prostate (Hannezo et al. 2017).

To study their model analytically, Hannezo et al. proposed the following system, which corresponds to the diffusive limit of the above stochastic dynamics. We restrict ourselves to the one-dimensional case. Due to a simple linear rescaling (“Appendix C”), we only need to consider the normalized reaction–diffusion systemHere, are the densities of active particles and inactive particles. The diffusion term describes the movement of the active particles, all other terms encode a growth process where the active particles eventually become inactive: the active particles branch with rate 1, produce inactive particles with rate , and become inactive upon collision with either an active or an inactive particle. The active particles grow logistically, which implies that the inactive particles grow at most exponentially. The resulting simple time-dependent bounds on A and I yield uniqueness and existence of smooth solutions, a suitable fixed-point theorem is presented in chapter 14 of Smoller (1994). More details about the underlying stochastic processes together with a non-rigorous derivation of this PDE can be found in Hannezo et al. (2017). Note that without the inactive particles, i.e. when , the remaining equation for A reduces to the well-known FKPP-equation (Fisher 1937; Kolmogorov et al. 1937).

The System (1.1) can be interpreted as a twofold degenerate Keller-Segel system (Keller and Segel 1971; Arumugam and Tyagi 2020): the active particles are not guided by a chemotactic gradient, but explore the space solely diffusively, and the inactive particles do not diffuse at all. Still, simulations of System (1.1) show that general solutions of (1.1) self-organize, a phenomenon which is typical for many different Keller-Segel systems (Painter 2019). The invading front of the system converges to a fixed shape: a pulse of active particles, that represents a layer of growth, is accompanied by a monotone wave of inactive particles, the resulting static tissue, as demonstrated in Fig. 1. This self-organization of the Reaction–Diffusion System (1.1) resembles that of the stochastic dynamics.

Fig. 1

Simulation of the Reaction–Diffusion System (1.1) for . Given a small initial heap of active particles and , two identical traveling fronts arise, the right one is shown. After the separation of the two fronts away from the origin, the density of the remaining inactive particles is given by and the front moves asymptotically with speed

For a wave speed , a right-traveling wave solves Eq. (1.1) via the Ansatz . We substitute , such that any traveling wave must be a solution ofThe occurrence of these seemingly stable traveling waves is quite surprising, since the System (1.1) features a continuum of steady state solutions:which is due to the fact that the inactive particles do not degrade. This continuum of steady states represents the difficulty when studying the system: we first need to find out which limiting states are chosen by the growth process.

Hannezo et al. presented a rich discussion of the Wave-Equation (1.2) along with numerics and several heuristics that show a deep connection with the original FKPP-equation, and predicted some of the following results. The goal of this paper is to give necessary and sufficient conditions for the existence of such traveling wave solutions and to analyze the shape of the wave form. Our main result characterizes a family of pulled traveling waves:

Theorem 1.1

Let and consider the System (1.1) and its traveling wave solutions given by (1.2). Set . For each pair such thatthere exists a unique bounded and positive traveling wave with speed c such thatThe function i(z) is decreasing, whereas a(z) has a unique local and global maximum. If , then convergence as is sub-exponentially fast and of order . If , then convergence as is exponentially fast. Convergence as is exponentially fast in all cases. The corresponding rates areMoreover, these are all bounded, non-negative, non-constant and twice differentiable solutions of Eq. (1.2).

Notice that this result is independent of the reproduction rate r, which affects the shape of the wave, but neither its limits nor the minimal speed of a positive solution. Hence, all non-negative and bounded traveling waves resemble the ones depicted in Fig. 2, consisting of a pulse of active particles and a monotone wave of inactive particles. These traveling wave solutions share many similarities with classical FKPP-waves of a single type of particles. Among other mathematical aspects, this will be discussed at the end of the paper, in Sect. 9. Notably, Theorem 1.1 analytically connects two continua of fixed points via a continuum of traveling waves. Our constructive approach is a novelty: we first prove the existence of almost constant solutions and then continuously deform these solution along the continuum of possible limits.

Fig. 2

Two different traveling waves with speed . The limits of the left wave are given by and . The limits of the right wave are given by and

Figure 1 shows a simulation of the System (1.1), starting with a small initial amount of active particles. After a short transition phase, we observe a front with fixed shape. Asymptotically, it equals the unique traveling wave with limits and speed , which is the minimal possible wave speed for this pair of limits. We observed this behavior for all compact initial data that we chose. Moreover, this wave seems to be stable against perturbations, as briefly discussed in the concluding Sect. 9. Even though it is only a first step into this direction, this paper sheds light at the ability of the Growth-Process (1.1) to self-organize and at the robustness of this mechanism, e.g. against errors of individual particles. Our theoretical analysis fortifies the numerical and biological findings of Hannezo et al., where a simple set of local rules organizes the growth of a complex epithelial structure. The underlying assumption of a logistic growth is quite natural, so similar rules might drive and regulate other growth processes as well, without the need for guiding gradients.

Outline of the paper

A sketch of the central ideas and techniques is presented in Sect. 2. The identity is proved in Sect. 3. The asymptotic behavior around the stable and unstable set of the traveling waves is analyzed in Sect. 4. A non-negative trapping region of a lower-dimensional sub-system is analyzed in Sect. 5. We use our knowledge about the sub-system to construct a suitable attractor of the full system in Sect. 6. Then, we connect the unstable manifold of the unstable set with this attractor, see Sect. 7. We complete the proof of Theorem 1.1 in Sect. 8. In Sect. 9, we highlight the similarity of the traveling waves with those of the original FKPP-equation and give a short outlook at their stability.

Identifying the correct limits

We reformulate the System (1.2) for a traveling wave as an equivalent system of first-order ODEs. Denoting differention with respect to z by a prime, we introduce the auxiliary variable , so that (1.2) becomeswith . We call a solution of Eq. (2.1) non-negative if . If a solution is in , then it is also in by a simple induction. This equation has a continuum of non-negative fixed points, similar to that of the PDE, cf. (1.3):Thus, in the first place, we need to find out which of these fixed points can be considered as limits of right-traveling waves. Any bounded and non-negative solution of System (2.1) can not be periodic and must converge since . It is now evident that the limits at must be fixed points of Eq. (2.1), thus we denote them as . Under mild assumptions regarding integrability, we can interrelate two different points on a given traveling wave, see Sect. 3. Most importantly, this leads to the correspondence of the limitsIn view of this, monotonicity of i implies that and .

The fixed points of the ODE System (2.1) are not isolated, hence its Jacobian D is degenerate there. It is easily verified that D is given byAt a fixed point , the eigenvalues of areHence, we can not apply the classical Theorem of Grobmann-Hartmann to linearize the asymptotic behavior around the fixed points. We apply center manifold theory to work out the higher moments of the approximation, see Sect. 4. The center manifold coincides with the continuum of fixed points. This implies that asymptotically, there is no flow along the direction of the eigenvector with zero eigenvalue. Hence, the asymptotic flow around any fixed point is two-dimensional and the stability of the fixed point (0, 0, K) is dictated by the two eigenvalues . When , the fixed point is unstable, while for , it is stable.

At the same time, the analysis of the asymptotic behavior also yields a necessary condition on the speed c of a non-negative wave. A traveling wave can only be non-negative if does not spiral while converging to 0. Therefore, the two eigenvalues at the limiting fixed point must be real-valued. In view of (2.5), for a fixed point (0, 0, K), this is given ifThus, if the stable fixed point is the limit of a non-negative traveling wave, where , it must by (2.6) further hold thatas in Theorem 1.1. In other words, is the minimal limiting density of inactive particles that is necessary for the existence of a non-negative traveling wave with speed c.

Construction of a traveling wave

We will explicitly construct a non-negative traveling wave such that the two necessary conditions and are fulfilled. Two key features of the model make it tractable: first, the monotonicity of i(z) allows us to investigate the convergence of the sub-system that arises for a fixed value of i, and then lift our result to almost constant solutions of the full system. Second, for extending this result to non-small solutions, we lean on an integral equation that allows us to interrelate two points on a given trajectory. However, the central Proposition 3.2 depends essentially on the logistic growth of the active particles. Apart from this, our general approach seems to be applicable to a broader class of systems.

Regarding the ODE System (2.1), our analysis of the flow around the fixed points in Sect. 4 reveals a suitable unstable setand a suitable stable setEach point has an unstable manifold of dimension one. Its restriction to is the only possible candidate for the tail of a non-negative traveling wave as . Each point is Lyapunov stable, which can also be seen in Fig. 3.

Fig. 3

Two-dimensional phase portrait of (a, i) of traveling waves (2.1) for and , omitting the coordinate . A unique trajectory emerges from each point in (where ) in positive direction of a and converges to (where ). Notice the correspondence of the limits

To begin with, we let , where , and follow its unstable manifold in positive direction of a. We prove that there exists a finite phase-time such that : the trajectory reaches a local maximum of active particles, again see Fig. 3. We denote it as . This is carried out in Sect. 7.

Thus, for finding a suitable attractor of , we analyze solutions that start in points of type . We first analyze the lower-dimensional subsystem in coordinates (a, b), imposing a fixed value of i, which is done Sect. 5. We construct a trapping region, wherein a converges and stays non-negative. The monotonicity of i(z) allows us to lift this result to the full system, see Sect. 6. Here, the continuum of fixed points comes at help: we first prove the existence of almost constant solutions, where and , that stay non-negative and converge. Then, we continuously deform these solutions: the Lyapunov-stability of the limits in implies continuity of the entire trajectory up to in initial data. We use this to derive sharp conditions regarding such that the trajectory stays non-negative and converges.

We show that the first local maximum along the instable manifold of does fulfill these conditions, see again Sect. 7. The technique is the same as for proving the identity , which is the starting point of our analysis and presented in the next section. The proof of Theorem 1.1 is completed in Sect. 8, where we bring together all the different pieces. The resulting continuous family of solutions is sketched in Fig. 3.

The mapping of the limits

We first verify global integrability of a non-negative solution:

Lemma 3.1

Let be a smooth, bounded and non-negative traveling wave that solves the ODE System (2.1). Then, as , vanishes and converges, and . Moreover, is decreasing and has a unique global and local maximum.

Proof

Let be a smooth, bounded and non-negative solution of Eq. (2.1). Since , it holds that is decreasing and by boundedness must converge as , so . The two limits must be fixed points and given by some and . Equality is only given if . If not , then there is at least one local maximum of active particles, which we denote as .

At this point, , so either and the wave is constant, or . In the second case case, assume that there is also a local minimum of , denoted as . Since vanishes as , we may assume that this be the first local minimum after passing through . As before, is already a fixed point or . Since is decreasing, must have been increasing, a contradiction to the assumption that this is the first local minimum after the maximum . Thus, there is only one local maximum of active particles, which is also the global one. Further, this implies . By , we know that is also in . We integrate over the finite interval , then send the boundaries to :We know that the right-hand is integrable since , and that both and vanish as . This impliesHence also , since . Finally, as a sum of integrable terms, also .

The following Proposition 3.2 will be used several times to interrelate two points on a traveling wave, where . By the previous Lemma, the necessary conditions regarding integrability are always verified for non-negative and bounded solutions.

Proposition 3.2

Let be a smooth and bounded solution of the ODE System (2.1) on some interval , where . Assume that . Further, assume that are integrable and define . The following three identities hold:

Any solution of the ODE System (2.1) also fulfills the original Wave Equations (1.2). We integrate these over , substitute and use that . This directly proves (3.3) and (3.4). Regarding Eq. (3.5), note that by integration by parts:

Remark

Equations (3.3) and (3.4) encode a mass transfer from the active to the inactive particles and are not specific for the chosen reactions. It is the quadratic Eq. (3.5) that relies on a logistic saturation mechanism, we do not see a (direct) way to generalize this result.

Given Proposition, the identity is a mere

Corollary 3.3

(Limits of traveling waves) Let be a non-negative and bounded traveling wave that solves the ODE System (2.1), and denote its limits as . Either implies , or the following identity holds:

We apply Proposition (3.2). Since , the Eqs. (3.3), (3.4) and (3.5) can be simplified toEither implies , or we divide the first equation by and solve the resulting linear system, which proves the claim.

Asymptotics around the fixed points

Let us recall that the Jacobian of the ODE System (2.1) at a fixed point (0, 0, K) has eigenvaluesThe existence of implies the existence of a center manifold. In the present case, it locally coincides with the set of fixed points . This implies that there is no flow along the center manifold, so the asymptotics are fully described by the remaining two linear terms. The calculations are standard and presented in “Appendix A”, along with a short review of the underlying theory. We only state the results here. First, regarding the unstable set , as defined in (2.8):

Theorem 4.1

(Unstable set) For , the point is an unstable fixed point of Dynamics (2.1). Locally, there exists a smooth unstable manifold of dimension one. Its restriction to is the unique trajectory that emerges from the fixed point such that as . It has the following properties:

By choice of , the eigenvalue is positive, whereas is negative. Denote by u, v, w the coordinates in the system of eigenvectors of the Jacobian at the fixed, where the fixed point is shifted to the origin. This transformation is done explicitly in Lemma A.11. By Theorem A.12, the dynamics in an open neighborhood around the fixed point are equivalent toHence, there is a stable and an unstable manifold, each of dimension one. The eigenvector describes the asymptotic direction of the unstable manifold, in coordinates a, b, i it is given bySince and , asymptotically along the branch of the unstable manifold in direction , where and , it also holds that and that .

Next, we prove Lyapunov-stability of the points in , defined in (2.9). Figure 4 shows how the phase lines converge to (0, 0) in the (a, b)-plane. For technical reasons, we require that . Later, we deal with this degenerate case via a continuity argument.

Fig. 4

Phase portrait of (a, b) of the Wave Eq. (2.1) if we impose a fixed value of , see also Sect. 5. The choices of c change the type of convergence towards the origin: spiraling for , one stable manifold with eigenvalue , which has algebraic multiplicity 2 and geometric multiplicity 1 for , two stable manifolds for

Theorem 4.2

(Stable set) For all and , such that , the point is Lyapunov stable under Dynamics (2.1). In a neighborhood of the fixed point, exponentially fast.

By choices of c and , both non-zero eigenvalues (4.1) of the Jacobian are real-valued and negative and it holds that . As before, denote by u, v, w the coordinates in the system of eigenvectors of the Jacobian at the fixed point, see Lemma A.11. By Theorem A.12, the dynamics of the system in a neighborhood around the fixed point are equivalent toTake some small enough initial data : in view of Eq. (4.5), does not vanish, but also does not propagate, whereas v and w converge to zero exponentially fast. Since a and b are represented in terms of v and w, see (A13), they vanish exponentially fast.

Proposition 4.3

Let and . There is no non-negative and non-constant traveling wave that converges to as .

As in the previous Theorem, the asymptotic behavior of Eq. (2.1) around the limiting fixed point is described by the linear System (4.5). But now, since , either or both eigenvalues have a non-vanishing imaginary part and thus, v and w spiral. Since a and b are represented in terms of v and w, see (A13), any trajectory that converges to can not stay non-negative in its a-component.

Attractor of a sub-system

Construction and result

We begin our search for a non-negative attractor of in an easier setting: we fix and investigate the two-dimensional sub-system in the remaining coordinates. To separate it from the full system, we write it as . For this system, we prove the existence of a suitable attractor. This set will be denoted as , to emphasize that it depends on the chosen value of i, which will be constant only in this section. The flow of the sub-system and the region are drawn in Fig. 5.

Fig. 5

The phase plot of following Eq. (5.1), displayed for several values of i and c. The only two fixed points are (0, 0) and . For , the orange triangles are invariant regions of Dynamics (5.1), see Prop. 5.3. They increase in : the point moves to the right and the two internal angles and increase. In the third case, implies that the system spirals around (0, 0) while converging

Definition 5.1

(Two-dimensional sub-system) For and , denote by the two-dimensional flow defined bywhich results from the Wave System (2.1) by fixing .

There are only two fixed points of (5.1), and . We denote the eigenvalues and eigenvectors of the Jacobian at (0, 0) asIt holds that . Moreover, for , it holds that . Note that are identical to the non-zero eigenvalues of the full system around the fixed point (0, 0, i), see (4.1). The eigenvectors are the projections of the corresponding three-dimensional eigenvectors into the (a, b)-plane.

The Jacobian at has eigenvalues and eigenvectorsand it holds that . These have no direct correspondence to the three-dimensional system.

We now define the region . It is a triangle, spanned by the two fixed points (0, 0) and and two adjacent eigenvectors:

Definition 5.2

(The triangle ) For and , let be the convex hull of the three points and C(i). Here, the point C(i) is the unique intersection of the two half-lineswith and defined in (5.2) and (5.3). We denote the internal angles of at (0, 0) and as and , respectively.

Visually, it can easily be seen in Fig. 5 that the set is invariant under Dynamics (5.1): the flow at the boundary of points inwards.

Proposition 5.3

(Invariant region of the reduced system) The set is an invariant region of Dynamics (5.1). If , thenIt holds that and within . Hence, if , then converges to 0 monotonically as .

For any non-negative solution of the full Wave System (2.1), it holds that . Thus, we are interested in how changes when i decreases:

Proposition 5.4

(Nested invariant regions) For a fixed , the set is increasing in . Thus, for all .

This proposition holds due to an easy geometric argument, again take a look at Fig. 5: when i decreases, the point moves to the right and the two internal angles and increase. The computations for both propositions are performed in the following Sect. 5.2. The reader might skip those and proceed with Sect. 6, where we investigate the full system.

Invariance and monotonicity of

We analyze the -system and the set in detail. We prove that the Flow (5.1) at the boundary of points inwards, and that the sets are increasing in . We formalize what is sketched Fig. 5, and begin by examining the eigenvector at the fixed point :

Lemma 5.5

Let , and let be defined as in (5.2). The quotient of the absolute values of the -component and -component of is increasing in .

Proof

The claim is equivalent toRecall that . A computation reveals thatthe last inequality holds since implies .

Now, since is the angle between the two vectors (1, 0) and , the previous Lemma directly implies

Corollary 5.6

Let and . It holds that , the angle is increasing in .

For the invariance of , we need

Lemma 5.7

For any , the Flow (5.1) at the point points inwards .

Let be orthogonal to and point inwards . The claim of the Lemma is now equivalent toLet . First compute the Flow (5.1):Its part in direction is given by

For the fixed point , we get similar results concerning its unstable eigenvector :

Lemma 5.8

Let , and let be defined as in (5.3). The quotient of the absolute values of the -component and -component of is increasing in :

Recall that . A computation reveals that

Now, since is the angle between the two vectors and , the previous Lemma implies

Corollary 5.9

Let and . It holds that , the angle is increasing in .

For the invariance of , we need

Lemma 5.10

For any , the Flow (5.1) at the point points inwards .

Let be orthogonal to and point inwards . The claim of the Lemma is now equivalent toLet . We compute the Flow (5.1) atIts part in direction is given bySince , it follows that . Since , the proof is complete if we can show thatAfter dividing by and rearranging, this is equivalent toThis is in fact an equality, since .

Considering the invariance of , we conclude the

Proof of Proposition 5.3

We need to show that the Flow (5.1) at the boundary of points inwards. The Lemmas 5.7 and 5.10 treat the left and right edge of , see again Fig. 5. For the third edge, we consider points of type , where . The derivative at is given by . Its -component is negative, hence it points inwards . The only points on the boundary of where the flow does not point strictly inwards are the two fixed points (0, 0) and .

Considering the monotonicity of , we conclude the

Proof of Proposition 5.4

Let . The point moves to the right as i decreases. Further, we have shown that the two internal angles increase in , and so does .

Attractor of the full system

We now analyze solutions of the full Wave System (2.1) under initial condition , such that as defined in the previous Section. In Sect. 6.1, we apply the results about the two-dimensional subsystem to the full system. Theorem 6.1 states that as long as , the -components of the full system stay within the triangle . Thus, it suffices to control , which we do in two steps.

In Sect. 6.2, we prove via some rough bounds that for sufficiently small initial values . This result is refined in Sect. 6.3: the Lyapunov-stability of the limiting point at implies that the entire trajectory including its limit is continuous in initial data. Carefully increasing , we increase the known attractor of the stable set , resulting in Theorem 6.10. This procedure is sketched in Fig. 6.

Fig. 6

Trajectories of of the Wave System (2.1) for . Initial values are , and such that . The upper bound is given in Definition 6.8. Trajectories with such initial data converge and stay non-negative, since

Assumption

If not explicitly stated otherwise, we will use the following setup over the entire Sect. 6: For , let and , which implies that . Let be the solution of the Wave Eq. (2.1) under initial values .

Invariant region of the full system

Theorem 6.1

(Invariant region of the full system) Assume that for all . We then can control the two remaining coordinates a(z), b(z) of the wave. It holds that

Within , and . Notice that while , it holds that . This directly implies the following

Corollary 6.2

Under the assumption that for all , the trajectory stays non-negative and converges as :

Proof of Theorem 6.1

In the full System (2.1) with coordinates (a, b, i), neither nor depend on , but only on a and i. Thus, we can easily compare the full system to the two-dimensional System (5.1) in coordinates . At a phase-time , the two vector fields and for fixed value are equal, compare (2.1) and (5.1).

By Proposition 5.4, this implies that points strictly inwards . There are two irrelevant exceptions: for , the system has already reached its limiting state. The point is a fixed point of the reduced, but not of the full system. In this case, since and , a Taylor-expansion reveals that for small times : changes much faster than and .

Importantly, the set is not decreasing as a function of . Hence, at each phase-time , the two components can not escape . In fact, since and is increasing in , the set is increasing in , at most up to . Thus, remain within for all .

With a similar argument, we can determine the rate of convergence:

Proposition 6.3

Assume that for all , such that as for some . If , then convergence is exponentially fast with rateFurther, if , which can only happen if , then the system converges sub-exponentially fast. As , the distance to the limit is of order .

In the case , all eigenvalues of the limit are simple, we refer to Sect. 4. The system converges exponentially fast, as shown in Theorem 4.2. It remains to determine the rate of convergence. The two candidates are . Corresponding to , the projections of the eigenvectors into the (a, b)-plane are given byWe know that for all . At (0, 0), the triangle is bounded by the line , see Definition 5.2 and Fig. 5. Since , the direction of is steeper than that of , such that the line lies outside for all . Thus, the two components a(z), b(z) cannot converge towards (0, 0) along . But since they converge exponentially fast, the only possible remaining rate of convergence is .

We do not have a complete description of the asymptotics around the fixed point for the degenerate case . However, under the assumption that the system stays non-negative and converges, the relation holds asymptotically. This will be proven in Theorem 6.10 (which does not rely on the type of convergence). The eigenvalue has algebraic multiplicity 2, but geometric multiplicity 1. It is well-known that this results in sub-exponential convergence, cf. chapter 9 in Boyce et al. (2017).

A small attractor

In view of the previous paragraph, convergence and non-negativeness follow if we can show that for all . If we choose small enough, some rough bounds do the trick. We control the total mass of active particles via

Lemma 6.4

Fix and let . Under the assumption that is true for all , there exists a finite constant , such that the following bound holds for all :

With the help of Theorem 6.1, we can use that for all . We integrate and use that :By monotonicity: and . It follows thatwhere we need to avoid a blow-up, which is true by our choice of . It remains to bound . Take a look at the flow in the (a, b)-plane in Fig. 5. It holds that stay within . Within the triangle , it holds for the left inner angle thatThus, also .

Since we can bound the total mass of active particles, we can also bound the change of i(z):

Proposition 6.5

(Small attractor of ) Fix and . There exists a constant , such that for all :Hence, also for all . The trajectory is non-negative and converges to as .

As long as , it must be that by Theorem 6.1. Assume there exists finite phase-time :where we used . For and sufficiently small, we can apply Lemma 6.4. This implies that there is a finite constant L, which does not depend on , such thatThe right-hand side is strictly larger than for sufficiently small , say , and the bound is independent of the phase-time. Thus, there is no such for .

Extending the attractor

The previous section ended with a condition of type , under which the system stays non-negative and converges. However, given and and under the assumption that the system converges, we can explicitly calculate its limit . Then, for fixed , we continuously deform the trajectory while increasing up to some upper bound , as sketched in Fig. 6. This results in Theorem 6.10.

We apply Proposition 3.2 to interrelate the limit of the trajectory to its initial data :

Lemma 6.6

If , such that the system stays non-negative and converges to as , then can be written as a function of and :The function is decreasing in , for .

We apply Proposition 3.2, and solve the resulting system for . In the present case, since , this results in a quadratic equation with two possible solutions. Since is decreasing, it must be that , which uniquely determines (6.12). A short computations proves that .

We look for values of that ensure . Thus, we rearrange (6.12) for , set , and choose the only possible positive solution of the resulting quadratic equation:

Lemma 6.7

Given and under the assumption that , the value of is uniquely determined by

Equation (6.13) can be restated as , but keep in mind that still have to prove convergence. It can easily be seen that . Since we require that , such that , this leads to our

Definition 6.8

(Upper bound for ) For fixed and, we define

This will hold as sharp upper bound for , such that the trajectory stays non-negative and converges. Before we state the corresponding theorem, we perform a last check that we are in the correct setup:

Lemma 6.9

Let and . If the system stays non-negative and converges to , thenwhere is given as in Lemma 6.6.

It holds that . The claim follows since is decreasing in and since , where .

After all these preparations, we are now ready to prove that the system converges and stays non-negative. We state

Theorem 6.10

(Attractor of ) For , let and . Let be the solution of Eq. (2.1) with initial data .

It holds that and for all . As , and converge to 0, and converges toThe type of convergence depends on and is given in Proposition 6.3.

Notation: We fix and change only . If , we must choose . For a compact notation, is the state of the system at phase-time , starting in . If the limit of a trajectory exists, we denoteStep 1: starting interval

For positive but small enough, Proposition 6.5 grants that for all :where is a bounded invariant region that contains only points such that . Thus, monotone. With the help of Lemma 6.6, we can explicitly calculate as stated in Eq. 6.16, and our claim holds on some small non-empty interval .

Step 2: neighborhood of existing trajectories

Pick some for which the statement is already proven. By choice of , it holds that . Thus, the limit is Lyapunov stable by our previous analysis of the asymptotics, see Theorem 4.2: for every , there exists a , such thatfor all . Choose and assure that . This grants after entering the -neighborhood. Within this attractor, also in view of Theorem 6.1, since .

Starting in , we follow the trajectory up to some finite time , where it has entered the -neighborhood:The derivative of the system is locally Lipschitz and all trajectories are within a bounded domain. Thus, the trajectories are uniformly continuous in initial data on any finite time interval [0, T], with respect to the maximum norm . This is a classic result and can easily be proven via a Grönwall’s inequality, we refer to chapter 2 of the textbook of Hsieh and Sibuya (1999). There exists some , s.t. for all with :This implies for all such trajectories :In particular, lies within the -neighborhood, so ultimatelyAgain, Theorem 6.1 implies for all . As before, the system is integrable and converges as , so we can explicitly calculate . Thus, our claim holds for all starting points a small open neighborhood of .

Step 3: limits of trajectories

Assume that the claim holds for all . For all trajectories starting in , where , it holds that is monotone and bounded from below by , such that stay within . Fix any finite time-horizon [0, T]. As mentioned before, the trajectories are uniformly continuous in initial data on finite time-intervals, and thus form a Cauchy-sequence on as . Since T can be chosen arbitrarily large and since the limits are also continuous in , our claim holds for the trajectory that starts in .

Step 4: conclusion

By step 1, the claim holds for in some small interval . By step 3, it then also holds for . If now , the claim holds for by step 2 for some . Iterating these two steps, the claim ultimately holds for all . In particular, we have proven that the trajectories are uniformly continuous with respect to initial data on . This continuity allows us to finish the proof of Proposition 6.3. In the non-critical cases where , the trajectories converge along a stable manifold with rate of convergence . As , the critical trajectory must converge along the limit of these manifolds.

The complete trajectory

We track the non-negative branch of the unstable manifold of , see (2.8), and show that it stays positive and enters the attractor of from the previous section, cf. Theorem 6.10.

We will use the following setup over the entire Sect. 7: Let and let a(z), b(z), i(z) be the unique solution of the ODE-System (2.1) that emerges from as and where asymptotically as .

For all , existence and uniqueness of these trajectories have been proven in Sect. 4. Moreover, we know their asymptotic behavior:

Lemma 7.1

The following holds as :

The first two inequalities are given by Theorem 4.1, which also yields asymptotically. Noticing that completes the proof.

The maximum of active particles

For connecting these trajectories with the attractor of , we will prove

Proposition 7.2

(The maximum of active particles) There exists a finite phase-time , such that for the first time.

We will prove that the sum decreases below 1. Given this, the term becomes negative, so b must eventually reach 0.

Lemma 7.3

As long as for all , it holds that

In view of the asymptotic behavior of the trajectory, described in Lemma 7.1, assume that there exists a finite time , such that for the first time , but still . The wave-equations and imply thatRearranging the third line yields . As long as , it holds that , hence alsoHowever, equality at implies thatwhich we can rewrite, using both (7.4) and (7.5) :But , so Eq. (7.4) implies that . This contradicts (7.8).

Lemma 7.4

As long as for all , it can not happen that converges to some finite .

By the previous lemma: while . Assume that converges to a finite value from above, which we denote as . This implies that also . By the Wave Equations (1.2), these two expressions are equivalent toThe first convergence indicates that after some time , since a is strictly increasing and hence positive. The second statement is equivalent to . Thus, also b is increasing. But for all and while , a contradiction.

We can now show that there exists a finite phase-time such that , finishing the

Proof of Proposition 7.2

By the previous lemma, decreases below every positive value as long as . In particular, for some : after some phase-time . Then for all , since :Either and the system has already passed a first local maximum of , or we may assume that . If now , then reaches zero after a finite time , which can not be larger than .

Reaching the attractor of

We now prove that lies in the attractor of the stable set . Therefore, we show that , then Theorem 6.10 ensures non-negativity and convergence as . We again use Proposition 3.2, now to interrelate and :

Lemma 7.5

The following holds at phase-time :In the case , then additionallyand the trajectory is non-negative on the interval .

As , all have exponential and hence integrable tails, cf. Theorem 4.1. We thus can apply Proposition 3.2 on the interval , for finite . This results in Eq. (7.13), we omit the intermediate steps. It holds that because for all .

In particular, implies that the second summand under the root in (7.13) must be strictly positive, which yields . Since for the first time, it must hold that . Given this, we use to bound . Since , this shows that .

If we assume additionally that , then implies that . Up to , is decreasing, which was proven in Lemma 7.3. Since a(z) is strictly increasing up to , i(z) must be strictly decreasing, but not below . Hence, the trajectory stays positive. The inequality implies that .

This allows us to connect the unstable manifold of with the attractor of :

Proposition 7.6

(Reaching the attractor of ) Let . The non-negative branch of the unstable manifold of reaches the point , where and . It then holds thatfor like in Definition 6.8. In view of Theorem 6.10, the trajectory that starts/continues in such a point converges to and stays non-negative.

We have just shown that and that . Recall Definition 6.8: , where is given byWe have already verified that , so proving suffices for proving . By (7.13), we know thatThe two expressions (7.16) and (7.17) are very similar. After some elementary steps, the claim is equivalent toThis is equivalent to , since and . But that is just how we have chosen .

Concluding the proof of the main result

With the results from the previous sections, we complete the

Proof of Theorem 1.1

Consider the ODE System (2.1) in coordinates a, b, i and let . The unstable manifold of the fixed point has dimension one. Its two branches are the only trajectories that leave the fixed point, which is stated in Theorem 4.1. There is one branch of the unstable manifold such that as , we follow this trajectory in positive direction of . There is a finite phase-time , such that for the first time , see Proposition 7.2. For all , it holds that . Denote the state of the system at as . Lemma 7.5 states that , Proposition 7.6 states that , for as in Definition 6.8. By Theorem 6.10, we then know that lies in a non-negative attractor of the set . Thus, as , where , and ultimately for all .

For any non-negative and bounded solution, the identity holds by Proposition 3.3. For and , the previous paragraph proves existence and uniqueness of the claimed wave. For , the constant solution can be the only non-negative and bounded one.

We then consider an arbitrary non-constant, bounded and non-negative solution. By monotonicity of i(z), it must converge as . If we assume that , it is one of the above solutions. If we assume that , then . In this case, the trajectory can not stay non-negative as , which is stated by Proposition 4.3, contradicting the assumption.

Discussion and outlook at stability

FKPP-waves

We have given a description of all bounded and non-negative traveling waves of the Reaction–Diffusion System (1.1). For the most related systems, the FKPP-equation (Kolmogorov et al. 1937; Fisher 1937), the FitzHugh-Nagumo-equation (FitzHugh 1961; Nagumo et al. 1962) and combustion equations (Berestycki et al. 1985), no such continuum of traveling waves has yet been constructed.

Still, the non-negative traveling waves of System (1.1) are closely related to pulled FKPP-waves with only a single type of particles. The equation for such a wave w(z) reads . For the purpose of a simple comparison, we let , where is the initial growth rate of the particles. In this case, Theorem 1.1 states that the convergence of System (1.1) as is identical to that of w, if , see e.g. (Uchiyama 1977). In words, the asymptotic growth speed of traveling waves of System (1.1) coincides with that of simple FKPP-waves in presence of a constant density of inhibiting particles. Moreover, Theorem 1.1 implies that for all . Thus, the minimal speed of an invasive front, where , is given by . Again, this coincides with the minimal wave speed of the associated FKPP-equation, i.e. in the absence of inactive particles. It is this critical front which can be interpreted as the most natural one, our simulations indicate that it always arises under compact initial data. If we assume convergence, a technique of Berestycki et al. (2018) yields an upper bound for the speed of the traveling front, just by ignoring the dampening influence of the inactive particles. For compact initial data, the system always chooses the smallest possible wave speed, as suggested.

The emergence of traveling fronts is known for many reaction–diffusion systems. We suggest the literature (Britton 2003; Volpert and Petrovskii 2009; Othmer et al. 2009) for more examples with a biological motivation. Rigorous proofs of these phenomena are rare. Often, only the form of the traveling waves is analyzed analytically. The FKPP-equation is one of the cases, where the convergence of the front of the PDE towards a traveling wave solution can be proved. The first rigorous proof was done by Kolmogorov, Petrovsky & Piscunov in 1937 Kolmogorov et al. (1937). Extensions of this result to more general initial data and a more precise description of the speed of the front have been provided by Uchiyama (1977) and Bramson (1983). The approach of Kolmogorov et al. and Uchiyama seems to be restricted to systems with only a single type of particles, as it relies on a maximum principle and monotonicity of the front. The approach of Bramson relies on a relationship between the FKPP-equation and branched Brownian motion, which can not be applied in the present case since the inactive particles do not diffuse. A singular perturbation of System (1.1), which again introduces a small diffusion to the inactive particles will be subject to future investigations. Despite the fact that this system would be biologically interesting, since no tissue or population is entirely static, this would also rule out some difficulties when analyzing the stability of the traveling waves against perturbations.

Stability of the traveling waves

We give a brief introduction to the stability of traveling waves against small perturbations, in the spirit of the introduction in Ghazaryan et al. (2013). A good overview, where the following concepts are presented in greater depth, has been written by Sandstede (2002).

Consider a reaction–diffusion systemwhere with , and R a smooth reaction. In the moving frame , the System readsA traveling wave w(z) with speed c is a constant solution of Eq. (9.2). The wave w is called non-linearly stable in a space , if any solution of the PDE (9.2) which starts in , where is a sufficiently small perturbation, converges to a shift of w. This type of stability is often encoded in the spectrum of the operator , that is obtained by linearizing the equation for the perturbation in (9.2) around to the constant part w:where JR is the Jacobian of the reaction R. Let be the operator given by . We say that the wave w is spectrally stable in if the spectrum of is contained in the half-plane , except maybe a simple eigenvalue at 0 (that corresponds to the traveling wave itself, if ). For diffusive systems, where , a quite general theory has been developed. If is appropriately chosen, spectral stability implies non-linear stability, we refer to the literature (Sandstede 2002; Ghazaryan et al. 2013). Classical results are e.g. given for subspaces of , the -Sobolev space.

As explained in the next paragraphs, we are not aware of any rigorous framework for studying the non-linear stability of System (1.1). Several technical problems arise, that so far have been treated only separately (Ghazaryan et al. 2013; Kirchgässner 1992).

Most importantly, the traveling waves of System (1.1) can not be stable against perturbations in the classical sense, since the inactive particles neither react nor diffuse. Any initial deviation remains for all times, as shown in Fig. 1. However, the actual front of the system does converge to a traveling wave. For capturing this idea, we introduce the weighted space with norm for some . Non-linear stability in is referred to as convective stability. Convergence of the PDE in the moving frame (9.2) in means that the front of the system approaches the traveling wave, whereas any initial finite and local deviation is convected towards and vanishes due to the weighting. A first rigorous result regarding convective stability was obtained by Ghazaryan et al. (2013). They could show that in some cases, spectral stability in implies convective stability against small perturbations in . For their approach, the authors require that the weight can be chosen such that all eigenvalues of except zero fulfill and such that the derivative of the traveling wave is an eigenfunction that corresponds to a simple eigenvalue at zero. Unfortunately, this setting is not suited for studying pulled FKPP-fronts: the assumption implies that the continuous spectrum of touches the origin, see e.g. chapter 6 in the work of Sattinger (1976).

Another difficulty arises when studying critical pulled fronts (with minimal possible speed) whose tail as converges sub-exponentially, as in Theorem 1.1. In this case, the requirement is only fulfilled for rather small values of , which do not suffice for shifting the continuous spectrum of to the left half-plane. For diffusive systems, non-linear stability of this more delicate case was first treated rigorously by Kirchgässner (1992), a recent overview is given in Faye and Holzer (2018). After introducing a small diffusion to the inactive particles, we could apply this theory to the critical front.

For the most natural traveling wave solution, the critical one with speed and , we performed a numerical analysis that strongly indicates that this wave is spectrally stable in , when we choose . The details are presented in “Appendix B”. Thus, based on the work of Ghazaryan et al. regarding convective stability (Ghazaryan et al. 2013) and the work of Kirchgässner regarding critical fronts (Kirchgässner 1992), we dare to make an educated guess: we expect that this traveling wave is convectively stable against small perturbations in .