Hysteresis-driven pattern formation in reaction-diffusion-ODE systems

The paper is devoted to analysis of far-from-equilibrium pattern formation in a system of a reaction-diffusion equation and an ordinary differential equation (ODE). Such systems arise in modeling of interactions between cellular processes and diffusing growth factors. Pattern formation results from hysteresis in the dependence of the quasi-stationary solution of the ODE on the diffusive component. Bistability alone, without hysteresis, does not result in stable patterns. We provide a systematic description of the hysteresis-driven stationary solutions, which may be monotone, periodic or irregular. We prove existence of infinitely many stationary solutions with jump discontinuity and their asymptotic stability for a certain class of reaction-diffusion-ODE systems. Nonlinear stability is proved using direct estimates of the model nonlinearities and properties of the strongly continuous diffusion semigroup.


1.
Introduction. This paper is devoted to mathematical analysis of hysteresisdriven pattern formation. Mathematically, a pattern is a stable spatially inhomogeneous stationary solution of an evolution equation describing the development of a system.
Most mathematical models of biological or ecological patterns are based on reaction-diffusion equations with the Turing mechanism (diffusion-driven instability) [30] which can explain de novo pattern formation [11]. Diffusion-driven instability is related to a local behavior of solutions in the neighborhood of a constant stationary solution that is destabilized by diffusion. It may lead to emergence of close to the equilibrium patterns that are stable continuous and spatially periodic structures around the destabilized constant equilibrium. The Turing concept became a paradigm for biological pattern formation and led to development of numerous theoretical models, though its biological validation has remained elusive. Another type of patterns are far-from-equilibrium patterns that arise in systems with multiple constant stationary solutions that usually exhibit bistability and hysteresis. Bistability (or more generally multistability) is a phenomenon of systems toggling among two (or more) stable equilibrium points. Usually, between the stable equilibria, there is an unstable intermediate equilibrium working as a threshold. Hysteresis refers to a system where the output does not only depend on the input, but also on the history of the system that may lead to different responses to the same input. Bistability and hysteresis are mechanisms important for biological systems, which generate oscillations and switching between discrete states, and integrate transient stimuli [1]. They can create an all-or-none response and transform a graded input signal into a discontinuous output. Moreover, hysteresis may be responsible for irreversibility, which is of particular importance in developmental processes [4].
To understand the role of bistability and hysteresis in pattern formation, we propose a prototype mathematical model based on a reaction-diffusion equation coupled to an ordinary differential equation (ODE) with hysteresis. Such models have recently been employed in biology [8,10,16,21,26,31] and in ecology [15]. They allow describing cell-to-cell communication coupled to intracellular dynamics as well as binding of a diffusing molecule to an immobile substrate [12] or receptor [16,17]. Reaction-diffusion-ODE 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 [23,18]. Furthermore, they may occur as reduced problems in the analysis of reaction-diffusion systems with small diffusion coefficients. For example, in [24] certain solutions of the reduced system have been used to construct large-amplitude solutions of a system of two reaction-diffusion equations having in interior transition layer, with one diffusion coefficient being sufficiently small and the other sufficiently large.
Reaction-diffusion-ODE models with hysteresis may additionally exhibit diffusion-driven instability what leads to an interesting dynamics as shown recently in Ref. [6,13]. In this paper, to streamline the presented analysis and explore the role of bistability and hysteresis, we focus on a model that does not exhibit the Turing mechanism. Instead, it has two spatially homogenous stationary solutions that are always locally asymptotically stable. Diffusion acts to average different states and is the cause of spatio-temporal patterns. We prove that bistability without the hysteresis effect is not sufficient for existence of stable spatially heterogenous patterns. Moreover, we provide a systematic description of stationary solutions that may differ from those of the usual reaction-diffusion systems [20,19,7]. In particular, we show a co-existence of infinite number of stationary solutions that exhibit jump-discontinuities in non-diffusing variables. In applications, discontinuities are useful for defining sharp boundaries of gene expression pattern [8].
This paper is organized as follows: In Section 2, we set up our problem and show that there is no spatially heterogeneous stable stationary solution if the nonlinearity does not exhibit hysteresis (Theorem 2.2). In the remainder of the paper we consider only nonlinearities with hysteresis. Section 3 is devoted to the existence (Theorem 3.2) and the asymptotic stability (Theorem 3.7 and Corollary 1) of monotone increasing (or decreasing) stationary solutions with jump discontinuity. In Section 4 we consider the location of the point of discontinuity (Theorem 4.3). Finally in Section 5 we consider the existence and stability of spatially periodic stationary solutions (Corollary 2) and then proceed to spatially irregular stationary solutions. We prove the uniqueness (Proposition 4) and the stability (Theorem 5.4) of spatially irregular stationary solutions, and discuss the existence (Example 5.3) of such solutions.
2. Problem setting. In this section we introduce a prototype reaction-diffusion-ODE model with hysteresis and discuss its basic properties. We consider a onedimensional problem supplemented with the homogeneous Neumann boundary condition for u and initial conditions To focus on a simple structure of nonlinearities exhibiting hysteresis, we consider a specific choice of the kinetic functions where α, β are positive constants and p(v) is a polynomial of degree three such that p(v) → +∞ as v → +∞. The parameters are chosen such that there exist three intersection points of f = 0 and g = 0 with non-negative coordinates. The above assumptions yield a bistable system with or without the hysteresis effect, see Fig. 1. We define: Case 1 (Bistability without hysteresis): Polynomial p(v) is assumed to be monotone increasing (cf. Fig. 1a). Case 2 (Hysteresis): Polynomial p(v) is assumed to be non-monotone. In this case, we denote by Moreover, we assume that the coordinates of H and T are positive and that lim v→+∞ p(v) = +∞ holds (cf. Fig. 1b). Here, by bistability it is meant that the kinetic system has two asymptotically stable equilibria. We check that the system is bistable in both cases but only a system with hysteresis allows emergence of stable patterns (spatially heterogenous stationary solutions).  Lemma 2.1. Assume that there exist three spatially homogeneous stationary solutions S 0 , S 1 and S 2 , 0 < u 1 < u 2 , 0 < v 1 < v 2 , of system (1)-(3) with nonlinearities given by (5). Then, S 0 and S 2 are asymptotically stable, while S 1 is unstable.
Proof. First, we check bistability of the kinetic system (6), i.e., asymptotic stability of the spatially homogenous stationary solutions S 0 and S 2 of system (1)-(3) under constant perturbations.
Let J(ū,v) be Jacobi-matrix of function (f (u, v), g(u, v)) at a stationary solution (ū,v). As α β is the slope of f (u, v) = 0 solved with respect to v, we obtain at the solutions S 0 and S 2 , we conclude about their asymptotic stability. For S 1 we have to distinguish between Case 1 and Case 2. In Case 2 it holds p (v 1 ) < 0, hence det J(S 1 ) < 0. In Case 1, since p (v 1 ) < α β , it also holds that det J(S 1 ) < 0. Therefore, linearization at S 1 has one positive and one negative eigenvalue.
The stability properties do not change for spatially heterogeneous perturbations, since diffusion causes a destabilization of stable spatially constant stationary solution (diffusion-driven instability) in a reaction-diffusion-ODE system if and only if the ODE component exhibits autocatalysis, i.e., g v (ū,v) > 0, see Theorem 2.11 in Ref. [19]. It is not the case for model (6), since p (0) > 0 and p (v 2 ) > 0 and hence g v (ū,v) < 0, independently of the choice of admissible model parameters. Consequently, S 0 and S 2 are asymptotically stable as spatially homogenous stationary solutions of system (1)-(3).
Next, we observe that bistability is not sufficient for emergence of stable patterns. Let (U, V ) be a stationary solution of system (1)-(3), Since, in Case 1, the polynomial p is monotone increasing, equation (8) can be uniquely solved with respect to V . We obtain an elliptic boundary value problem where h(U ) satisfies g(U, h(U )) = 0. Since f (U, h(U )) is a smooth function, equation (10) has classical solutions U ∈ C 2 [0, 1] , and hence also V = h(U ) ∈ C 2 [0, 1] ). Such solutions can be explicitly constructed by the phase plane analysis (cf. [2]). Their instability follows from analysis of the linearized problem using a similar argument as in case of a scalar reaction-diffusion equation. It results in Theorem 2.2. There exists no stable non-homogeneous stationary solution (U, V ) of system (1)-(4) in Case 1.
Details of the proof are provided in Appendix for the sake of completeness.
3. Patterns with jump-discontinuity in the model with hysteresis. In this section, we show that all stationary spatially heterogeneous solutions of problem (7)-(9) with parameters satisfying Case 2 have a jump-discontinuity in V (x). Then, we focus on construction and analysis of monotone solutions.
3.1. Construction of discontinuous steady-states. In Case 2, polynomial p is non-monotone and the stationary algebraic equation (8) cannot be solved uniquely with respect to V . We denote by V = h H (U ) the lower solution branch and by V = h T (U ) the upper one (cf. Fig. 1b). There is a third solution branch h 0 in the middle, but since the spatially constant steady state S 1 located on this branch is unstable, we do not expect to obtain stable stationary solutions using this branch. Indeed, we show in Section 3.2 that the stability condition is not fulfilled for V = h 0 (U ). The stationary problem reads where i = H or i = T . Using phase plane analysis, we conclude that there is no solution of this problem neither for i = H nor i = T , which means that there exists no stationary spatially heterogenous C 2 -solutions.
In the next step, we observe that the phase planes associated to the equation (11) for i = H and i = T overlap for u ∈ (u T , u H ) and search for solutions with jump discontinuity (cf. Fig. 2A). Heuristically, to construct a solution, we select a valueū ∈ (u T , u H ) and "glue" the phase planes together atū.
It is convenient to define the following functions: Let qū denote the function with discontinuity atū defined by qū(u) = q H (u) when u ≤ū q T (u) when u >ū.
We search for a weak solution of the boundary value problem i.e. a function U ∈ H 1 (0, 1) satisfying 1]) by virtue of the standard elliptic regularity theory.
Definition 3.1. A pair of functions U, V is called a solution with jump discontinuity atū of the stationary problem (7)-(9), if U ∈ H 1 (0, 1) is a weak solution of equation (13). Function V ∈ L ∞ (0, 1) is given for almost all x ∈ [0, 1] by The value 0 <x < 1 such that U (x) =ū is called the layer position of the solution U, V .
Proof. We apply the classical method of phase plane analysis [2] and the analysis of time-maps [28]. The first integral of equation (13) reads where the function is called the potential and corresponds in a physical context to the potential energy of the system. The constant E ∈ R is arbitrary and corresponds to the total energy of the system. Problem (13) can be written as a boundary value problem for the system of first order equations A monotone increasing solution of problem (13) is given by a solution of system (16) satisfying where u 0 , u e ∈ R are such that u 0 < u e and W (x) > 0 for all x ∈ (0, 1), see Fig.  2A. It follows from the first integral (14) where the constant is determined by the Neumann boundary condition at x = 0 and x = 1. Hence, because of (18), a necessary condition for the solution U (x) requires that Qū(U (x)) ≤ E holds for x ∈ [0, 1]. By definition Qū fulfils Q ū (u) = qū(u) and Qū(0) = 0. It is continuous for all u but not differentiable inū. Assuming thatū < u 2 , we conclude that the potential has local maxima at u = 0 and at u = u 2 , because qū(0) = f (0, h H (0)) = f (0, 0) = 0, and qū(u 2 ) = f (u 2 , h H (u 2 )) = f (u 2 , v 2 ) = 0, respectively. Furthermore, it has a local minimum at u =ū. Therefore, forū ∈ u T , min(u H , u 2 ) , it is possible to find values 0 < u 0 <ū < u e < u 2 fulfilling condition (19) and such that Qū(u) < Qū(u 0 ) for all u ∈ (u 0 , u e ).
Equation (18) yields We split this integral at the minimum of the potential and denote (21) T 1 u (u 0 ) and T 2 u (u e ) are called time-maps, since our variable x corresponds to time in classical models from mechanics. For a forward orbit in the phase plane representation of the system (16), with qū(u) = f (u, h H (u)) and starting at (u 0 , 0), T 1 u (u 0 ) is the "time" x required to reach the u =ū axis for the first time. In distinction to this, T 2 u (u e ) is the "time" x, required for a backward orbit in the phase plane representation system (16), with qū(u) = f (u, h T (u)) and starting at (u e , 0), to reach the u =ū axis for the first time.
Proof. The proofs of well-definedness and continuity are standard and can be found in [28]. We sketch the proof of surjectivity for completeness.
The integrand of T 1 u has a singularity at U = u 0 . Since Q ū (u 0 ) = q H (u 0 ) = 0, using the Taylor expansion we conclude that the integral is convergent.
Moreover, the Taylor expansion at 0 and the fact that Qū has a local maximum at 0 (since q H (0) = 0 and q H (0) < 0) allow to estimate For u 0 and u close toū, we use the Taylor expansion of Qū at u 0 and obtain with η ∈ (u 0 , u). Thus, we calculate the limit The corresponding results for T 2 u (u e ) are obtained in a similar fashion.
The maps T 1 u and T 2 u , resp., are defined for all u 0 ∈ (0,ū) and u e ∈ (ū, u 2 ). However, Tū(u 0 ) is defined only for u 0 such that there exists a u e satisfying (19). Thus, the domain of definition of Tū(u 0 ) depends on the local maxima of the potential. Therefore, if Qū(u 2 ) ≥ Qū(0) = 0, we denote u min = 0 andū < u max ≤ u 2 is the solution of Qū(u max ) = 0. If Qū(u 2 ) ≤ 0, we denote u max = u 2 and 0 ≤ u min <ū is the solution of Qū(u min ) = Qū(u 2 ). Proof. Clearly Tū is continuous as a sum and the composition of continuous functions. Because of the continuity of Qū we conclude that u 0 →ū implies u e →ū and, therefore, Finally, we know that that either u min = 0 or u max = u 2 holds true and therefore lim u0→umin because either the first or the second limit is infinite.
The previous lemma provides existence of a monotone increasing solution of problem (13) for any diffusion coefficient 1 γ > 0. Next, we show uniqueness of this solution, which is a consequence of the monotonicity of the time-map Tū. To prove this, we need the following representation of the derivatives of the time-maps.
Lemma 3.5. The time-maps are differentiable and the derivatives have the following form Proof. This formula has been shown in (2.18) and (2.19) of Ref. [14] under the assumption qū ∈ C 0 with piecewise continuous q ū and qū(ū) = 0. We adapt that proof to requirements satisfied by our problem. First, we notice that the derivative of T 2 u (u e ) can be written in the form where since the respective proof provided in Theorem 1 of [14] requires only existence of q T for u >ū. Now, we observe that it holds, for anyū

Using integration by parts of the second integral with
.
Finally, letting a →ū and b → u e , we obtain Together with representation (26), this integral yields the formula in the lemma. In a similar fashion we deduce the result for

Proposition 1.
Consider the stationary problem (7)-(9) in Case 2. Then the derivatives of T 1 u and T 2 u have the following signs dT 1 Therefore, for all jumpsū ∈ u T , min(u H , u 2 ) , the derivative of the time-map Tū is negative, i.e., Proof. To show the negativity of with a function l H defined by ) and the jumpū = 3.1. Here, it holds Qū(u2) < 0 and umin = 0.8957. We see that T 1 u is decreasing, whereas T 2 u is increasing, which leads to d du 0 Tū(u0) < 0. Furthermore, we observe that limu 0 →ū Tū(u0) = 0 and limu 0 →u min Tū(u0) = ∞ .
For this purpose, we multiply l H by the square of q H and calculate the derivative which leads together with To see that this expression is negative, we observe that by differentiating equation (19) with respect to u 0 . Therefore, we obtain 3.2. Stability analysis. A standard approach to showing stability of stationary solutions is based on linear stability analysis which in case of discontinuous solutions requires careful treatment, e.g., considering a special topology that allows excluding the discontinuity points [3,6]. In this section, we show the stability of stationary solutions by applying direct estimates. Numerical simulations of the basic model with hysteresis considered here suggest that the selected pattern depends on the initial conditions. More precisely, the ultimate stationary solution depends on the position of the initial functions relative to the separatrix of the saddle point of the kinetic system. For this reason it is suitable to consider the stability in L ∞ (0, 1) sense. A perturbation which is small in L ∞ (0, 1) cannot move a point on the stationary solution U, V lying on one side of the separatrix to the other side. In contrast, a perturbation which is small in L 2 (0, 1) may have high values on some small interval leading to another stationary solution. Thus, a stationary solutions which is stable in L ∞ (0, 1) sense might be unstable with respect to L 2 (0, 1)-perturbations, as illustrated in the following example.
We remark here that such strong dependence of the pattern on initial data is related to the specific simple nonlinearities we use to model hysteresis, i.e., the existence of two spatially homogeneous stationary solutions that are always stable. As shown recently in Ref. [6], hysteresis can be also exhibited by models with diffusion-driven instability of one of the spatially constant stationary solutions and, in such a case, pattern selection is more robust and in case of small perturbations depends on the diffusive scaling of the system (a phenomenon similar to the selection of Turing patterns).
Proof. First we choose a positive constant δ so small that the following inequalities are satisfied: By applying (34) to (32), we obtain Hence, Since p (V (x)) ≥ K in [0, 1] we obtain from (33) with the help of (31) that Since Ψ(t) is monotone increasing in t and R * ≤ R 0 , we see that Therefore, for any 0 ≤ t ≤ τ < T . Taking the supremum of the left-hand side over 0 ≤ t ≤ τ , and then changing τ to t, we conclude that Now consider the quadratic function Let D δ denote the discriminant of ν δ (y): By (38) we know that D δ > 0. Note also that the coefficient of the linear term in ν δ (y) is negative. Hence, the equation ν δ (y) = 0 has two distinct positive roots ρ 1,δ < ρ 2,δ . We observe that ν δ (Ψ(0)) > 0 whenever Moreover, we claim that Ψ(0) < ρ 1,δ . To see this, we note that 1 − α/[(β − δ)(K − δ)] < 2. Hence from (38) we obtain On the other hand, since ν δ (y) > 0 if and only if y < ρ 1,δ or y > ρ 2,δ , we see that Ψ(0) = ψ 0 ∞ < ρ 1,δ or Ψ(0) > ρ 2,δ , and the latter possibility is ruled out by the estimate above. From (41) we have ν δ (Ψ(t)) ≥ 0 for all 0 ≤ t < T , where T is the maximum existence time. The continuity of Ψ(t) and Ψ(0) < ρ 1,δ imply that Ψ(t) < ρ 1,δ for t ∈ [0, T ), hence (42) By virtue of (40) we obtain also that In view of (K − 2δ due to (36). Hence, (42) implies, in particular, ψ(t, ·) ∞ < 1 for all t ∈ [0, T ). It is now easy to verify that T = +∞, and this completes the proof of the theorem by choosing C * = max{C 1 , C 2 }. The stability condition seems to be difficult to check, because it depends on the stationary solution. However, we can translate it into a condition for the jump.  Proof. By definition of a stationary solution with jumpū, the function V (x) is given by h H (U (x)) for x ∈ [0, 1] fulfilling U (x) ≤ū. Therefore, for those x the requirement p (V (x)) > α β is equivalent to From this we see that q H (U (x)) < 0 holds exactly when U (x) < u 0 H . Therefore, u < u 0 H is necessary and sufficient for fulfilling the stability condition. Similarly, for x ∈ [0, 1] fulfilling U (x) >ū it holds V (x) = h T (U (x)) and the stability condition requires It holds q T (U (x)) < 0 if u 0 T < U (x), hence we deduce u 0 T <ū. Thus, the stability condition is fulfilled.
Definition 3.8. Kinetic functions such that u 0 T < u 0 H holds are called admissible. In this situation, valuesū ∈ (u 0 T , u 0 H ) are called admissible jumps. Consequently, Corollary 1 says that for an admissible kinetic function all monotone increasing stationary solutions having an admissible jumpū ∈ (u 0 T , u 0 H ) are stable.
The disadvantage of the above construction of patterns is the fact that we cannot directly set upū in a simulation of the time-dependent system. It is rather the choice of the initial condition by which we may influence the final pattern.   (15) with respect toū is given by Proof. Forū > u, it holds Q(ū, u) = u 0 q H (ũ)dũ which is independent ofū, whereas forū < u it holds Q(ū, u) = ū 0 q H (ũ)dũ− ū u q T (ũ)dũ. By calculating the derivative of these representations with respect toū the result follows immediately.
This can be solved with respect to du0 dū (ū) and leads together with iv) of this proposition to the assertion. vi) We apply the chain rule to the defining equation u e (ū) = u e (ū, u 0 (ū)).
We skip the details because we do not use them later on.

4.2.
Dependence of the layer position on the jump.
Definition 4.2. The layer positionx(ū) ∈ (0, 1) of a monotone increasing solution of equation (10) is defined byx This is the position where the monotone increasing solution of (13) switches from the phase plane corresponding to 1 γ U xx (x) + q H (U (x)) = 0 to that corresponding to 1 γ U xx (x) + q T (U (x)) = 0. Proposition 3. The layer positionx(ū) is continuously differentiable and its first derivative with respect toū is given by .
For admissible kinetic functions, the relation dx dū (ū) < 0 holds for all admissible jumpsū ∈ (u 0 T , u 0 H ). Proof. The functionx(ū) is continuously differentiable as a composition of such functions. Applying the chain rule to T 1 (ū, u 0 (ū)) yields dx dū where we used Proposition 2 i) & v) and equation (28). Also, for simplicity of the exposition, we have abbreviated For better readability, we write from now on only u 0 and u e and omit the dependencies onū. We multiply equation (45) by the denominator of the right-hand side q T (u e ) ∂T1 ∂u0 (ū, u 0 )+q H (u 0 ) ∂T2 ∂ue (ū, u e ) which is negative (cf. equation (28)). Therefore, showing the negativity of d dūx (ū) is equivalent to showing the positivity of (46) In Lemma 3.5, it is shown that the derivatives of the time-maps with respect to u 0 and u e , resp., have the form where we denote by int H and int T the functions Using these representations, we reformulate expression (46) Remarking that q H (u 0 )q T (u e ) is negative, the whole expression is positive if the term in the parentheses is negative. This is the case if it holds For showing the first estimate in (47), we investigate in detail the integral int H . Therefore, we split the integrand into a suitable sum of two summands.
We calculate the integral int 1 H using the fundamental theorem of calculus For the calculation of int 2 H , we recall thatū < u 0 H implies that q H (u) is negative for all u <ū (cf. Fig. 1B). Thus, the integrand is negative and we estimate We use here the estimate which holds true because Q(ū, u) is monotone decreasing in u for u ∈ (0,ū). Furthermore, it holds 0 > q H (u 0 ) > q H (ū), because q H is decreasing for u ≤ū. This yields q H (ū) q H (u0) > 1. Hence, we obtain the estimate (47) by putting the results (48) and (49) together as follows To show the second estimate in equation (47), we argue in the same way, which accomplishes the proof of this proposition.

Corollary 2.
For every jumpū ∈ u T , min(u H , u 2 ) and every diffusion coefficient 1 γ > 0, there are exactly two solutions U of problem (13), which are periodic of mode k. Restricted to the interval [0, 1 k ] one of them is monotone increasing, whereas the other one is monotone decreasing. If the kinetic functions are admissible and the jumpū is admissible, then these periodic stationary solutions are asymptotically stable.
Proof. The proof is standard and can be found, e.g., in [20]. It basically consists of using a monotone solutionŨ of problem (13) with jumpū for the diffusion coefficient which we continue periodically for x ∈ [0, 1] by formula (52), the function U is periodic of mode k by construction and a solution of (13) with jumpū for the diffusion coefficient 1 γ .
Additionally to periodic patterns, the model (1)-(3) admits another class of stationary solutions. This is a consequence from the observation that the phase planes of the equations 1 γ U xx + q H (U ) = 0 and 1 γ U xx + q T (U ) = 0 are overlapping. For a periodic solution the switch between these phase planes takes always place at the same valueū, as one can see for the blue trajectory in Fig. 6. However, a similar construction of discontinuous patterns can be performed with switches at different valuesū 1 ,ū 2 , . . . as one can see for the red trajectory in Fig. 6. Then, there exist sub-intervals of [0, 1] such that the solution restricted to each of them is a monotone stationary solution with jump atū 1 ,ū 2 , . . . , respectively. x 0 = 0 < x 1 < · · · < x i < x i+1 < · · · < x k−1 < 1 = x k and jumpsū i ∈ u T , min(u H , u 2 ) , for i = 1, . . . , k such that the restriction of U to the intervals [x i−1 , x i ] are alternately monotone increasing and monotone decreasing. Thereby, for is a monotone increasing (resp. decreasing) solution of the equation The Vcomponent of an irregular solution is given by We emphasise here that the function U is continuous, in particular at the values x i for i = 1, . . . , k − 1. Hence, this means that the ending value ofŨ i and the starting value ofŨ i+1 coincide. Using the notatioñ We recall that the time-map T (ū, u 0 ) is initially defined as the "time" x a monotone increasing stationary solution with jump atū needs to connect 0 < u 0 <ū withū < u e (u 0 ,ū) < u 2 . However, this can be extended for monotone decreasing solutions that connectū < u 0 < u 2 with 0 < u e (u 0 ,ū) <ū. Here, it holds T (ū, u 0 ) = T 2 (ū, u 0 ) + T 1 (ū, u e (u 0 ,ū)), because T 1 (ū, ·) is only defined for values smaller than u and T 2 (ū, ·) for values larger thanū. Thus, for proving the proposition, it will be enough to show that for a fixed set of jumpsū 1 , . . . ,ū k and prescribed monotonicity on the first sub-interval, there is only one possible partition of the interval [0, 1].

ALEXANDRA KÖTHE, ANNA MARCINIAK-CZOCHRA AND IZUMI TAKAGI
Here, i = 1, 3, 5, · · · ≤ k runs through all odd numbers. We start at i = 1 and use these relations together, which successively yields Finally, we calculate the derivative of the time-map ∂ ∂u 1 This proves the uniqueness of an irregular solution having jumps atū 1 ,ū 2 , . . . ,ū k which is monotone increasing on [0, x 1 ]. The second case is treated in a similar fashion by calculating the derivative of the time-map with respect to u 2 0 .
whereŨ 1 andŨ 2 are monotone solutions of equation (53). We require U (x) to be continuous and therefore the condition has to be fulfilled, because both values correspond to U (x 1 ).
Generalizing this, the existence of an irregular solutions for a prescribed set of jumps u 1 , . . . ,ū k depends essentially on the existence of values u 1 0 , . . . , u k 0 related to each other by equation (55) and such that equation (57) holds. Depending on the signs of Q(ū i , u 2 ) and Q(ū i+1 , u 2 ) there might be no solutions to equation (55) considered for i and i + 1 at the same time which fulfills u i+1 0 = u i e , when γ is larger than some γ max .
The value γ max depends on the relative position of the jumps. When allū i are equal, the solution is periodic and exists for all γ. Thus, the closer theū i are to each other, the larger is γ max . In particular, when all jumps are close to u * , the value γ max gets larger, because then Q(ū i , u 2 ) ≈ 0 for all i and the possible range for u i 0 is large. Proof. This theorem follows from Theorem 3.7, we only have to show that the stability condition (35) is fulfilled. As it holds essinf x∈[0,1] p (V (x)) = min i∈{1,...,k} it suffices to prove the stability condition essinf x∈[x i−1 ,x i ] p (V (x)) > α β for every subinterval [x i−1 , x i ]. Using thatū i ∈ (u 0 T , u 0 H ), it holds q H (U (x)) < 0 and q T (U (x)) < 0 for those x such that U (x) ≤ū and U (x) > 0, respectively (cf. Fig.  1). Remembering that h H and h T are local inverses of the polynomial p and using the chain rule yields q i (U (x)) = α 1/p (h i (U (x))) − α β for i = H, T . This yields the desired stability condition .
Finally, we would like to address the question if a certain initial function will lead to a stable irregular pattern. When an initial function is crossing the separatrix at several points, then these points set up the layer positions of the final pattern, provided that there is a stable pattern with such layer positions. So the problem can be reformulated as follows: Is there an irregular solution having a prescribed set of layer positions? Unfortunately, we are not able to deduce a formula for the partition x 1 , . . . , x k−1 , when we only know the layer positions, but not the jumps. Thus, we cannot apply Theorem 4.3, because the interval I 0 depends on the diffusion coefficient and it changes if we consider it on a subinterval [x i−1 , x i ] instead of [0, 1].
However, under the conditions which we have encountered at the end of section 4.3, the existence of a stable irregular solution having a prescribed set of layer positions is most likely to exists: the kinetics have to be admissible and there is a value u * ∈ (u 0 T , u 0 H ). In this situation, the range of layer positions in every subinterval [x i−1 , x i ] is large and most layer positions correspond to a jumpū i close to u * . This makes it more likely that a stable irregular solution exists.
To examine this, we perform simulations of model (1)-(3) with discontinuous initial conditions. To do so, we choose a number k ∈ N, k ≥ 2 and positions 0 <x 1 < · · · <x k−1 < 1 for the discontinuities, such thatx i −x i−1 ≥ 0.05 holds. Moreover, we select points P 1 , . . . , P k ∈ R 2 with positive coordinates, such that P i with even index i are attracted by the steady state S 0 , whereas those with odd index i are attracted by S 2 (or the other way around). Now we define the initial For the admissible kinetic functions (62) in Appendix 1 we observed that the broader region of high concentration persists, whereas the narrower region of high concentration disappears (Fig. 8C for (x 1 ,x 2 ,x 3 ,x 4 ) = (0.2, 0.3, 0.4, 0.8)) the narrow peak disappears and the broad peak persists). However, for the kinetic functions (61) in Appendix 1 we observed the opposite behaviour.
To sum up, the simulations show that a solution of the model strongly depends on the initial data. When there exists a stable stationary solution with layer positions where it has been prescribed by the initial function, we observe that the timedependent solution is quickly approaching this stationary solution. When there is no such stable stationary solution, we observe a moving front finally leading to a constant solution or a stationary solution with fewer layer positions than prescribed by the initial function. We never observed a moving front leading to a stable stationary solution with layer positions at positions different from those prescribed by the initial function.