# Numerical study of vanishing and spreading dynamics of chemotaxis systems with logistic source and a free boundary

The current paper is to investigate the numerical approximation of logistic type chemotaxis models in one space dimension with a free boundary. Such a model with a free boundary describes the spreading of a new or invasive species subject to the influence of some chemical substances in an environment with a free boundary representing the spreading front (see Bao and Shen [1], [2]). The main challenges in the numerical studies lie in tracking the moving free boundary and the nonlinear terms from chemical. To overcome them, a front fixing framework coupled with finite difference method is introduced. The accuracy of the proposed method, the positivity of the solution, and the stability of the scheme are discussed.The numerical simulations agree well with theoretical results such as the vanishing spreading dichotomy, local persistence, and stability. These simulations also validate some conjectures in our future theoretical studies such as the dependence of the vanishing-spreading dichotomy on the initial solution u0, initial habitat h0, the moving speed ν and the chemotactic sensitivity coefficients ḩi̧1,ḩi̧2.

## Authors

• 74 publications
• 2 publications
01/21/2022

### A convergent numerical scheme to a parabolic equation with a nonlocal boundary condition

In this paper, a numerical scheme for a nonlinear McKendrick-von Foerste...
07/09/2019

### Error analysis of finite difference/collocation method for the nonlinear coupled parabolic free boundary problem modeling plaque growth in the artery

The main target of this paper is to present a new and efficient method t...
06/03/2019

### Robust stability of moving horizon estimation for nonlinear systems with bounded disturbances using adaptive arrival cost

In this paper, the robust stability and convergence to the true state of...
11/30/2021

### A structure preserving front tracking finite element method for the Mullins–Sekerka problem

We introduce and analyse a fully discrete approximation for a mathematic...
07/04/2019

### Analysis and numerical simulation of the nonlinear beam equation with moving ends

The numerical analysis for the small amplitude motion of an elastic beam...
10/11/2020

### Convergence of a time-stepping scheme to the free boundary in the supercooled Stefan problem

The supercooled Stefan problem and its variants describe the freezing of...
10/12/2019

### Optimization of One-parameter Family of Integration Formulae for Solving Stiff Chemical-kinetic ODEs

A fast and robust Jacobian-free time-integration method - called Minimum...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The current paper is to study, in particular, numerically, the spreading and vanishing dynamiccs of the following attraction-repulsion chemotaxis system with a free boundary and logistic source,

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ut=uxx−χ1(uv1,x)x+χ2(uv2,x)x+u(a(t,x)−b(t,x)u),0

where in (1.1) is a positive constant, , , and () are nonnegative constants, and and satisfy the following assumption,

(H0) and are bounded functions on , and

 ainf:=inft∈R,x∈[0,∞)a(t,x)>0,binf:=inft∈R,x∈[0,∞)b(t,x)>0.

Biological backgrounds of (1.1) are discussed in the paper ([1]). The free boundary condition in (1.1) is also derived in [1] based on the consideration of “population loss” at the front which assumes that the expansion of the spreading front is evolved in a way that the average population density loss near the front is kept at a certain preferred level of the species, and for each given species in a given homogeneous environment, this preferred density level is a positive constant determined by their specific social and biological needs, and the environment.

One of the first mathematical models of chemotaxis were introduced by Keller and Segel ([13], [14]) to describe the aggregation of certain type of bacteria in 1970. Since their publications, considerable progress has been made in the analysis of various particular case of chemotaxis (Keller-Segel) model on both bounded and unbounded fixed domain (see [3], [5], [6], [9], [10], [12], [25], [30], [31], [33], [34], [35], [36], [37], [38], [39], [40], and the references therein). Among the fundamental problems in studying chemotaxis model are the existence of nonnegative solutions which are globally defined in time or blow up at a finite time and the asymptotic behavior of time global solutions.

Du and Lin studied the population invasion represented by Fisher-KPP free boundary problem in 2010 [7]. The breaking difference between the asymptotic behaviours of Fisher-KPP with a free boundary and on the fixed or fixed unbounded domain is the vanishing-spreading dichotomy, which is well supported by some empirical evidences, for example, the introduction of several bird species from Europe to North America in the 1900s was successful only after many initial attempts (see [23],[29]).

Compared to the studying chemotaxis model on fixed bounded or fixed unbounded domain and the asymptotic behaviour of Fisher-KPP equation with a free boundary, the central problems in studying system (1.1) are the existence of nonnegative solutions which are globally defined in time, the vanishing-spreading dichotomy, local persistence, local stability, and so on.

To state the main results of the current paper, we first recall some theoretical results proved in [2] which will all be validated in our numerical simulations.

Let

 Cbunif(R+)={u∈C(R+)|u(x)is% uniformly continuous and bounded onR+}

with norm , and

 Cbunif(R)={u∈C(R)|u(x)is % uniformly continuous and bounded onR}

with norm . Define

 M=min{ 1λ2((χ2μ2λ2−χ1μ1λ1)++χ1μ1(λ1−λ2)+), 1λ1((χ2μ2λ2−χ1μ1λ1)++χ2μ2(λ1−λ2)+)} (1.2)

and

 K=min{ 1λ2(|χ1μ1λ1−χ2μ2λ2|+χ1μ1|λ1−λ2|), 1λ1(|χ1μ1λ1−χ2μ2λ2|+χ2μ2|λ1−λ2|)}. (1.3)

Let (H1)- (H3) be the following standing assumptions.

(H1) .

(H2) .

(H3) .

Note that

 M≤χ2μ2.

Hence implies (H1). In the case , we can choose , and then and . Hence (H1) becomes , (H2) becomes , and (H3) becomes . In the case , we can also choose , and then and . Hence (H1) (resp.(H2), (H3)) becomes . Biologically, (H1), (H2), and (H3) indicate that the chemo-attraction sensitivity is relatively small with respect to logistic damping.

When (H1) holds, we put

 M0=asupbinf+χ2μ2−χ1μ1−M (1.4)

and

 m0=ainf(binf−(1+asupainf)χ1μ1+χ2μ2−M)(binf−χ1μ1+χ2μ2−M)(bsup−χ1μ1+χ2μ2). (1.5)

Note that if (H2) holds, then .

Let

 H(a,b)=cl{(a(t+⋅,⋅),b(t+⋅,⋅))|t∈R}

with open compact topology, where the closure is taken under the open compact topology.

The main results of the paper [2] are stated in the following. The first result is on the global existence of nonnegative solutions of (1.1).

Global existence [2, Theorem 1.2]: If (H1) holds, then for any , and any and any function on satisfying

 u0∈C2[0,h0],u0(x)≥0forx∈[0,h0],andu′0(0)=0,u0(h0)=0, (1.6)

(1.1) has a unique globally defined solution , , , with and . Moreover,

 0≤h′(t)≤2νM1C0, (1.7)
 0≤u(t,x;t0,u0,h0)≤max{∥u0∥∞,M0}∀t∈[t0,∞),x∈[0,h(t;t0,u0,h0)) (1.8)

and

 limsupt→∞supx∈[0,h(t;t0,u0,h0))u(t0+t,x;t0,u0,h0)≤M0, (1.9)

where is a big enough constant and .

Assume (H1). For any given , and any given and satisfying (1.6), by the nonnegativity of , for all . Hence exists. Put

 h∞(t0,u0,h0)=limt→∞h(t;t0,u0,h0).

We say vanishing occurs if and

 limt→∞∥u(t,⋅;t0,u0,h0)∥C([0,h(t;t0,u0,h0)])=0.

We say spreading occurs if and for any ,

 liminft→∞inf0≤x≤Lu(t,x;u0,h0)>0.

For given , consider the following linear equation,

 {vt=vxx+a(t,x)v,0

Let be the principal spectrum interval of (1.10) (see Definition 2.1 [2]). Let be such that for and (see Lemma 2.2, 2.3 [2] for the existence and uniqueness of ).

The second result is about the spreading and vanishing dichotomy scenario in (1.1).

Spreading-vanishing dichotomy [2, Theorem 1.3]: Assume that (H1) holds. For any given , and and satisfying (1.6), we have that either (i) vanishing occurs and ; or (ii) spreading occurs.

For given , and and satisfying (1.6), if spreading occurs, it is interesting to know whether local uniform persistence occurs in the sense that there is a positive constant independent of the initial data such that for any ,

 liminft→∞inf0≤x≤Lu(t,x;t0,u0,h0)≥~m0,

and whether local uniform convergence occurs in the sense that exists locally uniformly. We have the following result along this direction.

Persistence and convergence [2, Theorem 1.4]: Assume that (H1) holds and that and satisfy (1.6).

(i) (Local uniform persistence) For any given , if and (H2) holds, then for any ,

 liminft→∞inf0≤x≤Lu(t,x;t0,u0,h0)>m0,

where is as in (1.5).

(ii) (Local uniform convergence) Assume that (H3) holds, and that for any , there has a unique strictly positive entire solution , . Then for any given , if , there are such that to any , , for any ,

 limt→∞sup0≤x≤L|u(t,x;t0,u0,h0)−u∗(t,x;a,b)|=0. (1.11)

(iii) (Local uniform convergence) Assume that (H3) holds, and that and . Then for any given , if , then for any ,

 limt→∞sup0≤x≤L|u(t,x;u0,h0)−u∗(t)|=0,

where is the unique strictly positive entire solution of the ODE

 u′=u(a(t)−b(t)u) (1.12)

(see [11, Lemma 2.5] for the existence and uniqueness of strictly positive entire solutions of (1.12)).

###### Remark 1.1.

Biologically, the invasion or spreading of the population is depending on the initial solution, initial habitat, the moving speed ([7]). When the spreading happens, the local persistence and convergence can be guaranteed in Fisher-KPP equation with a free boundary, furthermore, there is a asymptotic spreading speed such that (see [7], [16]).

Compared to the vanishing-spreading dichotomy in Fisher-KPP equation with a free boundary, the chemotaxis system (1.1) do not have comparison principle which leads the following interesting open problems but has positive answers in Fisher-KPP case.

1. For given and , whether there is such that for , vanishing occurs, and for , spreading occurs.

2. For given , , and , whether there is such that for with , vanishing occurs, and for with , spreading occurs.

3. Whether there is a spreading speed such that as long as the spreading occurs.

###### Remark 1.2.

Vanishing-spreading result [2, Theorem 1.2] indicate there is a separating value , which is independent of the chemotactic sensitivity coefficients , such that in the vanishing scenario the limiting moving boundary and when the initial habitat the spreading guaranteed. The dependence of the dynamics of the system on the chemotactic sensitivity coefficients is another important and interesting questions [28], [38]. We also have the following question in this direction.

4.

If the asymptotic spreading speed exist, whether the limit depends on the chemotactic sensitivity coefficient and .

The objective of the current paper is to study the numerical effect of the parameters on the vanishing and spreading dynamics in the system (1.1) which will give us directions in the theoretical studies. For simplicity, we only consider the constant logistic coefficients in the system, where . In general, it is always difficult to handle the attraction term in chemotaxis system which may lead to convection dominant in the system [4], [19], [20], [26]. However, in the system (1.1), we have extra numerical challenges in efficiently and accurately handling the moving boundaries [24]. These two challenges require us to construct a new numerical algorithm in the numerical study.

Thanks to the maximum principle in the elliptic equation and death damping coefficient in the parabolic equation, the chemoattraction term can be controlled by the magnitude of the population density which has a global bound by the suppressing of the death rate. The front-fixing method has been successfully applied to solve one dimensional free boundary problem [15], [21], [22], [24] which changes the moving boundary to fixed domain, and is the main concern in our numerical studies of (1.1). Combined with the front-fixing method, finite difference in parabolic and finite volume method in the elliptic equations in the system (1.1), we construct a new algorithm in the numerical study, as a by-product we also obtain the consistency, monotonicity of the moving boundary, positivity of the solution and stability results.

Our numerical experiments validate the vanishing and spreading dichotomy in the numerical scheme of system (1.1) which is similar to Fisher-KPP equation with a free boundary and give evidences to our conjectures that:

1. For given and , there is such that for , vanishing occurs, and for , spreading occurs. Which means in order to spread to the half space , the moving speed should be large enough and otherwise the population will be extinct.

2. For given , , and , there is such that for with , vanishing occurs, and for with , spreading occurs. Biologically, large initial population density helps the establishment and spreading of invasion which is an indirect evidence to the early birds introduction problem in 1900s [23], [29].

3. There is a spreading speed such that as long as the spreading occurs, which is independent of chemotactic sensitivity coefficient and . Chemical and are produced by the species and the density is close to zero near the spreading front. In such case the decisive effect of the spreading speed should not depend on the chemotactic sensitivity coefficients and .

The rest of this paper is organized in the following way. In section 2, we first use Landau transformation to transfer the moving boundary to a fixed domain, then we use the finite difference, finite volume, and iteration method to approximate the continuous chemotaxis system. We also prove the monotonicity of the moving boundary, the positivity and stability of the discrete solutions. In section 3, we study the numerical spreading-vanishing dichotomy in (1.1) which validates our theoretical results (Vanishing-spreading dichotomy, local persistence and convergence). Our simulations also indicate the dependence or independence of the vanishing-spreading dichotomy on parameters and so on. In section 4, some future works are briefly discussed.

## 2 Numerical approximation of the free boundary problem

In this section, we study the numerical approximation of system (1.1) with constant logistic coefficients . First through the well-known Landau transformation (see [15]), we convert (1.1) into a fixed spatial domain problem. In such a way the length of the moving boundary is included as another variable to be solved apart from the population density. Then we solve the converted new problem on the basis of finite difference and finite volume method. There is a circulation that each time when time variable increase we solve the elliptic equations first and using the forward differential method to find the solution of the parabolic equation for the new time.

From the elliptic equations in (1.1), we know that

 ∂xxv1 = λ1v1−μ1u, (2.1) ∂xxv2 = λ2v2−μ2u. (2.2)

Combining (2.1), (2.2) and the first equation in (1.1), we have

 ut=uxx+(−χ1v1x+χ2v2x)ux+(−χ1λ1v1+χ2λ2v2+a)u+(χ1μ1−χ2μ2−b)u2. (2.3)

Then we introduce the Landau transformation,

 z(t,x)=xh(t),w(t,z)=u(t,x),V1(t,z)=v1(t,x),V2(t,z)=v2(t,x). (2.4)

Under this substitution the elliptic equations (2.1), (2.2) take the form:

 ∂2V1∂z2⋅1h2(t)−λ1V1+μ1w=0,0

The elliptic boundary conditions are

 V1,z(t,0)=V2,z(t,0)=0, (2.7) V1,z(t,1)=V2,z(t,1)=0. (2.8)

Equation (2.3) takes the form:

 ∂w∂t+∂w∂z(−h′(t)h(t)z)=∂2w∂z21h2(t)+(−χ1V1z+χ2V2z)∂w∂z1h2(t)+(−χ1λ1V1+χ2λ2V2+a)w+(χ1μ1−χ2μ2−b)w2,0

Let denote and multiply it on both sides of the above equation.

 G(t)∂w∂z+∂w∂z(−G′(t)⋅z2)=∂2w∂z2+(−χ1V1z+χ2V2z)∂w∂z+(−χ1λ1V1+χ2λ2V2+a)w⋅G(t)+(χ1μ1−χ2μ2−b)w2⋅G(t),0

Boundary conditions and Stefan condition take the form

 ∂w∂z(t,0)=0,w(t,1)=0,t>0 (2.11)

and

 G′(t)=−2ν∂w∂z(t,1),t>0. (2.12)

The initial conditions in (1.1) become:

 G(0)=h20,w(0,z)=w0(z)=U0(z⋅h0),0≤z≤1, (2.13)

and the initial function is changed into which maintains:

 w′0(0)=w0(1)=0,w0(z)>0,0≤z<1. (2.14)

Under the transformation, our aim is to solve the nonlinear parabolic partial differential system (2.5), (2.6), (2.10) in the fixed domain for the variables .

The following is the process according to the theory of the finite difference method. First we consider the time and space discretization , , which means the interval (0,1) is divided into equal cells

 0=z0

and the mesh points , with . For abbreviation, the approximate value of can be denoted by , the approximate values of and can be denoted by and . Besides we write for the value of . Let us consider the central approximation of the spatial derivatives,

 Vn1,j+1−Vn1,j−12h≈∂V1∂z(tn,zj),Vn2,j+1−Vn2,j−12h≈∂V2∂z(tn,zj), (2.15)
 wnj+1−wnj−12h≈∂w∂z(tn,zj),wnj−1−2wnj+wnj+1h2≈∂2w∂z2(tn,zj). (2.16)

By using the forward approximation of the time derivative, we get

 wn+1j−wnjτ≈∂w∂t(tn,zj),gn+1−gnτ≈G′(tn). (2.17)

Let us apply the approximation (2.15) on the elliptic equations (2.5), (2.6), then we get

 1G(t)1h2⋅Vn1,j−1+(1G(t)⋅−2h2−λ1)Vn1,j+1G(t)1h2⋅Vn1,j+1=−μ1wnj,1≤j≤M−1, (2.18) 1G(t)1h2⋅Vn2,j−1+(1G(t)⋅−2h2−λ2)Vn2,j+1G(t)1h2⋅Vn2,j+1=−μ2wnj,1≤j≤M−1. (2.19)

We mainly focus on solving the values of , the relevant results about can be obtained in a similar way. For (2.15), , we can get equations for and there are two other equations we need for the boundary. In order to achieve a higher accuracy (see [18]), we use the idea of finite volume method to handle the boundary conditions.

We get the equation for the left boundary:

 (h2d0+a1h)V1,0−a1hV1,1=h2ϕ0, (2.20)

where

 a1=(1h∫z1z0−G(tn)dz)−1,d0=2h∫z1/2z0−λ1dz,ϕ0=2h∫z1/2z0−μ1wndz,

and .

The equation for the right boundary is similar:

 −anhV1,M−1+(anh+h2dn)V1,M=h2ϕn, (2.21)

where

 an=(1h∫zMzM−1−G(tn)dz)−1,dn=2h∫zMzM−1/2−λ1dz,ϕn=2h∫zMzM−1/2−μ1wndz,

and .

By now we have equations which is enough to form a system of linear algebraic equations for . It is tridiagonal and there exists the unique solutions . Similarly, we can acquire the system of linear algebraic equations about and its corresponding solutions .

With the information of and , we can concentrate on solving the parabolic equation (2.10). From (2.15), (2.16), and (2.17), (2.10) is approximated by

 gnwn+1j−wnjτ−zj2wnj+1−wnj−12hgn+1−gnτ =wnj−1−2wnj+wnj+1h2+(−χ1Vn1,j+1−Vn1,j−12h+χ2Vn2,j+1−Vn2,j−12h)wnj+1−wnj−12h +(−χ1λ1Vn1,j+χ2λ2Vn2,j+a)wnjgn+(χ1μ1−χ2μ2−b)(wnj)2gn, (2.22)

for . Because of the initial conditions (2.14), we can assume a fictitious value at the point , and then

 wn1−wn−12h=0,wnM=0,n≥0. (2.23)

Considering the Stefan condition (2.12), according to (2.17) and three points backward spatial approximation of , we obtain

 gn+1−gnτ=−νh(3wnM−4wnM−1+wnM−2),n≥0. (2.24)

Because of (2.13), it can also be written as:

 gn+1=gn+τνh(4wnM−1−wnM−2),n≥0. (2.25)

Let us replace with (2.25) in (2.22), we get the explicit scheme:

 wn+1j = (−zj4hτν(4wnM−1−wnM−2)hgn+τgnh2−S1τ4h2gn)wnj−1 (2.26) +(1−2τgnh2+S2τ)wnj+S3τ(wnj)2 +(−zj4hτν(4wnM−1−wnM−2)hgn+τgnh2+S1τ4h2gn)wnj+1,

for where

 S1 = −χ1(Vn1,j+1−Vn1,j−1)+χ2(Vn2,j+1−Vn2,j−1), (2.27) S2 = −χ1λ1Vn1,j+χ2λ2Vn2,j+a, (2.28) S3 = χ1μ1−χ2μ2−b. (2.29)

The solution of (2.10) is classical (see [2, Lemma 3.1]

), we can obtain a optimal error estimates for the numerical approximation. Consider (

2.10), (2.11) and (2.12), we denote that

 L1(w,V1,V2,G) = ∂w∂t−z2G′(t)G(t)∂w∂z−1G(t)∂2w∂z2−1G(t)(−χ1∂V1∂z+χ2∂V2∂z)∂w∂z (2.30) −(−χ1λ1V1+χ2λ2V2+a)w−(χ1μ1−χ2μ2−b)w2=0, L2(w,V1,V2,G) =