Stable patterns with jump discontinuity in systems with Turing instability and hysteresis

Classical models of pattern formation are based on diffusion-driven instability (DDI) of constant stationary solutions of reaction-diffusion equations, which leads to emergence of stable, regular Turing patterns formed around that equilibrium. In this paper we show that coupling reaction-diffusion equations with ordinary differential equations (ODE) may lead to a novel pattern formation phenomenon in which DDI causes destabilization of both constant solutions and Turing patterns. Bistability and hysteresis effects in the null sets of model nonlinearities yield formation of far from the equilibrium patterns with jump discontinuity. We derive conditions for stability of stationary solutions with jump discontinuity in a suitable topology which allows disclosing the discontinuity points and leads to the definition of ({\epsilon}0 , A)-stability. Additionally, we provide conditions on stability of patterns in a quasi-stationary model reduction. The analysis is illustrated on the example of three-component model of receptor-ligand binding. The proposed model provides an example of a mechanism of de novo formation of far from the equilibrium patterns in reaction-diffusion-ODE models involving co-existence of DDI and hysteresis.


Introduction
Since the seminal paper of Alan Turing [26], mathematical models of biological pattern formation have been constructed by using systems of reaction-diffusion equations exhibiting diffusion-driven instability (DDI). DDI is related to a local behavior of solutions of a reaction-diffusion system in the neighborhood of a constant stationary solution that is destabilized through diffusion. It may lead to emergence of stable continuous and spatially periodic structures around the destabilized constant equilibrium. The Turing concept became a paradigm for pattern formation and led to development of numerous theoretical models, though its biological verification has remained elusive [1,6,21].
The models are based on the idea that cells differentiate according to positional information which is supplied to them by diffusing biochemical morphogens. Different morphogen concentrations are able to activate transcription of distinct target genes and thus lead to cell differentiation. However, both regulatory and signaling molecules (ligands) act by binding and activating receptor molecules, which are located in the cell membrane [11,22]. This observation leads to a hypothesis that the positional value of a cell may be determined by the density of bound receptors which do not diffuse. Taking dynamics of receptors into account leads to systems coupling reaction-diffusion equations with ordinary differential equations (ODE) [9,12,13], called also receptor-based models. Such systems can be obtained as a homogenization limit of the models describing coupling of cell-localized processes with cell-to-cell communication via diffusion in a cell assembly [14,19]. Nonlinear interactions of diffusive and non-diffusive components arise also from modeling of interactions between cellular or intracellular processes. Reactiondiffusion-ODE models have recently been employed in various biological contexts, see e.g., [8,9,12,17,23,27]. 1 Although receptor-based models may exhibit Turing instability as shown in [12] and discussed more recently on several examples from mathematical biology in [9], they are different from classical Turing-type models. Since DDI is induced by self-enhancement in the subsystem with smallest diffusion [2], in reaction-diffusion-ODE models instability is induced by growth properties of the ODE subsystem. Therefore, coupling of reactiondiffusion equations to ODEs may lead to DDI which does not arise if the non-diffusive components are neglected.
To understand the role of non-diffusive components in pattern formation process, we focus on systems involving a single reaction-diffusion equation coupled to ODEs on onedimensional domain I. In general, equations of such models can be represented by the following initial-boundary value problem, x ∈ I, t > 0, (u(0), v(0)) = (u 0 , v 0 ) ∈ C(I) × (C 2 (I) ∩ C(I)), with homogeneous Neumann boundary conditions for v.
It is an interesting case, since a scalar reaction-diffusion equation cannot exhibit stable spatially heterogeneous patterns [5]. Coupling it to an ODE fulfilling autocatalysis condition, i.e. ∂f /∂u > 0 at the constant equilibrium, leads to DDI. However, in this case all Turing patterns are unstable, i.e. the same mechanism which destabilizes constant solutions, destabilizes also all continuous spatially heterogeneous stationary solutions [15,16]. This instability result holds also for patterns with jump discontinuity in case of a special class of nonlinearities [15,16]. The question then arises as to which patterns, if any, can be exhibited in such models.
As shown in [15], it may happen that there exist no stable stationary patterns and the emerging spatially heterogeneous structures have a dynamical character. Simulations of different models of this form indicate formation of dynamical, multimodal and apparently irregular structures, the shape of which depends strongly on initial conditions [7,17,23]. On the other hand, reaction-diffusion-ODE models may give rise to discontinuous patterns due to the hysteresis effect in the model nonlinearities, i.e. when the equation f (u, v) = 0 has multiple quasi-stationary solutions v i = H i (u). Diffusion tries to average different states and may lead to formation of far from equilibrium structures. Hysteresis yields emergence of stationary solutions with a jump discontinuity, which may be monotone, periodic or irregular. The hysteresis-driven pattern formation has been investigated analytically in systems with two stable constant equilibria, i.e. without DDI [13,18,10]. In such a case patterns depend strongly on initial conditions and require large initial heterogeneity.
Consequently, the next question is whether a coupling of the two mechanisms, DDI and hysteresis, can provide a mechanism of de novo formation of stable far from equilibrium patterns. By de novo pattern formation we understand emergence of spatially heterogeneous structures starting from small perturbation of the constant stationary solution. According to our knowledge, systems with DDI and hysteresis have not been studied so far.
Motivated by these observations, we propose a new receptor-based model exhibiting DDI and hysteresis which lead to de novo formation of stable patterns with jump-type discontinuities. The model is a modification of the simplest model describing receptorligand binding dynamics developed in [19,12]. The modification accounts for saturation in production terms, which provides uniform boundedness of the model solutions. We consider two-and three-component systems, i.e., • a model consisting of one ordinary differential equation and one reaction-diffusion equation (1ODE+1RDE), • a model consisting of two ordinary differential equations and one reaction-diffusion equation (2ODE+1RDE), While the principle of DDI is similar to Turing's original idea, coupling to an ODE is crucial for stability of steady states in this paper. One interpretation of the observed phenomenon is that non-diffusive species, such as receptors, can facilitate formation of patterns with sharp transitions.
In this paper, we derive conditions for stability of such patterns. Since stationary solutions for nondiffusing variables have jump discontinuities, linearized stability analysis requires careful treatment. To cope with this difficulty, we consider a special topology which allows us to include the discontinuity points and yields a notion of (ε 0 , A)-stability. Our approach follows the one of Weinberger proposed in refs. [28,3] to study discontinuous solutions in a model of density dependent diffusion. A novelty of our paper lies in the stability analysis of patterns with jump discontinuity in systems involving several ODEs which may be reduced to a subsystem by using quasi-steady state approximation. We show that under suitable conditions on the nonlinearities, stability of stationary solutions is preserved under quasi-steady state reduction.
The paper is divided into five sections: In Section 2 the main results of this paper are formulated. We consider systems of type (1.1) with scalar-valued v and u having values in either R or R 2 and provide conditions under which steady states with jump discontinuity are stable in the topology introduced in ref. [28]. Then, we extend the stability results to systems with a small parameter and the corresponding quasi-stationary approximation. In Section 3, we present a specific receptor-based model of type (1.1) exhibiting both pattern formation mechanisms, i.e. DDI and hysteresis. We show that the model satisfies the theory developed in this paper. Model analysis is supplemented by numerical simulations showing de novo formation of far from equilibrium patterns exhibiting jump-discontinuities of which the location depends both on the size of diffusion and on the initial data. Section 4 presents a detailed proof of the stability theorems stated in Section 2. In Section 5, we prove all statements from Section 3 related to the analysis of DDI initiated formation of hysteresis-driven patterns with jump discontinuity.

Main results. Stability conditions for steady states with jump discontinuity
We state our results on the stability of steady states of (1.1) with finitely many jump discontinuities under the following assumptions. (1) v is a scalar, while u is either a scalar or a vector consisting of two components.

The topology.
In this section, we introduce (ε 0 , A)-stability following the idea of Weinberger in [28] and provide conditions for stability of steady states of system (1.1) in this topology.
For a bounded interval I = (0, l), let BV (I) denote the space of all functions of bounded variation defined on I. We define a neighborhood basis forũ ∈ BV (I) with values in [u, u], where u, u are given by Assumption 2.1, as L ∞ (R) < ε 2 and meas(I \ R) < ε 4 }. We are particularly interested in the situation where the initial function u 0 is close to a steady stateũ with finitely many jump-type discontinuities, but u 0 is still continuous on I. If we assume v 0 ∈ C 2 (I) ∩ C(I), then the solution (u(t, x), v(t, x)) of the initial-boundary value problem is continuous for all t > 0. If the steady stateũ is 'stable', then u(t, x) is expected to converge towards a function with jump discontinuity in space. Therefore, the uniform norm is not appropriate to measure the closeness, see Figure 2.1. It turns out that N ε,u,u (ũ) provides us a reasonable topology for this purpose. Following [28], (ε 0 , A)-stability in this topology is defined as Definition 2.2 ((ε 0 , A)-stability). A stationary solution (ũ,ṽ) of system (1.1) is said to be (ε 0 , A)-stable for positive constants ε 0 and A if the initial functions (u 0 , v 0 ) satisfy for some R ⊂ I with meas(I \ R) < ε 4 , and for some ε ∈ (0, ε 0 ), then for all t > 0. Here H 1 (I) = {u ∈ L 2 (I) | u ′ ∈ L 2 (I)} and u H 1 (I) = u L 2 (I) + u ′ L 2 (I) .
2.2. Stability conditions. Definition 2.2 allows us to give conditions for the (ε 0 , A)-stability of steady states with jump discontinuity. First, we analyze a system of one ordinary differential equation coupled to one reaction-diffusion equation.
Theorem 2.3 (Stability for a scalar-valued u). Let I = (0, l) be a bounded interval in R. Under Assumption 2.1, consider the system of equations Let (ũ,ṽ) be a steady state with finitely many discontinuities ofũ. Denote the Jacobian matrix of the kinetic system at the steady state by and assume that the following inequalities hold for all x ∈ I: where c 1 is a positive constant independent of x. Then (ũ,ṽ) is (ε 0 , A)-stable for a pair (ε 0 , A) with 0 < ε 0 , A < ∞.
The proof of Theorem 2.3 is similar to that of stability for a particular model with cross-diffusion in [28]. We shall present it as a corollary to more general Theorem 2.4 in Section 4. Theorem 2.3 shows that the steady states of the model are stable, if both species are self-inhibitory and cross-inhibition/cross-proliferation is asymmetric. Even if each species inhibits all species, a steady state can be stable if the cross-inhibition is sufficiently small compared to the self-inhibition. However, in this case, external influx is necessary.
Next, we consider a model comprising of two ordinary differential equations and one reaction-diffusion equation, i.e., u = (u 1 , u 2 ): We are interested in the stability question on a given steady state (ũ 1 (x),ũ 2 (x),ṽ(x)) of (2.8) with finitely many jump discontinuities. It is convenient to put Theorem 2.4 (Stability for vector-valued u). Let I ⊂ R be a bounded interval. Under Assumption 2.1, consider a system of type (2.8). Let (ũ 1 ,ũ 2 ,ṽ) be a steady state with finitely many discontinuities of (ũ 1 ,ũ 2 ) in x. Denote the Jacobian matrix of the kinetic system at the steady state by A(x) as in (2.10). Let A ij denote the matrix resulting from omitting the i-th row and the j-th column, i.e. A ij = (a kl ) k =i,l =j . Assume that A(x) and for all x ∈ I, where κ is a positive constant.Then the steady state (ũ 1 ,ũ 2 ,ṽ) is (ε 0 , A)stable for some positive constants ε 0 and A.
Remark. Condition (2.11) is the Routh-Hurwitz condition for the matrix A(x), i.e., all eigenvalues of A(x) have negative real part if and only if (2.11) holds. Note that it implies 3 j=1 det A jj (x) > 0. The above theorem deals with general three-component systems, and it is non-trivial to check the set of conditions (2.11)-(2.14), practically speaking. However, in the case where one species, say u 2 , reacts very rapidly, the situation becomes more manageable. Therefore, we consider the following system with a small positive parameter δ: (2.15) Corollary 2.5 (Stability for vector-valued u with a small parameter). Let I = (0, l) ⊂ R be a bounded interval. Under Assumption 2.1, consider a system of type (2.15). Let (ũ 1 ,ũ 2 ,ṽ) be a steady state with finitely many discontinuities of (ũ 1 ,ũ 2 ) in x. Let A(x) be the matrix defined by (2.10). Let A ij denote the matrix resulting from omitting the i-th row and the j-th column, i.e., A ij = (a kl ) k =i,l =j . Assume that the following inequalities are satisfied for all x ∈ I: where k 1 , k 2 , k 3 are positive constants independent of x. Then there exists a positive constant δ * such that the steady state (ũ 1 ,ũ 2 ,ṽ) is (ε 0 , A)-stable for some positive constants ε 0 and A, provided that 0 < δ < δ * .

Quasi-steady state approximation.
By the quasi-steady state approximation of (2.15), we understand system (2.15) for δ = 0, i.e., the following system of differential-algebraic equations: Assuming that f 2 (u 0 1 , u 0 2 , v 0 ) = 0 has an isolated solution u 0 2 = u * 2 (u 0 1 , v 0 ), we can re-write (2.21) as It follows immediately that any steady state of the reduced system (2.22) gives rise to a steady state of (2.21). Hence, we want to investigate the following questions: (1) Does system (2.15) exhibit DDI if its quasi-steady state reduction (2.22) does?

Application to a system exhibiting DDI and hysteresis
In this section, we present a model exhibiting the coexistence of DDI and hysteresis. First, we propose a model consisting of two ordinary differential equations coupled to one reaction-diffusion equation. We reduce this model using quasi-steady state approximation. Then, we rescale the reduced model and show that it satisfies the conditions of Theorem 2.3. Finally, based on Theorem 2.6, we show that the stability of steady states of the reduced model implies stability of corresponding patterns in the original model.
The model is a modification of the receptor-based model proposed in [12,19]. Numerical investigations for κ = 0 show spike solutions blowing up in infinite time, see [7]. Introducing κ > 0 reflects saturation of de novo production. The modification assures uniform boundedness of the model solutions.
Variables u 1 and u 2 describe cell surface receptors at free and occupied states, respectively. The state of receptors is being changed through binding to diffusive ligands, denoted by v, and a reverse process of dissociation of bound receptors. Free receptors and ligands are produced de novo, what is described by a function with saturation effect for κ > 0.
For δ = 0, (3.1) reduces to a system of differential-algebraic equations. In this case, the second equation of (3.1) can be solved uniquely for u 0 2 leading to a system of one ordinary differential equation coupled to one reaction-diffusion equation. This reduction leads to simplification of analytic investigation of the model. We show that the limit δ → 0 is regular in the sense of stability-preservation of steady states.
To simplify notation, we define We show that the model has a spatially homogeneous steady state that is stable to spatially homogeneous perturbation and is destabilized due to autocatalysis in dynamics of nondiffusive component. Due to the nonlinear nature of the reaction terms, there exists a 'far from equilibrium' regime, where the species become self-inhibitory, stabilizing steady states with jump discontinuity.
Our first step is to show that models (3.3)-(3.5) and (3.1) satisfy Assumption 2.1, i.e., the solutions are uniformly bounded for nonnegative initial functions. Numerical simulations of the model solutions are presented in Figure 3.2. We formulate the results in the remainder of this section and defer their proofs to Section 5.

Existence of steady states.
By definition of quasi-stationary reduction and linear rescaling, there exists a one-toone mapping between the sets of steady states of model (3.3)-(3.5) and of model (3.1). Therefore, we may focus on existence of steady states of model may consist of more than one element. We characterize different possible branches of solutions for given v: We show in Lemma 3.7 that the steady states have different stability properties depending on the branch.
The existence of steady states can be summarized in the following lemma, illustrated in Figure 3.1:  (2) Let m 1 < m 2 and Then, there exists a positive k * 1 such that for all k < k * 1 , system (3.3)-(3.5) has exactly two positive homogeneous steady states. Moreover, if k < min(k * 1 , ((m 2 − m 1 )/(m 1 µ 3 )) 2 ), both nontrivial positive homogeneous steady states are of type u − .
We show in Lemma 3.7 that the conditions of Theorem 2.3 are never satisfied if u is of type u − at some x ∈ I, but always satisfied if u is of type u 0 or u + for all x ∈ I. Hence, the latter type of solutions is constructed. Lemma 3.6 (Existence of steady states with jump discontinuities). Assume 2 √ k < m 1 < m 2 . Then, there exist infinitely many steady states (ũ,ṽ) ∈ (L ∞ (I) ∩ BV (I)) × (H 2 (I) ∩ C 1 (I)) which are of type u + or u 0 for all x ∈ I.
3.3. Stability of steady states for the quasi-steady state approximation. Following the same principle as in Section 2, we prove stability of steady states of the quasi-steady state approximation (3.3)-(3.5). In the next subsection, we conclude stability of steady states of model (3.1). In order to apply Theorem 2.3, we investigate how the entries of Jacobian matrix at a steady state depend on the branch: and, in particular, for u(x) = u + (v(x)), Applying Theorem 2.3 yields stability of the steady states constructed in Lemma 3.6. Assume D > 0 and 2 √ k < m 1 < m 2 . Then, there exist infinitely many (ε 0 , A)-stable steady states with jump discontinuity of model (3.3)-(3.5), which are for every x ∈ I either of type u 0 or u + .
To obtain DDI, it is left to prove that at least one of the spatially homogeneous steady states of type u − is stable under spatially constant perturbations, i.e. a stable steady state of the kinetic system. Instability of all spatially homogeneous steady states of type u − follows from [16] due to autocatalysis of the non-diffusive species. Existence of exactly two homogeneous steady states of type u − results from Lemma 3.5. If there exist two homogeneous steady states Furthermore, all spatially homogeneous steady states of type u − are unstable for k ≥ 0 and D > 0.
A numerically obtained solution, using the finite element library deal.ii [4], is presented in Figure 3.2. We observe destabilization of the spatially homogeneous steady state and convergence towards a pattern with jump discontinuity. Compartment u is discretized with cell-wise constant finite elements and v is discretized with cell-wise linear, globally continuous finite elements. The formulation is the weak formulation. The time-stepping scheme is Crank-Nicholson.

Stability of steady states for the unreduced model.
In this section we show that model (3.1) satisfies the conditions of the 'stability transfer' Theorem 2.6.
there exists a constant κ * in the interval and a positive constant δ * such that for all 0 ≤ δ < δ * and for all 0 < κ < κ * the following statements hold true: (1) system (3.1) has exactly two strictly positive spatially homogeneous steady states, is an unstable steady state of the kinetic system for (3.1) and of the original system (3.1).
is a stable steady state of the kinetic system of (3.1) and an unstable steady state of (3.1). (4) (0, 0, 0) is a stable steady state of the kinetic system for (3.1) and of the original system (3.1) . (5) system (3.1) has infinitely many (ε 0 , A)-stable steady states with jump discontinuity, with component u 1 being at x ∈ I of type u + or u 0 .

Numerical results: Effects of initial conditions and diffusion coefficient.
In the previous subsection, we showed that a spatially homogeneous steady state is destabilized and that infinitely many, discontinuous steady states are stable, both due to introduction of diffusion. In this subsection, we show some numerical results showing that the arising pattern depends on both, the initial conditions and the diffusion coefficient. In figure 3.3, numerical approximation of solutions to problem (3.3)-(3.5) are shown. We observe that for the same initial conditions and parameters of the kinetic system, the number of jump-type discontinuities varies according to the size of diffusion coefficient D.
A trend towards a higher number for smaller diffusion-coefficient can be observed and the shape of the arising pattern is different from the shape of the initial conditions. For large diffusion coefficient, a 'plateau' arises close to the spatial position on which the initial condition of u assumes its maximal value. Here, initial conditions seem to determine the position of 'plateaus' as can be observed in  To handle Theorem 2.3 we introduce the following notation for simplicity: . Then −L 1 is a sectorial operator on (L 2 (I)) 2 , i.e., there exist positive constants M 1 , κ 0 and ω ∈ (0, π/2) such that where B((L 2 (I)) 2 ) denotes the Banach space of all bounded linear operators on (L 2 (I)) 2 equipped with the operator norm. Hence L 1 generates the analytic semigroup e tL 1 on (L 2 (I)) 2 and satisfies the estimate Here M and κ 1 are positive constants independent of (φ, ψ).
Since Lemma 4.1 is obtained as a corollary to Lemma 4.2, we prove only Lemma 4.2. (Indeed, choosing α > 0 sufficiently large, we see that the matrix Proof of Lemma 4.2. Instead of −L 2 we study L 2 ; the estimate (4.5) is obtained by replacing λ with −λ in (4.41) below. Hence, we consider the nonhomogeneous equation for r 1 , r 2 , s ∈ L 2 (I) under homogeneous Neumann boundary conditions on ψ.
Proof of Sublemma 4.3. First we prove the boundedness of B. It is convenient to introduce the following notation: for all λ satisfying Re λ + 2κ ≥ 0 or |Im λ| ≥ Γ 0 , which implies the desired inequality: Second, we prove the coerciveness of B. Let us start with the case where |λ| is sufficiently large. From as |λ| → +∞, uniformly in x ∈ I. On the other hand, we have Therefore, we see that as |λ| → +∞, uniformly in x ∈ I. Since Re λ − a 33 (x) ≥ Re λ + 3κ and T 3 (x) < 0 by assumption (2.14), we conclude that there exists a positive constant Γ 1 ≥ Γ 0 such that whenever Re λ ≥ −2κ and |Im λ| ≥ Γ 1 .
We now turn to the proof of Theorem 2.4. In order to control the nonlinear terms, we first estimate the H 1 (I) norm of ψ, which gives us a bound on the L ∞ (I) norm of ψ by virtue of the Sobolev embedding theorem.
Put X(t, x) = |φ 1 (t, x)| 2 + |φ 2 (t, x)| 2 , Y (t, x) = |ψ(t, x)| 2 and Z(t) = ψ(t, ·) 2 H 1 (I) . So far we have established the following inequalities: Also, by the Sobolev embedding theorem, Here and in what follows, we use C to denote various positive constants independent of t, X(t), Y (t) and Z(t). From (4.58) we get It is straightforward to see that the right-hand side is equal to Hence, from (4.59) we conclude that (4.61) For any measurable set R ⊂ I, let Then from (4.57) and (4.60), we have Decomposing the integral over I into the sum of those over R and I \ R, and denoting the Lebesgue measure of R by µ(R), we see (4.64) From (4.61) and (4.64) it follows that From (4.63) we obtain that is, Therefore, m(t) = max 0≤τ ≤t (ξ R (τ ) + Z(τ )) satisfies m(t) ≤ C m(0) + m(t) 2 + µ(I \ R) . Now it is easy to show that, for any A satisfying A > max{C, 1}, if we choose ε 0 > 0 so small that then for any 0 < ε < ε 0 , as long as m(0) ≤ ε 2 and µ(I \R) < ε 4 , the inequality m(t) < Aε 2 holds for all t ≥ 0. This completes the proof of Theorem 2.4. Theorem 2.3 on two-component systems is proved in exactly the same way and hence we omit it.

Quasi-steady state reduction of the ODE subsystem.
In this subsection we address the question of whether Diffusion Driven Instability can be investigated or not based on the quasi-steady state approximation. For δ > 0, consider a system of type Then by the implicit function theorem there exists a smooth function u * Hence, (u 1 , u 2 , v) is a constant steady state of (4.71)-(4.73), and (u 1 , v) is a constant steady state of its quasi-steady state reduction:  (u 1 , v), then there exists a positive number δ * such that the unreduced system (4.71)-(4.74) exhibits DDI at (u 1 , u 2 , v) as long as 0 < δ < δ * .
Let B = (b ij ) 1≤i,j≤2 denote the Jacobian matrix around the equilibrium (u 1 , v) of the kinetic system for the reduced system, i.e. (4.77)-(4.78) with D = 0. Let A δ = (a δ ij ) 1≤i,j≤3 denote the Jacobian matrix around (u 1 , u 2 , v) of the kinetic system for the unreduced system, i.e., (4.71)-(4.73) with D = 0. We prove this proposition by demonstrating the following two lemmas, both of which are algebraic in nature. First, in Lemma 4.7 we prove that if δ is sufficiently small then two of the three eigenvalues of A δ remain in the neighborhood of the two eigenvalues of B, while one of the eigenvalues of A δ tends to −∞ as δ → 0, which in particular shows that (u 1 , u 2 , v) is stable under spatially homogeneous disturbance, provided that δ is sufficiently small. Second, in Lemma 4.9 we show that DDI of the reduced system (4.77)-(4.79) implies instabilty of (u 1 , u 2 , v) for all D > 0. Combining these lemma proves Proposition 4.6.
Remark. If we assume, instead of (4.77)-(4.78), that both species diffuse and consider the system , v 0 ) then the constant steady state (u 1 , v) is destabilized only for D > D cr with D cr being a positive number. This is quite different from the case D 1 = 0 in which (u 1 , v) is always unstable for any D > 0. The same observation applies to the three component system (4.71)-(4.73). See [2] for the case where all three species diffuse.

Proof of statements concerning the example system
First, we prove uniform boundedness of solutions and classification of different branches.
Proof of Lemma 3.1. To obtain bounds on u and v, we see that for v ≥ 0 The first inequality of (5.1) proves that u remains positive if v ≥ 0.
On the other hand, since m 2 > m 1 , it holds that Therefore, if both (5.33) and (5.34) are satisfied, and hence (u, l(u)) and (u, v f,g (u)) intersect twice in R. This establishes the existence of two homogeneous steady states for sufficiently small k.
To prove that both steady states are of type u − , we use the fact that the steady states component u is less than (m 1 µ 3 )/(m 2 − m 1 ). Since v fr=0 (u) = −1 + m 1 u/(1 + ku 2 ) is strictly increasing on (0, 1/ √ k), it holds that at any steady state u * . By the characterization Lemma 5.1, both spatially homogeneous steady states are of type ( Existence of spatially inhomogeneous steady states can be proved similarly to [20] or [18]. To apply this work, we solve f r (u, v) = 0 in u and substitute this into The shape of g r (u(v), v) for different branches of the solution u(v) of f r (u, v) = 0 is shown in Fig. 5  g r (u ± (0), 0) > 0, (5.44) If there exist two positive spatially homogeneous steady states of type u − , it holds that uniformly on any finite interval v ∈ (0, c) as k ց 0.
Proof. Inserting v = 0 into the right-hand side (5.43) leads to We differentiate g r (u + (v), v) with respect to v and obtain for v ≥ 0, hence (5.45) and (5.46) hold true.
If there exist two spatially homogeneous steady states of type u − , then g r (u ± (v r ), v r ) must be positive because g r (u − (0), 0) > 0 and g r (u − (v), v) has exactly two roots of order one in the interval 0 < v < v r . Hence (5.47) is proved. Inequality (5.48) follows immediately from u + (0) ≤ m 1 /k and (5.50).
To see (5.49), we differentiate g r (u − (v), v) with respect to v: The limits for k → 0 of the first and the last term are clear. To obtain the limit of the second term, we use l'Hôpital's rule and obtain which proves (5.49).
Using this knowledge, especially about the sign of g r , we can prove existence of steady states with jump discontinuity. The proof is oriented on the shooting method and follows the principle presented in [20], [3] and [18]. However, for completeness we state Proof of Lemma 3.6. The concept of the proof is illustrated in Figure 5.3. Step 1 . We construct steady-states of (3.3)-(3.5) by using two building blocks. The first one is the solution of the initial value problem which is solved explicitly: We denote this solution by v(x; a). The second one is the solution of the initial value problem We denote this solution by w(x; b), where the initial data w 0 is taken from the interval (0, v r ). Denote −g r (u + (w), w) by h(w), which is defined for w ∈ [0, v r ]. Then we know from Lemma 5.2 that Clearly w(x) is positive and monotone decreasing in a certain interval 0 ≤ x < x b where w(x b ) = 0. Note that x b < 2b/|h(b)| since w(x) = b + h(w(θx))x 2 /2 with 0 < θ < 1.
Step 2 . We would like to connect w(x; b) with v(x; a), or v(x; a) with w(x; b), so that the boundary conditions are satisfied. Fix a number c in the interval (0, b) and let x c = x c (c, b) ∈ (0, x b ) be such that w(x c ; b) = c. Now we connect w(·; b) to v(·; a) at w = c as follows: We would like to choose an appropriate a ∈ (0, c) and a y c = y c (c, a) > 0 so that the matching conditions is satisfied. Indeed, if (5.60) holds, then there is a unique ξ c < 0 such that tanh ξ c = w ′ (x c ; b)/( √ µ 3 c). Then we obtain y c = −ξ c / √ µ 3 and a = c/ cosh ξ c . We say that w(x; b) is switchable to v(x) at w = c if (5.59) is satisfied. We remark that for each b there exists a unique c * = c * (b) ∈ (0, b) such that w(x; b) is switchable to v(x) at c if and only if c * (b) < c < b. This is proved easily by using the following three facts: (w ′ /w) ′ = h(w)/w − (w ′ /w) 2 < 0, w ′ (0)/w(0) = 0 and lim x↓xw 0 w ′ (x)/w(x) = −∞.
Next, fix a number γ ∈ (a, v r ) and let y c = y c (γ, a) > 0 be such that v(y c ; a) = γ. We extend w(x; b) to −x b ≤ x ≤ 0 by defining w(x; b) = w(−x; b) for x ∈ [−x b , 0]. This time we would like to find b ∈ (a, v r ) and x c = x c (γ, b) ∈ (0, x b ) such that which is equivalent to the following conditions: We say that v(x; a) is switchable to w(x) at v = γ if there exist b and x c (γ, b) satisfying the matching condition (5.61). Likewise for w(x; b), we shall prove that there exists a unique c * = c * (a) such that v(x; a) is switchable to w(x) at v = γ if and only if a < γ < c * (a).
Step 3 . Assuming the existence of c * (a), we describe how to construct a steady state of (3.3)-(3.5). First we choose b ∈ (0, v r ) arbitrarily and then choose any c ∈ (c * (b), b). We now obtain a and y c (c, a) which satisfy (5.59). Let us put Clearly W 1 (x) = W ((x c (c, b) + y c (c, a))x) is a weak solution of (3.3)-(3.5) for D = 1/(x c + y c ) 2 . Next, choose γ ∈ (a, c * (a)) arbitrarily and find a pair (b 1 , x c (γ, b 1 )) satisfying (5.62). We extend W (x) to the interval [0, x c (c, b) + y c (c, a) + y c (γ, a) + y c (γ, b 1 )] by putting , then this gives rise to a weak solution of (3.3)-(3.5) for D = 1/L 2 3 , which has two switching points. We can repeat these procedures to obtain steady states of (3.3)-(3.5) with an arbitrary number of switching points. In the same manner, starting with v(x; a), one obtains another type of steady states with any number of switching points.
Step 5 . Finally we examine the range of D which has monotone decreasing solutions. Note that c * is characterized by the property that |w ′ (x c (c * , b); b)| = √ µ 3 c * . Therefore, tanh ξ c → 1, i.e., ξ c → −∞ as c ↓ c * (b), so that y c → +∞ and a → 0. This means that D = 1/(x c + y c ) 2 → 0 as c ↓ c * (b). On the other hand, as c ↑ b, we have (i) x c → 0 and (ii) w ′ (x c ; b) → 0 and hence ξ c → 0. This implies that y c → 0 and a → c as c ↑ 0. In particular, D = 1/(x c + y c ) 2 → +∞ as c ↑ b. Observing that x c and y c are continuous function of c, we conclude that for each D > 0 and for each b ∈ (0, v r ) there exists a monotone decreasing solution W 1 (x; b) of (3.3)-(3.5) as stated in Step 3.
we hence obtain infinitely many steady states for each D > 0. (The arguments above also prove the existence of infinitely many steady states which are not monotone in x for each D > 0.) This completes the proof of Lemma 3.6.
x w v Now that we have established the existence of spatially homogeneous steady states and spatially inhomogeneous steady states, we turn to the proof of the lemmas stating that the conditions of Theorem 2.3 are satisfied: Proof of Lemma 3.7. First, we calculate the Jacobian matrix of the kinetic system for given arbitrary (u, v): Combining (5.67) with u − (v) ≤ 1/ √ k ≤ u + (v) yields the result for b 11 . We investigate b 21 : The continuity of b 21 and b 11 (v r ) = 0 imply the result for b 21 .