Stability of nonconstant stationary solutions in a reaction-diffusion equation coupled to the system of ordinary differential equations

In this paper we study pattern formation arising in a system of a single reaction-diffusion equation coupled with subsystem of ordinary differential equations, describing spatially-distributed growth of clonal populations of precancerous cells, whose proliferation is controled by growth factors diffusing in the extracellular medium and binding to the cell surface. We extend the results on the existence of nonhomogenous stationary solutions obtained in the previous work of one of the authors to a general Hill-type production function and full parameter set. Using spectral analysis and perturbation theory we derive conditions for the linearized stability of such spatial patterns.


Introduction
Partial differential equations of diffusion type have long served to model regulatory feedbacks and pattern formation in aggregates of living cells. Classical mathematical models of pattern formation in cell assemblies have been constructed using reaction-diffusion equations. They have been applied to describe pattern formation of animal coat markings, bacterial and cellular growth patterns, tumor growth and tissue development, see e.g., [10] and [11] and references therein. One of the mechanisms of pattern formation in reaction-diffusion systems, prevalent in the modeling literature since the seminal paper of Allan Turing [14], is diffusion driven instability (Turing-type instability).
Diffusion-driven instability arises in a reaction-diffusion system, when there exists a spatially homogeneous solution, which is asymptotically stable in the sense of linearized stability in the space of constant functions, but it is unstable with respect to spatially inhomogeneous perturbation. The majority of theoretical studies in theory of pattern formation focus on the analysis of the systems of two or more reaction-diffusion equations. In many biological applications it is relevant to consider systems consisting of a single reaction-diffusion equation coupled with a system of ordinary differential equations. Such models can also exhibit diffusion-driven instability. However, they are very different from classical Turing-type models and the spatial structure of the pattern emerging from the destabilisation of the spatially homogeneous steady state cannot be concluded from a linear stability analysis [8,10]. The models exhibit qualitatively new patterns of behavior of solutions, including, in some cases, a strong dependence of the emerging pattern on initial conditions and quasi-stability followed by rapid growth of solutions [8].
Mathematical theory exists only for special cases of such solutions arising in the models, which can be simplified to the single reaction-diffusion equation with nonlocal terms and with fast growing kinetics. The example of such systems are Gierer-Meinhard and Gray-Scott models. The existence and stability of the solutions of such systems were intensively studied using singular perturbation analysis and spectral analysis of the eigenvalue problem associated with the linearization around the "spike-like" solution, e.g. [3,4,15]. The structure of the kinetics involved in the cell proliferation model considered by us is different, in particular the autocatalysis is excluded by the biological assumptions and the solutions of our model are uniformly bounded. As shown in [7,9] the models have stationary solutions of periodic type, the maxima and minima of which may be of the spike or plateau type. Numerical simulations show that in some cases solutions of the model converge to a spatially heterogeneous pattern, which persists for an arbitrary long time while in other cases the transient growth of nonconstant pattern is observed and ultimately the solution converges to a stable spatially homogeneous state.
In this paper we approach the issue of the stability of nonconstant stationary solutions and investigate linear stability of spatially heterogeneous solutions of the model of the reaction-diffusion equation coupled with two ordinary differential equations, which was proposed in [9]. The paper is organized as follow: In Section 2 the model is introduced and the results on existence, regularity and boundedness of the model solutions are presented. In Section 3 existence of a spatially nonconstant steady state is shown. Section 4 is devoted to the linearized stability analysis.

Model description
We consider a model of a cell population controlled by a diffusive growth factor, with initial conditions where c denotes the concentration of precancerous cells, whose proliferation rate is reduced by cell crowding but enhanced in a paracrine manner by a hypothetical biomolecular growth factor b bound to cells. Free growth factor g is secreted by the cells, then it diffuses among cells with diffusion constant 1/γ, and binds to cell membrane receptors becoming the bound factor b. Then, it dissociates at a rate ν, returning to the b-pool. Parameter µ denotes a small influx of new precancerous cells due to mutation. All coefficients in the system (1) are assumed to be constant and positive, and k ∈ N. For k larger than 1, production of growth factor molecules is given by a Hill function and models a process with fast the saturation effect. Proof. Using the existence and regularity theory for systems of parabolic and ordinary differential equations (see [5], [12]), for (c in , b in , g in ) ∈ C 2+α ([0, 1]) 3 , we obtain, due to local Lipschitz continuity of reaction terms in the system, the existence of a local solution of (1), (c, b, g) ∈ C 1+α/2,2+α ([0, The theory of bounded invariant rectangles (see [2], [13]) and the properties of 1], t ≥ 0} is positive invariant and solutions are nonnegative for nonnegative initial data. Then, from the first equation in (1) follows the estimate and, by Gronwall inequality, we obtain, that c is bounded for all finite time, i.e., sup The boundedness of c implies also the boundedness of b and g, and existence of a global solution.
3 Existence of positive non-constant steady states Stationary solutions of (1) can be found from the system We look for a positive solution (c(x), b(x), g(x)) defined in (0, 1) that has at least one non-constant component. It is clear that for any solution of (2) either all functions c, b and g are constant or all depend on x effectively. The next observation is that the function c doesn't change the sign. In fact, if c(x 0 ) = 0 for some x 0 ∈ (0, 1), then it follows from the first equation that µ = 0, which is impossible. Let us solve the system with respect to c and b.
Lemma 3.1. If θ < d c , then there exists a unique solution (c(g), b(g)) of (3), which is continuous and positive for all g ≥ 0. In the case θ = d c there is only one positive solution and for θ > d c there are two positive solutions that exist only in some interval (0, g * ), where g * is a finite number.
Proof. Let us introduce the temporary notation (g) = αµ d b +ν g. First observe that and hence b is positive as soon as g is positive. The main question is whether there in fact exists a positive c. Substituting expression (4) into the first equation of (3) yields This quadratic equation with respect to c admits real roots provided Note that (p + q) 2 − 4εpq ≥ 0 for all positive p, q if and only if ε ≤ 1. Hence, under the assumption θ < d c , equation (5) has the solution which is positive for all g > 0. Furthermore, In the the critical case θ = d c there exists a unique solution of (5) where c is positive for g ∈ 0, dc(d b +ν) αµ . In the case θ > d c there are two positive solutions of (5), αµ . But only one of them is bounded at g = 0, and this solution is given by (6). Moreover the solution is also bounded at the point g * : Note that g * is the smallest positive point in which the discriminant D = D(g) vanishes. Observe also that (g) < d c for all g ≤ g * .
We next come back to system (2). Let ω = αd b d b +ν . In view of Lemma 3.1 the function g must be a positive solution of the nonlinear boundary value problem where c is given by (6), if θ = d c , and by (8) otherwise. We must keep in mind that the function c depends on parameters α, µ, ν, θ, d b and d c .

Nonlinear boundary value problem
Let us consider the boundary value problem where and c is given by (6) or (8). The equation of this type, so called system with one degree of freedom, was intensively studied in classical mechanics. For a deeper analysis of the system we refer to [1,Sec.12].
Theorem 3.2. If function h = h(g) has no less than three positive roots, then there exists a set Γ ⊂ R + of diffusion constants γ for which boundary value problem (11) admits a positive solution.
Proof. We assume for a while that γ = 1. Suppose that the potential energy has a positive critical point g 0 that is a local minimum of U . The trajectories of dynamic system g = z, z = h(g) resemble ellipses near the point (g 0 , 0) of the phase space. Let us consider the trajectory ζ that intersects the axis z = 0 at positive points q 1 and q 2 (see Fig. 1). Therefore there exists a periodic solution g = g(x), z = z(x) of the dynamic system subject to initial conditions g(0) = q 1 , z(0) = 0, and its period is given by ) . Remark that the energy U without local minimum points is either monotonic or a function with a single maximum point. In both cases no trajectory starting at a point (g 1 , 0) can return again to the line z = 0, with the exception of the equilibrium positions. Next, g is a positive solution of the initial problem g = h(g), g(0) = q 1 , g (0) = 0. Moreover, g ( nT 2 ) = 0 for any natural n, because g ( nT 2 ) = z(0) = 0 for even n and g ( nT 2 ) = z( T 2 ) = 0 for odd n. Given n, we consider the function u(x) = g( x √ γ ) that is obviously the solution to equation u = γh(u), and u (0) = 0. We set γ n = 4 n 2 T 2 , so that Hence, u is a positive solution to (11). It follows that any close trajectory ζ lying in the half-plane g > 0 produces a countable set of positive solutions u n (ζ, ·) to (11) with γ = γ n (ζ). All these sequences {γ n (ζ)} n∈N form the set Γ, which is in general uncountable.
The first and primary question is whether the energy U has a local minimum. We consider first the case θ < d c . The energy has a local minimum at point g = g 0 if only h(g 0 ) = 0, h > 0 in a left-side neighborhood of g 0 and h < 0 in a right-side one. Let .
We look for a root g 0 of the equation h 1 (g) = h 2 (g) for which h 1 (g) > h 2 (g) the left of g 0 and h 1 (g) < h 2 (g) on its right. A trivial verification shows that c = c(g) (given by (6)) is a steadily increasing function and Set c min = µ dc and c max = µ dc−θ . It is convenient to consider the c-representation of functions h j on interval (c min , c max ): α(dc−θ) . Both functions are strictly increasing in the region under study. In regard to h 1 , its derivative can be written as and is obviously positive for c ∈ (c min , c max ). The function h 1 has two vertical asymptotes c = 0 and c = c max , while h 2 is bounded. In addition, h 2 has an inflection point c = k k−1 k+1 , whereas the inflection point of h 1 depends on parameters of the model. Fig. 2 shows the typical plots of h 1 and h 2 . Evidently, the energy U has a local minimum if plots of the functions intersect at three points.
Similarly, for θ = d c we obtain that c(g) is monoton increasing and The function h 1 is nearly linear for small g and growth quadratically as g → dc(d b +ν) αµ , whereas Thus, for k ≥ 2 there exist sets of parameters d c , d b , ν, α, µ, β, such that the function h(g) has three roots in (0, dc(d b +ν) αµ ). For θ > d c and c(g) defined by (6), we obtain and again h 1 is nearly linear, h 2 (g) ∼ cg k for small g, h 1 (0) = 0 < h 2 (0), and both functions are bounded on (0, g * ). For k ≥ 2 there exist two roots, as functions of parameters, of h(g) on the interval (0, g * ). The third root of h(g) lies on the interval (g * * , ∞), where g * * is the largest positive point in which D = D(g) vanishes.
Remark 1. Let d c < θ, the solution c(g) of (5) is defined by parameters d c , d b , µ, ν, α satisfy the assumption 0 < g * < (d b + ν)d c /(αµ), and β is choosen such that h(g * ) < 0. Then the function h is continuous in (0, g * ], and h(g) → +∞ as g → +0. This properties of h and the assumption h(g * ) < 0 provide existence of a local minimum of the energy U at a point g 0 ∈ (0, g * ). Thus, there exists a set Γ ⊂ R + of diffusion constants γ for which the boundary value problem (11), with c(g) defined by (13), admits a positive solution. This situation was considered in [9] with k = 1.
If U has a local minimum g 0 , then there are a local maximum g − in the interval (0, g 0 ) and a local maximum g + in the interval (g 0 , ∞) as shown in Fig. 1. Set ε 0 = min{|g 0 − g − |, |g 0 − g + |}.

Stability and destabilisation of steady states
Let X be a complex Banach space with norm · and let A be a closed operator on X with a dense domain D(A). Suppose that −A is a sectorial operator and Re λ < 0 for all λ ∈ σ(A). Hereafter, σ(T ) stands for the spectrum of an operator T . For each s ≥ 0 we introduce the interpolation space X s = D((−A) s ) equipped with the norm x s = (−A) s x . Clearly X 0 = X. We consider the equation ∂ t u = Au + f (u) (14) in the Banach space X. Let u * ∈ D(A) be a stationary solution of (14) .
To study stability of the stationary solutions we apply the following proposition about stability and instability by the linear approximation that is a version of Theorems 5.1.1 and 5.1.3 in [5] adapted for our purposes.
Proposition 1. Let f : U → X be a locally Lipschitz continuous map in a neighborhood U ⊂ X s of a steady state u * , for some s ∈ (0, 1). Suppose that for z ∈ X s the map f admits the representation provided z s is small enough. If • B is a linear bounded operator from X s to X, • ||p(u * , z)|| = o(||z|| s ) as ||z|| s → 0, • the spectrum of A + B lies in the set {λ ∈ C : Re λ < h} for some h < 0, then the equilibrium solution u * of (14) is asymptotically stable in X s . Moreover, if the spectrum σ(A + B) and the half-plane {λ ∈ C : Re λ > 0} have a non-empty intersection, then u * is unstable.

Linearised analysis
For simplicity of notation, we denote the vector (c, b, g) by u = (u 1 , u 2 , u 3 ) so that system (1) reads and N 0 is the Sturm-Liouville operator given by on the interval [0, 1], subject to the Neumann boundary conditions, We have that A is a closed densely defined operator on X and the spectrum of A is given by Also for λ / ∈ σ(A) we have the estimate .
The function f is smooth in R 3 + = {y ∈ R 3 : y k > 0, k = 1, 2, 3}. Therefore, f admits the representation f (y + z) = f (y) + B(y)z + p(y, z), where the remainder satisfies the estimate in a neighborhood of any point y ∈ R 3 + with a continuous function ϑ. Here Suppose now that u * is a positive steady state, that is Au * + f (u * ) = 0 and u * (x) ∈ R 3 + for all x ∈ (0, 1). Assume also that u * j are C 2 -functions, j = 1, 2, 3. Obviously, B(u * ) is a linear bounded operator from X s to X, for each s ∈ (0, 1).
Next, using (17) we obtain f (u * + z) = f (u * ) + B(u * )z + p(u * , z) with the estimate of the remainder Consequently, for any z ∈ X s using the fact that v C[0,1] ≤ C v H 1 (0,1) we conclude as z s → 0, with constants c j , where j = 1, 2, being independent on t, and s ∈ [1/2, 1). For nonlinear model (1) the linearized system at the steady state u * can be written as ∂ t z = (A + B(u * ))z. Then turning back to the previous notation, we obtain The first and primary question is whether the spectrum of A + B(u * ) lies in the half-plane {λ ∈ C : Re λ < h} for some h < 0.
Moreover, the first equation in (2) yields for all x ∈ [0, 1]. Hence h * is a negative number.
Next, we will consider the operator A as the perturbation of A 0 by the operator A 1 . Let T and S be operators with the same domain space H such that D(T ) ⊂ D(S) and ||Su|| ≤ a||u|| + b||T u||, u ∈ D(T ), where a, b are nonnegative constants. Then, we say that S is relatively bounded with respect to T or simply T -bounded. Assume that T is closed and there exists a bounded operator T −1 , and S is T -bounded with constants a, b satisfying the inequality Then, T + S is a closed and bounded invertible operator [6, Th.
Remark 2. It is enough to verify (25) in a local minimum g 0 of the potential U from Theorem 3.2 and for the sets of parameters such that the local minimum exists. Then, Corollary 1 and the strict inequality in (25) provide stability of a stationary solution with the small amplitude.

Conclusions
In this paper we performed linearized stability analysis of the nonhomogeneous stationary solutions of a single reaction-diffusion equation coupled to ordinary differential equations. We focused on a specific model of pattern formation in a system of cells with the proliferation regulated by a diffusive growth factor. We extended results on the stationary problem obtained in [9] to the larger space of parameters and showed existence of inifinitely many stationary solutions. Numerical calculations indicate that the relationship between d c and θ is significant for the existence of nonconstant stationary solutions. It is observed that for d c = θ the system (1) admits the nonconstant stationary solutions when the difference |d c − θ| is small enough. For d c = θ the value of the parameter β plays an important role in existence of the stationary solutions. Performing spectral analysis using sectorial operators and the perturbation theory we obtained conditions for the linearized stability of the spatial patterns. Formulation of the conditions in terms of the model parameters independent of the values of stationary solutions is difficult to obtain due to the complex structure of the stationary problem. For some parameter values, due to existence of multiple solutions of the ODEs subsystem, the model exhibits also discontinuous stationary solutions and it is difficult to decide which patterns we see in numerical simulations. The mechanism of pattern selection is a topic of further analytical and numerical studies.