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

## Abstract

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

yu_holovaty@franko.lviv.ua

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

anna.marciniak@iwr.uni-heidelberg.de

Mariya Ptashnyk

Department of Mathematics

University of Dundee

DD1 4HN Dundee, Scotland, UK

m.ptashnyk@dundee.ac.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,

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂tc=θbcb+c−dcc+μ in (0,∞)×(0,1),∂tb=αc2g−dbb−νb in (0,∞)×(0,1),∂tg=1γ∂2xg−αc2g−dgg+βck1+ck+νb in (0,∞)×(0,1),∂xg(t,0)=∂xg(t,1)=0 in (0,∞), (1)

with initial conditions

 c(0,x)=cin(x),b(0,x)=bin(x),g(0,x)=gin(x) for x∈(0,1),

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), .

###### Proof.

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

 ∂tc≤(θ−dc)c+μ,

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

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩θbcb+c−dcc+μ=0,αc2g−dbb−νb=0,1γg′′−αc2g−dgg+βck1+ck+νb=0,g′(0)=g′(1)=0. (2)

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

 ⎧⎨⎩θbcb+c−dcc+μ=0,αc2g−dbb−νb=0 (3)

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.

###### Proof.

Let us introduce the temporary notation . First observe that

 b(g)=1μℓ(g)c2(g), (4)

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

 (θ−dc)ℓc2+μ(ℓ−dc)c+μ2=0. (5)

 D=μ2((ℓ−dc)2−4(θ−dc)ℓ)=μ2((ℓ+dc)2−4θdcdcℓ)≥0.

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

 c(g)=μ⋅ℓ(g)−dc+√(ℓ(g)−dc)2+4(dc−θ)ℓ(g)2(dc−θ)ℓ(g), (6)

which is positive for all . Furthermore,

 c(0+)=μdc,c(+∞)=μdc−θ. (7)

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

 c(g)=μdc−ℓ(g), (8)

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 :

 c(g∗)=μ(dc−ℓ(g∗))2(θ−dc)ℓ(g∗). (9)

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

 ⎧⎪⎨⎪⎩1γg′′=g(dg+ωc2(g))−βck(g)1+ck(g) in (0,1),g′(0)=0,g′(1)=0, (10)

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

 g′′=γh(g),x∈(0,1),g′(0)=0,g′(1)=0, (11)

where

 h(g)=dgg(1+ωc2(g))−βck(g)1+ck(g),ω=αdbdg(db+ν),

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

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

###### Proof.

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

 U(g)=−∫g0h(ξ)dξ

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

 T=2∫q2q1dη√2(U(q1)−U(η)).

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

 u′(1)=1√γng′(1√γn)=1√γng′(nT2)=0.

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

 h1(g)=dgg(1+ωc2(g)),h2(g)=βck(g)1+ck(g).

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

 c(+0)=μdc,c(g)→μdc−θas g→+∞. (12)

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

 h1(c)=a(c−cmin)(1+ωc2)c(cmax−c),h2(c)=βck1+ck,

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

 h′1(c)=a⋅(1+ωc2)(cmin(cmax−c)+c(c−cmin))+2ωc2(c−cmin)(cmax−c)c2(cmax−c)2

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

 c(+0)=μdc,c(g)→∞ as g→dc(db+ν)αμ.

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 .

For and defined by (6), we obtain

 c(+0)=μdc,c(g∗)=μ(dc−l(g∗))2(θ−dc)l(g∗),

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

 c(g)=μ⋅ℓ(g)−dc−√(ℓ(g)−dc)2+4(dc−θ)ℓ(g)2(dc−θ)ℓ(g), (13)

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

 |gn(x)−g0|≤ε

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

## 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

 ∂tu=Au+f(u) (14)

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

 f(u∗+z)=f(u∗)+Bz+p(u∗,z), (15)

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

 ∂tu=Au+f(u) (16)

in the Banach space , where

 A=⎛⎜⎝−dc000−(db+ν)000N0⎞⎟⎠,f(u)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝θu1u2u1+u2+μαu21u3βuk11+uk1+νu2−αu21u3⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠

and is the Sturm-Liouville operator given by

 N0v=1γd2vdx2−dgv

on the interval , subject to the Neumann boundary conditions,

 D(N0)={v∈H2(0,1):v′(0)=0,v′(1)=0}.

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

 σ(A)={−dc,−db−ν}∪{−dg−γ−1π2j2}∞j=0.

Also for we have the estimate

 ∥(A−λE)∥X≤1dist(λ,σ(A)).

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

 f(y+z)=f(y)+B(y)z+p(y,z),

where the remainder satisfies the estimate

 ∥p(y,z)∥R3≤ϑ(y)∥z∥2R3 (17)

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

 B(y)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝θy22(y1+y2)2θy21(y1+y2)202αy1y30αy21kβyk−11(1+yk1)2−2αy1y3ν−αy21⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

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

 ∥p(u∗,z)∥R3≤ϑ(u∗(x))∥z(t,x)∥2R3≤maxx∈[0,1]ϑ(u∗(x))∥z(t,x)∥2R3≤ϑmax∥z∥2R3.

Consequently, for any using the fact that we conclude

 ∥p(u∗,z)∥X≤c1(∥z1∥2C[0,1]+||z2||2C[0,1]+||z3||2H1(0,1))≤c1(∥z1∥2C[0,1]+||z2||2C[0,1]+||z3||2H2s(0,1))≤c2∥z∥2s=o(∥z∥s),

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

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂tz1=(θb2∗(c∗+b∗)2−dc)z1+θc2∗(c∗+b∗)2z2,∂tz2=2αc∗g∗z1−(db+ν)z2+αc2∗z3,∂tz3=(kβck−1∗(1+ck∗)2−2αc∗g∗)z1+νz2+(1γ∂2x−αc2∗−dg)z3. (18)

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

 A=⎛⎜⎝l(x)m(x)0n(x)−ϰq(x)r(x)νN⎞⎟⎠

with

 ϰ=db+ν,n=2αc∗g∗,q=αc2∗, l=θb2∗(b∗+c∗)2−dc,m=θc2∗(b∗+c∗)2,r=kβck−1∗(1+ck∗)2−2αc∗g∗,

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

 A0=⎛⎜⎝lm00−ϰq00N⎞⎟⎠andA1=⎛⎜⎝000n00rν0⎞⎟⎠.

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 .

###### Proof.

To analyse the spectrum of we consider

 ⎧⎨⎩(l−λ)ϕ1+mϕ2=f1,−(ϰ+λ)ϕ2+qϕ3=f2,Nϕ3−λϕ3=f3, (19)

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

 ϕ1(x,λ) =f1(x)l(x)−λ+m(x)(f2(x)−q(x)ϕ3(x))(ϰ+λ)(l(x)−λ), (20) ϕ2(x,λ) =−f2(x)−q(x)ϕ3(x)ϰ+λ, ϕ3(x,λ) =(N−λ)−1f3.

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

 {−ϰ}∪{λj}∞j=1∪R(l),

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

Since