Instability of turing patterns in reaction-diffusion-ODE systems

The aim of this paper is to contribute to the understanding of the pattern formation phenomenon in reaction-diffusion equations coupled with ordinary differential equations. Such systems of equations arise, for example, from modeling of interactions between cellular processes such as cell growth, differentiation or transformation and diffusing signaling factors. We focus on stability analysis of solutions of a prototype model consisting of a single reaction-diffusion equation coupled to an ordinary differential equation. We show that such systems are very different from classical reaction-diffusion models. They exhibit diffusion-driven instability (turing instability) under a condition of autocatalysis of non-diffusing component. However, the same mechanism which destabilizes constant solutions of such models, destabilizes also all continuous spatially heterogeneous stationary solutions, and consequently, there exist no stable Turing patterns in such reaction-diffusion-ODE systems. We provide a rigorous result on the nonlinear instability, which involves the analysis of a continuous spectrum of a linear operator induced by the lack of diffusion in the destabilizing equation. These results are extended to discontinuous patterns for a class of nonlinearities.


Introduction
In this paper we focus on diffusion-driven instability (DDI) in systems of equations consisting of a single reaction-diffusion equation coupled with an ordinary differential equation system. Such systems are important for systems biology applications; they arise for example in modeling of interactions between processes in cells and diffusing growth factors, such as in refs. Hock et al. (2013), Klika et al. (2012), Marciniak-Czochra (2003), Kimmel (2006, 2008), Pham et al. (2011), Umulis et al. (2006). In some cases they can be obtained as homogenization limits of models describing coupling of cell-localized processes with cell-to-cell communication via diffusion in a cell assembly (Marciniak-Czochra and Ptashnyk 2008;Marciniak-Czochra 2012). Other examples are discussed e.g. in refs. Chuan et al. (2006), Evans (1975), Marciniak-Czochra et al. (2015), Mulone and Solonnikov (2009), Wang et al. (2013) and in the references therein. A detailed discussion of the DDI phenomena in the three-component systems with some diffusion coefficients equal to zero is found in the recent work (Anma et al. 2012).
Diffusion-driven instability, also called the Turing instability, is a mechanism of de novo pattern formation, which has been often used to explain self-organization observed in nature. DDI is a bifurcation that arises in a reaction-diffusion system, when there exists a spatially homogeneous stationary solution which is asymptotically stable with respect to spatially homogeneous perturbations but unstable to spatially heterogeneous perturbations. Models with DDI describe a process of a destabilization of stationary spatially homogeneous steady states and evolution of the system towards spatially heterogeneous steady states. DDI has inspired a vast number of mathematical models since the seminal paper of Turing (1952), providing explanations of symmetry breaking and de novo pattern formation, shapes of animal coat markings, and oscillating chemical reactions. We refer the reader to the monographs by Murray (2002Murray ( , 2003 and to the review article (Suzuki 2011) for references on DDI in the two component reaction-diffusion systems and to the paper Satnoianu et al. (2000) in the several component systems.
However, in many applications there are components which are localized in space, which leads to systems of ordinary differential equations coupled with reactiondiffusion equations. Our main goal is to clarify in what manner such models are different from the classical Turing-type models and to demonstrate that the spatial structure of the pattern emerging via DDI cannot be determined based on linear stability analysis.
To understand the role of non-diffusive components in the pattern formation process, we focus on systems involving a single reaction-diffusion equation coupled to ODEs. It is an interesting case, since a scalar reaction-diffusion equation cannot exhibit stable spatially heterogenous patterns (Casten and Holland 1978) and hence in such models it is the ODE component that yields the patterning process. As shown in ref. Marciniak-Czochra et al. (2013), it may happen that there exist no stable stationary patterns and the emerging spatially heterogeneous structures are of a dynamical nature. In numerical simulations of such models, solutions having the form of unbounded periodic or irregular spikes have been observed (Härting and Marciniak-Czochra 2014).
Thus, the aim of this work is to investigate to which extent the results obtained in Marciniak-Czochra et al. (2013), concerning the instability of all stationary structures, apply to a general class of reaction-diffusion-ODE models with a single diffusion operator.
We focus on the following two-equation system in a bounded domain ⊂ R N for N ≥ 1, with a C 2 -boundary ∂ , supplemented with the Neumann boundary condition where ∂ ν = ∂ ∂ν and ν denotes the unit outer normal vector to ∂ , and with initial data (1.4) The nonlinearities f = f (u, v) and g = g (u, v) are arbitrary C 3 -functions. Notice that Eq. (1.2) may contain an arbitrary diffusion coefficient which, however, can be rescaled and assumed to be equal to one.
In this paper we investigate stability properties of stationary solutions of the problem (1.1)-(1.3). Our main results are Theorems 2.1 and 2.11 which assert that, under a natural assumption satisfied by a wide variety of systems, stationary solutions are unstable. We call this assumption the autocatalysis condition (see Theorem 2.1) following its physical motivation in the model. We show in Section 3 that this condition is satisfied for all stationary solutions of a wide class of systems from mathematical biology. Our results are different in continuous and discontinuous stationary solutions. In the latter case, additional assumptions on the structure of nonlinearities are required.
As a complementary result to the instability theorems, we prove Theorem 2.9 which states that each non-constant regular stationary solution intersecting (in a sense to be defined) constant steady states with the DDI property, has to satisfy the autocatalysis condition. It is a classical idea by Turing that stable patterns appear around the constant steady state in systems of reaction-diffusion equations with DDI. Mathematical results on stability of such patterns can be found, e.g., in refs. Iron et al. (2004), Wei (2008), Wei and Winter (2007, 2008 and in the references therein. In the current work, combining Theorems 2.1 and 2.9, we show that this is not the case in the reaction-diffusion-ODE problems (1.1)-(1.4). In other words, the same mechanism which destabilizes constant solutions of such models, destabilizes also non-constant stationary solutions, a behavior that does not fit the usual paradigm of the reactiondiffusion-type equations. See Remark 2.10 for more details.
Mathematically, in the proof of our main result we need to consider a nonempty continuous spectrum of the linearized operator. This seems to be a novelty in the study of reaction-diffusion equations, and is caused by the absence of diffusion in one of the equations. In Section 4, we provide a rigorous proof of the nonlinear instability of steady states by using ideas from fluid dynamics equations.
The paper is organized as follows. In Section 2, we state the main results. Section 3 provides relevant mathematical biology-related examples of reaction-diffusion-ODE systems. Proofs are postponed to Sects. 4 and 5. Section 4 is devoted to showing instability of the stationary solutions under the autocatalysis condition. A proof of the instability of discontinuous solutions requires additional conditions on the model nonlinearities. In Section 5, the continuous stationary solutions are characterized and it is shown that the autocatalysis condition is satisfied in the class of reaction-diffusion-ODE problems (1.1)-(1.4) exhibiting DDI. Appendix contains additional information on the model of early carcinogenesis which was the main motivation for our research.

Results and comments
First we formulate a condition which leads to instability of regular stationary solutions of the problem (1.1)-(1.4). Then, we show that it is the necessary condition for DDI in reaction-diffusion-ODE systems. Finally, we extend the instability results to a class of discontinuous stationary solutions satisfying additional assumptions.

Instability of regular steady states
First, we focus on regular stationary solutions (U, V ) of problem (1.1)-(1.3). For this, we assume that there exists a solution (not necessarily unique) of the equation satisfies the elliptic problem We show that all regular stationary solutions of the problem (1.1)-(1.4) are unstable under a simple assumption on the first equation.
Inequality (2.7) can be interpreted as autocatalysis in the dynamics of u at the steady state (U, V ) at some point x 0 ∈ . Stability of the stationary solution is understood in the Lyapunov sense. Moreover, we prove nonlinear instability of the stationary solutions of the problem (1.1)-(1.3) and not only their linear instability, i.e. the instability of zero solution of the corresponding linearized problem, see Section 4 for more explanations.
Each constant solution (ū,v) ∈ R 2 of the problem (2.1)-(2.3) is a particular case of a regular solution. Thus, Theorem 2.1 provides a simple criterion for the diffusiondriven instability (DDI) of (ū,v).
then it has the DDI property.
This corollary follows directly from Theorem 2.1, because the second and the third inequality in (2.8) imply that (ū,v) is stable under homogeneous perturbations; see Remark 2.4 for more details.

Assumption 2.3 (Non-degeneracy of the stationary solutions) Let all stationary solu
(2.9) Remark 2.4 Let us note that the first two conditions in (2.9) allow us to study the asymptotic stability of (ū,v) treated as a solution of the corresponding system of ordinary differential equations by analyzing eigenvalues of the corresponding linearization matrix. Indeed, the conditions for linearized stability read 1. If then the Jacobi matrix has all eigenvalues with negative real parts, and hence (ū,v) is an asymptotically stable solution of system (2.10). 2. On the other hand, if (2.13) then the linearization matrix (2.12) has an eigenvalue with a positive real part, and consequently, the pair (ū,v) is an unstable solution of (2.10). Now, we state a simple but fundamental property of the stationary solutions of the problem (1.1)-(1.4). To prove Proposition 2.5, it suffices to integrate equation (2.2) over and to use the Neumann boundary condition (2.3) to obtain In the case described by Proposition 2.5, we say that a non-constant solution (U, V ) intersects a constant solution (ū,v). Now, we prove an important property of the constant solutions that are intersected by non-constant regular solutions.
Proposition 2.6 Let U (x), V (x) be a regular non-constant stationary solution of problem (1.1)-(1.3) and assume that all constant stationary solutions that are intersected by (U, V ) are non-degenerate, i.e. relations (2.9) are satisfied. Then, at least at one of those constant solutions, denoted here by (ū,v), the following inequality holds (2.14) The proof of Proposition 2.6 is based on the properties of the solutions of the elliptic Neumann problem (2.4)-(2.5) (see Theorem 5.1, below), which we prove in Sect. 5.
Remark 2.7 Every non-degenerate constant solution (ū,v) of the problem (1.1)-(1.4) satisfying inequality (2.14) is unstable. If both factors on the left-hand side of inequality (2.14) are positive, then, in particular, the autocatalysis condition f u (ū,v) > 0 is satisfied. Hence, the constant solution (ū,v) is an unstable solution of the reactiondiffusion-ODE system (1.1)-(1.4) by Theorem 2.1. On the other hand, if both factors on the left-hand side of inequality (2.14) are negative, then, in particular, the determinant in inequality (2.14) is negative and the constant vector (ū,v) is an unstable solution of the corresponding kinetic system (2.10), see the alternative (2.13) in Remark 2.4.
Remark 2.8 It is worth to emphasize the following particular case of the phenomenon described in Remark 2.7, because we shall encounter it in our examples, further on. Suppose that the problem (1.1)-(1.4) has a non-constant regular stationary solution (U, V ) intersecting only one constant and non-degenerate steady state (ū,v) which is asymptotically stable as a solution of the kinetic system (2.10). In such case, inequality (2.14) together with the second inequality in (2.11) directly imply the autocatalysis condition f u (ū,v) > 0. Thus, by Theorem 2.1, (ū,v) is an unstable solution of the reaction-diffusion-ODE problem (1.1)-(1.4), i.e. the constant steady state (ū,v) has the DDI property. Below, in Theorem 2.9, we show that the non-constant stationary solution (U, V ) also satisfies the autocatalysis condition (2.7), and hence, it is unstable. Now, we are in the position to show that the autocatalysis condition (2.7) has to be satisfied in reaction-diffusion-ODE systems (1.1)-(1.3) with non-constant regular stationary solutions which intersect constant steady states with the DDI property.
Theorem 2.9 Let (U, V ) be a non-constant regular stationary solution of problem (1.1)-(1.4). Denote by (ū,v) a non-degenerate constant solution which intersects (U, V ), and satisfies inequality (2.14). Assume that (ū,v) is an asymptotically stable solution of the kinetic system (2.10). Then, there exists x 0 ∈ such that The following remark emphasizes importance of the above results.
Remark 2.10 The instability results from Theorem 2.1 and Corollary 2.2 combined with Theorem 2.9 can be summarized in the following way. This is a classical idea that, in a system of reaction-diffusion equations with a constant solution having the DDI property, one expects stable patterns to appear around that constant steady state. Such stationary solutions are called the Turing patterns. For the initial-boundary value problem for a reaction-diffusion-ODE system with a single diffusion equation (1.1)-(1.3), such stationary solutions can be constructed in the case of several models of interest (see Section 3). However, the same mechanism that destabilizes constant solutions of such models, also destabilizes the non-constant solutions. In other words, all Turing patterns in the reaction-diffusion-ODE problems (1.1)-(1.3) are unstable.

Instability of non-regular steady states
The initial-boundary value problem (1.1)-(1.4) may also have non-regular steady states in the case when the equation f (U, V ) = 0 is not uniquely solvable. Choosing different branches of solutions of the equation f U (x), V (x) = 0, we obtain the relation U (x) = k V (x) with a discontinuous, piecewise C 1 -function k. Here, we recall that a pair U, holds for all test functions ϕ ∈ W 1,2 ( ).
In this work, we do not prove the existence of such discontinuous solutions and we refer the reader to classical works (Aronson et al. 1988;Mimura et al. 1980;Sakamoto 1990) as well as to our recent paper Marciniak-Czochra et al. (2013, Thm. 2.9) for information about how to construct such solutions to one dimensional problems using phase portrait analysis. Our goal is to formulate a counterpart of the autocatalysis condition (2.7), which leads to instability of the weak (including discontinuous) stationary solutions.

Theorem 2.11 Assume that (U, V ) is a weak bounded solution of the problem (2.1)-(2.3) satisfying the following counterpart of the autocatalysis condition
is an unstable solution of the initial-boundary value problem (1.1)-(1.4).
We prove Theorem 2.11 in Subsect. 4.6 by applying ideas developed for the analysis of the Euler equation and other fluid dynamics models. In that approach, it suffices to show that the spectrum of the linearized operator with the Neumann boundary condition ∂ ν v = 0, has so-called spectral gap, namely, there exists a subset of the spectrum σ (L), which has a positive real part, separated from zero. Here, we prove that σ x ∈ } and of isolated eigenvalues of L, see Section 4 and, in particular, Fig. 1 for more detail. One should emphasize that the instability of steady states from Theorems 2.1 and 2.11 is caused not by an eigenvalue with a positive real part, but rather by positive numbers from the set Range f u (U, V ) which is contained in the continuous spectrum of the operator (L, D(L)), see Theorem 4.5 below for more details.
In fact, in the case of particular nonlinearities, we do not need to assume that condition (2.15) holds true for almost all x ∈ . Indeed, if f (0, v) = 0, one may have stationary solutions U = U (x) such that U (x) = 0 on a subset of and U (x) > 0 on a complement. Such stationary solutions can be, for example, constructed for the carcinogenezis model (3.5)-(3.7) presented below (see Marciniak-Czochra et al. (2013)), and for several other one-dimensional equations discussed in ref. Mimura et al. (1980). In the following corollary, we show instability of the discontinuous stationary solutions, under the autocatalysis condition only for x ∈ such that U (x) = 0.
Corollary 2.12 (Instability of weak solutions) Assume that the nonlinear term in the Remark 2.13 A typical nonlinearity satisfying the assumptions of Corollary 2.12 has the form f(u, v) = r(u, v)u. It can be found in the models, where the unknown variable u evolves according to the Malthusian law with a density dependent growth rate r .
We defer the proofs of Theorems 2.1 and 2.11 as well as of Corollary 2.12 to Subsection 4.6. Theorem 2.9 is somewhat independent of Theorems 2.1 and 2.11 and it is proven in Section 5.

Model examples
In this section, our results are illustrated by applying them to some models from mathematical biology.

Gray-Scott model
First we consider a reaction-diffusion-ODE model with nonlinearities as in the celebrated Gray-Scott system describing pattern formation in chemical reactions (Gray and Scott 1983). The system with a non-diffusing activator takes the form with the zero-flux boundary condition for v and with nonnegative initial conditions. The constants B and k are assumed to be positive. The system exhibits the instability phenomenon described in Sect. 2.
Here, every regular positive stationary solution (U, V ) of the Neumann boundaryinitial value problem for equations (3.1)-(3.2) has to satisfy the relation U = (B + k)/V , hence, All continuous positive solutions of such boundary value problem in one dimensional case have been constructed in our recent paper Marciniak-Czochra et al. (2013, Sec. 5). A construction of discontinuous stationary solutions of the reaction-diffusion-ODE problem for (3.1)-(3.2) can be also found in Marciniak-Czochra et al. (2013, Thm. 2.9).
Instability results in Theorems 2.1 and 2.9 imply that all stationary solutions (constant, regular as well as discontinuous) of the reaction-diffusion-ODE problem (3.1)-(3.2) are unstable under heterogeneous perturbations. For the proof, it suffices to notice that the autocatalysis assumptions (2.7) and (2.15) are satisfied, since, for

Model of early carcinogenesis
The main motivation for the research reported in this work has been the study of the reaction-diffusion system of three ordinary/partial differential equations modeling the diffusion-regulated growth of a cell population of the following form supplemented with zero-flux boundary conditions for the function w and with nonnegative initial conditions, Marciniak-Czochra et al. (2013). Here, the letters a, d c , d b , d g , d, D, κ 0 denote positive constants. The theory developed in this paper applies to a reduced two-equation version of the model (3.5)-(3.7), obtained using a quasi-steady state approximation of the dynamics of v. Applying the quasi-steady state approximation in Eq. (3.6) (i.e., setting v t ≡ 0), we obtain the relation v = u 2 w/(d b + d), which after substituting into the remaining Eqs. (3.5) and (3.7) yields the initial-boundary value problem for the following reaction-diffusion-ODE system A rigorous derivation of the two equation model (3.8)-(3.9) from the model (3.5)-(3.7) as well as other properties of the solutions to (3.8)-(3.9) are presented in Appendix A of this work. Moreover, numerical simulations suggest that the two-equation model exhibits qualitatively the same dynamics as system (3.5)-(3.7). The autocatalysis assumptions (2.7) and (2.15) are satisfied by simple calculations, similar to those in the previous example (see Marciniak-Czochra et al. (2013) for more details). As a consequence, all nonnegative stationary solutions of the system (3.8)-(3.9) (regular and non-regular) are unstable due to Theorems 2.1 and 2.11. This corresponds to our results on the three-equation model (3.5)-(3.7) proved in ref. Marciniak-Czochra et al. (2013).
Stability analysis of the space homogeneous solutions of the two equation model (3.8)-(3.9) is reported in Appendix B. In particular, by Remark 2.7, constant steady states of (3.8)-(3.9) are either unstable solutions of the corresponding kinetic system or they have the DDI property.

Model of glioma invasion
Our results can be also applied to the "go-or-grow" model introduced in ref. Pham et al. (2011) to investigate the dynamics of a population of glioma cells switching between a migratory and a proliferating phenotype in dependence on the local cell density. The model consists of two reaction-diffusion equations where tumor cells are decomposed into two sub-populations: a migrating population with density v(x, t) and a proliferating population with density u(x, t) (Caution: we changed the notation from Pham et al. (2011), where ρ 1 = v and ρ 2 = u). In this model, the constant μ > 0 is the rate at which cells change their phenotype and the constant r ≥ 0 is the proliferation rate. The function = (ρ) has the following explicit form with constant α > 0 and ρ * > 0. It describes two complementary mechanisms for the phenotypic transmissions.
Go-or-rest model. Let us first look at a particular version of model (3.10)-(3.11) with no proliferation rate (namely r = 0) which is called in Pham et al. (2011) as the "go-or-rest model": One can check by a simple calculation that this system supplemented with the Neumann boundary condition for v(x, t) has a one parameter family of constant stationary solutions: (3.14) These vectors are degenerate (i.e. they do not satisfy Assumption 2.3) because the determinant in (2.9) vanishes in this case. However, by an elementary analysis of the phase portrait of the system of the ODEs, one can show that vectors (3.14) are stable solutions of system (3.15). The constant steady state (3.14) satisfies the autocatalysis condition (2.7) if see Pham et al. (2011, Ch. 3.1) for further discussion. Thus, by our Theorem 2.1, constant stationary solutions (3.14) are unstable solutions of the reaction-diffusion-ODE system (3.12)-(3.13). System (3.12)-(3.13) has no heterogeneous stationary solutions, because the counterpart of the boundary-value problem (2.1)-(2.3) reduces in this case to the problem which has constant solutions, only.
It is beyond the scope of this work to study positive heterogeneous stationary solutions of the go-or-grow model with r > 0. However, if there exist regular and strictly positive stationary solutions, then under the assumption ( (1) + (1)) + (1 − (1)) < 0, they must be unstable by Theorems 2.1 and 2.9, see also Remark 2.10. In conclusion, the structures shown in simulations of the models in ref. Pham et al. (2011) are not Turing patterns.

Instability of the stationary solutions 4.1 Existence of solutions
We begin our study of properties of solutions to the initial-boundary value problem (1.1)-(1.4) by recalling results on local-in-time existence and uniqueness of solutions for all bounded initial conditions.
We recall that a mild solution of problem (1.1)-(1.4) is a pair of measurable functions u, v : [0, T ] × → R satisfying the following system of integral equations where e t is the semigroup of linear operators generated by Laplacian with the Neumann boundary condition. Since our nonlinearities f = f (u, v) and g = g (u, v) are Remark 4.2 If u 0 and v 0 are more regular, i.e. if for some α ∈ (0, 1) we have u 0 ∈ C α ( ), v 0 ∈ C 2+α ( ) and ∂ ν v 0 = 0 on ∂ , then the mild solution of problem (1.1)-(1.4) is smooth and satisfies u ∈ C 1,α [0, T ] × and v ∈ C 1+α/2, 2+α [0, T ] × . We refer the reader to Rothe (1984, Thm. 1, p. 112) as well as to Garroni et al. (2009) for studies of general reaction-diffusion-ODE systems in the Hölder spaces.

Linearization of reaction-diffusion-ODE problems
with the Neumann boundary condition, ∂ ν v = 0, where the linear operator L and the nonlinearity N are defined by formulas (4.4) and (4.6), resp.  1)-(1.4). We consider the following linear system Proof Here, we use the Sobolev space Notice that L is a bounded perturbation of the operator which generates an analytic semigroup on L p ( ) × L p ( ) for each p ∈ (1, ∞). Thus, it is well-known (see e.g. Engel and Nagel (2000, Ch. III.1.3) and Yagi (2010, Theorems 2.15 and 2.19)) that the same property holds true for the operator (L, D(L)).
Next, we show certain elementary estimate of the nonlinearity in equation (4.3).
for all u, v ∈ L ∞ such that u L ∞ < ρ and v L ∞ < ρ, where ρ > 0 is an arbitrary constant. If, moreover, U, W ∈ W 1, p ( ) then Proof The proofs of both inequalities consist in using the Taylor formula applied to the C 3 -nonlinearities f = f (u, v) and g = g(u, v) in problem (1.1)-(1.2).

Continuous spectrum of the linear operator
Now, we are in a position to study the spectrum σ (L) of the linear operator L, given by the formula (4.4) when we linearize the reaction-diffusion-ODE problem (1.1)-(1.4) at a regular stationary solution. Proof We show that for each λ ∈ [λ 0 , 0 ] the operator , etc., cannot have a bounded inverse. Suppose, a contrario, that (L−λI ) −1 exists and is bounded. Then, for a constant K = (L−λI ) −1 , we have for all (ϕ, ψ) ∈ L p ( ) × W 2, p ( ) or, equivalently, using the usual norms in L p ( )× W 2, p ( ) and in L p ( ) × L p ( ): (4.10) A contradiction will be obtained by showing that inequality (4.10) cannot be true for all (ϕ, ψ) ∈ L p ( ) × W 2, p ( ).
To prove this claim, first we observe that, for each λ ∈ [λ 0 , 0 ], there exists Next, for arbitrary ψ ∈ C ∞ c ( ) satisfying supp ψ ⊂ B ε , we choose ϕ ∈ L p ( ) such that supp ϕ ⊂ B ε and in such a way that ψ + g u ϕ + (g v − λ)ψ = ζ , where the function ζ ∈ L p ( ) satisfies ζ L p ( ) ≤ ε ϕ L p ( ) . Let us explain that such a choice of ϕ, ζ ∈ L p ( ) is always possible. We cut g u at the level ε in the following way Thus, we obtain Thus, substituting functions ϕ, ψ, and ζ into inequality (4.10), we obtain the estimate (4.11) Hence, choosing ε > 0 in such a way that 2K ε ≤ 1 and compensating the term 2K ε ϕ L p ( ) on the right-hand side of inequality (4.11) by its counterpart on the left-hand side, we obtain the estimates which, obviously, cannot be true for all ψ ∈ C ∞ c ( ) such that supp ψ ⊂ B ε . We have completed the proof that each λ ∈ [λ 0 , 0 ] belongs to σ (L).

Eigenvalues
In the Hilbert case D(L) = L 2 ( ) × W 2,2 N ( ), the remainder of the spectrum of L, D(L) consists of a discrete set of eigenvalues {λ n } ∞ n=1 ⊂ C \ [λ 0 , 0 ]. Here, we sketch the proof of this result, however, it does not play any role in our instability results.
As the usual practice, we analyze the corresponding resolvent equations with arbitrary F, G ∈ L 2 ( ). Here, one should notice that for every λ ∈ C \ [λ 0 , 0 ], one can solve equation (4.12) with respect to ϕ. Thus, after substituting the resulting (4.13), we obtain the boundary value problem (4.17) For a fixed λ ∈ C \ [λ 0 , 0 ], by the Fredholm alternative, either the inhomogeneous problem (4.15)-(4.16) has a unique solution (so, λ is not an element of σ (L)) or else the homogeneous boundary value problem has a nontrivial solution ψ. Hence, it suffices to consider those λ ∈ C \ [λ 0 , 0 ], for which problem (4.18)-(4.19) has nontrivial solution. Now, we are in a position to prove that the set σ (L) \ [λ 0 , 0 ] consists of isolated eigenvalues of L, only. Here, it suffices to use the following general result on a family of compact operators, which we state for the reader's convenience. The proof can be found in the Reed and Simon book Reed and Simon (1980, Thm. VI.14).

Theorem 4.6 (Analytic Fredholm theorem) Assume that H is a Hilbert space and denote by L(H ) the Banach space of all bounded linear operators acting on H . For an open connected set D ⊂ C, let f : D → L(H ) be an analytic operator-valued function such that f (z) is compact for each z ∈ D. Then, either
(a) (I − f (z)) −1 exists for no z ∈ D, or (b) (I − f (z)) −1 exists for all z ∈ D \ S, where S is a discrete subset of D (i.e. a set which has no limit points in D).
Let us rewrite problem (4.18)-(4.19) in the form where the operator G = "( − I ) −1 " supplemented with the Neumann boundary conditions is defined in the usual way. Here, ∈ R is a fixed number different from each eigenvalue of Laplacian with the Neumann boundary condition. Recall that, for each λ ∈ C \ [λ 0 , 0 ], the operator R(λ) : L 2 ( ) → L 2 ( ) is compact as the superposition of the compact operator G and of the continuous multiplication operator with the function q(λ) + ∈ L ∞ ( ). Moreover, the mapping λ → R(λ) from the open set C \ [λ 0 , 0 ] into the Banach space of linear compact operators is analytic, which can be easily seen using the explicit form of q(λ) in (4.17). Thus, the set σ (L) \ [λ 0 , 0 ] consists of isolated points due to the analytic Fredholm Theorem 4.6. Here, to exclude the case (a) in Theorem 4.6, we have to show that the operator I − R(λ) is invertible for some λ ∈ C \ [λ 0 , 0 ]. This is, however, the consequence of the fact that the inhomogeneous boundary value problem (4.15)-(4.16) has a unique solution if λ > 0 is chosen so large that q(x, λ) < 0.

Linearization principle
The next goal in this section is to recall that, under appropriate conditions, the linear instability of the stationary solutions of a reaction-diffusion-ODE problem implies their nonlinear instability. Such a theorem is well-known for ordinary differential equations. Furthermore, in the case of reaction-diffusion equations where the spectrum of a linearized problem is discrete, one my apply the abstract result from the book by Henry (1981, Thm.5.1.3). However, in the case of reaction-diffusion-ODE problems, the linearized operator at a stationary solution (either smooth or discontinuous) may have a non-empty continuous spectrum (cf. Theorem 4.5). Hence, checking the assumptions of general results from Henry (1981) does not seem to be straightforward. Therefore, here, we propose a different approach.
Let us consider a general evolution equation where L is a generator of a C 0 -semigroup of linear operators {e tL } t≥0 on a Banach space X and N is a nonlinear operator such that N (0) = 0. First, we recall an idea introduced by Shatah and Strauss (2000) which asserts that, under relatively strong assumption on a nonlinearity in equation (4.20), the existence of a positive part of the spectrum of the linear operator L is sufficient to show that the zero solution of equation (4.3) is unstable. This is the precise statement of that result.

Theorem 4.7 (Shatah and Strauss (2000, Thm 1)) Consider an abstract problem (4.20), where
(1) the linear operator L generates a strongly continuous semigroup of linear operators on a Banach space X , (2) the intersection of the spectrum of L with the right half-plane {λ ∈ C : : Re λ > 0} is nonempty.

Then the zero solution of this equation is (nonlinearly) unstable.
We apply Theorem 4.7 to show an instability of regular steady states. In the case of discontinuous stationary solutions, we are unable to show that the nonlinearity in equation (4.3) satisfies the the condition (3) of Theorem 4.7. One may overcome this obstacle by assuming that the the spectrum σ (L) has so-called spectral gap. This classical method has been recently used by Mulone and Solonnikov Mulone and Solonnikov (2009) to show the instability of regular stationary solutions to certain reaction-diffusion-ODE problems, however, assumptions imposed in Mulone and Solonnikov (2009) are not satisfied in our case.
The crucial idea underlying this approach is to use two Banach spaces: a "large" space Z where the spectrum of a linearized operator is studied and a "small" space X ⊂ Z where an existence of solutions can be proved. More precisely, let (X, Z ) be a pair of Banach spaces such that X ⊂ Z with a dense and continuous embedding. A solution w ≡ 0 of the Cauchy problem (4.20) is called (X, Z )-nonlinearly stable if for every ε > 0, there exists δ > 0 so that if w(0) ∈ X and w(0) Z < δ, then (1) there exists a global in time solution of (4.20) such that w ∈ C([0, ∞); X ); (2) w(t) Z < ε for all t ∈ [0, ∞).
An equilibrium w ≡ 0 that is not stable (in the above sense) is called Lyapunov unstable.
In this work, we drop the reference to the pair (X, Z ). Let us also note that, under this definition of stability, a loss of the existence of a solution of (4.20) is a particular case of instability. Now, we recall a result linking the existence of the so-called spectral gap to the nonlinear instability of a trivial solution to problem (4.20).
(2) The nonlinear term N satisfies the inequality for some constants C 0 > 0 and ρ > 0.
Then, the trivial solution w 0 ≡ 0 of the Cauchy problem (4.20) is nonlinearly unstable.
The proof of this theorem can be found in the work by Friedlander et al. Friedlander et al. (1997, Thm. 2.1).

Remark 4.9
The operator L considered in this work satisfies the "spectral mapping theorem": σ (e tL ) \{0} = e tσ (L) , see Lemma 4.3. Thus, due to the relation |e z | = e Re z for every z ∈ C, the spectral gap condition required in Theorem 4.8 holds true if for every λ ∈ σ (L), either Re λ ∈ (κ, μ) or Re λ ∈ (M, ).
Remark 4.10 The authors of the reference Friedlander et al. (1997, Thm. 2.1) formulated their instability result under the spectral gap condition for a group of linear operators {e tL } t∈R and in the case of a finite constant κ (caution: in Friedlander et al. (1997), the letter λ is used instead of κ). However, the proof of Friedlander et al. (1997, Thm. 2.1) holds true (with a minor and obvious modification) in the case of a semigroup {e tL } t≥0 as well as κ = −∞ is allowed, as stated in Theorem 4.8. This extension is important to deal with the operator L introduced in Lemma 4.3, which generates a semigroup of linear operators, only, and which may have an unbounded sequence of eigenvalues.

Proofs of instability results
Proof of Theorem 2.1 Here, it suffices to apply Theorem 4.7 to the semi-linear equation (4.3) with the Banach space Recall the well-known embedding W 1, p ( ) ⊂ L ∞ ( ) for every p ∈ (n, ∞].. We refer the reader to Yagi (2010, Ch. 2) for the proof that the operator L discussed in Lemma 4.3 generates a semigroup of linear operators on X . The autocatalysis condition (2.7) combined with Theorem 4.5 imply that σ (L) meets the right-hand plane of C. Due to the embedding X ⊂ L ∞ ( ) × L ∞ ( ), inequalities (4.7) and (4.8) imply that the nonlinear mapping N in (4.6) satisfies the condition (3) of Theorem 4.7 with η = 1.
Hence, the regular stationary solution (U, V ) is unstable.

Proof of Theorem 2.11
To show an instability of non-regular stationary solution, we begin as in the proof of Theorem 2.1. First, we linearize our problem at a weak bounded stationary solution (U, V ) and we notice that assumptions of Lemmas 4.3 and 4.4 are satisfied. Next, following the arguments from the proof of Theorem 4.5 we show that the number f u (U (x 0 ), V (x 0 )) belongs to σ (L), where f u (U (x), V (x)) is positive at x 0 and continuous in its neighborhood. Notice that we do not need to show that all numbers from Range f u (U, V ) are in σ (L) to show the spectral gap condition required by Theorem 4.8. The reasoning from Subsection 4.4 concerning eigenvalues can be copied here without any change because q(λ, x) defined in (4.17) is a bounded function for every λ ∈ C \ [λ 0 , 0 ]. Now, let us show that the operator L has a spectral gap as required in assumption (1) of Theorem 4.8.
Proof of Corollary 2.12 Here, the analysis is similar to the case of regular stationary solutions discussed in Theorem 2.1, hence, we only emphasize the most important steps.
Let (U, V ) be a weak solution of problem (2.1)-(2.3) and denote by I ⊂ its null set, namely, a measurable set such that U (x) = 0 for all x ∈ I and U (x) = 0 for all x ∈ \ I. For a null set I, we define the associate L 2 -space supplemented with the usual L 2 -scalar product, which is a Hilbert space as the closed subspace of L 2 ( ). In the same way, we define the subspace L ∞ Obviously, when the measure of I equals zero, we have L 2 I ( ) = L 2 ( ). The imposed assumptions imply that I is different from the whole interval. Now, observe that if u 0 (x) = 0 for some x ∈ then by equations (1.1) with the nonlinearity f (u, v) = r (u, v)u we have u(x, t) = 0 for all t ≥ 0. Hence, the spaces The spectrum σ (L) is marked by thick dots and by the interval [λ 0 , 0 ] in the sector δ,ω 0 . The spectral gap is represented by the strip {λ ∈ C : μ ≤ Re λ ≤ M} without elements of σ (L) are invariant for the flow generated by problem (1.1)-(1.4) (notice that there is no "I" in the second coordinates of X I and Z I ). The crucial part of our analysis is based on the fact that, as long as we work in the space X I and Z I , we can linearize problem (1.1)-(1.4) at the weak solution (U, V ). Moreover, for each x ∈ \ I, the corresponding linearized operator agrees with L defined in Lemma 4.3. Hence, the analysis from the proof of Theorem 2.1 can be directly adapted to discontinuous steady states in the following way.
We fix a weak stationary solution (U I , V I ) with a null set I ⊂ . The Fréchet derivative of the nonlinear mapping F : Z I → Z I defined by the mappings (U, V ) → ( f (U, V ), g(U, V )) at the point (U I , V I ) ∈ Z I has the form This results immediately from the definition of the Fréchet derivative. Next, we study spectral properties of the linear operator in the Hilbert space Z I (see (4.22)) with the domain D(L I ) = L 2 I ( ) × W 2,2 ( ). Here, the reasoning from the proof of Theorem 2.1 can be directly adapted with modifications as in the proof of Theorem 2.11.
Finally, we may study the discrete spectrum of L I in the same way as in Subsection 4.4 because the corresponding function q(λ, x) is bounded for λ ∈ C \ [λ 0 , 0 ]. The proof of instability of the stationary solution (U I , V I ) is completed by Theorem 4.8 and Lemmas 4.3-4.4.

Constant steady states which are intersected by non-constant stationary solutions
First, we prove a certain property of stationary solutions to a general elliptic Neumann problem. This result will imply immediately Proposition 2.6.
Then, there exists x 0 ∈ and a 0 ∈ R such that Proof First, as in the proof of Proposition 2.5, we integrate the equation in (5.1) and we use the Neumann boundary condition to obtain h(V (x)) dx = 0. Hence, there exists x 0 ∈ and a 0 ∈ R such that V (x 0 ) = a 0 and h(a 0 ) = 0. Now, we suppose that h (a 0 ) < 0, and consider two cases: x 0 ∈ and x 0 ∈ ∂ , separately.
Let x 0 ∈ . Since h(a 0 ) = 0, we have Using the well-known formula Hence, there exists an open neighbourhood U ⊂ of x 0 such that r (x, a 0 ) < 0 for all x ∈ U. Suppose that r (x, a 0 ) < 0 for all x ∈ . Multiplying both sides of equation (5.3) by V (x) − a 0 and integrating over , we obtain This implies that V (x) ≡ a 0 , which is a contradiction, because we assume that V = V (x) is a non-constant solution. Therefore, there exists x 1 ∈ ∂U ∩ such that r (x 1 , a 0 ) = 0. It follows from equation (5.3) that V (x 1 ) = 0, and consequently, from equation (5.1) we have h(V (x 1 )) = 0. Hence, there exists a 1 ∈ R such that V (x 1 ) = a 1 and h(a 1 ) = 0. Note that a 1 = a 0 . Thus, if the equation h(V ) = 0 has only one root, the proof of Theorem 5.1 is completed. Now, we consider two cases. Case I: The equation h(V ) = 0 has no solution between a 0 and a 1 . Thus, we define the function and, without loss of generality, we can assume that a 0 < ψ(s) < a 1 for all 0 < s < 1.
Since h(a 0 ) = 0 and h (a 0 ) < 0, we have h(a 0 + θ 0 ) < 0 for small θ 0 > 0. If we also suppose that h (a 1 ) < 0, then, we can find small θ 1 > 0 such that h(a 1 − θ 1 ) > 0. Noting that ψ(s) is continuous function, we see that there exist s * , s * * ∈ (0, 1) such that ψ(s * ) = a 0 + θ 0 and ψ(s * * ) = a 1 − θ 1 . This implies that there existsŝ ∈ (s * , s * * ) for which h(V (x 0 +ŝ(x 1 − x 0 ))) = 0, and from the assumption for ψ(s), This is a contradiction with the hypothesis that the equation h(V ) = 0 has no roots between a 0 and a 1 . Hence, h (a 1 ) ≥ 0. Case II: The equation h(V ) = 0 has a solution a m between a 0 and a 1 . It is clear that V (x) has to intersect a m , too. Choosing a m the closest root of h(V ) = 0 to a 0 , we repeat the argument from Case I to show that h (a m ) ≥ 0.
Next, let x 0 ∈ ∂ . Following the previous reasoning and using the hypothesis h (a 0 ) < 0, we find a ball B ⊆ such that x 0 ∈ ∂ and r (x, a 0 ) < 0 for all x ∈ B. Moreover, we can assume that either V (x) > a 0 or V (x) < a 0 for all x ∈ B, because, if there exists x 1 ∈ B such that V (x 1 ) = a 0 , then we can apply the same argument as in the first part of this proof to obtain h (a 0 ) ≥ 0.
Let V (x) < a 0 for all x ∈ B, and we apply the Hopf boundary lemma to equation (5.3) in the ball B. If V is a non-constant solution satisfying V (x) − a 0 < 0 and V (x 0 ) − a 0 = 0, then necessarily ∂ V (x 0 )/∂ν > 0, which contradicts the Neumann boundary condition satisfied by V at x 0 ∈ ∂ . Now, we consider the case V (x) > a 0 for all x ∈ B. Here, the function Moreover, we differentiate both sides of the equation Finally, choosing v =v and u =ū = k(v) and substituting equation (5.8) into (5.7), we obtain This completes the proof of inequality (2.14).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A. Model of early carcinogenesis -quasi-stationary approximation
In Appedix, we discuss in detail the model example from Subsection 3.2. Initial-boundary value problems for reaction-diffusion-ODE systems arise in the modeling of the growth of a spatially-distributed cell population, where proliferation is controlled by endogenous or exogenous growth factors diffusing in the extracellular medium and binding to cell surface as proposed by Marciniak and Kimmel in the series of recent papers Marciniak-Czochra and Kimmel (2006, 2008. Here, we consider the following particular case of such systems, which was studied by us in Marciniak-Czochra et al. (2013) on (0, ∞) × , supplemented with zero-flux boundary conditions for w ε ∂ ν w ε (x, t) = 0 for all t > 0, x ∈ ∂ (6.4) and with nonnegative, smooth and bounded initial data In equations (6.1)-(6.3), the letters a, d c , d b , d g , d, κ 0 , D denote positive constants.
As it was shown in Marciniak-Czochra et al. (2013, Sec. 3), solutions of this problem are nonnegative and stay bounded on every interval [0, T ]. Here, we study the behavior of these solutions as ε → 0.
First, we notice that choosing ε = 0 in equation (6.2), we obtain the identity Substituting it to the two remaining equations (6.1), (6.3) yields the system u(x, 0) = u 0 (x) and w(x, 0) = w 0 (x). (6.10) Repeating our reasoning from Marciniak-Czochra et al. (2013) one can show that this new system has also a unique, global-in-time nonnegative, smooth solution for nonnegative and smooth initial data.
In this part of Appendix, we show that solutions of the quasi-stationary system (6.6)-(6.10) are good approximation of solutions of the original three-equation model (6.1)-(6.5). Our main result concerns uniform estimates for an approximation error for u and w on each finite time interval [0, T ].
First, we show that solutions of ε-problem (6.1)-(6.5) are uniformly bounded with respect to small ε > 0, locally in time.
It will be clear from the proof of Lemma 6.1 that the constant C(T ) growths exponentially in T > 0.
Proof of Lemma 6.1 Since v ε /(u ε + v ε ) ≤ 1 for nonnegative solutions, equation (6.1) yields the inequality, (6.11) Applying the comparison principle to the parabolic equation (6.3) with the Neumann boundary condition we obtain the estimate where the function C w = C w (t) satisfies the Cauchy problem (6.13) and is given by the formula Next, we use equation (6.2) to obtain (6.15) Thus, using inequality (6.11) and plugging the above estimate into (6.14) yields Changing the order of integration, we can simplify the last term on the right-hand side hence, for 0 < ε < d b +d d g , using inequalities (6.12) we obtain where C = C( w 0 L ∞ ( ) , v 0 L ∞ ( ) ) is independent of T and of ε. Finally, the Gronwall inequality applied to (6.17) implies the estimate (6.18) and for all sufficiently small ε > 0 (e.g. ε ∈ 0, (d b + d)/(2d g ) ), where positive constants C and C 1 (T ) are independent of ε. Finally, estimate (6.18) applied to inequality (6.15) implies an analogous bound for v ε .
Proof Letting α = u ε − u, β = v ε − v and δ = w ε − w, we obtain by the Taylor expansion the following system (6.24) supplemented with the initial conditions with v 0 obtained from u 0 and w 0 via formula (6.6), and with the Neumann boundary condition for δ(x, t). In equations (6.22)-(6.24) the following coefficients are calculated in certain intermediate points and are bounded independently of ε due to Lemma 6.1. The proof is divided into three steps.

Remark 6.3
To obtain a better estimate of β, one should construct an initial value layer, since β| t=0 = 0.

Appendix B. Model of early carcinogenesis -constant steady states
Here, we consider the space homogeneous solutions of the two-equation model (6.6)-(6.10) which satisfy the corresponding kinetic system where a, d c , d b , d, d g , κ 0 are positive constants and we have always assumed that a > d c . The structure of constant steady states of this system is the same as of the original three-equation model and can be characterized by the following lemma (for the proof see Marciniak-Czochra et al. (2013)).