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

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 [9] 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.

Yuriy Golovaty

Department of Mechanics and Mathematics

Ivan Franko National University of L’viv

Universytetska str. 1, L’viv 79000, Ukraine

Anna Marciniak-Czochra

University of Heidelberg,

Interdisciplinary Center for Scientific Computing (IWR)

Institute of Applied Mathematics and BIOQUANT

Im Neuenheimer Feld 267, 69120 Heidelberg, Germany

Mariya Ptashnyk

Department of Mathematics

University of Dundee

DD1 4HN Dundee, Scotland, UK

2000 MSC Primary: 35K57, 35J57; Secondary: 92B99.

Keywords: Pattern formation, reaction-diffusion equations, linearized stability, spectral analysis

1 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.

2 Model description

We consider a model of a cell population controlled by a diffusive growth factor,


with initial conditions

where 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 bound to cells. Free growth factor is secreted by the cells, then it diffuses among cells with diffusion constant , and binds to cell membrane receptors becoming the bound factor . Then, it dissociates at a rate , returning to the -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 . For larger than 1, production of growth factor molecules is given by a Hill function and models a process with fast the saturation effect.

Theorem 2.1.

For nonnegative initial data there exists a nonnegative global solution of system (1), .


Using the existence and regularity theory for systems of parabolic and ordinary differential equations (see [5], [12]), for , we obtain, due to local Lipschitz continuity of reaction terms in the system, the existence of a local solution of (1), for some .

The theory of bounded invariant rectangles (see [2], [13]) and the properties of functions , , , i.e. are smooth for and in , for all , , for all , , and for all , , imply that the set 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 is bounded for all finite time, i.e., . The boundedness of implies also the boundedness of and , 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 defined in that has at least one non-constant component. It is clear that for any solution of (2) either all functions , and are constant or all depend on effectively. The next observation is that the function doesn’t change the sign. In fact, if for some , then it follows from the first equation that , which is impossible. Let us solve the system


with respect to and .

Lemma 3.1.

If , then there exists a unique solution of (3), which is continuous and positive for all . In the case there is only one positive solution and for there are two positive solutions that exist only in some interval , where is a finite number.


Let us introduce the temporary notation . First observe that


and hence is positive as soon as is positive. The main question is whether there in fact exists a positive . Substituting expression (4) into the first equation of (3) yields


This quadratic equation with respect to admits real roots provided

Note that for all positive , if and only if . Hence, under the assumption , equation (5) has the solution


which is positive for all . Furthermore,


In the the critical case there exists a unique solution of (5)


where is positive for . In the case there are two positive solutions of (5), defined for with , . But only one of them is bounded at , and this solution is given by (6). Moreover the solution is also bounded at the point :


Note that is the smallest positive point in which the discriminant vanishes. Observe also that for all . ∎

We next come back to system (2). Let . In view of Lemma 3.1 the function must be a positive solution of the nonlinear boundary value problem


where is given by (6), if , and by (8) otherwise. We must keep in mind that the function depends on parameters , , , , and .

3.1 Nonlinear boundary value problem

Let us consider the boundary value problem



and 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].

Figure 1: Plots of the energy and the trajectories of a dynamical system.
Theorem 3.2.

If function has no less than three positive roots, then there exists a set of diffusion constants for which boundary value problem (11) admits a positive solution.


We assume for a while that . Suppose that the potential energy

has a positive critical point that is a local minimum of . The trajectories of dynamic system , resemble ellipses near the point of the phase space. Let us consider the trajectory that intersects the axis at positive points and (see Fig. 1). Therefore there exists a periodic solution , of the dynamic system subject to initial conditions , , and its period is given by

Remark that the energy without local minimum points is either monotonic or a function with a single maximum point. In both cases no trajectory starting at a point can return again to the line , with the exception of the equilibrium positions.

Next, is a positive solution of the initial problem , , . Moreover, for any natural , because for even and for odd . Given , we consider the function that is obviously the solution to equation , and . We set , so that

Hence, is a positive solution to (11).

It follows that any close trajectory lying in the half-plane produces a countable set of positive solutions to (11) with . All these sequences form the set , which is in general uncountable.

The first and primary question is whether the energy has a local minimum. We consider first the case . The energy has a local minimum at point if only , in a left-side neighborhood of and in a right-side one. Let

We look for a root of the equation for which the left of and on its right. A trivial verification shows that (given by (6)) is a steadily increasing function and


Set and . It is convenient to consider the -representation of functions on interval :

where . Both functions are strictly increasing in the region under study. In regard to , its derivative can be written as

and is obviously positive for . The function has two vertical asymptotes and , while is bounded. In addition, has an inflection point , whereas the inflection point of depends on parameters of the model. Fig. 2 shows the typical plots of and . Evidently, the energy has a local minimum if plots of the functions intersect at three points.

Similarly, for we obtain that is monoton increasing and

The function is nearly linear for small and growth quadratically as , whereas for such that , and is bounded for ; additionally, . Thus, for there exist sets of parameters , such that the function has three roots in .

Figure 2: Sample plots of and corresponding to nonexistence and existence of a nonhomogeneous stationary solution in the left and right panels, respectively.

For and defined by (6), we obtain

and again is nearly linear, for small , , and both functions are bounded on . For there exist two roots, as functions of parameters, of on the interval . The third root of lies on the interval , where is the largest positive point in which vanishes. ∎

Remark 1.

Let , the solution of (5) is defined by


parameters , , , , satisfy the assumption , and is choosen such that . Then the function is continuous in , and as . This properties of and the assumption provide existence of a local minimum of the energy at a point . Thus, there exists a set of diffusion constants for which the boundary value problem (11), with defined by (13), admits a positive solution. This situation was considered in [9] with .

If has a local minimum , then there are a local maximum in the interval and a local maximum in the interval as shown in Fig. 1. Set .

Corollary 1.

For any there exists a countable set of solutions to boundary value problem (11) with and as , such that

for all . Moreover is rapidly oscillating function as goes to .

Figure 3: Spatially nonhomogeneous solutions of the stationary problem for the different values of . Parameters: , , , , , , , and and in the left and right panels, respectively.

4 Stability and destabilisation of steady states

Let be a complex Banach space with norm and let be a closed operator on with a dense domain . Suppose that is a sectorial operator and for all . Hereafter, stands for the spectrum of an operator . For each we introduce the interpolation space equipped with the norm . Clearly . We consider the equation


in the Banach space . Let 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 be a locally Lipschitz continuous map in a neighborhood of a steady state , for some . Suppose that for the map admits the representation


provided is small enough. If

  • is a linear bounded operator from to ,

  • as ,

  • the spectrum of lies in the set for some ,

then the equilibrium solution of (14) is asymptotically stable in . Moreover, if the spectrum and the half-plane have a non-empty intersection, then is unstable.

4.1 Linearised analysis

For simplicity of notation, we denote the vector by so that system (1) reads


in the Banach space , where

and is the Sturm-Liouville operator given by

on the interval , subject to the Neumann boundary conditions,

We have that is a closed densely defined operator on and the spectrum of is given by

Also for we have the estimate

Hence is a sectorial operator, see e.g. [5], and for all .

For a sectorial operator we can consider interpolation spaces , for , each of which is a Banach subspace of , see e.g. [5]. Set .

The function is smooth in . Therefore, admits the representation

where the remainder satisfies the estimate


in a neighborhood of any point with a continuous function . Here

Suppose now that is a positive steady state, that is and for all . Assume also that are -functions, . Obviously, is a linear bounded operator from to , for each .

Next, using (17) we obtain with the estimate of the remainder

Consequently, for any using the fact that we conclude

as , with constants , where , being independent on , and .

For nonlinear model (1) the linearized system at the steady state can be written as . Then turning back to the previous notation, we obtain


The first and primary question is whether the spectrum of lies in the half-plane for some .

4.2 Spectral analysis of linearised problem

Let us consider the Sturm-Liouville operator in the domain . For abbreviation, we write instead of , and consider the operator


where . Of course, constant and functions , , are positive in , since the stationary solution is positive.

The matrix is regarded as an unbounded operator in the space with the domain . It can be written in the form , where

Notice that is a bounded operator in .

Theorem 4.1.

Let , where is the largest eigenvalue of . The spectrum of lies in the set and .


To analyse the spectrum of we consider


for each . Since the last equation in (19) is independent of and , for all and all we have a unique solution . Notice that due to the Sobolev embedding theorem .

If , , and , then we have a unique solutions of (19) in , given by


Here denotes the range of the continuous function on . Thus is a point of the resolvent set of .

If either , , or , then the operator does not have a solution for some or and , defined in (20), are not continuous for some .

Hence, the spectrum of coincides with the set

where are eigenvalues of the Sturm-Liouville operator subject to the Neumann boundary conditions.