Positivity-preserving high-order compact difference method for the Keller-Segel chemotaxis model

1.   Introduction

Chemotaxis, which is the phenomenon of directional movement for cells or microorganisms along the concentration gradient of chemical stimuli in the tissue or living environment, plays a key role in many complex biological processes such as tumor growth and development, disease state transfer, human tissue development, pattern formation, bacterial aggregation and so on ^[1]. As early as the 1950s–1980s, many researchers studied the phenomena of chemotaxis in biology and proposed a series of partial differential equations (PDEs) models related to chemotaxis ^{[2,3,4,5,6,7]}. The most classic is the chemotaxis model proposed by Keller and Segel ^[3,4] (namely the Keller-Segel chemotaxis model) to simulate the aggregation phenomena of amoeba and dictyostelium discoideum. In this paper, we consider the following common formulation of the Keller-Segel chemotaxis model, which is written in the dimensionless form as ^[3,7,8]

with initial and homogeneous Neumann boundary conditions

The meaning of each symbol is as follows

The total mass $\int_\Omega u(x, t)dx = \int_\Omega u(x, 0)dx$ is a constant as the temporal evolution under the boundary conditions, so the system (1.1)–(1.3) is isolated.

Since the Keller-Segel chemotaxis models were proposed, it had fascinated a large number of researchers in the mathematics and biology circles. So far, many types of chemotaxis models have been developed. Such as, the minimal model ^[9], the quasilinear Keller-Segel model ^[10], the Keller-Segel model with logistic source ^[11,12] and so on (More detailed introduction can be seen in ^[13,14,15] and references therein). A common property of these chemotaxis models is their ability to mathematically cause the solutions to grow rapidly in the small neighborhood of the concentration point/curve. The solutions may blow up or exhibit very strange spike behavior ^[16,17]. This blow up represents a mathematical description of the cell aggregate phenomenon that occurs in real biological system ^{[7,18,19,20,21,22]}. In order to employ the Keller-Segel model to explore the phenomenon of cell aggregation, usually there are two theoretical research focuses: one is to prove that the Keller-Segel parabolic equation blowns up in a finite time, and the $\delta$-function generated by the blown-up solution is used to describe it ^[7]. The other is to prove that the global solution of the system (1.1)–(1.3) is uniformly bounded with respect to time and tends to the stationary solution (such as spikes and patterns, etc.) under a large time ^[8,23,24]. For example, in one-dimensional (1D) space, the homogeneous Neumann initial-boundary value problem (IBVP) for the Keller-Segel chemotaxis model has a global solution and the solution is bounded ^[25,26].

Capturing the numerical solutions of such blowup or spikes is still a very challenging task ^[17], and it is very important for the modeling and analysis of chemotaxis system to establish accurate and efficient numerical methods. The research of numerical solutions can not only be intuitively observed by the dynamic behavior of the chemotaxis system, but also verify and supplement the theoretical results. In recent years, the research on the numerical methods of the Keller-Segel chemotaxis model mainly employs finite element methods (FEMs) ^{[27,28,29,30,31,32,33]} and finite volume methods (FVMs) ^{[16,34,35,36]}, as well as some other fractional step (operator splitting) methods ^[37,38], etc. But finite difference methods (FDMs) ^[39,40] are rarely used. For example, (ⅰ) In terms of FEMs, Epshteyn and Kurganov ^[31] proposed a series of new internal penalty discontinuous Galerkin (DG) methods for solving the Keller-Segel chemotaxis model, and simulated the pre-blow-up time, non-negativity and non-oscillation of the solutions. Li et al.^[28] improved the method in ^[31] by using the locally discontinuous Galerkin (LDG) method, and gave the optimal convergence rate of the Keller-Segel chemotaxis model in the special finite element space before blow-up. (ⅰ) In terms of FVMs, Filbet ^[34] employed the FVM to study the Patlak-Keller-Segel chemotaxis model, and simulated the blow-up characteristics of the chemotaxis model under the large mass of cells. Chertock and Kurganov ^[16] proposed a second-order finite-volume center-upward scheme based on the work in ^[35], and simulated the blow-up problem of the Keller-Segel model in the corners and the center area, as well as the pattern formation of bacteria in the semi-solid state and liquid medium, respectively. (ⅲ) In terms of FDMs, Saito and Suzuki ^[39] developed a finite difference scheme for the Keller-Segel chemotaxis model with the parabolic-elliptic coupling, which satisfies the conservation of a discrete $L_1$ norm error by applying the upwind technique. In order to develop a numerical scheme to satisfy conservation of positivity and total mass, Saito ^[40] proposed a conservative FDM for the Keller-Segel chemotaxis model with nonlinear diffusion.

Unfortunately, most of the numerical methods of the Keller-Segel chemotaxis models (1.1)–(1.3) developed above are low accuracy ^{[16,27,28,30,31,32,34,35,37,38,39,40]}. Although some numerical methods can achieve third-order ^[33,36] or fourth-order ^[17,29] accuracy in space, but the highest accuracy is third-order in time. In addition, many numerical schemes are non-compact on spatial grid, and conditionally stable, that is, a small time step size needs to be taken to satisfy its stability condition and to obtain a high-resolution numerical solution. As we all know, a high-order compact (HOC) difference scheme has the characteristics of small computational error, better numerical stability, and small grid-point stencil, which has aroused extensive concerns of many researchers. Backward differencing multi-step schemes have been proved to be the suitable method with highly accuracy and relatively large stability area ^[41,42]. Stiffly-stable schemes, which are backward differencing multi-step schemes in nature, appear first in ^[43]. In order to obtain the high accurate numerical solutions of the Keller-Segel chemotaxis model (1.1)–(1.3), in this paper, we attempt to develop a compact difference scheme, which is fourth-order accuracy in both time and space for the interior grid points and the boundaries, and stiffly-stable. To this end, first, we rewrite the Keller-Segel chemotaxis models (1.1)–(1.3) as its equivalent form: nonlinear advection diffusion reaction (ADR) system as follows

where

$[a, b]$ and $(0, T]$ represent space and time domains, respectively. $D$ is constant diffusion coefficient. Then, a fully implicit fourth-order compact difference scheme for the interior points is proposed for solving the above-mentioned system. And at the same time, we refer to the method in ^[44] to deal with homogeneous Neumann boundary condition (1.6), and develop a boundary scheme with fourth-order accuracy. However, since the unknown functions $u(x, t)$ and $v(x, t)$ in the Keller-Segel chemotaxis model have actual biological significance, their values cannot be negative, otherwise their biological significance will be lost. Therefore, we try to employ the methods in ^[33,45,46] to establish a suitable positivity-preserving algorithm for the developed HOC difference scheme, so that the value of the obtained cell density $u(x, t)$ is non-negative, and the fourth-order accuracy is still maintained for the positivity-preserving algorithm.

The outline of this paper is arranged as follows. In Section 2, the definition of the basic symbols and lemmas required in the subsequent sections are given. An HOC difference scheme to solve the system (1.1)–(1.3) is proposed, and the computational strategy for the nonlinear chemotaxis term is provided in Section 3. We design a positive-preserving numerical algorithm and a time advancement algorithm in Section 4. In Section 5, we choose several types of numerical examples to simulate the non-constant stationary solution of the Keller-Segel chemotaxis model, and demonstrate the accuracy, stability and positivity-preserving of the proposed method, respectively. The conclusion is given at the end of this paper.

2.   Preliminaries

First, we divide the domain $[a, b]\times(0, T]$ with positive integers $M$ and $N$. The space step size $h = (b - a)/N$, and the time step size $\tau = T/M$, respectively. ($x_i, t_n$) denotes the mesh points, where ${x_i} = a + ih, \; \; {t_n} = n\tau$, $0\leq i\leq N\;$, $\; 0\leq n \leq M$. Figure 1 represents the one-dimensional space-time discrete subdomain.

Then, we define $\Omega_h = \{x_i|0\leq i\leq N\}$, $\Omega_{\tau} = \{t_n|0\leq n\leq M\}$, $\Omega_{h\tau} = \Omega_h\times \Omega_{\tau}$. For any grid function $w \in \mathcal{W}_\tau = \{w^n|0\leq n\leq M\}$ defined on $\Omega_{\tau}$, denote

For any grid function $w \in \mathcal{W}_h = \{w_i|0\leq i\leq N\}$ defined on $\Omega_{h}$, denote

Next, for simplicity, we denote $\frac{\partial^pw}{\partial \varsigma^p}(x, t): = w_{\varsigma^p}(x, t)$, $p\in {\bf{Z}}^{+}$, $\varsigma$ may represents $x$ or $t$. In order to obtain the boundary schemes with the fourth-order accuracy, we use the Taylor expansion with Peano remainder to process the homogeneous Neumann boundary condition (1.3), and obtain the following theorem.

Theorem 2.1. For any function $w(x): \Omega_h\rightarrow \mathcal{W}_h$, we have

(1) If $w(x)\in C^5[x_0, x_4]$, and $w_x(x_0) = 0$, then

where $O_0(h^4) = -\frac{2}{25}h^4\int_0^1(1-s)^4[w_{x^5}(x_0+sh)-24w_{x^5}(x_0+2sh)+81w_{x^5}(x_0+3sh)-64w_{x^5}(x_0+4sh)]ds.$

(2) If $w(x)\in C^5[x_{N-4}, x_N]$, and $w_x(x_N) = 0$, then

where $O_N(h^4) = \frac{2}{25}h^4\int_0^1(1-s)^4[w_{x^5}(x_{N}-sh)-24_{x^5}(x_{N}-2sh)+81w_{x^5}(x_{N}-3sh)-64w_{x^5}(x_{N}-4sh)]ds.$

Proof. First, we expand $w(x_1)$, $w(x_2)$, $w(x_3)$, $w(x_4)$ at $w(x_0)$ by using Taylor expansion with Peano remainder ^[44]

That is

Second, we multiply the left and right sides of the above Eqs (2.4)–(2.6) by $\tilde{\alpha}$, $\tilde{\beta}$ and $\tilde{\gamma}$, and respectively add with the left and right sides of Eq (2.3) to obtain the following linear form,

where

We let the coefficients of the second, third, and fourth-order derivatives in Eq (2.7) equal to $0$, that is,

It is easy to solve the above linear Eqs (2.8)–(2.10) to obtain, $\tilde{\alpha} = -\frac{3}{4}$, $\tilde{\beta} = \frac{1}{3}$, $\tilde{\gamma} = -\frac{1}{16}$. Substitute them into Eq (2.7), and notice $w_x(x_0) = 0$, we have

then

where $O_0(h^4) = -\frac{2}{25}h^4\int_0^1(1-s)^4[w_{x^5}(x_0+sh)-24w_{x^5}(x_0+2sh)+81w_{x^5}(x_0+3sh)-64w_{x^5}(x_0+4sh)]ds.$ In the similar way, Eq (2.2) also can be proved.

For simplicity, for any grid function $w\in \mathcal{W}_h$, define the following average operator $\mathcal{A}$, i.e.,

In addition, we are inspired by Theorem 2.1, suppose $w(x, t)\in C^{6, 5} ([a, b]\times(0, T])$, for $1\leq i\leq N-1$, define

3.   Derivation for the HOC difference scheme

In this section, we derive an HOC difference scheme for solving the Keller-Segel chemotaxis models (1.1)–(1.3). Define a grid function $\mathbb{U} = \{\mathbb{U}_i^n|0\leq i\leq N, 0\leq n \leq M\}$ on $\Omega_{h\tau}$ to be the solution of (1.4)–(1.6), where $\mathbb{U}_i^n = \begin{bmatrix}u(x_i, t_n)\\ v(x_i, t_n)\end{bmatrix}$. Considering the semi-discrete difference form of Eq (1.4) at point $x_i$, we have

For the second derivative term on the right-hand side (RHS) of Eq (3.1), by using the fourth-order compact difference formula ^[47], i.e., ${\mathbb{U}_{xxi}} = {\mathcal{A}^{ - 1}}{\delta _x^2}{\mathbb{U}_i} + \mathcal{A}^{-1}O_{xxi}({h^4})$, $1\leq i\leq N-1$, and rearranging it, we have

3.1.   Fully implicit compact difference scheme

We consider Eq (3.2) at the $n$th time step, i.e.,

And discrete the time derivative term on the left-hand side (LHS) of Eq (3.3) by using the fourth-order backward differencing formula (BDF)^[43], i.e.,

Using the notations above, combining the initial and boundary conditions (1.5)–(1.6), and Lemma 2.1, we can obtain

where $Q_i^n = O_{ti}^n(\tau^4)+O_{xxi}^n(h^4)$, $Q_0^n = O_0^n(h^4)$, and $Q_N^n = O_N^n(h^4)$. Then, there exist constants $c_1$, $c_2$ and $c_3$ such that

Omitting the small terms ${\bf{Q}}_i^n$ in Eqs (3.4)–(3.7), replacing $\mathbb{U}_i^n$ with ${\bf{U}}_i^n$, we can get the difference scheme for solving the Keller-Segel models (1.1)–(1.3) when $n\geq 4$

Considering the ${\bf{F}}_x$ in the difference schemes (3.11)–(3.14) are nonlinear, we move them to the RHS of the equation, and use a nonlinear cycle iterative method to compute them. To obtain the fourth-order accuracy in the spatial direction, we employ the following fourth-order $Pad\acute{e}$ compact difference scheme ^[48]

in which, for the values of ${\bf{F}}_x$ at the boundary points, we employ the consistent fourth-order boundary scheme in ^[49], which is very stable for advection dominant problems, to compute them as follows

For the computation of the $v_x$ in ${\bf{F}}_x$, we also use the fourth-order $Pad\acute{e}$ scheme similar to Eq (3.4) and the fourth-order boundary difference schemes similar to Eqs (3.13) and (3.14) to compute them.

Remark 3.1. In the computation, the nonlinear cycle iterative process is shown as follows

● At the $n$th time step, using ${({\bf{F}}_x)_{i}^{n-1}}$ to approximate ${({\bf{F}}_x)_{i}^{n}}$, and then compute ${\bf{U}}_{i}^{n}$.

● Updating the values of ${({\bf{F}}_x)_{i}^{n}}$, then compute ${\bf{U}}_{i}^{n}$ again.

Repeat computations in turn till $\|{\bf{U}}_{i}^{n, (k)}-{\bf{U}}_{i}^{n, (k-1)}\|$ less than a small amount $\sigma$, set ${\bf{U}}_{i}^{n, (*)}\leftarrow {\bf{U}}_{i}^{n, (k)}$, $1\leq i\leq N-1{, 1\leq n \leq M}$, where $k$ is iteration number and $*$ is the approximate convergent solution. Then proceed the computation to the next time step.

Remark 3.2. Noticing that the schemes (3.11)–(3.14) requires to solve the system of nonlinear algebraic equations at each time step, and the difference equation is tri-diagonal at each time step, we can use the Thomas algorithm ^[50] to solve them. See the time advancement algorithm in Subsection 4.2 for details.

Remark 3.3. It is not difficult to find that Eqs (3.11)–(3.14) is a five-level fully implicit compact difference scheme for solving the system (1.4)–(1.6). Nonetheless, the new difference scheme not only has fourth-order accuracy in the spatial direction, but also has fourth-order accuracy in the temporal direction. And the compactness above only refers to the spatial direction with three grid points, however, it is not compact in the temporal direction since five time steps are involved.

Remark 3.4. Noticing the schemes (3.11)–(3.14) is five-level, when we use it to solve the system (1.4)–(1.6) step by step, we need to compute the first, second, and third start-up time steps. We give the computational method for them in the next part.

3.2.   Computation for the start-up time steps

We consider Eq (3.2) at the $(n-\frac{1}{2})$th time step, i.e.,

The Crank-Nicolson method is used to discretize the time derivative term on the LHS of Eq (3.7), i.e.,

We take the weighted average of the $n$th and $(n-1)$th time steps for other items. That is, ${\bf{\Theta}}_{i}^{n-\frac{1}{2}} = \frac{1}{2}(\bf{\Theta}_{i}^{n}+\bf{\Theta}_{i}^{n-1})+\tilde{O}_i^{n-\frac{1}{2}}(\tau^2)$, $1\leq i\leq N-1$, $1\leq n\leq 3$. (${\bf{\Theta}}$ could be $\mathbb{U}, {\bf{F}}, {\bf{R}}$ in Eq (3.7)). Similarly, combining the initial and boundary conditions (1.5), (1.6) and Lemma 2.1, we have

where ${\bf{P}}_i^{n-\frac{1}{2}} = O_{ti}^{n-\frac{1}{2}}({\tau ^2})+O_{xxi}^{n-\frac{1}{2}}({h^4})$, ${\bf{P}}_0^{n-\frac{1}{2}} = O_0^{n-\frac{1}{2}}(h^4)$, and ${\bf{P}}_N^{n-\frac{1}{2}} = O_N^{n-\frac{1}{2}}$. Then, there exist constants $c_4$, $c_5$ and $c_6$ such that

Omitting the small terms ${\bf{P}}_i^{n-\frac{1}{2}}$ in Eqs (3.19)–(3.22), replacing $\mathbb{U}_i^{n-\frac{1}{2}}$ with ${\bf{U}}_i^{n-\frac{1}{2}}$, we can get the start-up time steps difference scheme for solving the Keller-Segel models (1.1)–(1.3) as follows

Remark 3.5. For the computational strategy for chemotaxis term $({\bf{F}}_{x})_i^{n-\frac{1}{2}}$ in the difference schemes (3.26)–(3.29), similarly, we can employ Eqs (3.4)–(3.6) to compute them. At the same time, we employ Eq (3.4) and the boundary difference schemes (3.28) and (3.29) to compute the $v_x$ in ${\bf{F}}_{x}$.

Remark 3.6. Equations (3.26)–(3.29) is two-level implicit difference scheme. It has the second-order accuracy in the temporal direction and the fourth-order accuracy in the spatial direction.

Remark 3.7. According to the derivation process above, for the start-up time steps schemes (3.26)–(3.29), the Crank-Nicolson method is employed to discretize the temporal derivative term, and the fourth-order compact difference formula is used to discretize the spatial derivative terms, so the start-up time steps difference schemes (3.26)–(3.29) is unconditionally stable. For the fully implicit four-order compact difference schemes (3.11)–(3.14), the BDF-4 formula is employed to discretize the temporal derivative term, and the space derivative term is computed by using the fourth-order compact difference formulas and the fourth-order Padé formulas respectively. Therefore, according to ^[41,42,43], the fully implicit four-order compact difference schemes (3.11)–(3.14) is stiffly-stable.

3.3.   Richardson's extrapolation in time dimension

We can employ the difference schemes (3.26)–(3.29) to compute the start-up time steps (i.e., $n = 1, 2, 3$) from the initial value ${\bf{U}}^{0}$. Notice that the HOC difference schemes (3.11)–(3.14) have the fourth-order accuracy in the temporal direction, however, the difference schemes (3.26)–(3.29) have only second-order accuracy. To obtain the same accuracy in the temporal direction as the schemes (3.11)–(3.14), the accuracy of the start-up time steps is improved to the fourth-order by using the Richardson extrapolation technique as follows

where ${\bf{U}}^n(\tau, h)$ in the LHS of Eq (3.30) is the extrapolated numerical solution of the start-up time steps. ${\bf{U}}^n(\tau, h)$ in the RHS of Eq (3.30) is the numerical solution when time step length is taken as $\tau$. ${\bf{U}}^{2n}(\tau/2, h)$ is the numerical solution when time step length is taken as $\tau/2$.

Remark 3.8. We only use Eq (3.30) to extrapolate the start-up time steps (i.e., $n = 1, 2, 3$) step by step on the basis of the difference schemes (3.26)–(3.29), and the extrapolated solutions for the start-up time steps has the fourth-order accuracy in the temporal direction.

4.   Numerical algorithm

In order to obtain the biologically non-negative positivity-preserving of the cell density $u(x_i, t_n)$ and the concentration of the chemoattractant $v(x_i, t_n)$, we establish the positivity-preserving algorithm and the time advancement algorithm based on the developed schemes in this section.

4.1.   Positivity-preserving algorithm

The solution obtained by using the difference schemes (3.11)–(3.14) and (3.26)–(3.29) may become negative, especially in the area where the value changes greatly in the physical domain, the solution exhibits a large gradient. To force the positivity of the cell density $u(x_i, t_n)$ and chemoattractant concentration $v(x_i, t_n)$ at all time discrete steps, similar to the ^[33] and using Lemma 2.1, combining with the boundary discrete methods (3.6) and (3.7) and Eqs (3.21) and (3.22) in Section 3, we define the following solutions average

We use the following positivity-preserving limiter (PPL) ^[33,45,46] to filter the solution ${\bf{U}}_i^n, \; i = 0, 1, \cdots, N$, $0 < n\leq M$, at each time step and obtain the positivity-preserving solution $\widetilde{\bf{U}}_i^n\geq 0, \; i = 0, 1, \cdots, N$, $0 < n\leq M$, without losing the fourth-order accuracy.

where

Then, we establish the following positivity-preserving algorithm.

4.2.   Time advancement algorithm

Noticing that the implicit schemes (3.11)–(3.14) and Eqs (3.26)–(3.29) require to solve the system of nonlinear algebraic equations at each time step. Because the difference equation is tri-diagonal and only three grid points are used at each time step, we can solve them by using the Thomas algorithm ^[50]. And since the schemes (3.11)–(3.14) and Eqs (3.26)–(3.29) are nonlinear, linearized iterative strategy is used to solve them. The time advancement algorithm with linearized iteration process is described next page.

5.   Numerical experiments

In this section, we employ several examples to verify the various performances of the proposed method for solving the Keller-Segel chemotaxis models (1.1)–(1.3). First, we simulate the non-constant stationary solution of the Keller-Segel chemotaxis model in Subsection 5.1, and separately analyze the evolution of the solution over time and the asymptotic behavior of the solution under the strong chemotaxis sensitivity coefficient $\chi$. Second, we test the accuracy, convergence, stability and positivity-preserving of the proposed method for solving the Keller-Segel chemotaxis models (1.1)–(1.3) in Subsection 5.2.

We use homogeneous Neumann boundary conditions for all numerical examples. The computations are programmed in Fortran 90 language with double precision arithmetic and are implemented on a Intel(R)Core(TM) i7-8550U CPU@1.80 GHz 2.00 GHz, 8 GB RAM Windows workstation. The $L_\infty$, $L_2$ norm errors and convergence rate are defined as

where ${\bf{U}}(x_i, t_M)$ and ${\bf{U}}_{i}^M$ represent the exact and numerical solutions at point $(x_i, t_M)$, respectively. $L_\rho (h_1)$ and $L_\rho (h_2)$ denote the relevant $L_\rho$ norm errors when the grid sizes are $h_1$ and $h_2$, respectively, in which, $\rho$ could be $\infty$ or $2$.

Note: $k$ is iteration number and $*$ is convergent approximation solution.

5.1.   Asymptotic behavior of the solution

In this subsection, we employ the proposed method to solve two different types of Keller-Segel's minimal models with different initial values, and analyze the evolution of the solutions over time and the asymptotic behavior under strong chemotaxis sensitivity coefficients $\chi$.

5.1.1.   Keller-Segel's minimal model Ⅰ

First, we consider the following Keller-Segel' minimal chemotaxis model^[8,23]

where $\Omega = [0, \pi]$, and the initial conditions as follows

This system has a constant solution $(\bar{u}, \bar{v})$, and the initial values of the system are a small disturbance of the $(\bar{u}, \bar{v})$. When $\chi > \pi$, the constant solution is unstable. According to ^[23], this minimal model has a monotonically decreasing non-constant stationary solution. In particular, the asymptotic expression for the steady problem of the minimal model is given in ^[8] as follows

where

and $\varphi_1(y): = \tanh y$, $\varphi_2(y): = 1-y\tanh y$, $O(1)$ denotes a generic quantity bounded by a constant depending only on $l = \pi$.

Case 1. The evolution process of the system (5.1)–(5.3) over time. Figure 2 shows the formation process of the monotonic edge peak solution of the minimal model in $[0, \pi]$, with $N = 27, \tau = 0.1h$, $\chi = 5$. From left to right, they represent the initial values, intermediate changes and stationary solution of $(u, v)$, respectively. Figure 3 shows the space-time distributions of the stationary solution formed by the evolution of $u$ and $v$ over time, respectively. From these figures, it is not difficult to find that $(u, v)$ stays near the above initial value for a long period of time, and starts to approach the left end of the $x-$direction rapidly at $t\approx10$, and a sharp peak solution is formed at $t\approx18$ in the left boundary of $x-$direction, and finally a stable edge peak solution is formed at $t\approx20$.

Case 2. The asymptotic behavior of the solution for the system (5.1)–(5.3) under the strong chemotaxis coefficient $\chi$. Figures 4 and 5 show the $\delta-$function structure formed when the chemotaxis coefficient $\chi$ is large. In each subfigure, the dotted dashed line (blue in the color map) represents the stationary solution of the minimal models (5.1)–(5.3), and the solid line (purple in the color map) represents the asymptotic expansion (5.4) and (5.5) of the stationary solution in ^[8]. As can be seen from these subfigures, when the chemotaxis coefficient $\chi$ increases, the distribution of the cell density $u(x, t)$ tends to a $\delta-$function, and the distribution of the chemical concentration $v(x, t)$ tends to a Green's function. And these monotonic stationary solutions are stable, which is in good agreement with the asymptotic expansion (5.4) and (5.5).

5.1.2.   Keller-Segel's minimal model Ⅰ

Next, we consider another Keller-Segel's minimal chemotaxis model^[9]

where $\Omega = [0, 1]$, and the initial conditions as follows

In the computation below, we take the parameters ${d} = 0.1$, $\chi = 5$, $N = 50$ and $\tau = 0.5\; h$. In Figure 6, numerical simulations of the minimal model is plotted on a fixed domain $[0, 1]$ with homogeneous Neumann boundary, at different moments $t = 0, 0.5, 1, 1.5, 2, 10$. The cell density $u(x, t)$ and chemical concentration $v(x, t)$ are plotted at distinct times, showing the growth of the minimal models (5.6)–(5.8) solution as cells accumulate into a sharp boundary peak. Figure 7 shows the stationary solution and space-time distribution of minimal models (5.6)–(5.8) evolving over time.

5.2.   Performance testing

In this subsection, we employ three different types of examples to test the performances of the developed difference schemes in this paper for solving the Keller-Segel chemotaxis models (1.1)–(1.3), which include the accuracy, stability and positivity-preserving.

5.2.1.   Accuracy and stability test without PPL

First, we study the accuracy and stability of the difference scheme without PPL for solving the system (1.1)–(1.3). To this end, we consider the Keller-Segel chemotaxis model with two additional fluxes $r_{1}(x, t)$ and ${r_2}(x, t)$ on the RHS of system (1.1)–(1.3) as follows

with homogeneous Neumann boundary conditions, where $\chi = d = 1$ and the fluxes are given by

The initial conditions $u(x, 0)$ and $v(x, 0)$ are given by the exact solution

In this computation, we take the time step length $\tau = 0.5\; h$, and $T = 1$. Figure 8 shows that the $L_{\infty}$, $L_{2}$ norm errors and convergence rates of the cell density $u(x, t)$ and the chemoattractant concentration $v(x, t)$, respectively. Where, we take space step length $h = \pi/512, \pi/256, \pi/128, \pi/64, \pi/32$ from left to right. It is not difficult to find from Figure 8 that in the two cases, when we take a large time step length $\tau = O(h)$, the difference scheme can still converge with the fourth-order. In Table 1, we set grid ratio $\lambda = \tau/h$, and take numbers of grids $N = 128,256,512$, to test the stability of the difference schemes when the grid ratio $\lambda$ gradually increases, respectively. According to the Table 1, we find the difference schemes have better stability. Through this numerical experiment, the theoretical results of this paper are well verified.

5.2.2.   Accuracy and stability test with PPL

Next, we test the accuracy and stability of the proposed method with positivity or bound-preserving limiter (BPL) for solving Keller-Segel chemotaxis models (1.1)–(1.3). Meanwhile, we will compare the computed results with those without PPL. To this end, we continue to consider the Keller-Segel chemotaxis models (5.9) and (5.10) with additional fluxes in Subsection 5.2.1. We take $\chi = d = 1$, but the fluxes are

The initial conditions $u(x, 0)$ and $v(x, 0)$ are given by the exact solution $u(x, t) = v(x, t) = e^{-t}\cos(x).$

In Table 2, we list the $L_{\infty}$, $L_{2}$ norm errors and convergence rate of the cell density $u(x, t)$ with/without BPL when $\tau = 0.1h, T = 0.2$, for Problem 5.2.2, and compare the computed results with those in ^[33]. It can be observed from this table that the results obtained by the proposed method can not only converge with the fourth-order accuracy before adopting the BPL, but also the fourth-order accuracy is maintained after adopting the BPL. However, the computed result in ^[33] has only third-order accuracy, and the errors obtained by the proposed method is 2-3 orders of magnitude lower than the errors obtained in ^[33]. Tables 3 and 4 show that the $L_{\infty}$, $L_{2}$ norm errors and convergence rate of the cell density $u(x, t)$ and chemoattractant concentration $v(x, t)$ with/without BPL when $\tau = 0.5h, T = 1$ for Problem 5.2.2, respectively. From these two tables, it can still be found that when we take a large time step $\tau = O(h)$, the results obtained before and after using the proposed method with the BPL can converge well at the fourth-order rate, which is consistent with the theoretical analysis. In Table 5, we run the code to terminal time $T = 10$, and test the accuracy and stability of the proposed method for solving the example before and after the BPL is applied under the large time. In Table 6, we continue to set the grid ratio $\lambda$ gradually increases, take the numbers of spatial grids to be $N = 128,256,512$, and test the stability of the proposed method before and after the application of the BPL. According to the results in Tables 5 and 6, it is not difficult to find that the proposed method has better stability. At the same time, it can still maintain its better stability and accuracy after using the BPL.

5.2.3.   Mass conservation and non-negativity of cell density with PPL

In this example, we test the mass conservation and the capability of positivity-preserving by using proposed HOC difference schemes for the cell density and chemoattractant concentration problem. We choose the following initial conditions for the Keller-Segel chemotaxis models (1.1)–(1.3)

where $\chi = d = 1$. Noticing that the solution of the system does not blow up under 1D setting, but the cell density $u(x, t)$ will appear high concentration aggregation phenomenon in the central area of the domain $[-0.5, 0.5]$. Next we take $\tau = 10^{-4}h$ to solve the Keller-Segel chemotaxis model by using the proposed method with/without PPL applied, respectively.

Figure 9 shows the comparison of cell density $u(x, t)$ and chemoattractant concentration $v(x, t)$ with $T = 6\times10^{-5}$, $N = 50$ in the whole domain $[-0.5, 0.5]$ and near the left and right boundary points before and after applying the PPL by using the proposed method. From these two figures, it can be seen that before applying the PPL, the cell density $u(x, t)$ appears high concentration aggregation and strong oscillations around the center $x = 0$ of the domain $[-0.5, 0.5]$. Meanwhile, it can be seen that there are many small negative values around $x = 0$ and the left and right boundaries $x = -0.5$ and $0.5$. These negative values makes the cell density $u(x, t)$ lose its biological significance. However, the chemoattractant concentration $v(x, t)$ dose not appear negative in the whole domain $[-0.5, 0.5]$ and near the left and right boundary points, and the concentration curve is smooth and non-oscillating. After applying the PPL, it can be observed that all the negative values of cell density $u(x, t)$ have been removed, and the oscillation around $x = 0$ is effectively reduced. The chemoattractant concentration $v(x, t)$ is highly consistent with that before applying the PPL. Figure 10 shows the three-dimensional space-time density distribution of cell density $u(x, t)$ and the projections on the $xou$ plane and the $xot$ plane when $T = 6\times10^{-5}$ and $N = 50$, before and after applying the PPL, respectively. From these two figures, it is obvious that the PPL maintains the non-negativity of the cell density $u(x, t)$.

Figure 11 shows the comparison of the cell density curve and its projection on the $xou$ plane when $T = 1\times10^{-4}$, $N = 50$ and $N = 100$ by using the proposed method before and after the application of the PPL, respectively. According to these two figures and combined with Figure 9, it can be observed that under different spatial grid numbers $N$ and terminal time $T$, the cell density $u(x, t)$ has high concentration aggregation and strong oscillation around $x = 0$, as well as varying degrees of negative values. The proposed method by applying the PPL can remove the corresponding negative values and reduce its oscillation, so that the density curve becomes smoother. In addition, we can also observe that when the same spatial grid number $N$ is used, as the terminal time $T$ gradually increases, the values of the cell density $u(x, t)$ around $x = 0$ increases, and the oscillation range becomes larger. When the terminal time $T$ remains unchanged, with the increase of the spatial grid number $N$, the cell density $u(x, t)$ also shows an increasing trend around $x = 0$, but the oscillation range shrinks. This shows that the increase of the spatial grid number $N$ and the different terminal time $T$ appear different degrees of cell aggregation phenomenon for the Keller-Segel chemotaxis model.

Table 7 shows that the discrete mass of the proposed scheme at different time $T$ with $\tau = 2\times10^{-6}, h = 2\times10^{-2}$ for Problem 5.2.3. According to this table, we can observe that the discrete mass of the proposed scheme is conserved without PPL. It is still conserved by using PPL when the function appears smooth and less oscillation, but the PPL will affect the mass conservation if the oscillation is large.

6.   Conclusions

In this paper, an HOC difference scheme with positivity-preserving for solving the Keller-Segel chemotaxis model is proposed, and the computational strategy for the nonlinear chemotaxis term in the difference scheme is provided. The truncation error is $O(\tau^4+h^4)$, that is, it has fourth-order accuracy in both spatial and temporal directions. A positivity-preserving numerical algorithm that ensures the non-negativity of the cell density at all time and maintains the fourth-order without accuracy loss, and a time advancement algorithm is established. Finally, through several numerical experiments, the non-constant stationary solutions of the Keller-Segel chemotaxis model are simulated, and the accuracy, stability and positivity-preserving of the present difference scheme are validated.

The proposed high accurate numerical method can be directly extended to solve the two-dimensional Keller-Segel chemotaxis model. However, because the chemotaxis model is a nonlinear coupled system in two-dimensional cases, it may encounter a higher computational time cost. At the same time, the solution of the chemotaxis model maybe blow up in finite time^[15], so capturing the blow-up solution will be a challenging task ^[16]. We are trying to solve the two-dimensional Keller-Segel chemotaxis model by using adaptive moving grid and multi-grid acceleration techniques. This work is in progress, and the research results will be reported in the near future.

Acknowledgments

We would like to thank the editors and anonymous reviewers for their constructive comments and suggestions that are helpful to improve the quality of the paper. This work is partially supported by National Natural Science Foundation of China (12161067, 11772165, 11961054, 11902170), National Natural Science Foundation of Ningxia (2022AAC02023), the Key Research and Development Program of Ningxia (2018BEE03007), National Youth Top-notch Talent Support Program of Ningxia, and the First Class Discipline Construction Project in Ningxia Universities: Mathematics.

Conflict of interest

The authors declare there is no conflict of interest.