Loading [MathJax]/jax/element/mml/optable/BasicLatin.js
Research article Special Issues

Identification of drug side effects with a path-based method


  • Received: 10 January 2022 Revised: 20 February 2022 Accepted: 07 March 2022 Published: 06 April 2022
  • The study of drug side effects is a significant task in drug discovery. Candidate drugs with unaccepted side effects must be eliminated to prevent risks for both patients and pharmaceutical companies. Thus, all side effects for any candidate drug should be determined. However, this task, which is carried out through traditional experiments, is time-consuming and expensive. Building computational methods has been increasingly used for the identification of drug side effects. In the present study, a new path-based method was proposed to determine drug side effects. A heterogeneous network was built to perform such method, which defined drugs and side effects as nodes. For any drug and side effect, the proposed path-based method determined all paths with limited length that connects them and further evaluated the association between them based on these paths. The strong association indicates that the drug has a side effect with a high probability. By using two types of jackknife test, the method yielded good performance and was superior to some other network-based methods. Furthermore, the effects of one parameter in the method and heterogeneous network was analyzed.

    Citation: Meng Jiang, Bo Zhou, Lei Chen. Identification of drug side effects with a path-based method[J]. Mathematical Biosciences and Engineering, 2022, 19(6): 5754-5771. doi: 10.3934/mbe.2022269

    Related Papers:

    [1] Lin Zhang, Yongbin Ge, Zhi Wang . Positivity-preserving high-order compact difference method for the Keller-Segel chemotaxis model. Mathematical Biosciences and Engineering, 2022, 19(7): 6764-6794. doi: 10.3934/mbe.2022319
    [2] Sunwoo Hwang, Seongwon Lee, Hyung Ju Hwang . Neural network approach to data-driven estimation of chemotactic sensitivity in the Keller-Segel model. Mathematical Biosciences and Engineering, 2021, 18(6): 8524-8534. doi: 10.3934/mbe.2021421
    [3] Thierry Colin, Marie-Christine Durrieu, Julie Joie, Yifeng Lei, Youcef Mammeri, Clair Poignard, Olivier Saut . Modeling of the migration of endothelial cells on bioactive micropatterned polymers. Mathematical Biosciences and Engineering, 2013, 10(4): 997-1015. doi: 10.3934/mbe.2013.10.997
    [4] Shangbing Ai, Zhian Wang . Traveling bands for the Keller-Segel model with population growth. Mathematical Biosciences and Engineering, 2015, 12(4): 717-737. doi: 10.3934/mbe.2015.12.717
    [5] Tong Li, Zhi-An Wang . Traveling wave solutions of a singular Keller-Segel system with logistic source. Mathematical Biosciences and Engineering, 2022, 19(8): 8107-8131. doi: 10.3934/mbe.2022379
    [6] Maghnia Hamou Maamar, Matthias Ehrhardt, Louiza Tabharit . A nonstandard finite difference scheme for a time-fractional model of Zika virus transmission. Mathematical Biosciences and Engineering, 2024, 21(1): 924-962. doi: 10.3934/mbe.2024039
    [7] Qigang Deng, Fugeng Zeng, Dongxiu Wang . Global existence and blow up of solutions for a class of coupled parabolic systems with logarithmic nonlinearity. Mathematical Biosciences and Engineering, 2022, 19(8): 8580-8600. doi: 10.3934/mbe.2022398
    [8] Z. Jackiewicz, B. Zubik-Kowal, B. Basse . Finite-difference and pseudo-spectral methods for the numerical simulations of in vitro human tumor cell population kinetics. Mathematical Biosciences and Engineering, 2009, 6(3): 561-572. doi: 10.3934/mbe.2009.6.561
    [9] Benjamin Wacker, Jan Christian Schlüter . A non-standard finite-difference-method for a non-autonomous epidemiological model: analysis, parameter identification and applications. Mathematical Biosciences and Engineering, 2023, 20(7): 12923-12954. doi: 10.3934/mbe.2023577
    [10] Fadoua El Moustaid, Amina Eladdadi, Lafras Uys . Modeling bacterial attachment to surfaces as an early stage of biofilm development. Mathematical Biosciences and Engineering, 2013, 10(3): 821-842. doi: 10.3934/mbe.2013.10.821
  • The study of drug side effects is a significant task in drug discovery. Candidate drugs with unaccepted side effects must be eliminated to prevent risks for both patients and pharmaceutical companies. Thus, all side effects for any candidate drug should be determined. However, this task, which is carried out through traditional experiments, is time-consuming and expensive. Building computational methods has been increasingly used for the identification of drug side effects. In the present study, a new path-based method was proposed to determine drug side effects. A heterogeneous network was built to perform such method, which defined drugs and side effects as nodes. For any drug and side effect, the proposed path-based method determined all paths with limited length that connects them and further evaluated the association between them based on these paths. The strong association indicates that the drug has a side effect with a high probability. By using two types of jackknife test, the method yielded good performance and was superior to some other network-based methods. Furthermore, the effects of one parameter in the method and heterogeneous network was analyzed.



    In this paper, the finite-difference method (FDM) is applied to solve the following two-dimensional (2D) Keller-Segel model [1,2]:

    {ut+(χuv)=dΔu,vt=Δvv+u,(x,y,t)Ω×(0,T]. (1.1)

    The model (1.1) describes the evolutionary process of cell density u(x,y,t) and a chemical stimulus (chemoattractant) concentration v(x,y,t) over a time t and location (x,y), where Ω={(x,y)|ax,yb}R2 is a convex bounded domain and a and b are constants. and Δ are gradient and Laplacian operators, respectively. χ>0 represents the chemotactic sensitivity constant; d>0 denotes the diffusion rate of cells. In addition, the initial conditions related with (1.1) are given as

    u(x,y,0)=u0(x,y),v(x,y,0)=v0(x,y),(x,y)Ω, (1.2)

    and the boundary conditions are assumed to be homogeneous Neumann boundary conditions, that is,

    un=vn=0,(x,y,t)Ω×(0,T], (1.3)

    where Ω is the boundary of Ω and n is the outward normal of Ω. With this condition (1.3), the total mass

    Massu=Ωu(x,y,t)dxdy=Ωu(x,y,0)dxdy

    is conserved as temporal evolution. Besides the above, the model (1.1) has the following form of free energy:

    E(u,v)(t)=Ω(uln(u)+χ2|v|2+χ2v2χuv)dxdy,t>0. (1.4)

    From mathematical analysis, the following equation can be verified by direct calculation of Eq. (1.4),

    E(u,v)t=Ω(u|(ln(u)χv)|2+χ(vt)2)dxdy0,t>0.

    Many early works show that the free energy (1.4) is decreasing over time and is mainly employed to demonstrate the existence of solutions for the chemotaxis system; see [3,4] and the references therein.

    Chemotaxis refers to the directional movement, which includes toward or away from the higher concentrations of cells or microorganisms that are stimulated by chemical substances in the external environment along the gradient directions of the concentration for stimuli. Chemotaxis phenomena play a crucial role in numerous intricate biological evolutions, such as bacterial aggregation, angiogenesis, pattern formation, embryonic development and so on[5]. From as early the 1950 to 1980, the chemotaxis phenomena have been extensively studied by many applied mathematicians and biologists, and a class of models of partial differential systems closely associated with taxis have been proposed; see [1,2,6,7,8,9,10]. Among them the most classical is the above-named system (1.1) (first proposed by Keller and Segel [1,2] in 1970), which simulated the aggregation phenomenon for Amoebae and Dictyostelium, as well as the traveling wave migration phenomenon for Escherichia coli in a capillary filled with nutrients.

    Since the model was proposed, many researchers have systematically analyzed the properties of its solution, including its global existence, asymptotic profile [11], global boundedness [12,13], finite-time blow-up [14,15,16], etc. Particularly, if the initial mass Ωu(x,y)dxdy of the cells in the 2D case satisfies a critical threshold value, its solution will blow up in finite time [7,8,9,14,15,16]. This blow-up denotes a mathematical concept of the bacterial aggregation arising in real biotic environments [7,17,18]. However, it is arduous if we want to obtain the analytical solutions of the model due to the strong nonlinear characteristics of the model itself. Meanwhile, it is still difficult to better numerically capture the blow-up or spike solutions. Therefore, it is desperately needed to establish a high accuracy and more stable numerical method to investigate the properties of the solutions for the chemotaxis system given by (1.1)–(1.3). At the same time, the numerical investigation of the chemotaxis system also facilitates the theoretical exploration of its dynamic behavior.

    In recent years, some numerical methods have been involved via the investigation of the chemotaxis systems [3,4,16,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35]. For example, Saito and Suzuki [19] proposed a conservative numerical scheme by using the FDM to solve a parabolic-elliptic coupling chemotaxis system. Meanwhile, in order to obtain a positivity-preserving scheme under total mass that is conservative, Saito [20] proposed a conservative scheme by using the FDM for the system given by (1.1)–(1.3) with nonlinear diffusion. Xiao et al. [21] derived a semi-implicit scheme by using a characteristic finite element method (FEM) to simulate the blow-up solutions of the chemotaxis system on surfaces, as well as pattern formulation and aggregation phenomenon of bacteria. And, their method has second-order accuracy in L2-norm and H1-norm errors. Epshteyn and Kurganov [22] introduced chemotaxis concentration gradient variables to rewrite the original Keller-Segel model given by (1.1)–(1.3), and they designed an internal penalty discontinuous Galerkin (DG) method for the rewritten system to border on the solution of the original system. Li et al. [23] used a local DG (i.e., LDG) method to modify that in [22], and they obtained the optimal convergence rates based on a specific finite element space before blow-up for the chemotaxis system. They also deduced a positivity-preserving P1 LDG scheme and proved that it is L1-stable. The LDG method was also used to solve the system given by (1.1)–(1.3) in [3], and the energy dissipation with the LDG discretization was proved. Sulman and Nguyen [24] proposed an adaptive moving mesh method by applying an implicit-explicit FEM to solve the system given by (1.1)–(1.3). The method has second-order accuracy in the spatiotemporal directions, and it is noteworthy that the obtained solutions are positive at all time steps if the initial values of the system (1.1) are positive. Qiu et al. [25] proposed a new scheme by using a interface-corrected direct DG method to solve this model (1.1), and their method satisfies the positivity-preserving requiremet without losing third-order accuracy. Based on the gradient flow structure, Shen and Xu [4] proposed a class of numerical schemes for solving the Keller-Segel model given by (1.1)–(1.3). Among them, the first-order accuracy scheme satisfies the mass conservation, bounded positivity, unique solvability and energy dissipation of the original differential equation, while the second-order accuracy scheme satisfies the first three properties. Chen et al. [16] analyzed the error of the numerical scheme proposed in [4], and they deduced the finite-time blow-up of non-radial numerical solutions under certain assumptions. Based on the generalized smoothed particle hydrodynamics meshless method, Dehghan and Abbaszadeh [26] established a second-order-accuracy numerical scheme for some chemotaxis models with the blow-up phenomena during tumor growth, and they numerically simulated the blow-up problems of the original chemotaxis model. Filbet [27] employed the finite volume method (FVM) to approximate the solution of the system given by (1.1)–(1.3), and they simulated the blow-up problems. Chertock and Kurganov [28] developed a center-upward scheme with second-order accuracy by using the FVM based on their previous work in [29] and the blow-up problems of the system (1.1)–(1.3), the pattern formation of bacteria were also simulated. Epshteyn [30] proposed an upwind-difference potentials scheme to solve the problems in complex geometries for the system given by (1.1)–(1.3). In addition, some other numerical methods [31,32,33,34,35], such as the fractional step (or operator splitting) method [31,32], hybrid finite-volume-finite-difference methods [33], generalized FDM [35], etc., were also used to solve the chemotaxis system given by (1.1)–(1.3).

    Although some numerical methods mentioned above for solving the system given by (1.1)–(1.3) can achieve high accuracy in the spatial direction, such as those in [25,30] (third-order accuracy) and [3,33] (fourth-order accuracy), most numerical methods have low accuracy, especially in the temporal direction; see [4,16,19,20,21,22,23,24,26,27,28,29,31,32,34,35]. Meanwhile, many high-accuracy numerical methods are non-compact in space, and the stability conditions are relatively harsh. In other words, if we want to employ these methods to solve real problems, we must take a small time step length to satisfy their stability conditions, which will expend expensive computational time. However, the high-order compact (HOC) difference scheme has attracted many researchers because of its strong advantages, such as fewer computational nodes, small computational errors, better numerical stability and non-complicated boundaries. Meanwhile, the backward differentiation formula with fourth-order accuracy (BDF-4), which is an A-stabilized method and appears first in [36], has been verified to be a feasible method to obtain high accuracy, and it has a relatively large stability condition range [37,38]. In addition, the major difficulty in solving the system comes from the nonlinearities, such as the chemotaxis term (χuv). One issue is that the coupling form will increase the discrete difficulty when we want to obtain the high accuracy and satisfy mass conservative schemes. Another difficulty is positivity preservation since u and v in the system (1.1)-(1.3) have real biological significance in complex chemotaxis phenomenon and cannot have negative values. To acquire the high-accuracy, positivity-preserving numerical solutions for the system given by (1.1)–(1.3), in this work, our purpose was to derive a compact difference scheme to approximate the solutions of the original chemotaxis system, and the scheme has space-time fourth-order accuracy and is stable, as well as positivity-preserving.

    The remainder of this paper is organized as follows. Some preliminary preparations with basic symbols, definitions and theorems are provided in Section 2. In Section 3, we deduce an HOC scheme for the system given by (1.1)–(1.3), and give the computational strategies for the initial time steps and the nonlinear terms. In Section 4, a time advancement algorithm combined with a multigrid method and a positivity-preserving algorithm are proposed. In Section 5, some numerical examples are employed to verify the accuracy, stability, positivity-preserving property, mass conservation and energy dissipation. The finitetime blow-up problems for the chemotaxis system given by (1.1)–(1.3) are simulated by using the proposed method. Finally the conclusion is provided in the end.

    First, the domain {(x,y,t)|ax,yb,0tT} is divided by uniform meshes N2×M, M,NZ+. Denote h=(ba)/N to represent the spatial step size, and τ=T/M stands for the temporal step length. Let xi=a+ih,yj=a+jh and tn=nτ, 0i,jN, 0nM, and mark the mesh points (xi,yj,tn). Figure 1 shows the 2D spatial mesh point stencil.

    Figure 1.  Spatial mesh point stencil.

    Second, we define Ωh={(xi,yj)|0i,jN}; its discrete boundary Ωh={(0,j),(N,j)|0jN}{(i,0),(i,N)|1iN1}; let Ωτ={tn|0nM} and Ωhτ=Ωh×Ωτ. For any mesh function wWhτ={wni,j|0i,jN,0nM} defined on Ωhτ, we have

    wn12i,j=12(wni,j+wn1i,j),δtwn12i,j=1τ(wni,jwn1i,j),δxwni,j=12h(wni+1,jwni1,j),δywni,j=12h(wni,j+1wni,j1),δ2xwni,j=1h2(wni+1,j2wni,j+wni1,j),δ2ywni,j=1h2(wni,j+12wni,j+wni,j1),Δtwni,j=112τ(25wni,j48wn1i,j+36wn2i,j16wn3i,j+3wn4i,j).

    Next, for simplicity, we let pwςp(x,y,t):≜wςp(x,y,t), pZ+, where ς represents x, y or t. And, denote A=[A1,A2,A3,A4,A5]=[25,48,36,16,3], B=[B1,B2,B3]=[1,10,1], S=[S1,S2,S3,S4]=[1,34,13,116]. To obtain a higher accuracy on the boundary described by Eq.(1.3), we employ the Taylor expansion with the Peano remainder to cope with Eq. (1.3). Thus, the following theorem holds.

    Theorem 2.1. Denote Wh={w(xi,yj)|0i,jN}; for mapping w:ΩhWh, we have the following:

    (1) If w(x,y)C5,0([x0,x4]×[y0,yN]), then wx(x0,yj)=112h5k=1Akw(xk1,yj)+O0,j(h4);

    (2) If w(x,y)C5,0([xN4,xN]×[y0,yN]), then wx(xN,yj)=112hN1k=N5ANkw(xk+1,yj)+ON,j(h4);

    (3) If w(x,y)C0,5([x0,xN]×[y0,y4]), then wy(xi,y0)=112h5k=1Akw(xi,yk1)+Oi,0(h4);

    (4) If w(x,y)C0,5([x0,xN]×[yN4,yN]), then wy(xi,yN)=112hN1k=N5ANkw(xi,yk+1)+Oi,N(h4), where the local truncation errors are

    O0,j(h4)=h464k=110k5Skwx5(x0+ksh,yj)(1s)4ds,0jN,ON,j(h4)=h464k=110k5Skwx5(xNksh,yj)(1s)4ds,0jN,Oi,0(h4)=h464k=110k5Skwy5(xi,y0+ksh)(1s)4ds,0iN,Oi,N(h4)=h464k=110k5Skwy5(xi,yNksh)(1s)4ds,0iN.

    Proof. First, we prove (1). We assume that w(x,y)Ck+1,0([xi1,xi+1]×[y0,yN]), and, according to the Taylor expansion with the Peano remainder, we have

    w(xi±h,yj)=kl=0(1)lhll!wxl(xi,yj)+hk+1k!10wxk+1(xi±sh,yj)(1s)kds, (2.1)

    where 1iN1 and 0jN. We suppose that w(x,y)C5,0([x0,x4]×[y0,yN]) and expand w(x1,yj), w(x2,yj), w(x3,yj) and w(x4,yj) at (x0,yj) by using Eq. (2.1) above, that is, we respectively take k=1,2,3,4, and have

    w(xk,yj)=4l=0(kh)ll!wxl(x0,yj)+(kh)52410wx5(x0+ksh,yj)(1s)4ds. (2.2)

    Then, we perform the following the operation: Eq. (2.2)k=1+α× Eq. (2.2)k=2+β× Eq. (2.2)k=3+γ× Eq. (2.2)k=4 for Eq. (2.2), and denote E=[E1,E2,E3,E4]=[1,α,β,γ], where α,β and γ, are constants; then, we can obtain

    4k=1Ekw(xk,yj)=5m=14k=1(kh)m1(m1)!Ekwxm1(x0,yj)+h5244k=1k5Ek10wx5(x0+ksh,yj)(1s)4ds. (2.3)

    Suppose that the coefficients of wx2(x0,yj), wx3(x0,yj) and wx4(x0,yj) in Eq. (2.3) equal to 0; we obtain

    {4k=1k2Ek=0,4k=1k3Ek=0,4k=1k4Ek=0.{α=34,β=13,γ=116.

    Substituting them into Eq. (2.3), denote A=[25,48,36,16,3] and S=[1,34,13,116], and we have

    wx(x0,yj)=112h5k=1Akw(xk1,yj)+O0,j(h4),0jN,

    where O0,j(h4)=h464k=110k5Skwx5(x0+ksh,yj)(1s)4ds. Similarly, (2) and (3) are also easily obtained. The proof is complete.

    According to Theorem 2.1, and by considering wx(x0,yj)=0, wx(xN,yj)=0, wy(xi,y0)=0 and wy(xi,yN)=0, we can easily derive the following fourth-order-accuracy boundary approximation formulas:

    w(x0,yj)125[48w(x1,yj)36w(x2,yj)+16w(x3,yj)3w(x4,yj)],w(xN,yj)125[48w(xN1,yj)36w(xN2,yj)+16w(xN3,yj)3w(xN4,yj)],w(xi,y0)125[48w(xi,y1)36w(xi,y2)+16w(xi,y3)3w(xi,y4)],w(xi,yN)125[48w(xi,yN1)36w(xi,yN2)+16w(xi,yN3)3w(xi,yN4)].

    In addition, we suppose that w(x,y,t)C6,6,5([xi1,xi+1]×[yj1,yj+1]×[tn,tn1]) for 1i,jN1, 1nM, and we denote I=[1,24,81,64] to derive the following truncation errors based on the Taylor expansion with the Peano remainder, that is,

    (Oxx)ni,j(h4)=h424102k=1[13(1s)315(1s)5]wx6(xi+(1)k1sh,yj,tn)ds,(Ot)ni,j(τ4)=τ4610(1μ)44k=1Ikwt5(xi,yj,tnkμτ)dμ,(Ox)ni,j(h4)=h412102k=1[13(1s)314(1s)4]wx5(xi+(1)k1sh,yj,tn)ds,(Ot)n12i,j(τ2)=τ21610(1μ)22k=1wt3(xi,yj,tn12+(1)k1μτ2)dμ,˜On12i,j(τ2)=τ2410(1μ)2k=1wt2(xi,yj,tn12+(1)k1μτ2)dμ.

    Similarly, we can easily derive the expressions of (Oyy)ni,j(h4), (Oy)ni,j(h4), (Ox)n1i,j(h4), (Oy)n1i,j(h4), (Oxx)n12i,j(h4) and (Oyy)n12i,j(h4). To facilitate mathematical analysis below, the following definitions are given.

    Definition 2.2. For any mesh function {wi,j|0i,jN}, the average operators A and B are defined as

    Awi,j={1254k=1Ak+1wk,j,i=0,0jN,1123k=1Bkwi+k2,j,1i,jN1,125N1k=N4ANk+1wk,j,i=N,0jN,Bwi,j={1254k=1Ak+1wi,k,j=0,0iN,1123k=1Bkwi,j+k2,1i,jN1,125N1k=N4ANk+1wi,k,j=N,0iN.

    Definition 2.3. The maximum norm and 2-norm errors are defined as

    ||||=max

    where w(x_i, y_j, t_M) and w_{{i, j}}^M stand for the exact and numerical solutions at the discrete mesh point (x_i, y_j, t_M) , respectively.

    Definition 2.4. Denote w_{\max} (t):\triangleq\max\limits_{(x, y)}w(x, y, t) ; the relatively maximum norm and 2 -norm errors in the temporal dimension are defined as

    \begin{align} &Rel_{{\infty}} = \frac{\|w_{\max} (t)-w_{\max}^* (t)||_\infty}{||w_{\max}^*\|_\infty} = \frac{\max\limits_{{0\leqslant n\leqslant M}}|w_{\max} (t_n)-w_{\max} ^{*}(t_n)|}{\max\limits_{{0\leqslant n\leqslant M}}|w_{\max} ^{*}(t_n)|}, \\ &Rel_{{2}} = \frac{\|w_{\max} (t)-w_{\max}^* (t)||_2}{||w_{\max}^*\|_2} = \frac{\sqrt{\tau\sum\limits_{n = 0}^{M}[w_{\max} (t_n)-w_{\max} ^{*}(t_n)]^2}}{\sqrt{\tau\sum\limits_{n = 0}^{M}[w_{\max} ^{*}(t_n)]^2}}, \end{align}

    where w_{\max}^* (t) represents the reference solution.

    Definition 2.5. The convergence rate is defined as

    \begin{align} &\text{Rate} = \frac{\log[L_\nu (h_1)/L_\nu (h_2)]}{\log(h_1/h_2)}, \end{align}

    where \nu stands for \infty or 2 and L_\nu (h_1) and L_\nu (h_2) are corresponding ||\cdot||_\nu norm errors which are closely related to h_1 and h_2 , respectively.

    In this part, an HOC scheme is derived to border on the solution of the Keller-Segel system given by (1.1)–(1.3). To facilitate the numerical analysis, we rewrite an equivalent form for the system given by (1.1)–(1.3), that is,

    \left\{ \begin{array}{l} \textbf{U}_{t}+\textbf{F}_{x}+\textbf{G}_{y} = \textbf{D}\Delta\textbf{U}+\textbf{R}, & a < x, y < b, 0 < t\leqslant T, & (3.1)\\ \textbf{U}(x, y, 0) = \textbf{U}_0(x, y), & a\leqslant x, y\leqslant b, & (3.2)\\ \textbf{U}_x(a, y, t) = \textbf{U}_x(b, y, t) = 0, & 0 < t\leqslant T, & (3.3)\\ \textbf{U}_y(x, a, t) = \textbf{U}_y(x, b, t) = 0, & 0 < t\leqslant T, & (3.4) \end{array} \right.

    where

    \textbf{U}:\triangleq\begin{bmatrix}u \\ v\end{bmatrix}, \quad \textbf{F}:\triangleq\begin{bmatrix}\chi uv_x \\ 0 \end{bmatrix}, \quad \textbf{G}:\triangleq\begin{bmatrix}\chi uv_y \\ 0 \end{bmatrix}, \quad \textbf{D}:\triangleq\begin{bmatrix}d\; \; \; \; \; \; \; \; \\\; \; \; \; \; \; \; 1 \end{bmatrix}, \quad \textbf{R}:\triangleq\begin{bmatrix}0\\-v+u \end{bmatrix}, \quad \textbf{U}_0:\triangleq\begin{bmatrix}u_0 \\ v_0\end{bmatrix}.

    Equation (3.1) is a nonlinear advection-diffusion-reaction equation. Then, we define

    \begin{align} &\mathbb{U} = \{[u(x_i, y_j, t_n), v(x_i, y_j, t_n)]^\intercal:\triangleq\mathbb{U}_{i, j}^n|0\leqslant i, j\leqslant N, 0\leqslant n \leqslant M\}, \\ &\mathbb{F} = \{[\chi u(x_i, y_j, t_n)v_x(x_i, y_j, t_n), 0]^\intercal:\triangleq\mathbb{F}_{i, j}^n|0\leqslant i, j\leqslant N, 0\leqslant n \leqslant M\}, \\ &\mathbb{G} = \{[\chi u(x_i, y_j, t_n)v_y(x_i, y_j, t_n), 0]^\intercal:\triangleq\mathbb{G}_{i, j}^n|0\leqslant i, j\leqslant N, 0\leqslant n \leqslant M\}, \\ &\mathbb{R} = \{[0, u(x_i, y_j, t_n)-v(x_i, y_j, t_n)]^\intercal:\triangleq\mathbb{R}_{i, j}^n|0\leqslant i, j\leqslant N, 0\leqslant n \leqslant M\} \end{align}

    on \Omega_{h\tau} , respectively, and employ the FDM to discretize Eqs. (3.1)–(3.4) for the interior points.

    Now, we focus on Eq. (3.1) at (x_i, y_j, t_n) for 1\leqslant i, j\leqslant N-1 and 0 < n \leqslant M ; we have

    \begin{align} \textbf{U}_{t}(x_i, y_j, t_n) &+ \textbf{F}_{x}(x_i, y_j, t_n)+\textbf{G}_{y}(x_i, y_j, t_n) = \textbf{D}[\textbf{U}_{xx}(x_i, y_j, t_n)+\textbf{U}_{yy}(x_i, y_j, t_n)] + \textbf{R}(x_i, y_j, t_n). \end{align} (3.5)

    First, the following formulas are applied to discretize the second derivative terms on the right-hand side of Eq. (3.5), i.e.,

    \begin{align} &\textbf{U}_{xx}(x_i, y_j, t_n) = \mathcal{A}^{-1}\delta _x^2\mathbb{U}_{i, j}^n + \mathcal{A}^{-1}(O_{xx})_{i, j}^n({h^4}), \; \; \; \; \textbf{U}_{yy}(x_i, y_j, t_n) = \mathcal{B}^{-1}\delta _y^2\mathbb{U}_{i, j}^n + \mathcal{B}^{-1}(O_{yy})_{i, j}^n({h^4}), \end{align} (3.6)

    where 1\leqslant i, j\leqslant N-1 , and \delta_x^2 , and \delta_y^2 are the central difference operators. Then, for the nonlinear advection terms \textbf{F}_{x}(x_i, y_j, t_n) and \textbf{G}_{y}(x_i, y_j, t_n) on the left-hand side of Eq. (3.5), the following Padé compact schemes [39] are applied to compute them, that is,

    \begin{align} &(1+\frac{h^2}{6}\delta_x^2)\textbf{F}_{x}(x_i, y_j, t_n) = \delta_x\mathbb{F}_{i, j}^{n}+(O_x)_{i, j}^n(h^4), \; \; \; \; \; \; \; (1+\frac{h^2}{6}\delta_y^2)\textbf{G}_{y}(x_i, y_j, t_n) = \delta_y\mathbb{G}_{i, j}^{n}+(O_y)_{i, j}^n(h^4), \end{align} (3.7)

    where 1\leqslant i, j\leqslant N-1 and 0\leqslant n\leqslant M . To be consistent with the accuracy of the interior points above, the following fourth-order boundary schemes [40] are employed, i.e.,

    \begin{align} &\textbf{F}_x(x_0, y_j, t_n)+\frac{14}{15}\textbf{F}_x(x_0+h, y_j, t_n) = \frac{1}{h}(\vec{\textbf{A}}\vec{ \mathbb{F}}_0)+(O_x)_{0, j}^n(h^4), \; \; \; \; \; \; \; 0\leqslant j\leqslant N, 0 < n\leqslant M, \end{align} (3.8)
    \begin{align} &\textbf{F}_x(x_N, y_j, t_n)-\frac{14}{15}\textbf{F}_x(x_N-h, y_j, t_n) = \frac{1}{h}(\vec{\textbf{B}}\vec{ \mathbb{F}}_N)+(O_x)_{N, j}^n(h^4), \; \; \; \; \; 0\leqslant j\leqslant N, 0 < n\leqslant M, \end{align} (3.9)
    \begin{align} &\textbf{G}_y(x_i, y_0, t_n)+\frac{14}{15}\textbf{G}_y(x_i, y_0+h, t_n) = \frac{1}{h}(\vec{\textbf{A}}\vec{ \mathbb{G}}_0)+(O_y)_{i, 0}^n(h^4), \; \; \; \; \; \; \; 0\leqslant i\leqslant N, 0 < n\leqslant M, \end{align} (3.10)
    \begin{align} &\textbf{G}_y(x_i, y_N, t_n)-\frac{14}{15}\textbf{G}_y(x_i, y_N-h, t_n) = \frac{1}{h}(\vec{\textbf{B}}\vec{ \mathbb{G}}_N)+(O_y)_{i, N}^n(h^4), \; \; \; \; \; 0\leqslant i\leqslant N, 0 < n\leqslant M, \end{align} (3.11)

    where

    \begin{align} &\vec{\textbf{A}} = \Big[-\frac{184}{75}, \frac{703}{180}, -\frac{89}{30}, \frac{67}{30}, -\frac{77}{90}, \frac{41}{300}\Big], \; \; \; \; \; \; \vec{\textbf{B}} = \Big[\frac{52}{25}, -\frac{1067}{180}, \frac{67}{10}, -\frac{41}{10}, \frac{133}{90}, -\frac{69}{300}\Big], \\ &\vec{ \mathbb{F}}_0 = [\mathbb{F}_{0, j}^n, \mathbb{F}_{1, j}^n, \mathbb{F}_{2, j}^n, \mathbb{F}_{3, j}^n, \mathbb{F}_{4, j}^n, \mathbb{F}_{5, j}^n]^\intercal, \; \; \; \; \; \; \; \; \vec{ \mathbb{F}}_N = [\mathbb{F}_{N, j}^n, \mathbb{F}_{N-1, j}^n, \mathbb{F}_{N-2, j}^n, \mathbb{F}_{N-3, j}^n, \mathbb{F}_{N-4, j}^n, \mathbb{F}_{N-5, j}^n]^\intercal, \\ &\vec{ \mathbb{G}}_0 = [\mathbb{G}_{i, 0}^n, \mathbb{G}_{i, 1}^n, \mathbb{G}_{i, 2}^n, \mathbb{G}_{i, 3}^n, \mathbb{G}_{i, 4}^n, \mathbb{G}_{i, 5}^n]^\intercal, \; \; \; \; \; \; \vec{ \mathbb{G}}_N = [\mathbb{G}_{i, N}^n, \mathbb{G}_{i, N-1}^n, \mathbb{G}_{i, N-2}^n, \mathbb{G}_{i, N-3}^n, \mathbb{G}_{i, N-4}^n, \mathbb{G}_{i, N-5}^n]^\intercal, \\ &(O_x)_{0, j}^n(h^4) = -\frac{1}{60}h^4\textbf{F}_{x^5}(x_0, y_j, t_n)+\frac{169}{1800}h^5\textbf{F}_{x^6}(x_0, y_j, t_n)+O(h^6), \\ &(O_x)_{N, j}^n(h^4) = -\frac{1}{60}h^4\textbf{F}_{x^5}(x_N, y_j, t_n)-\frac{281}{1800}h^5\textbf{F}_{x^6}(x_N, y_j, t_n)+O(h^6), \\ &(O_y)_{i, 0}^n(h^4) = -\frac{1}{60}h^4\textbf{F}_{y^5}(x_i, y_0, t_n)+\frac{169}{1800}h^5\textbf{F}_{y^6}(x_i, y_0, t_n)+O(h^6), \\ &(O_y)_{i, N}^n(h^4) = -\frac{1}{60}h^4\textbf{F}_{y^5}(x_i, y_N, t_n)-\frac{281}{1800}h^5\textbf{F}_{y^6}(x_i, y_N, t_n)+O(h^6). \end{align}

    Next, we employ the BDF-4 [36] to cope with the temporal derivative on the left-hand side of Eq. (3.5), i.e., for 1\leqslant i, j\leqslant N-1 and 4\leqslant n\leqslant M,

    \begin{align} {\textbf{U}_{t}(x_i, y_j, t_n)}& = \frac{1}{12\tau}\Big(25\mathbb{U}_{{i, j}}^{n}-\sum\limits_{k = 1}^4A_{k+1}\mathbb{U}_{{i, j}}^{n-k}\Big)+(O_t)_{i, j}^n(\tau^4) {:\triangleq\Delta_t\mathbb{U}_{i, j}^n+(O_t)_{i, j}^n(\tau^4)}. \end{align} (3.12)

    Substituting Eqs. (3.6)–(3.12) into Eq. (3.5), using the definitions above, combining Eqs. (3.2)–(3.4) and applying Theorem 2.1 for 4\leqslant n\leqslant M , we have

    \left\{ \begin{array}{l} \mathcal{AB}[\Delta_t\mathbb{U}_{i, j}^n + (\mathbb{F}_{x})_{i, j}^{n} +(\mathbb{G}_{y})_{i, j}^{n}] = \textbf{D}(\mathcal{B}\delta_x^2 + \mathcal{A}\delta _y^2)\mathbb{U}_{i, j}^{n}+\mathcal{AB}\mathbb{R}_{i, j}^{n} +\textbf{Q}_{i, j}^n, 1\leqslant i, j\leqslant N-1, & (3.13) \\ \mathbb{U}_{i, j}^0 = \mathbb{U}_0(x_{i}, y_j), \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0\leqslant i, j\leqslant N, & (3.14) \\ \mathbb{U}_{0, j}^n = \mathcal{A}\mathbb{U}_{0, j}^n+\textbf{Q}_{0, j}^n, \ \ \ \ \ \ \ \ \ \ \mathbb{U}_{N, j}^n = \mathcal{A}\mathbb{U}_{N, j}^n+\textbf{Q}_{N, j}^n, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0\leqslant j\leqslant N, & (3.15)\\ \mathbb{U}_{i, 0}^n = \mathcal{B}\mathbb{U}_{i, 0}^n+\textbf{Q}_{i, 0}^n, \ \ \ \ \ \ \ \ \ \ \mathbb{U}_{i, N}^n = \mathcal{B}\mathbb{U}_{i, N}^n+\textbf{Q}_{i, N}^n, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0\leqslant i\leqslant N, & (3.16) \end{array} \right.

    where \textbf{Q}_{{i, j}}^{n} = \mathcal{AB}(O_t)_{i, j}^n(\tau^4) +\{[\mathcal{B}(O_{xx})_{i, j}^n+\mathcal{A}(O_{yy})_{i, j}^n] +\mathcal{AB}[(O_x)_{i, j}^n+(O_x)_{0, j}^n+(O_x)_{N, j}^n+(O_y)_{i, j}^n+(O_y)_{i, 0}^n+(O_y)_{i, N}^n]\}(h^4) , \textbf{Q}_{0, j}^{n} = O_{0, j}^n(h^4) , \textbf{Q}_{N, j}^{n} = O_{N, j}^n(h^4) , \textbf{Q}_{i, 0}^{n} = O_{i, 0}^n(h^4) and \textbf{Q}_{i, N}^{n} = O_{i, N}^n(h^4) . Then, there exist positive constants c_1 , c_2 , c_3 , c_4 and c_5 ; it holds that

    {|\textbf{Q}_{i, j}^n|\leqslant} \left\{ \begin{array}{l}c_1h^4, & i = 0, 0\leqslant j\leqslant N, 4\leqslant n \leqslant M, \nonumber\\ c_2h^4, &j = 0, 0\leqslant i\leqslant N, 4\leqslant n \leqslant M, \nonumber\\ c_3(\tau^4+h^4), &1\leqslant i, j\leqslant N-1, 4\leqslant n \leqslant M, \nonumber\\ c_4h^4, & i = N, 0\leqslant j\leqslant N, 4\leqslant n \leqslant M, \nonumber\\ c_5h^4, &j = N, 0\leqslant i\leqslant N, 4\leqslant n \leqslant M.\nonumber\end{array} \right.

    Omitting \textbf{Q}_{i, j}^n in Eqs. (3.13)–(3.16) and replacing \mathbb{U} , \mathbb{F}_x , \mathbb{G}_y and \mathbb{R} with \textbf{U} , \textbf{F}_x , \textbf{G}_y and \textbf{R} , respectively, the following HOC scheme for solving the chemotaxis system given by (1.1)–(1.3) for 4\leqslant n\leqslant M can be obtained:

    \left\{ \begin{array}{l} \mathcal{AB}[\Delta_t\textbf{U}_{i, j}^n + (\textbf{F}_{x})_{i, j}^{n} +(\textbf{G}_{y})_{i, j}^{n}] = \textbf{D}(\mathcal{B}\delta_x^2 + \mathcal{A}\delta _y^2)\textbf{U}_{i, j}^{n}+\mathcal{AB}\textbf{R}_{i, j}^{n}, \ \ \ 1\leqslant i, j\leqslant N-1, & (3.17) \\ \textbf{U}_{i, j}^0 = \textbf{U}_0(x_{i}, y_j), \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0\leqslant i, j\leqslant N, & (3.18)\\ \textbf{U}_{0, j}^n = \mathcal{A}\textbf{U}_{0, j}^n, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \textbf{U}_{N, j}^n = \mathcal{A}\textbf{U}_{N, j}^n, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0\leqslant j\leqslant N, & (3.19)\\ \textbf{U}_{i, 0}^n = \mathcal{B}\textbf{U}_{i, 0}^n, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \textbf{U}_{i, N}^n = \mathcal{B}\textbf{U}_{i, N}^n, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0\leqslant i\leqslant N. & (3.20) \end{array} \right.

    Remark 3.1. For the function values containing the unknown time steps \textbf{U}_{i, j}^{n} in Eqs. (3.17)–(3.20), the following method is used to compute them, where k = 1, 2, 3, \cdots , which stands for the nonlinear cyclic iterative parameter, \textbf{U}^{n, (*)} denotes the approximate convergent solution at the n th time step and \sigma represents a small amount. That is, for the n th time step, the following is performed:

    Approximate \textbf{U}^{n, (k)} using \textbf{U}^{n, (k-1)} and solve v_x^{n, (k)} and v_y^{n, (k)} using Eqs. (3.3), (3.4) and (3.7);

    Compute the chemotaxis terms \textbf{F}_x^{n, (k)} and \textbf{G}_y^{n, (k)} using Eqs. (3.7)–(3.11);

    Solve the linear system given by (3.17)–(3.20) for \textbf{U}^{n, (k)} ;

    Set \textbf{U}^{n, (*)}\leftarrow \textbf{U}^{n, (k)} ; if \|\textbf{U}^{n, (k)}-\textbf{U}^{n, (k-1)}\| < \sigma , then proceed to the (n+1) th time step.

    Remark 3.2. Equations (3.17)–(3.20) demonstrate fully implicit compact scheme including five time steps, and its truncation error is O(\tau^4+h^4) , i.e., it has space-time fourth-order accuracy. It is noteworthy that the compactness here only refers to the spatial direction, since no more than nine mesh points are required for the 2D spatial mesh subdomain. However, it is non-compact in the temporal direction.

    Remark 3.3. For the HOC scheme given by (3.17)–(3.20), besides the values of \textbf{U}^{0} being known, the values of \textbf{U}^{1} , \textbf{U}^{2} and \textbf{U}^{3} are also needed; then, we can use the scheme given by (3.17)–(3.20) to compute the values of \textbf{U}^{n}(n\geqslant 4) in turn. Thus, in the following part, we discuss the computational strategies for solving \textbf{U}^{1} , \textbf{U}^{2} and \textbf{U}^{3} .

    Next, we consider Eq. (3.1) at (x_i, y_j, t_{n-\frac{1}{2}}) for 1\leqslant i, j\leqslant N-1 and 1\leqslant n \leqslant 3 ; we have

    \begin{align} \textbf{U}_{t}(x_i, y_j, t_{n-\frac{1}{2}}) + \textbf{F}_{x}(x_i, y_j, t_{n-\frac{1}{2}})&+\textbf{G}_{y}(x_i, y_j, t_{n-\frac{1}{2}})\\ & = \textbf{D}[\textbf{U}_{xx}(x_i, y_j, t_{n-\frac{1}{2}})+\textbf{U}_{yy}(x_i, y_j, t_{n-\frac{1}{2}})] + \textbf{R}(x_i, y_j, t_{n-\frac{1}{2}}). \end{align} (3.21)

    First, we employ the Crank-Nicolson (C-N) method to cope with the temporal derivative term in Eq. (3.21), that is,

    \begin{align} \textbf{U}_{t}(x_i, y_j, t_{n-\frac{1}{2}})& = \delta_t\mathbb{U}_{i, j}^{n-\frac{1}{2}}+ (O_t)_{i, j}^{n-\frac{1}{2}}({\tau ^2}) , \; \; \; \; 1\leqslant i, j\leqslant N-1, 1\leqslant n\leqslant 3. \end{align}

    Then, we employ the following method for the remaining parts, i.e.,

    \begin{align} \mathbf{\Theta}(x_i, y_j, t_{n-\frac{1}{2}}) = \frac{1}{2}[\mathbf{\Theta}(x_i, y_j, t_{n}) +\mathbf{\Theta}(x_i, y_j, t_{n-1})]+(\tilde{O}_\mathbf{\Theta})_{i, j}^{n-\frac{1}{2}}(\tau^2), 1\leqslant i, j\leqslant N-1, 1\leqslant n\leqslant 3, \end{align}

    where \mathbf{\Theta} could be \textbf{U}, \textbf{F}_x, \textbf{G}_y and \textbf{R} in Eq. (3.21). Then, Eq. (3.21) becomes

    \begin{align} &\delta_t\mathbb{U}_{i, j}^{n-\frac{1}{2}} + \frac{1}{2}[\textbf{F}_{x}(x_i, y_j, t_n)+\textbf{F}_{x}(x_i, y_j, t_{n-1}) +\textbf{G}_{y}(x_i, y_j, t_{n})+\textbf{G}_{y}(x_i, y_j, t_{n-1})] = \frac{1}{2}\textbf{D}[\textbf{U}_{xx}(x_i, y_j, t_{n})\\ &+\textbf{U}_{xx}(x_i, y_j, t_{n-1}) +\textbf{U}_{yy}(x_i, y_j, t_{n})+\textbf{U}_{yy}(x_i, y_j, t_{n-1})] +\frac{1}{2}[\textbf{R}(x_i, y_j, t_{n})+\textbf{R}(x_i, y_j, t_{n-1})] +[(O_t)_{i, j}^{n-\frac{1}{2}}\\ &+(\tilde{O}_{xx})_{i, j}^{n-\frac{1}{2}}+(\tilde{O}_{yy})_{i, j}^{n-\frac{1}{2}} +(\tilde{O}_{\textbf{F}_x})_{i, j}^{n-\frac{1}{2}} +(\tilde{O}_{\textbf{G}_x})_{i, j}^{n-\frac{1}{2}} +(\tilde{O}_{\textbf{R}})_{i, j}^{n-\frac{1}{2}}](\tau^2), \; 1\leqslant i, j\leqslant N-1, 1\leqslant n \leqslant 3. \end{align} (3.22)

    We use Eq. (3.6) to discretize \textbf{U}_{xx}(x_i, y_j, t_{n}) , \textbf{U}_{xx}(x_i, y_j, t_{n-1}) , \textbf{U}_{yy}(x_i, y_j, t_{n}) and \textbf{U}_{yy}(x_i, y_j, t_{n-1}) in the spatial direction of Eq. (3.22), and we use Eqs. (3.7)–(3.11) to compute \textbf{F}_{x}(x_i, y_j, t_n) , \textbf{F}_{x}(x_i, y_j, t_{n-1}) , \textbf{G}_{y}(x_i, y_j, t_{n}) and \textbf{G}_{y}(x_i, y_j, t_{n-1}) in the spatial direction. According to Eqs. (3.2)–(3.4) and Theorem 2.1, after rearrangement, the following form can be obtained:

    \left\{ \begin{array}{l} [\frac{1}{\tau}\mathcal{AB}-\frac{1}{2}\textbf{D}(\mathcal{B}\delta _x^2+\mathcal{A}\delta _y^2)]\mathbb{U}_{i, j}^{n} + \frac{1}{2}\mathcal{AB}[(\mathbb{F}_{x})_{i, j}^{n}+(\mathbb{F}_{x})_{i, j}^{n-1} + (\mathbb{G}_{y})_{i, j}^{n}+(\mathbb{G}_{y})_{i, j}^{n-1}] = [\frac{1}{\tau}\mathcal{AB}\nonumber\\ + \frac{1}{2}\textbf{D}(\mathcal{B}\delta _x^2+\mathcal{A}\delta _y^2)]\mathbb{U}_{i, j}^{n-1} + \frac{1}{2}\mathcal{AB}(\mathbb{R}_{i, j}^{n}+\mathbb{R}_{i, j}^{n-1}) + \textbf{P}_{i, j}^{n}, 1\leqslant i, j\leqslant N-1, 1\leqslant n\leqslant 3, & (3.23)\\ \mathbb{U}_{i, j}^0 = \mathbb{U}_0(x_{i}, y_{j}), \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad 0\leqslant i, j\leqslant N, & (3.24)\\ \mathbb{U}_{0, j}^{n} = \mathcal{A}\mathbb{U}_{0, j}^{n}+\textbf{P}_{0, j}^{n}, \qquad \mathbb{U}_{N, j}^{n} = \mathcal{A}\mathbb{U}_{N, j}^{n}+\textbf{P}_{N, j}^{n}, \qquad 0\leqslant j\leqslant N, 1\leqslant n\leqslant 3, & (3.25)\\ \mathbb{U}_{i, 0}^{n} = \mathcal{B}\mathbb{U}_{i, 0}^{n}+\textbf{P}_{i, 0}^{n}, \qquad \mathbb{U}_{i, N}^{n} = \mathcal{B}\mathbb{U}_{i, N}^{n}+\textbf{P}_{i, N}^{n}, \qquad 0\leqslant i\leqslant N, 1\leqslant n\leqslant 3, & (3.26) \end{array} \right.

    where \textbf{P}_{{i, j}}^{n} = \mathcal{AB}[(O_t)_{i, j}^{n-\frac{1}{2}} +(\tilde{O}_{xx})_{i, j}^{n-\frac{1}{2}}+(\tilde{O}_{yy})_{i, j}^{n-\frac{1}{2}} +(\tilde{O}_{\textbf{F}_x})_{i, j}^{n-\frac{1}{2}} +(\tilde{O}_{\textbf{G}_x})_{i, j}^{n-\frac{1}{2}} +(\tilde{O}_{\textbf{R}})_{i, j}^{n-\frac{1}{2}}](\tau^2) +\frac{1}{2}\{[\mathcal{B}((O_{xx})_{i, j}^n+(O_{xx})_{i, j}^{n-1}) +\mathcal{A}(O_{yy})_{i, j}^n+\mathcal{A}(O_{yy})_{i, j}^{n-1}] +\mathcal{AB}[(O_x)_{i, j}^n+(O_x)_{0, j}^n+(O_x)_{N, j}^n+(O_y)_{i, j}^n+(O_y)_{i, 0}^n+(O_y)_{i, N}^n +(O_x)_{i, j}^{n-1}+(O_x)_{0, j}^{n-1}+(O_x)_{N, j}^{n-1}+(O_y)_{i, j}^{n-1}+(O_y)_{i, 0}^{n-1}+(O_y)_{i, N}^{n-1}]\}(h^4) . \textbf{P}_{0, j}^{n} = O_{0, j}^n(h^4) , \textbf{P}_{N, j}^{n} = O_{N, j}^n(h^4) , \textbf{P}_{i, 0}^{n} = O_{i, 0}^n(h^4) , and \textbf{P}_{i, N}^{n} = O_{i, N}^n(h^4) . Then, there exist positive constants c_6 , c_7 , c_8 , c_9 and c_{10} such that

    {|\textbf{P}_{i, j}^{n}|\leqslant} \left\{ \begin{array}{l} c_6h^4, & i = 0, 0\leqslant j\leqslant N, 1\leqslant n \leqslant 3, \nonumber\\ c_7h^4), & j = 0, 0\leqslant i\leqslant N, 1\leqslant n \leqslant 3, \nonumber\\ c_8(\tau^2+h^4), & 1\leqslant i, j\leqslant N-1, 1\leqslant n \leqslant 3, \nonumber\\ c_9h^4, & i = N, 0\leqslant j\leqslant N, 1\leqslant n \leqslant 3, \nonumber\\ c_{10}h^4, & j = N, 0\leqslant i\leqslant N, 1\leqslant n \leqslant 3. \end{array} \right.

    Omitting \textbf{P}_{i, j}^n in Eqs. (3.23)–(3.26) and replacing \mathbb{U} , \mathbb{F}_x , \mathbb{G}_y and \mathbb{R} with \textbf{U} , \textbf{F}_x , \textbf{G}_y and \textbf{R} , respectively, we can get the following scheme for the initial time steps to solve the chemotaxis system given by (1.1)–(1.3) as follows:

    \left\{ \begin{array}{l} [\frac{1}{\tau}\mathcal{AB}-\frac{1}{2}\textbf{D}(\mathcal{B}\delta _x^2+\mathcal{A}\delta _y^2)]\textbf{U}_{i, j}^{n} + \frac{1}{2}\mathcal{AB}[(\textbf{F}_{x})_{i, j}^{n}+(\textbf{F}_{x})_{i, j}^{n-1} + (\textbf{G}_{y})_{i, j}^{n}+(\textbf{G}_{y})_{i, j}^{n-1}]\nonumber\\ = [\frac{1}{\tau}\mathcal{AB}+ \frac{1}{2}\textbf{D}(\mathcal{B}\delta _x^2+\mathcal{A}\delta _y^2)]\textbf{U}_{i, j}^{n-1} + \frac{1}{2}\mathcal{AB}(\textbf{R}_{i, j}^{n}+\textbf{R}_{i, j}^{n-1}), 1\leqslant i, j\leqslant N-1, 1\leqslant n\leqslant 3, & (3.27) \\ \textbf{U}_{i, j}^0 = \textbf{U}_0(x_{i}, y_{j}), \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad 0\leqslant i, j\leqslant N, & (3.28)\\ \textbf{U}_{0, j}^{n} = \mathcal{A}\textbf{U}_{0, j}^{n}, \qquad \textbf{U}_{N, j}^{n} = \mathcal{A}\textbf{U}_{N, j}^{n}, \qquad\qquad\qquad\qquad 0\leqslant j\leqslant N, 1\leqslant n\leqslant 3, & (3.29)\\ \textbf{U}_{i, 0}^{n} = \mathcal{B}\textbf{U}_{i, 0}^{n}, \qquad \textbf{U}_{i, N}^{n} = \mathcal{B}\textbf{U}_{i, N}^{n}, \qquad\qquad\qquad\qquad 0\leqslant i\leqslant N, 1\leqslant n\leqslant 3. & (3.30) \end{array} \right.

    Remark 3.4. The same method in the previous section can be used to compute (\textbf{F}_{x})_{i, j}^{n} , (\textbf{F}_{x})_{i, j}^{n-1} , (\textbf{G}_{y})_{i, j}^{n} , (\textbf{G}_{y})_{i, j}^{n-1} , \textbf{R}_{i, j}^{n} and \textbf{R}_{i, j}^{n-1} in Eq. (3.27), as described in Remark 3.1.

    Remark 3.5. Equations (3.27)–(3.30) constitute a two-level implicit scheme, and its truncation error is O(\tau^2+h^4) .

    Remark 3.6. Since the C-N method is used to discretize the time derivative term, the scheme given by (3.27)–(3.30) for the initial time steps is unconditionally stable. On the other hand, the formula (3.12) is a k -step method for time integration. In this work, only the fourth-order BDF scheme ( k = 4 ) is used, which is widely employed to solve stiff differential equations [37]. Therefore, according to [36,37,38], the HOC scheme given by (3.17)–(3.20) is A -stable.

    As described in Subsections 3.1 and 3.2 above, we use the two-level implicit scheme given by (3.27)–(3.30) to compute \textbf{U}_{i, j}^{1} , \textbf{U}_{i, j}^{2} and \textbf{U}_{i, j}^{3} based on a known \textbf{U}_{i, j}^{0} . The Richardson extrapolation technique is applied to enhance the accuracy of the scheme given by (3.27)–(3.30) in time and ensure that it is consistent with that of the scheme given by (3.17)–(3.20); the following extrapolation formula is used, that is,

    \begin{align} &\widehat{\textbf{U}}_{i, j}^n(h, \tau) = \frac{1}{3}[4\textbf{U}_{i, j}^{2n}(h, \tau/2)-\textbf{U}_{i, j}^n(h, \tau)], \; \; \; \; \; \; {0\leqslant i, j \leqslant N, 1\leqslant n \leqslant 3.} \end{align} (3.31)

    Here, \widehat{\textbf{U}}_{i, j}^n(h, \tau) stands for the extrapolated solution at the n th time step. \textbf{U}_{i, j}^n(h, \tau) represents the computed solution obtained by using the scheme given by (3.27)–(3.30) with the temporal step length \tau , whereas \textbf{U}_{i, j}^{2n}(h, \tau/2) denotes the computed solution obtained by using the scheme given by (3.27)–(3.30) with \tau/2 .

    Remark 3.7. Equation (3.31) can extrapolate the accuracy in the temporal direction from second to fourth order.

    In this part, based on the compact difference schemes proposed above, we will establish the corresponding numerical algorithm. On one hand, because the slow convergence speed of the classical iterative method (for instance, Gauss–Seidel, Jacobi or successive over-relaxation iterations) used to solve the algebraic systems arise out of the fully implicit scheme at each time step, we aim to employ a multigrid method to accelerate the convergence speed. On the other hand, we want to establish a positivity-preserving algorithm to obtain the non-negative solutions of u and v without losing the fourth-order accuracy. Finally, we structure a time advancement algorithm based on these algorithms above to solve the system given by (1.1)–(1.3).

    Since the idea of the multigrid method was proposed in the 1930s, until Professor Brandt published his pioneering work [41] in 1977, this method was increasingly applied to solve engineering and technical problems. Up to now, the multigrid method has been widely applied in numerous disciplines and engineering fields, such as computational fluid dynamics and computational biology. The multigrid method has been theoretically proved to be a first-rank numerical computational method for linear elliptic problems [41,42]. Its convergence rate is independent of the mesh scale, and the computational speed of the existing computational program using the classical iterative method can be increased by employing the multigrid method with 1 to 2 orders of magnitude, which is especially suitable for application in large-scale engineering numerical computational problems. Nowadays, the multigrid method is also used more and more to solve non-elliptic problems to speed up the convergence of each time step and satisfy the needs of practical problems [43,44].

    The multigrid method is implemented by using a cycling algorithm. A multigrid cycle V includes three elements: relaxation, restriction and interpolation operators. We use alternating direction line Gauss–Seidel (ALGS) [45] iteration for the relaxation operator, as well as half-weighted and fully weighted restriction operators and a bilinear interpolation operator [42]. As the schemes given by (3.27)–(3.30) and (3.17)–(3.20) are both nonlinear, we employ a multigrid full approximation scheme (FAS) [42,44] to accommodate nonlinearities. For simplicity, we formally denote the algebraic equations arising from Eqs. (3.27)–(3.30) and (3.17)–(3.20) at each time step as

    \begin{align} \mathcal{L}^{h}\textbf{U}^{h} = \mathcal{F}^{h}. \end{align} (4.1)

    The multigrid FAS algorithm is implemented recursively. Here, we only give the following two-level FAS V(\nu_1, \nu_2) cycle algorithm based on Eq. (4.1):

    Algorithm 4.1: Two-level FAS V(\nu_1, \nu_2) cycle algorithm
    \textbf{ Step 1: } \ Using \ ALGS \ iteration \ \nu_{1} \ times \ {at} \ the \ fine \ mesh \ level \ to \ solve \ \mathcal{L}^{h}\textbf{U}^{h} \ = \ \mathcal{F}^{h} \ and \ computing
    \quad\quad\quad its \ residual, \ i.e., \ Res^{h} \ = \ \mathcal{F}^{h}-\mathcal{L}^{h}\textbf{U}^{h} ;
    \textbf{ Step 2: } \ Restricting \ \textbf{U}^{h} \ and \ Res^{h} \ to \ the \ next \ coarse \ mesh: \ \overline{\textbf{U}}^{2h}\leftarrow \ I_{h}^{2h}\textbf{U}^{h}, \ \overline{Res}^{2h}\leftarrow \ \overline{\mathcal{I}}_{h}^{2h}Res^{h};
    \textbf{ Step 3: } Using \ ALGS \ iteration \ \nu_{1} \ times \ {at} \ the \ coarse \ mesh \ level \ to \ solve \ \mathcal{L}^{2h}\textbf{U}^{2h} \ = \ \mathcal{L}^{2h}\overline{\textbf{U}}^{2h}+\overline{Res}^{2h};
    \textbf{ Step 4: } \ Interpolating \ the \ correction \ errors \ from \ the \ coarse \ to \ the \ fine \ mesh \ levels: \ \Delta \ \textbf{U}^{h}\leftarrow \ \mathcal{I}_{2h}^{h}(\textbf{U}^{2h}-\overline{\textbf{U}}^{2h});
    \textbf{ Step 5: } \ Updating \ \textbf{U}^{h} \ {at} \ the \ fine \ mesh \ level: \ \textbf{U}^{h}\leftarrow \ \textbf{U}^{h}+\Delta \ \textbf{U}^{h};
    \textbf{ Step 6: } \ Using \ ALGS \ iteration \ \nu_{2} \ times \ {at} \ the \ fine \ mesh \ level \ to \ solve \ \mathcal{L}^{h}\textbf{U}^{h} \ = \ \mathcal{F}^{h}.

    Remark 4.1. Here, \mathcal{L}^{h} and \mathcal{L}^{2h} are the difference operators; \mathcal{I}_h^{2h} (half-weighted) and \overline{\mathcal{I}}_h^{2h} (fully weighted) are the restriction operators, which are used to restrict the approximate value and residual from the fine mesh level h to coarse mesh level 2h ; and, \mathcal{I}_{2h}^{h} is the interpolation operator, which is used to transfer the error correction from the coarse to fine mesh levels. Res^h is the residual at the fine mesh level, and Res^{2h} is the residual at the adjacent coarse mesh level; \nu_1 and \nu_2 represent the pre-smooth and post-smooth numbers, respectively.

    Remark 4.2. In this work, we employ fully multi-level algorithms for all calculations. For instance, if the finest grid is 256\times256 , an eight-level algorithm is used. Under such conditions, the coarsest grid level only has 2\times2 grids.

    When the schemes given by (3.27)–(3.30) and (3.17)–(3.20) are employed to approximate the solution of the chemotaxis system, it may generate negative values where blow-up occurs. To obtain the non-negative values of u and v for all time levels, we define the following average solution for each time step:

    \begin{align} &\overline{\textbf{U}}_{i, j}^n = \frac{1}{\Delta x_i\Delta y_j}\iint\limits_{\Omega_{i, j}}\textbf{U}(x, y, t_n)dxdy, \end{align} (4.2)

    where \Delta x_i = x_{i+1}-x_{i-1} , \Delta y_j = y_{j+1}-y_{j-1} , 1\leqslant i, j\leqslant N-1 and \Omega_{i, j} = \{[x_{i-1}, x_{i+1}]\times[y_{j-1}, y_{j+1}]|1\leqslant i, j\leqslant N-1\} . Similar to the numerical integration formula of the one-variable function (i.e., Simpson formula), we structure the numerical integration of the two-variable function \textbf{U}(x, y) on \Omega_{i, j} for each time step, i.e., for 1\leqslant i, j\leqslant N-1 ,

    \begin{align} &\iint\limits_{\Omega_{i, j}}\textbf{U}(x, y, t_n)dxdy\thickapprox \frac{1}{6}\Delta x_i\Delta y_j\sum\limits_{k = 1}^3[\textbf{U}(x_{i+k-2}, y_{j}, t_n)+\textbf{U}(x_{i}, y_{j+k-2}, t_n)]. \end{align} (4.3)

    Substituting Eq. (4.3) into Eq. (4.2) at each time step, we obtain the following average solutions \overline{\textbf{U}} , that is, for 1\leqslant i, j\leqslant N-1 and 0 < n\leqslant M , we have

    \begin{align} \overline{\textbf{U}}_{i, j}^n\thickapprox& \frac{1}{6}\sum\limits_{k = 1}^3[\textbf{U}_{i+k-2, j}^n+\textbf{U}_{i, j+k-2}^n]. \end{align} (4.4)

    Next, a positivity-preserving limiter (PPL) [25,46] is employed to eliminate the negative values of \textbf{U}_{i, j}^n , 0\leqslant i, j\leqslant N-1 , 0 < n\leqslant M , that is,

    \begin{align} &\widetilde{\textbf{U}}_{i, j}^n = \min\left\{\Big|\frac{\overline{\textbf{U}}_{i, j}^n}{\overline{\textbf{U}}_{i, j}^n-\varpi_{i, j}^n}\Big|, 1\right\} (\textbf{U}_{i, j}^n-\overline{\textbf{U}}_{i, j}^n)+\overline{\textbf{U}}_{i, j}^n, \; \; \; \; \; 0\leqslant i, j\leqslant N, 0 < n\leqslant M, \end{align} (4.5)

    where \varpi_{i, j}^n = \min\{\textbf{U}_{i, j}^n, \textbf{U}_{i-1, j}^n, \textbf{U}_{i+1, j}^n, \textbf{U}_{i, j+1}^n, \textbf{U}_{i, j-1}^n\}, \; \; 1\leqslant i, j\leqslant N-1, \; 0 < n\leqslant M. For 1\leqslant i, j\leqslant N-1 and 0 < n\leqslant M , we can get the positivity-preserving solution \widetilde{\textbf{U}}_{i, j}^n\geqslant 0 without losing the fourth-order accuracy obtained by the proposed scheme. Finally, we establish the following algorithm for the n th time step:

    Algorithm 4.2: Positivity-preserving algorithm
    \textbf{ 1. } \ Compute \ \textbf{U}_{i, \ j}^n \ using \ the \ {schemes \ given \ by} \ (3.27) \ - \ (3.30) \ and \ (3.17) \ - \ (3.20);
    \textbf{ 2. } \ Compute \ the \ averaged \ solution \ \overline{\textbf{U}}_{i, \ j}^n \ using \ Eq{.} \ (4.4);
    \textbf{ 3. } Compute \ the \ positive \ solution \ \widetilde{\textbf{U}}_{i, \ j}^n \ using \ the \ PPL \ (4.5).

    On the basis of the derivation process described in Section 3 and the proposed algorithms in Subsections 4.1 and 4.2, we employ the difference schemes given by (3.27)–(3.30) and (3.17)–(3.20) to approximate the solution of the original system given by (1.1)–(1.3), which mainly includes the following three steps: first, provide the initial value of \textbf{U}_{i, j}^{0} ; second, we employ Eqs. (3.7)–(3.11) and (3.27)–(3.31) to compute the values of \textbf{U}_{i, j}^{1} , \textbf{U}_{i, j}^{2} and \textbf{U}_{i, j}^{3} ; finally, we employ Eqs. (3.7)–(3.11) and (3.17)–(3.20) to compute the value of \textbf{U}_{i, j}^{n} (n\geqslant 4) . As (3.27)–(3.30) and (3.17)–(3.20) involve nonlinear parts, we establish the following algorithm with a nonlinear iteration tactic to solve them.

    Algorithm 4.3: Time advancement algorithm
    \ {Establish} \ the \ initial \ values: \ \textbf{U}_{i, j}^0 = \textbf{U}_0(x_i, y_j) ;
    \textbf{ for } n = 1, 2, \cdots, M \textbf{ do }
       \textbf{ if } 1\leqslant n \leqslant 3 \textbf{ then }
         \textbf{R}_{i, j}^{2(n-1)} = { \textbf{R}}_{i, j}^{2(n-1), (*)}(h, \tau/2), \textbf{R}_{i, j}^{(n-1)} = { \textbf{R}}_{i, j}^{(n-1), (*)}(h, \tau) ;
         Computing \ the \ following \ {chemotaxis \ parts \ at \ {the} \ (n-1)th \ time \ step} \ using \ {Eqs.} \ (3.7)-(3.11) ;
         (\textbf{F}_x)_{i, j}^{2(n-1)} = { (\textbf{F}_x)}_{i, j}^{2(n-1), (*)}(h, \tau/2), (\textbf{F}_x)_{i, j}^{(n-1)} = { (\textbf{F}_x)}_{i, j}^{(n-1), (*)}(h, \tau) ;
         (\textbf{G}_y)_{i, j}^{2(n-1)} = { (\textbf{G}_y)}_{i, j}^{2(n-1), (*)}(h, \tau/2), (\textbf{G}_y)_{i, j}^{(n-1)} = { (\textbf{G}_y)}_{i, j}^{(n-1), (*)}(h, \tau) ;
    \textbf{ for } k = 1, 2, 3, \cdots \textbf{ do }
         {Computing} \ \ \textbf{U}_{0, j}^n, \ \textbf{U}_{N, j}^n, \ \textbf{U}_{i, 0}^n \ {and} \ \textbf{U}_{i, N}^n \ using \ Eqs. \ (3.29)-(3.30), \ i.e.,
         { \textbf{U}_{0, j}^{2n, k} = \mathcal{A}\textbf{U}_{0, j}^{2n, k-1}(h, \tau/2), \textbf{U}_{N, j}^{2n, k} = \mathcal{A}\textbf{U}_{N, j}^{2n, k-1}(h, \tau/2);} {\textbf{U}_{0, j}^{n, k} = \mathcal{A}\textbf{U}_{0, j}^{n, k-1}(h, \tau), \textbf{U}_{N, j}^{n, k} = \mathcal{A}\textbf{U}_{N, j}^{n, k-1}(h, \tau); }
         { \textbf{U}_{i, 0}^{2n, k} = \mathcal{A}\textbf{U}_{i, 0}^{2n, k-1}(h, \tau/2), \textbf{U}_{i, N}^{2n, k} = \mathcal{A}\textbf{U}_{i, N}^{2n, k-1}(h, \tau/2);} {\textbf{U}_{i, 0}^{n, k} = \mathcal{A}\textbf{U}_{i, 0}^{n, k-1}(h, \tau), \textbf{U}_{i, N}^{n, k} = \mathcal{A}\textbf{U}_{i, N}^{n, k-1}(h, \tau); }
         \textbf{R}_{i, j}^{2n, (k)} = { \textbf{R}}_{i, j}^{2n, (k-1)}(h, \tau/2), \textbf{R}_{i, j}^{n, (k)} = { \textbf{R}}_{i, j}^{n, (k-1)}(h, \tau) ;
         Computing \ the \ following \ {chemotaxis \ parts \ at \ {the} \ (n)th \ time \ step} \ using \ {Eqs.} \ (3.7)-(3.11) ;
         (\textbf{F}_x)_{i, j}^{2n, (k)} = { (\textbf{F}_x)}_{i, j}^{2n, (k-1)}(h, \tau/2), (\textbf{F}_x)_{i, j}^{n, (k)} = { (\textbf{F}_x)}_{i, j}^{n, (k-1)}(h, \tau) ;
         (\textbf{G}_y)_{i, j}^{2n, (k)} = { (\textbf{G}_y)}_{i, j}^{2n, (k-1)}(h, \tau/2), (\textbf{G}_y)_{i, j}^{n, (k)} = { (\textbf{G}_y)}_{i, j}^{n, (k-1)}(h, \tau) ;
         {Computing} \ \ \textbf{U}_{i, j}^{2n, (k)}(h, \tau/2) \ and \ \textbf{U}_{i, j}^{n, (k)}(h, \tau) \ \ using \ {Eqs.} \ \ (3.27)-(3.30) \ {and} \ \ Algorithm \ 4.1; \
         \textbf{ if } ||\textbf{U}_{{i, j}}^{2n, (k)}(h, \tau/2)-\textbf{U}_{{i, j}}^{2n, (k-1)}(h, \tau/2)|| < \sigma and ||\textbf{U}_{{i, j}}^{n, (k)}(h, \tau)-\textbf{U}_{{i, j}}^{n, (k-1)}(h, \tau)|| < \sigma \textbf{ then }
           \textbf{U}_{{i, j}}^{2n, (*)}(h, \tau/2) = \textbf{U}_{{i, j}}^{2n, (k)}(h, \tau/2), \textbf{U}_{{i, j}}^{n, (*)}(h, \tau) = \textbf{U}_{{i, j}}^{n, (k)}(h, \tau) ;
         \textbf{ end if }
       \textbf{ end for }
         Computing \ the \ extrapolated \ solutions \ \widehat{\textbf{U}}_{{i, \ j}}^{n, \ (*)} \ using \ {Eq.} (3.31);
    \textbf{ else }
       \textbf{ for } \ k = 1, 2, 3, \cdots, \textbf{ do }
         {Computing} \ \ \textbf{U}_{0, j}^n, \ \textbf{U}_{N, j}^n, \ \textbf{U}_{i, 0}^n \ {and} \ \textbf{U}_{i, N}^n \ using \ Eqs. \ (3.15)-(3.16) \, \ i.e., \
         { \textbf{U}_{0, j}^{n, k} = \mathcal{A}\textbf{U}_{0, j}^{n, k-1}(h, \tau), \textbf{U}_{N, j}^{n, k} = \mathcal{A}\textbf{U}_{N, j}^{n, k-1}(h, \tau); } { \textbf{U}_{i, 0}^{n, k} = \mathcal{A}\textbf{U}_{i, 0}^{n, k-1}(h, \tau), \textbf{U}_{i, N}^{n, k} = \mathcal{A}\textbf{U}_{i, N}^{n, k-1}(h, \tau); }
         \textbf{R}_{i, j}^{n, (k)} = { \textbf{R}}_{i, j}^{n, (k-1)}(h, \tau) ;
         Computing \ the \ following \ {chemotaxis \ parts \ at \ {the} \ \ (n)th \ time \ step} \ using \ {Eqs.} \ (3.7)-(3.11) ;
         (\textbf{F}_x)_{i, j}^{n, (k)} = { (\textbf{F}_x)}_{i, j}^{n, (k-1)}(h, \tau), (\textbf{G}_y)_{i, j}^{n, (k)} = { (\textbf{G}_y)}_{i, j}^{n, (k-1)}(h, \tau) ;
         {Computing} \ \ \ \textbf{U}_{i, j}^{n, (k)}(h, \tau) \ using \ Eqs. \ (3.17)-(3.20) \ and \ Algorithm \ 4.1;
         \textbf{ if } ||\textbf{U}_{{i, j}}^{n, (k)}(h, \tau)-\textbf{U}_{{i, j}}^{n, (k-1)}(h, \tau)|| < \sigma \textbf{ then } \textbf{U}_{{i, j}}^{n, (*)}(h, \tau) = \textbf{U}_{{i, j}}^{n, (k)}(h, \tau) ;
         \textbf{ end if }
        \textbf{ end for }
       \textbf{ end if }
        Computing \ the \ positive \ solutions \ \widetilde{\textbf{U}}_{i, j}^n \ {using} \ Algorithm \ 4.2.
    \textbf{ end for }

    In this part, we give several numerical experiments to test the various properties of the proposed scheme and algorithm when solving the system given by (1.1)–(1.3). First, in Subsection 5.1 below, the accuracy and stability of the proposed scheme and algorithm in the absence of PPL and in the presence of a PPL when solving the chemotaxis system are tested. Second, in Subsection 5.2 below, we simulate the blow-up phenomena for the system given by (1.1)–(1.3), compare the obtained results from the proposed scheme and algorithm before and after using the positivity-preserving algorithm and verify the mass conservation and energy dissipation of the proposed method.

    We consider the following general type of Eq. (1.1) with source terms to test the accuracy, that is,

    \begin{align} \left\{ \begin{array}{ll} u_t = \triangle u-\nabla\cdot(u\nabla v)-\frac{4u}{2+\sin(x+y)}+\frac{2u^2[\cos(2x+2y)-2\sin(x+y)]}{[2+\sin(x+y)]^2}, &0\leqslant x, y\leqslant 2\pi, \; t > 0, \\ v_t = \triangle v-v+u-\frac{4v}{2+\sin(x+y)}, &0\leqslant x, y\leqslant 2\pi, \; t > 0. \end{array} \right. \end{align}

    We use the periodic boundary conditions, and its analytical solution is u(x, y, t) = v(x, y, t) = (2+\sin(x+y))e^{-2t}. By observing the expression of this analytical solution, we can see that the global solutions of this problem are positive values in the entire physical domain. Therefore, in this computation, we will use the difference scheme proposed in this work to approximate the solution for this problem in the absence of a PPL.

    In Table 1, we take \tau = 0.1h and test the convergence of the HOC scheme at the final time T = 0.1 for Problem 5.1.1. From this table, the fourth-order accuracy is obtained by using the proposed scheme. The obtained results for the LDG method in [3] (the time step size \tau = 0.001h^2 is used to solve this problem) represent third-order accuracy, which shows the advantage of our method, that is, we can obtain higher accuracy with a larger time step. Table 2 displays the convergence of the HOC scheme at T = 1 and T = 5 for Problem 5.1.1, where \tau = 0.5h . Table 3 shows the ||\cdot||_2 errors of u(x, y, t) and v(x, y, t) for different step ratios \lambda = \tau/h for Problem 5.1.1, where T = 1 and \tau = \lambda h . From these tables, we obtain that the proposed method still converges at a fourth-order rate when computing long times and large step ratios, which indicates that our method has better stability.

    Table 1.  Convergence of the HOC scheme for Problem 5.1.1 at T = 0.1 .
    N LDG scheme [3] HOC scheme
    ||u||_\infty Rate ||u||_2 Rate ||u||_\infty Rate ||u||_2 Rate
    8 5.79E-03 2.53E-02 9.415E-03 3.974E-02
    16 8.79E-04 2.69 3.87E-03 2.71 3.359E-04 4.81 1.379E-03 4.85
    32 1.19E-04 2.91 5.18E-04 2.90 1.912E-05 4.13 7.711E-05 4.16
    64 1.49E-05 3.00 6.57E-05 2.98 1.112E-06 4.10 4.431E-06 4.12

     | Show Table
    DownLoad: CSV
    Table 2.  Convergence of the HOC scheme for Problem 5.1.1 with \tau = 0.5h .
    N T=1 T=5
    ||\cdot||_\infty Rate ||\cdot||_2 Rate ||\cdot||_\infty Rate ||\cdot||_2 Rate
    Cell density u(x, y, t)
    20 1.472E-03 6.477E-03 5.428E-06 2.415E-05
    40 1.046E-04 3.81 4.621E-04 3.81 2.289E-07 4.57 1.018E-06 4.57
    80 6.404E-06 4.03 2.832E-05 4.03 1.840E-08 3.64 5.259E-08 4.27
    160 3.902E-07 4.04 1.726E-06 4.04 6.967E-10 4.72 3.043E-09 4.11
    Chemoattractant concentration v(x, y, t)
    20 1.459E-03 6.488E-03 5.428E-06 2.415E-05
    40 1.041E-04 3.81 4.629E-04 3.81 2.289E-07 4.57 1.018E-06 4.57
    80 6.382E-06 4.03 2.837E-05 4.03 1.840E-08 3.64 5.259E-08 4.27
    160 3.889E-07 4.04 1.729E-06 4.04 6.967E-10 4.72 3.043E-09 4.11

     | Show Table
    DownLoad: CSV
    Table 3.  ||\cdot||_{2} errors of the HOC scheme for Problem 5.1.1 with T = 1, \tau = \lambda h .
    N \lambda ||u||_{2} ||v||_2
    0.4 2.8366512125E-05 2.8485324431E-05
    64 0.8 4.2917900191E-04 4.2928242249E-04
    1.6 5.7178346316E-03 5.7179298620E-03
    0.8 2.8033654005E-05 2.8041090035E-05
    128 1.6 4.1876675539E-04 4.1877318296E-04
    3.2 5.5523910128E-03 5.5523969211E-03
    1.6 2.7710430090E-05 2.7710912582E-05
    256 3.2 4.1352388720E-04 4.1352430255E-04
    6.4 5.4716502752E-03 5.4716506553E-03

     | Show Table
    DownLoad: CSV

    Next, we consider the following general type of Eq. (1.1) with two additional fluxes[24]:

    \begin{align} \left\{ \begin{array}{ll} u_t = \triangle u-\nabla\cdot(\chi u\nabla v)+r_1(x, y, t), &0\leqslant x, y\leqslant 2\pi, \; t > 0, \\ v_t = \triangle v-v+u+r_2(x, y, t), &0\leqslant x, y\leqslant 2\pi, \; t > 0, \end{array} \right. \end{align}

    where \chi = 1 , and the homogeneous Neumann boundary condition (1.3) is applied.

    Case 1: The additional fluxes are r_1(x, y, t) = -2 \exp(-2t)(\cos^2(x)+\cos(x)\cos(y)+\cos^2(y)-1) , r_2(x, y, t) = 0 , with the initial condition u(x, y, 0) = v(x, y, 0) = \cos(x)+\cos(y), 0\leqslant x, y\leqslant 2\pi , and the analytical solution[23,24] {u(x, y, t) = v(x, y, t) = \exp(-t)(\cos(x)+\cos(y))}.

    Case 2: The additional fluxes are r_1(x, y, t) = -2e^{-4t}\cos(2x+2y) , r_2(x, y, t) = 0 , with the initial condition u(x, y, 0) = v(x, y, 0) = \cos(x+y), 0\leqslant x, y\leqslant 2\pi , and the analytical solution u(x, y, t) = v(x, y, t) = e^{-2t}\cos(x+y) .

    Table 4 shows the convergence of the HOC scheme when we take the time step size \tau = 0.001 at T = 0.1 when using the proposed method to solve Problem 5.1.2 Case 1: the obtained results are compared with those in [24]. According to Table 4, the results obtained via the proposed method in the absence of a PPL can converge at the fourth-order rate. After using the positivity-preserving algorithm, although the errors obtained via the proposed scheme in the presence of a PPL decreases slightly, it still converges at the fourth-order rate and is better than the results in [24], as it is only second-order accuracy in [24]. In Tables 5 and 6, the results computed by using the proposed method for Problem 5.1.2 are given respectively. And, the convergence of the HOC scheme in the presence and absence of a PPL are compared when \tau = 0.5h , T = 1 and T = 5 . By observing these tables, it can be found that the computed results for the proposed method can achieve the fourth-order accuracy in both cases, including with and without a PPL. Finally, we take the spatial meshes N = 64,128 and 256 , respectively, at the final time T = 1 . By using the proposed method to solve Problem 5.1.2 Case 1, we obtain the ||\cdot||_2 norm errors computed with and without a PPL for different step ratios \lambda . The computed results are listed in Table 7. On the basis of this table, with the increase of the step ratio \lambda , the proposed method can converge well regardless of PPL existence, which also reflects that the proposed HOC difference scheme has better stability.

    Table 4.  Convergence of the HOC scheme for Problem 5.1.2 Case 1 with \tau = 0.001, T = 0.1 .
    Ref. [24] HOC scheme
    Without PPL With PPL
    N ||u||_{2} Rate ||u||_{2} Rate ||u||_{2} Rate
    10 1.02E-01 9.275E-03 2.170E-01
    20 2.74E-02 1.90 6.000E-04 3.95 2.870E-02 2.92
    40 6.30E-03 2.12 2.624E-05 4.52 2.469E-03 3.54
    80 1.45E-03 2.12 1.593E-06 4.04 1.476E-04 4.06

     | Show Table
    DownLoad: CSV
    Table 5.  Convergence of the HOC scheme for Problem 5.1.2 Case 1 when \tau = 0.5h, T = 1 .
    N Without PPL With PPL
    ||\cdot||_\infty Rate ||\cdot||_2 Rate ||\cdot||_\infty Rate ||\cdot||_2 Rate
    Cell density u(x, y, t)
    16 4.670E-03 1.022E-02 4.846E-03 1.424E-02
    32 1.009E-04 5.53 2.447E-04 5.38 2.684E-04 4.17 7.271E-04 4.29
    64 1.686E-06 5.90 5.902E-06 5.37 2.048E-05 3.71 4.974E-05 3.87
    128 7.683E-08 4.46 2.329E-07 4.66 1.508E-06 3.76 3.539E-06 3.81
    Chemoattractant concentration v(x, y, t)
    16 3.154E-03 8.025E-03 3.457E-03 1.217E-02
    32 9.167E-05 5.10 2.391E-04 5.07 2.769E-04 3.64 7.280E-04 4.06
    64 2.233E-06 5.36 6.687E-06 5.16 2.129E-05 3.70 5.074E-05 3.84
    128 4.873E-08 5.52 1.986E-07 5.07 1.558E-06 3.77 3.574E-06 3.83

     | Show Table
    DownLoad: CSV
    Table 6.  Convergence of the HOC scheme for Problem 5.1.2 Case 1 when \tau = 0.5h, T = 5 .
    N Without PPL With PPL
    ||u||_\infty Rate ||u||_2 Rate ||u||_\infty Rate ||u||_2 Rate
    16 2.693E-03 1.418E-02 4.424E-03 2.430E-02
    32 6.517E-05 5.37 3.530E-04 5.33 2.106E-04 4.39 1.172E-03 4.37
    64 1.631E-06 5.32 9.106E-06 5.28 1.331E-05 3.98 7.517E-05 3.96
    128 4.703E-08 5.12 2.672E-07 5.09 8.955E-07 3.89 5.098E-06 3.88

     | Show Table
    DownLoad: CSV
    Table 7.  ||\cdot||_2 errors of the HOC scheme for Problem 5.1.2 Case 1 with T = 1, \tau = \lambda h .
    N \lambda Without PPL With PPL
    ||u||_2 ||v||_2 ||u||_2 ||v||_2
    0.4 6.4023871303E-06 7.2616742789E-06 6.1464555035E-05 6.2581655559E-05
    64 0.8 1.2105436765E-05 1.1596839477E-05 3.5050251455E-05 3.5272598246E-05
    1.6 1.4457697771E-04 1.4260648214E-04 1.4620731799E-04 1.4406453705E-04
    0.8 9.4507876168E-07 8.4039918076E-07 2.4504429125E-06 2.4129793104E-06
    128 1.6 1.4859666027E-05 1.4234907841E-05 1.4967935575E-05 1.4314292650E-05
    3.2 1.5091392012E-04 1.4863144808E-04 1.5093408941E-04 1.4864245029E-04
    1.6 1.0545853150E-06 9.9391781457E-07 1.0618094946E-06 9.9866861508E-07
    256 3.2 1.5143113254E-05 1.4499368302E-05 1.5145042129E-05 1.4500157168E-05
    6.4 1.5233575478E-04 1.4973374665E-04 1.5233626588E-04 1.4973395499E-04

     | Show Table
    DownLoad: CSV

    Table 8 lists the numerical convergence of HOC scheme Problem 5.1.2 Case 2 in the absence and presence of a PPL, where T = 0.2 and \tau = 0.1h , and it compares the computed results with those obtained via the method in [25]. The results in [25] have third-order accuracy before and after using a PPL, while the results computed via the method in this work converge at the fourth-order rate for both cases, including without and with a PPL, and the ||\cdot||_{\infty} , ||\cdot||_{2} errors are better than those in [25]. This also reflects the superiority of the HOC difference scheme proposed in this work.

    Table 8.  Convergence of the HOC scheme for Problem 5.1.2 Case 2 at T = 0.2 .
    N Ref. [25] HOC scheme
    ||u||_\infty Rate ||u||_2 Rate ||u||_\infty Rate ||u||_2 Rate
    Without PPL
    20 3.65E-03 3.04E-03 6.189E-04 9.217E-04
    40 4.55E-04 3.00 3.67E-04 3.05 3.554E-05 4.12 7.089E-05 3.70
    80 5.68E-05 3.00 4.58E-05 3.00 2.535E-06 3.81 4.677E-06 3.92
    160 7.02E-06 3.01 5.77E-06 2.99 1.741E-07 3.86 3.077E-07 3.93
    With PPL
    20 3.55E-02 3.62E-02 6.337E-04 1.019E-03
    40 4.37E-03 3.02 4.54E-03 3.00 3.788E-05 4.06 7.260E-05 3.81
    80 5.40E-04 3.02 5.61E-04 3.02 2.562E-06 3.89 4.706E-06 3.95
    160 6.74E-05 3.01 6.92E-05 3.02 1.745E-07 3.88 3.082E-07 3.93

     | Show Table
    DownLoad: CSV

    In the 2D case of the system given by (1.1)–(1.3), when the initial mass satisfies certain critical values, the solutions will blow up in a finite time. To capture the blow-up time of the solutions, referring to [47], the following adaptive technology is used to obtain the optimal time step sizes, i.e.,

    \begin{align} &\tau_n = h\parallel\mathbb{U}^n\parallel_{\infty}^{-1}, \; \; \; \; \; n\geqslant0. \end{align} (5.1)

    The computation is terminated if ||\mathbb{U}||_{\infty}\geqslant 10^{\nu} , where \nu denotes the maximum order of magnitude of cell density when the actual problem experiences a blow up. And, t_n = \sum\limits_{m = 0}^{n-1}\tau_m is regarded as the approximation of the blow-up time T .

    In a rectangular region [-0.5, 0.5]\times[-0.5, 0.5] , we first take \chi = d = 1 and test the initial-boundary value problem (IBVP) [23,24,25] for the chemotaxis system given by (1.1)–(1.3) with the homogeneous Neumann boundary conditions given by (1.3), where its initial conditions are as follow:

    \begin{align} \left\{ \begin{array}{ll} u(x, y, 0) = 840e^{-84(x^2+y^2)}, &-0.5\leqslant x, y\leqslant 0.5, \; t > 0, \\ v(x, y, 0) = 420e^{-42(x^2+y^2)}, &-0.5\leqslant x, y\leqslant0.5, \; t > 0. \end{array} \right. \end{align}

    The solution of this problem will blow up in finite time at the center of this rectangular region.

    (1) Finite-time blow up. Figure 2 plots the numerical solutions of cell density u(x, y, t) as obtained by using the proposed method in the presence of a PPL, to solve the IBVP 5.2.1 at t = 0 , t = 5\times10^{-5} , t = 10^{-4} and t = 1.2\times10^{-4} , where the initial time step size \tau_0 = 5\times10^{-6} and the fixed spatial mesh numbers are 100\times 100 . In [22], the blow-up time t is approximately equal to 1.21\times10^{-4} , which is verified by the authors. Based on observation, we can expect to see that Problem 5.2.1 displays blow-up phenomena around t\thickapprox 1.2\times10^{-4} at the center of the rectangular region [-0.5, 0.5]\times[-0.5, 0.5] . It shows that the maximum peak values of the cell density u(x, y, t) gradually increase over time and present a very sharp aiguille structure at the central zone. We find that u(x, y, t) is strictly positive during time evolution by using the proposed method in the presence of a PPL. Figure 3 displays the projections of u(x, y, t) on the xou plane for Problem 5.2.1 at two pre-blow-up times t = 10^{-4} and t = 1.2\times10^{-4} , with \tau_0 = 5\times10^{-6} and N = 100 . They intuitively show that the proposed method in this work can achieve positivity-preserving capability.

    Figure 2.  Cell density u(x, y, t) with PPL for Problem 5.2.1.
    Figure 3.  Projection for u(x, y, t) on the xou plane for Problem 5.2.1.

    (2) Non-negativity. To highlight the performance of positivity-preserving capability for our method in the presence of a PPL, we employ the proposed scheme in the absence and presence of a PPL to solve the IBVP 5.2.1 in two different space grids N = 50 and N = 100 at pre-blow-up times t = 10^{-4} and t = 1.2\times10^{-4} , where \tau_0 = 5\times10^{-6} . In Figure 4, we can see that u(x, y, t) has negative values in the absence of a PPL, and that all negative values are eliminated after the PPL is used. At the same time, the negative values decrease with the increase of spatial grid numbers, and the solutions become smoother. However, for the same number of grids, the oscillation of the blow-up area increases with time.

    Figure 4.  One-dimensional profile of u(x, y, t) along y = 0 for Problem 5.2.1.

    (3) Influence of chemotaxis sensitivity coefficient \chi . We employ the proposed method to compute the blow-up solution in a finite time for this IBVP 5.2.1 under different chemotaxis sensitivity coefficients \chi = 2 and \chi = 3 at the same pre-blow-up time t = 5\times10^{-5} . The results are shown in Figure 5, where \tau_0 = 10^{-6} and N = 100 . We find that, for different values of \chi , the cell density u(x, y, t) achieves negative values in the absence of a PPL. After using the PPL, the algorithm effectively eliminates the negative values. At the same time, the maximum peak value when \chi = 3 is one order of magnitude higher than that when \chi = 2 . It can be seen that the peak value is gradually increased over \chi . Therefore, the sensitivity intensity \chi also affects the value of u(x, y, t) and the blow-up degree of Problem 5.2.1.

    Figure 5.  One-dimensional profile of u(x, y, t) along y = 0 for Problem 5.2.1 with t = 5\times10^{-5} .

    (4) Chemoattractant concentration. In Figure 6, we plot the numerical solutions of chemoattractant concentration v(x, y, t) for Problem 5.2.1 when t = 10^{-4} , \tau_0 = 5\times10^{-6} and N = 100 , and we have compared the computed results in the absence and presence of a PPL. Observing Figures 4b and 6, we find that, although we applied the same parameters in the absence of a PPL, u(x, y, t) achieves negative values, but v(x, y, t) does not have negative values, and that the numerical solutions in the presence of a PPL is highly consistent with those in the absence of a PPL. In summary, our method can effectively eliminate these negative values in this computation. It is consistent with the results obtained in the absence of a PPL if the obtained solutions do not have negative values. Meanwhile, Problem 5.2.1 always experiences a finite-time blow-up at the center of the rectangular region under this initial condition.

    Figure 6.  Chemoattractant concentration v(x, y, t) for Problem 5.2.1 with t = 10^{-4} .

    (5) Accuracy and mass conservation. Since the analytical solution for this IBVP 5.2.1 is very difficult to obtain, to verify the accuracy and mass conservation, we take the fine mesh 800\times800 as the reference solution. The convergence of the HOC scheme in the presence of a PPL are tested by using the relative errors between the maximum value \textbf{U}_{\max} (t) of u(x, y, t) and the reference solutions \textbf{U}_{\max}^*(t) ; the results are shown in Table 9. We find that u(x, y, t) can converge at a fourth-order relative error convergence rate. Meanwhile, we plot u_{max}(t) , which is the time evolution of the maximum value of u(x, y, t) , on the left of Figure 7, and the numerical energy with time on the right of Figure 7, using the proposed scheme in the presence of a PPL to solve this IBVP 5.2.1 with PPL at \tau_0 = 5\times 10^{-6} , N = 50, 80 , 100 . We find that the energy is strictly positive and monotonically decreasing with the evolution of time, which verifies the energy dissipation nature of the original problem. In addition, in Table 10, we give the discrete masses of u(x, y, t) and v(x, y, t) as obtained by using the proposed scheme in the absence and presence of a PPL to solve this IBVP 5.2.1 at different times T with \tau_0 = 5\times10^{-6}, N = 100 . According to this table, we can find that the proposed scheme in the absence of a PPL well verifies the mass conservation. After using the PPL, it is still conserved before blow up, and the mass conservation is slightly affected near the blow-up time.

    Table 9.  Relative ||\cdot||_{\infty} , ||\cdot||_{2} errors and convergence rates of u_{\max}(t) for Problem 5.2.1 with \tau_0 = 10^{-6}, {t = 5\times10^{-5}} .
    N Rel_{{\infty}} Rate Rel_{{2}} Rate
    100 3.204E-02 2.104E-02
    200 1.864E-03 4.10 1.197E-03 4.14
    300 3.290E-04 4.28 2.134E-04 4.25
    400 7.872E-05 4.97 5.902E-05 4.47
    500 3.332E-05 3.85 2.154E-05 4.52

     | Show Table
    DownLoad: CSV
    Table 10.  Mass_u and Mass_v of the HOC scheme for Problem 5.2.1 at time T with \tau_0 = 5\times10^{-6}, N = 100 .
    T Without PPL With PPL
    Mass_u Mass_v Mass_u Mass_v
    0 31.4159265323 31.4156968041 31.4159265323 31.4156968041
    1.0\times10^{-5} 31.4159265336 31.4157188735 31.4159265336 31.4157188735
    2.0\times10^{-5} 31.4159265339 31.4157222761 31.4159265339 31.4157222761
    3.0\times10^{-5} 31.4159265341 31.4157255274 31.4159265341 31.4157255274
    4.0\times10^{-5} 31.4159265344 31.4157287900 31.4159265344 31.4157287900
    5.0\times10^{-5} 31.4159265348 31.4157319661 31.4159265348 31.4157319661
    7.0\times10^{-5} 31.4159265354 31.4157381680 31.4292462953 31.4157381680
    1.0\times10^{-4} 31.4159265366 31.4157472806 37.7315206418 31.4157952866
    1.2\times10^{-4} 31.4159265376 31.4157533114 48.2548641463 31.4159976499

     | Show Table
    DownLoad: CSV
    Figure 7.  Left: Plots of u_{max}(t) : the time evolution of the maximum value for u(x, y, t) ; Right: Plots of the free energy with time for Problem 5.2.1 with the PPL at N = 50, 80 and 100 .

    In the section, we consider the following IBVP [24] for the chemotaxis system given by (1.1)–(1.3) with the boundary conditions given by (1.3) for the rectangular region [-0.5, 0.5]\times[-0.5, 0.5] . Its initial condition is given as follows:

    \begin{align} &u(x, y, 0) = 1000e^{-100[(x-0.15)^2+(y-0.15)^2]}, \; \; \; \; v(x, y, 0) = 0. \end{align}

    In the computation, we take \chi = d = 1 and employ the proposed scheme in the presence of a PPL to simulate the temporal evolution of u(x, y, t) and v(x, y, t) . The solution for this Problem 5.2.2 will blow up in finite time [8,9,24] at the corner of the rectangular region [-0.5, 0.5]\times[-0.5, 0.5] . Next, we take the initial time step size \tau_0 = 10^{-3} and fix a mesh of 200\times 200 to compute them.

    Figures 8 and 9 respectively display the evolutionary processes of u(x, y, t) and v(x, y, t) over time t . It can be seen in these figures that the maximum value of u(x, y, t) gradually moves toward the corner of the rectangular region [-0.5, 0.5]\times[-0.5, 0.5] in the form of a Gaussian profile. When it is close to the corner of the rectangular region, the value of u(x, y, t) increases rapidly, finally blowing up in a finite time. At the same time, v(x, y, t) also increases rapidly at the corner of the rectangular region.

    Figure 8.  Evolution of u(x, y, t) over time t for Problem 5.2.2.
    Figure 9.  Evolution of v(x, y, t) over time t for Problem 5.2.2.

    In this work, first, a fourth-order compact difference scheme for solving a 2D Keller-Segel model was derived and the computing strategies for the initial time steps and the nonlinear chemotaxis terms has been presented. Then, a multigrid algorithm was established to improve the computational efficiency and the binary function integration method has been used to establish the positivity-preserving algorithm. Third, the accuracy and reliability of the proposed method in the absence and presence of a PPL have been respectively verified by several numerical examples with flux source terms. And, we simulated the blow-up phenomena for the chemotaxis system in the center and corner of the rectangular region, respectively. Finally, the energy dissipation and mass conservation were also verified by numerical experiments. In summary, compared with the classical high-order numerical methods for solving the chemotaxis system in the literature, this proposed method has the following advantages:

    ● The proposed HOC scheme is compact in the spatial directions and fully implicit, as no more than nine mesh points are required for the 2D spatial mesh subdomain.

    ● The truncation error of the proposed HOC scheme is O(\tau^4+h^4) , i.e., it has space-time fourth-order accuracy. Thus, we can take large time step sizes to approximate the solution of this chemotaxis system.

    ● The proposed numerical algorithm can effectively filter the negative values in the solution process to ensure that the values of cell density are not negative at all time steps, while maintaining the original fourth-order accuracy without loss.

    ● The computed results show that the proposed method is more accurate than most such schemes reported in the literature.

    Nevertheless, our method also has some disadvantages. First, the new scheme has five layers in time. When it is used to compute the actual problem, three start-up time steps are required, which is slightly more troublesome for programming. Second, the proposed positivity-preserving algorithm would slightly affect its mass conservation when the actual problem blows up. The above deficiencies are also the driving force for our next work. At the same time, it is also possible to extend our approach to solve more complex chemotaxis and haptotaxis systems, which will be incorporated into our next work.

    Support of the study was received from the National Natural Science Foundation of China (12161067, 12001015, 12261067), National Natural Science Foundation of Ningxia, China (2022AAC02023) and National Youth Top-Notch Talent Support Program of Ningxia, China.

    All authors assert no conflict of interest.

    The data that support the findings of this study are available from the authors upon reasonable request.



    [1] S. Shabani-Mashcool, S. A. Marashi, S. Gharaghani, NDDSA: A network- and domain-based method for predicting drug-side effect associations, Inform. Process. Manag., 57 (2020), 102357. https://doi.org/10.1016/j.ipm.2020.102357 doi: 10.1016/j.ipm.2020.102357
    [2] Y. J. Ding, J. J. Tang, F. Guo, Identification of drug-side effect association via multiple information integration with centered kernel alignment, Neurocomputing, 325 (2019), 211–224. https://doi.org/10.1016/j.neucom.2018.10.028 doi: 10.1016/j.neucom.2018.10.028
    [3] A. Lakizadeh, S. M. H. Mir-Ashrafi, Drug repurposing improvement using a novel data integration framework based on the drug side effect, Inform. Med. Unlocked, 23 (2021), 100523. https://doi.org/10.1016/j.imu.2021.100523 doi: 10.1016/j.imu.2021.100523
    [4] E. Pauwels, V. Stoven, Y. Yamanishi, Predicting drug side-effect profiles: a chemical fragment-based approach, BMC Bioinformatics, 12 (2011), 169. https://doi.org/10.1186/1471-2105-12-169 doi: 10.1186/1471-2105-12-169
    [5] S. Jamal, S. Goyal, A. Shanker, A. Grover, Predicting neurological adverse drug reactions based on biological, chemical and phenotypic properties of drugs using machine learning models, Sci. Rep., 7 (2017), 872. https://doi.org/10.1038/s41598-017-00908-z doi: 10.1038/s41598-017-00908-z
    [6] Y. Zheng, H. Peng, S. Ghosh, C. Lan, J. Li, Inverse similarity and reliable negative samples for drug side-effect prediction, BMC Bioinformatics, 19 (2019), 554. https://doi.org/10.1186/s12859-018-2563-x doi: 10.1186/s12859-018-2563-x
    [7] M. Liu, Y. Wu, Y. Chen, J. Sun, Z. Zhao, X. W. Chen, et al., Large-scale prediction of adverse drug reactions using chemical, biological, and phenotypic properties of drugs, J. Am. Med. Inform. Assoc., 19 (2012), e28–35. https://doi.org/10.1136/amiajnl-2011-000699 doi: 10.1136/amiajnl-2011-000699
    [8] S. Dey, H. Luo, A. Fokoue, J. Hu, P. Zhang, Predicting adverse drug reactions through interpretable deep learning framework, BMC Bioinformatics, 19 (2018), 476. https://doi.org/10.1186/s12859-018-2544-0 doi: 10.1186/s12859-018-2544-0
    [9] L. Chen, T. Huang, J. Zhang, M. Y. Zheng, K. Y. Feng, Y. D. Cai, et al., Predicting drugs side effects based on chemical-chemical interactions and protein-chemical interactions, BioMed Res. Int., 2013 (2013), 485034. https://doi.org/10.1155/2013/485034 doi: 10.1155/2013/485034
    [10] W. Zhang, F. Liu, L. Luo, J. Zhang, Predicting drug side effects by multi-label learning and ensemble learning, BMC Bioinformatics, 16 (2015), 365. https://doi.org/10.1186/s12859-015-0774-y doi: 10.1186/s12859-015-0774-y
    [11] N. Atias, R. Sharan, An algorithmic framework for predicting side effects of drugs, J. Comput. Biol., 18 (2011), 207–218. https://doi.org/10.1089/cmb.2010.0255 doi: 10.1089/cmb.2010.0255
    [12] E. Muñoz, V. Novácek, P. Y. Vandenbussche, Facilitating prediction of adverse drug reactions by using knowledge graphs and multi-label learning models, Brief. Bioinform., 20 (2017), 190–202. https://doi.org/10.1093/bib/bbx099 doi: 10.1093/bib/bbx099
    [13] W. Zhang, Y. Chen, S. Tu, F. Liu, Q. Qu, Drug side effect prediction through linear neighborhoods and multiple data source integration, in IEEE International Conference on Bioinformatics and Biomedicine, (2016), 427–434. https://doi.org/10.1109/BIBM.2016.7822555
    [14] E. Munoz, V. Novacek, P. Y. Vandenbussche, Using drug similarities for discovery of possible adverse reactions, AMIA Annu. Symp. Proc., 2016 (2016), 924–933.
    [15] X. Zhao, L. Chen, J. Lu, A similarity-based method for prediction of drug side effects with heterogeneous information, Math. Biosci., 306 (2018), 136–144. https://doi.org/10.1016/j.mbs.2018.09.010 doi: 10.1016/j.mbs.2018.09.010
    [16] H. Liang, L. Chen, X. Zhao, X. Zhang, Prediction of drug side effects with a refined negative sample selection strategy, Comput. Math. Method. M., 2020 (2020), 1573543. https://doi.org/10.1155/2020/1573543 doi: 10.1155/2020/1573543
    [17] X. Zhao, L. Chen, Z. H. Guo, T. Liu, Predicting drug side effects with compact integration of heterogeneous networks, Curr. Bioinform., 14 (2019), 709–720. https://doi.org/10.2174/1574893614666190220114644 doi: 10.2174/1574893614666190220114644
    [18] X. Guo, W. Zhou, Y. Yu, Y. Ding, J. Tang, F. Guo, A novel triple matrix factorization method for detecting drug-side effect association based on kernel target alignment, BioMed Res. Int., 2020 (2020), 4675395. https://doi.org/10.1155/2020/4675395 doi: 10.1155/2020/4675395
    [19] Y. Ding, J. Tang, F. Guo, Identification of drug-side effect association via semi-supervised model and multiple kernel learning, IEEE J. Biomed. Health, 23 (2019), 2619–2632. https://doi.org/10.1109/JBHI.2018.2883834 doi: 10.1109/JBHI.2018.2883834
    [20] H. Tong, C. Faloutsos, J. Pan, Fast random walk with restart and its applications, in Sixth International Conference on Data Mining, (2006), 613–622. https://doi.org/10.1109/ICDM.2006.70
    [21] D. E. Carlin, B. Demchak, D. Pratt, E. Sage, T. Ideker, Network propagation in the cytoscape cyberinfrastructure, PLoS Comput. Biol., 13 (2017), e1005598. https://doi.org/10.1371/journal.pcbi.1005598 doi: 10.1371/journal.pcbi.1005598
    [22] M. Kuhn, M. Campillos, I. Letunic, L. J. Jensen, P. Bork, A side effect resource to capture phenotypic effects of drugs, Mol. Syst. Biol., 6 (2010), 343. https://doi.org/10.1038/msb.2009.98 doi: 10.1038/msb.2009.98
    [23] M. Kuhn, D. Szklarczyk, S. Pletscher-Frankild, T. H. Blicher, C. von Mering, L. J. Jensen, et al., STITCH 4: integration of protein–chemical interactions with user data, Nucleic Acids Res., 42 (2014), D401–D407. https://doi.org/10.1093/nar/gkt1207 doi: 10.1093/nar/gkt1207
    [24] D. Weininger, SMILES, a chemical language and information system. 1. Introduction to methodology and encoding rules, J. Chem. Inf. Comput. Sci., 28 (1988), 31–36. https://doi.org/10.1021/ci00057a005 doi: 10.1021/ci00057a005
    [25] X. Xiao, W. Zhu, B. Liao, J. Xu, C. Gu, B. Ji, et al., BPLLDA: predicting lncRNA-Disease associations based on simple paths with limited lengths in a heterogeneous network, Front. Genet., 9 (2018), 411. https://doi.org/10.3389/fgene.2018.00411 doi: 10.3389/fgene.2018.00411
    [26] W. Ba-Alawi, O. Soufan, M. Essack, P. Kalnis, V. B. Bajic, DASPfind: new efficient method to predict drug-target interactions, J. Cheminformatics, 8 (2016), 15. https://doi.org/10.1186/s13321-016-0128-4 doi: 10.1186/s13321-016-0128-4
    [27] Z. H. You, Z. A. Huang, Z. Zhu, G. Y. Yan, Z. W. Li, Z. Wen, et al., PBMDA: a novel and effective path-based computational model for miRNA-disease association prediction, PLoS Comput. Biol., 13 (2017), e1005455. https://doi.org/10.1371/journal.pcbi.1005455 doi: 10.1371/journal.pcbi.1005455
    [28] J. Gao, B. Hu, L. Chen, A path-based method for identification of protein phenotypic annotations, Curr. Bioinform., 16 (2021), 1214–1222. https://doi.org/10.2174/1574893616666210531100035 doi: 10.2174/1574893616666210531100035
    [29] S. Kohler, S. Bauer, D. Horn, P. N. Robinson, Walking the interactome for prioritization of candidate disease genes, Am. J. Hum. Genet., 82 (2008), 949–958. https://doi.org/10.1016/j.ajhg.2008.02.013 doi: 10.1016/j.ajhg.2008.02.013
    [30] Y.J. Li, J. C. Patra, Genome-wide inferring gene-phenotype relationship by walking on the heterogeneous network, Bioinformatics, 26 (2010), 1219–1224. https://doi.org/10.1093/bioinformatics/btq108 doi: 10.1093/bioinformatics/btq108
    [31] X. Chen, M. X. Liu, G. Y. Yan, Drug-target interaction prediction by random walk on the heterogeneous network, Mol. BioSyst., 8 (2012), 1970–1978. https://doi.org/10.1039/C2MB00002D doi: 10.1039/C2MB00002D
    [32] L. Chen, T. Liu, X. Zhao, Inferring anatomical therapeutic chemical (ATC) class of drugs using shortest path and random walk with restart algorithms, BBA Mol. Basis Dis., 1864 (2017), 2228–2240. https://doi.org/10.1016/j.bbadis.2017.12.019 doi: 10.1016/j.bbadis.2017.12.019
    [33] L. Chen, Y. H. Zhang, Z. Zhang, T. Huang, Y. D. Cai, Inferring novel tumor suppressor genes with a protein-protein interaction network and network diffusion algorithms, Mol. Ther. Methods Clin. Dev., 10 (2018), 57–67. https://doi.org/10.1016/j.omtm.2018.06.007 doi: 10.1016/j.omtm.2018.06.007
    [34] S. Lu, K. Zhao, X. Wang, H. Liu, X. Ainiwaer, Y. Xu, et al., Use of laplacian heat diffusion algorithm to infer novel genes with functions related to uveitis, Front. Genet., 9 (2018), 425. https://doi.org/10.3389/fgene.2018.00425 doi: 10.3389/fgene.2018.00425
    [35] H. Y. Liang, B. Hu, L. Chen, S. Q. Wang, Aorigele, Recognizing novel chemicals/drugs for anatomical therapeutic chemical classes with a heat diffusion algorithm, BBA Mol. Basis. Dis., 1866 (2020), 165910. https://doi.org/10.1016/j.bbadis.2020.165910 doi: 10.1016/j.bbadis.2020.165910
    [36] M. Imanishi, Y. Hori, M. Nagaoka, Y. Sugiura, Design of novel zinc finger proteins: towards artificial control of specific gene expression, Eur. J. Pharm. Sci., 13 (2001), 91–97. https://doi.org/10.1016/S0928-0987(00)00212-8 doi: 10.1016/S0928-0987(00)00212-8
    [37] M. Alirezaei, E. Mordelet, N. Rouach, A. C. Nairn, J. Glowinski, J. Premont, Zinc-induced inhibition of protein synthesis and reduction of connexin-43 expression and intercellular communication in mouse cortical astrocytes, Eur. J. Neurosci., 16 (2002), 1037–1044. https://doi.org/10.1046/j.1460-9568.2002.02180.x doi: 10.1046/j.1460-9568.2002.02180.x
    [38] K. H. Ibs, L. Rink, Zinc-altered immune function, J. Nutr., 133 (2003), 1452s–1456s. https://doi.org/10.1093/jn/133.5.1452S doi: 10.1093/jn/133.5.1452S
    [39] Z. A. Bhutta, R. E. Black, K. H. Brown, J. M. Gardner, S. Gore, A. Hidayat, et al., Prevention of diarrhea and pneumonia by zinc supplementation in children in developing countries: Pooled analysis of randomized controlled trials, J. Pediatr., 135 (1999), 689–697. https://doi.org/10.1016/S0022-3476(99)70086-7 doi: 10.1016/S0022-3476(99)70086-7
    [40] D. E. Roth, S. A. Richard, R. E. Black, Zinc supplementation for the prevention of acute lower respiratory infection in children in developing countries: meta-analysis and meta-regression of randomized trials, Int. J. Epidemiol., 39 (2010), 795–808. https://doi.org/10.1093/ije/dyp391 doi: 10.1093/ije/dyp391
    [41] D. Hulisz, Efficacy of zinc against common cold viruses: an overview, J. Am. Pharm. Assoc., 44 (2004), 594–603. https://doi.org/10.1331/1544-3191.44.5.594.Hulisz doi: 10.1331/1544-3191.44.5.594.Hulisz
    [42] R. O. Suara, J. E. Crowe, Effect of zinc salts on respiratory syncytial virus replication, Antimicrob. Agents Ch., 48 (2004), 783–790. https://doi.org/10.1128/AAC.48.3.783-790.2004 doi: 10.1128/AAC.48.3.783-790.2004
    [43] D. Li, L. Z. Wen, H. Yuan, Observation on clinical efficacy of combined therapy of zinc supplement and jinye baidu granule in treating human cytomegalovirus infection, Zhongguo Zhong xi yi jie he za zhi, 25 (2005), 449–451.
    [44] F. Femiano, F. Gombos, C. Scully, Recurrent herpes labialis: a pilot study of the efficacy of zinc therapy, J. Oral Pathol. Med., 34 (2005), 423–425. https://doi.org/10.1111/j.1600-0714.2005.00327.x doi: 10.1111/j.1600-0714.2005.00327.x
    [45] M. Singh, R. R. Das, Zinc for the common cold, Cochrane Database Syst. Rev., 6 (2013), CD001364. https://doi.org/10.1002/14651858.CD001364.pub4 doi: 10.1002/14651858.CD001364.pub4
    [46] M. Lazzerini, H. Wanzira, Oral zinc for treating diarrhoea in children, Cochrane Database Syst. Rev., 12 (2016), CD005436. https://doi.org/10.1002/14651858.CD005436.pub5 doi: 10.1002/14651858.CD005436.pub5
    [47] F. Sakai, S. Yoshida, S. Endo, H. Tomita, Double-blind, placebo-controlled trial of zinc picolinate for taste disorders, Acta oto-laryngol., 122 (2002), 129–133. https://doi.org/10.1080/00016480260046517 doi: 10.1080/00016480260046517
    [48] A. R. Watson, A. Stuart, F. E. Wells, I. B. Houston, G. M. Addison, Zinc supplementation and its effect on taste acuity in children with chronic renal failure, Hum. Nutr. Clin. Nutr., 37 (1983), 219–225.
    [49] J. Cervantes, A. E. Eber, M. Perper, V. M. Nascimento, K. Nouri, J. E. Keri, The role of zinc in the treatment of acne: A review of the literature, Dermatol. Ther., 31 (2018), e12576. https://doi.org/10.1111/dth.12576 doi: 10.1111/dth.12576
    [50] A. Y. Bedikian, M. Valdivieso, L. K. Heilbrun, R. H. Withers, G. P. Bodey, E. J. Freireich, Glycerol: an alternative to dexamethasone for patients receiving brain irradiation for metastatic disease, South. Med. J., 73 (1980), 1210–1214.
    [51] M. S. Frank, M. C. Nahata, M. D. Hilty, Glycerol: a review of its pharmacology, pharmacokinetics, adverse reactions, and clinical use, Pharmacotherapy, 1 (1981), 147–160. https://doi.org/10.1002/j.1875-9114.1981.tb03562.x doi: 10.1002/j.1875-9114.1981.tb03562.x
    [52] J. Wang, Y. Ren, S. F. Wang, L. D. Kan, L. J. Zhou, H. M. Fang, et al., Comparative efficacy and safety of glycerol versus mannitol in patients with cerebral oedema and elevated intracranial pressure: A systematic review and meta-analysis, J. Clin. Pharm. Ther., 46 (2021), 504–514. https://doi.org/10.1111/jcpt.13314 doi: 10.1111/jcpt.13314
    [53] J. Wang, Y. Ren, L. J. Zhou, L. D. Kan, H. Fan, H. M. Fang, Glycerol Infusion Versus Mannitol for Cerebral Edema: A Systematic Review and Meta-analysis, Clin. Ther., 43 (2021), 637–649. https://doi.org/10.1016/j.clinthera.2021.01.010 doi: 10.1016/j.clinthera.2021.01.010
    [54] E. Righetti, M. G. Celani, T. A. Cantisani, R. Sterzi, G. Boysen, S. Ricci, Glycerol for acute stroke, Cochrane Database Syst. Rev., 2 (2004), CD000096. https://doi.org/10.1002/14651858.CD000096.pub2 doi: 10.1002/14651858.CD000096.pub2
    [55] A. Frei, C. Cottier, P. Wunderlich, E. Lüdin, Glycerol and dextran combined in the therapy of acute stroke. A placebo-controlled, double-blind trial with a planned interim analysis, Stroke, 18 (1987), 373–379. https://doi.org/10.1161/01.STR.18.2.373 doi: 10.1161/01.STR.18.2.373
    [56] E. Lin, Glycerol utilization and its regulation in mammals, Annu. Rev. Biochem., 46 (1977), 765–795. https://doi.org/10.1146/annurev.bi.46.070177.004001 doi: 10.1146/annurev.bi.46.070177.004001
    [57] Y. Yu, C. Kumana, I. Lauder, Y. Cheung, F. Chan, M. Kou, et al., Treatment of acute cortical infarct with intravenous glycerol. A double-blind, placebo-controlled randomized trial, Stroke, 24 (1993), 1119–1124. https://doi.org/10.1161/01.STR.24.8.1119 doi: 10.1161/01.STR.24.8.1119
    [58] B. á Rogvi-Hansen, G. Boysen, Intravenous Glycerol Treatment of Acute Stroke – A Statistical Review, Cerebrovasc. Dis., 2 (1992), 11–13. https://doi.org/10.1159/000108981 doi: 10.1159/000108981
    [59] H. L. Philpott, S. Nandurkar, J. Lubel, P. R. Gibson, Drug-induced gastrointestinal disorders, Frontline Gastroente., 5 (2014), 49–57. http://dx.doi.org/10.1136/flgastro-2013-100316 doi: 10.1136/flgastro-2013-100316
    [60] S. Saleem, How to induce arrhythmias with dopamine, in Arrhythmia Induction in the EP Lab, Springer, (2019), 81–89. https://doi.org/10.1007/978-3-319-92729-9_9
    [61] R. Ceravolo, C. Rossi, E. Del Prete, U. Bonuccelli, A review of adverse events linked to dopamine agonists in the treatment of Parkinson's disease, Expert Opin. Drug Saf., 15 (2016), 181–198. https://doi.org/10.1517/14740338.2016.1130128 doi: 10.1517/14740338.2016.1130128
  • mbe-19-06-269-supplementary.pdf
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2652) PDF downloads(99) Cited by(6)

Figures and Tables

Figures(6)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog