Processing math: 58%
Research article

Global behavior of a multi-group SEIR epidemic model with age structure and spatial diffusion

  • Different epidemic models with one or two characteristics of multi-group, age structure and spatial diffusion have been proposed, but few models take all three into consideration. In this paper, a novel multi-group SEIR epidemic model with both age structure and spatial diffusion is constructed for the first time ever to study the transmission dynamics of infectious diseases. We first analytically study the positivity, boundedness, existence and uniqueness of solution and the existence of compact global attractor of the associated solution semiflow. Based on some assumptions for parameters, we then show that the disease-free steady state is globally asymptotically stable by utilizing appropriate Lyapunov functionals and the LaSalle's invariance principle. By means of Perron-Frobenius theorem and graph-theoretical results, the existence and global stability of endemic steady state are ensured under appropriate conditions. Finally, feasibility of main theoretical results is showed with the aid of numerical examples for model with two groups which is important from the viewpoint of applications.

    Citation: Pengyan Liu, Hong-Xu Li. Global behavior of a multi-group SEIR epidemic model with age structure and spatial diffusion[J]. Mathematical Biosciences and Engineering, 2020, 17(6): 7248-7273. doi: 10.3934/mbe.2020372

    Related Papers:

    [1] Bertrand Haut, Georges Bastin . A second order model of road junctions in fluid models of traffic networks. Networks and Heterogeneous Media, 2007, 2(2): 227-253. doi: 10.3934/nhm.2007.2.227
    [2] Mohamed Benyahia, Massimiliano D. Rosini . A macroscopic traffic model with phase transitions and local point constraints on the flow. Networks and Heterogeneous Media, 2017, 12(2): 297-317. doi: 10.3934/nhm.2017013
    [3] Caterina Balzotti, Maya Briani . Estimate of traffic emissions through multiscale second order models with heterogeneous data. Networks and Heterogeneous Media, 2022, 17(6): 863-892. doi: 10.3934/nhm.2022030
    [4] Emiliano Cristiani, Smita Sahu . On the micro-to-macro limit for first-order traffic flow models on networks. Networks and Heterogeneous Media, 2016, 11(3): 395-413. doi: 10.3934/nhm.2016002
    [5] Maya Briani, Emiliano Cristiani . An easy-to-use algorithm for simulating traffic flow on networks: Theoretical study. Networks and Heterogeneous Media, 2014, 9(3): 519-552. doi: 10.3934/nhm.2014.9.519
    [6] Alberto Bressan, Khai T. Nguyen . Conservation law models for traffic flow on a network of roads. Networks and Heterogeneous Media, 2015, 10(2): 255-293. doi: 10.3934/nhm.2015.10.255
    [7] Simone Göttlich, Camill Harter . A weakly coupled model of differential equations for thief tracking. Networks and Heterogeneous Media, 2016, 11(3): 447-469. doi: 10.3934/nhm.2016004
    [8] Paola Goatin . Traffic flow models with phase transitions on road networks. Networks and Heterogeneous Media, 2009, 4(2): 287-301. doi: 10.3934/nhm.2009.4.287
    [9] Michael Herty, Adrian Fazekas, Giuseppe Visconti . A two-dimensional data-driven model for traffic flow on highways. Networks and Heterogeneous Media, 2018, 13(2): 217-240. doi: 10.3934/nhm.2018010
    [10] Alberto Bressan, Ke Han . Existence of optima and equilibria for traffic flow on networks. Networks and Heterogeneous Media, 2013, 8(3): 627-648. doi: 10.3934/nhm.2013.8.627
  • Different epidemic models with one or two characteristics of multi-group, age structure and spatial diffusion have been proposed, but few models take all three into consideration. In this paper, a novel multi-group SEIR epidemic model with both age structure and spatial diffusion is constructed for the first time ever to study the transmission dynamics of infectious diseases. We first analytically study the positivity, boundedness, existence and uniqueness of solution and the existence of compact global attractor of the associated solution semiflow. Based on some assumptions for parameters, we then show that the disease-free steady state is globally asymptotically stable by utilizing appropriate Lyapunov functionals and the LaSalle's invariance principle. By means of Perron-Frobenius theorem and graph-theoretical results, the existence and global stability of endemic steady state are ensured under appropriate conditions. Finally, feasibility of main theoretical results is showed with the aid of numerical examples for model with two groups which is important from the viewpoint of applications.


    A taxis is the movement of an organism in response to a stimulus such as chemical signal or the presence of food. Taxes can be classified based on the types of stimulus, such as chemotaxis, prey-taxis, galvanotaxis, phototaxis and so on. According to the direction of movements, the taxis is said to be attractive (resp. repulsive) if the organism moves towards (resp. away from) the stimulus. In the ecosystem, a widespread phenomenon is the prey-taxis, where predators move up the prey density gradient, which is often referred to as the direct prey-taxis. However some predators may approach the prey by tracking the chemical signals released by the prey, such as the smell of blood or specific odo, and such movement is called indirect prey-taxis (cf. [1]). Since the pioneering modeling work by Kareiva and Odell [2], prey-taxis models have been widely studied in recent years (cf. [3,4,5,6,7,8,9,10,11,12]), followed by numerous extensions, such as three-species prey-taxis models (cf. [13,14,15]) and predator-taxis models (cf. [16,17]). The indirect prey-taxis models have also been well studied (cf. [18,19,20]).

    Recently, a predator-prey model with attraction-repulsion taxis mechanisms was proposed by Bell and Haskell in [21] to describe the interaction between direct prey-taxis and indirect chemotaxis, where the direct prey-taxis describes the predator's directional movement towards the prey density gradient, while the indirect chemotaxis models a defense mechanism in which the prey repels the predator by releasing odour chemicals (like a fox breaking wind in order to escape from hunting dogs). The model reads as

    {ut=dΔu+u(a1a2ua3v),xΩ, t>0,vt=(v+χvwξvu)+ρv(1v)+ea3uv,xΩ, t>0,wt=ηΔw+ruγw,xΩ, t>0,uν=vν=wν=0,xΩ,t>0,(u,v,w)(x,0)=(u0,v0,w0)(x),xΩ, (1.1)

    where the unknown functions u(x,t), v(x,t) and w(x,t) denote the densities of the prey, predator and prey-derived chemical repellent, respectively, at position xΩ and time t>0. Here, ΩRn is a bounded domain (habitat of species) with smooth boundary Ω, and ν is the unit outer normal vector of Ω. The parameters d, η, χ, ξ, a1, a2, a3, e, ρ, r, γ are all positive, where χ>0 and ξ>0 denote the (attractive) prey-taxis and (repulsive) chemotaxis coefficients, respectively. The predator v is assumed to be a generalist, so that it has a logistic growth term ρv(1v) with intrinsic growth rate ρ>0. More modeling details with biological interpretations are referred to in [21]. We remark that the predator-prey model with attraction-repulsion taxes has some similar structures to the so-called attraction-repulsion chemotaxis model proposed originally in [22], where the species elicit both attractive and repulsive chemicals (see [23,24,25,26] and references therein for some mathematical studies).

    The initial data satisfy the following conditions:

    v0C0(¯Ω),u0,w0W1,(Ω), and u0, v0, w00 in ¯Ω. (1.2)

    In [21], the global existence of strong solutions to (1.1) was established in one dimension (n=1), and the existence of nontrivial steady state solutions alongside pattern formation was studied by the bifurcation theory. The main purpose of this paper is to study the global dynamics of (1.1) in higher dimensional spaces, which are usually more physical in the real world. Specifically, we shall show the existence of global classical solutions in all dimensions and explore the global stability of constant steady states, by which we may see how parameter values play roles in determining these dynamical properties of solutions.

    The first main result is concerned with the global existence and boundedness of solutions to (1.1). For the convenience of presentation, we let

    K1=max{a1a2,u0L(Ω)},  K2=max{a1K1+a2K21,a3K1} (1.3)

    and

    K3(z)=23z12zdz(n+2(z1)K22z+1)z+12((z1)(4z2+n)K21)z12+23zz2d1z((z1)ξ2z+1)z+1z((4z2+n)K21)1z. (1.4)

    Then, the result on the global boundedness of solutions to (1.1) is stated as follows.

    Theorem 1.1 (Global existence). Let ΩRn(n1) be a bounded domain with smooth boundary and parameters d, η, χ, ξ, a1, a2, a3, e, ρ, r, γ be positive. If

    ρ{>0,n2,2K3([n2]+1)[n2]+1,n>2,

    where K3(p) is defined in (1.4), then for any initial data (u0,v0,w0) satisfying (1.2), the system (1.1) admits a unique classicalsolution (u,v,w) satisfying

    u, v,wC0(¯Ω×[0,+))C2,1(¯Ω×(0,+)),

    and u,v,w>0 in Ω×(0,+). Moreover, there exists a constant C>0 independent of t such that

    u(,t)W1,(Ω)+v(,t)L(Ω)+w(,t)W1,(Ω)Cfor all t>0.

    Our next goal is to explore the large-time behavior of solutions to (1.1). Simple calculations show the system (1.1) has four possible homogeneous equilibria as classified below:

    {(0,0,0), (0,1,0), (a1a2,0,ra1γa2),if a1a3,(0,0,0), (0,1,0), (a1a2,0,ra1γa2),(u,v,w),if a1>a3,

    with

    u=ρ(a1a3)ρa2+ea23,v=ea1a3+ρa2ρa2+ea23,w=rρ(a1a3)γ(ρa2+ea23) (1.5)

    where the trivial equilibrium (0,0,0) is called the extinction steady state, (0,1,0) is the predator-only steady state, and (u,v,w) is the coexistence steady state. We shall show that if a1>a3, then the coexistence steady state is globally asymptotically stable with exponential convergence rate, provided that ξ and χ are suitably small, while if a1a3, the predator-only steady state is globally asymptotically stable with exponential or algebraic convergence rate when ξ and χ are suitably small. To state our results, we denote

    Γ=4dρ(a1a3)K21(ea1a3+ρa2),Φ=2a2ρa23+e,Ψ=γηa23K21(ρa2+ea23)dρ2r2(a1a3) (1.6)

    and

    A=ξ24d,B=ea2a1,D=16ηγa1r2, (1.7)

    where K1 is defined in (1.3). Then, the global stability result is stated in the following theorem.

    Theorem 1.2 (Global stability). Let the assumptions in Theorem 1.1 hold. Then, the following results hold.

    (1) Let a1>a3. If ξ and χ satisfy

    ξ2<Γ(Φ+Φ2e2) and χ2<Ψmaxy[a,b](Γyξ2)(y2+2Φye2)y,

    where a=max{ξ2Γ,ΦΦ2e2},b=Φ+Φ2e2, then there exist some constants T, C, α>0 such that the solution (u,v,w) obtained in Theorem 1.1 satisfies for all tT

    u(,t)uL(Ω)+v(,t)vL(Ω)+w(,t)wL(Ω)Ceαt.

    (2) Let a1a3, If ξ and χ satisfy

    ξ2<4dea2a1andχ2<D(A+B2AB),

    then there exist some constants T, C, β>0 such that the solution (u,v,w) obtained in Theorem 1.1 satisfies, for all tT,

    u(,t)L(Ω)+v(,t)1L(Ω)+w(,t)L(Ω){Ceβt if a1<a3,C(t+1)1 if a1=a3.

    Remark 1.1. In the biological view, the relative sizes of a1 and a2 determine the coexistence of the system. The results indicated that a large a1a2 facilitates the coexistence of the species.

    The rest of this paper is organized as follows. In Section 2, we state the local existence of solutions to (1.1) with extensibility conditions. Then, we deduce some a priori estimates and prove Theorem 1.1 in Section 3. Finally, we show the global convergence to the constant steady states and prove Theorem 1.2 in Section 4.

    For convenience, in what follows we shall use Ci(i=1,2,) to denote a generic positive constant which may vary from line to line. For simplicity, we abbreviate t0Ωf(,s)dxds and Ωf(,t)dx as t0Ωf and Ωf, respectively. The local existence and extensibility result of problem (1.1) can be directly established by the well-known Amman's theory for triangular parabolic systems (cf. [27,28]). Below, we shall present the local existence theorem without proof for brevity, and we refer to [21] for the proof in one dimension as a reference.

    Lemma 2.1 (Local existence and extensibility). Let ΩRn be a bounded domain with smooth boundary. The parameters d, η, χ, ξ, a1, a2, a3, e, ρ, r, γ are positive. Then, for the initial data (u0,v0,w0) satisfying (1.2), there exists Tmax(0,] such that the system (1.1) admits a unique classicalsolution (u,v,w) satisfying

    u, v, wC0(¯Ω×[0,Tmax))C2,1(¯Ω×(0,Tmax)),

    and u,v,w>0 in Ω×(0,Tmax). Moreover, we have

    either Tmax=+ or lim suptTmax(u(,t)W1,(Ω)+v(,t)L(Ω)+w(,t)W1,(Ω))=+.

    We recall some well-known results which will be used later frequently. The first one is an uniform Grönwall inequality [29].

    Lemma 2.2. Let Tmax>0, τ(0,Tmax). Suppose that c1, c2, y are three positive locally integrable functions on (0,Tmax) such that y is locally integrable on (0,Tmax) and satisfies

    y(t)c1(t)y(t)+c2(t)for all t(0,Tmax).

    If

    t+τtc1C1,t+τtc2C2,  t+τtyC3for all t[0,Tmaxτ),

    where Ci(i=1,2,3) are positive constants, then

    y(t)(C3τ+C2)eC1for all t[τ,Tmax).

    Next, we recall a basic inequality [30].

    Lemma 2.3. Let p[1,). Then, the following inequality holds:

    Ω|u|2(p+1)2(4p2+n)u2L(Ω)Ω|u|2(p1)|D2u|2

    for any uC2(ˉΩ) satisfying uν=0 on Ω, where D2u denotes the Hessian of u.

    The last one is a Gagliardo-Nirenberg type inequality shown in [31,Lemma 2.5].

    Lemma 2.4. Let Ω be a bounded domain in R2 with smooth boundary. Then, for any φW2,2(Ω) satisfying φν|Ω=0, there exists a positive constant C depending only on Ω such that

    φL4(Ω)C(Δφ12L2(Ω)φ12L2(Ω)+φL2(Ω)). (2.1)

    In this section, we establish the global boundedness of solutions to the system (1.1). To this end, we will proceed with several steps to derive a priori estimates for the solution of the system (1.1). The first one is the uniform-in-time L(Ω) boundedness of u.

    Lemma 3.1. Let (u,v,w) be the solution of (1.1) and K1 be as defined in (1.3). Then, we have

    uL(Ω)K1for all t(0,Tmax).

    Furthermore, there is a constant C>0 such that for any 0<τ<min{Tmax,1}, it follows that

    t+τt|u|2C  for all  t(0,Tmaxτ).

    Proof. The result is a direct consequence of the maximum principle applied to the first equation in (1.1). Indeed, if we let ˉu=max{a1a2,u0L(Ω)}, then ˉu satisfies

    {ˉutdΔˉu+ˉu(a1a2ˉua3v),xΩ,t>0,ˉuν=0,xΩ,t>0,ˉu(x,0)u0(x),xΩ.

    Apparently, the comparison principle of parabolic equations gives uˉu on Ω×(0,Tmax).

    Next, we multiply the first equation of (1.1) by u and integrate the result to get

    ddtΩu2+dΩ|u|2=a1Ωu2Ωu(a2u+a3v)a1K21|Ω|.

    Then, the integration of the above inequality with respect to t over (t,t+τ) completes the proof by noting that Ωu20 is bounded.

    Having at hand the uniform-in-time L(Ω) boundedness of u, the a priori estimate of w follows immediately.

    Lemma 3.2. Let (u,v,w) be the solution of (1.1). We can find a constant C>0 satisfying

    wW1,(Ω)Cfor all t(0,Tmax).

    Proof. Noting the boundedness of uL(Ω) from Lemma 3.1, we get the desired result from the third equation of (1.1) and the regularity theorem [32,Lemma 1].

    Now, the a priori estimate of v can be obtained as below.

    Lemma 3.3. Let (u,v,w) be the solution of (1.1). Then, there exists a constant C>0 such that

    ΩvCfor all t(0,Tmax), (3.1)

    and

    t+τtΩv2Cfor all t(0,Tmaxτ), (3.2)

    where τ is a constant such that 0<τ<min{Tmax,1}.

    Proof. Integrating the second equation of (1.1) over Ω by parts, using Young's inequality and Lemma 3.1, we find some constant C1>0 such that

    ddtΩv=ρΩvρΩv2+ea3Ωuv(ρ+ea3supt(0,Tmax)uL(Ω))ΩvρΩv2Ωvρ2Ωv2+C1for all t(0,Tmax). (3.3)

    Hence, (3.1) is obtained by the Grönwall inequality. Integrating (3.3) over (t,t+τ), we get (3.2) immediately.

    Due to the estimates of u and v obtained in Lemmas 3.1 and 3.3 respectively, we have the following improved uniform-in-time L2(Ω) boundedness of u and the space-time L2 boundedness of Δu when n=2.

    Lemma 3.4. Let (u,v,w) be the solution of (1.1). If n=2, then we can find a constant C>0 such that

    Ω|u|2Cfor all t(0,Tmax) (3.4)

    and

    t+τtΩ|Δu|2Cfor all t(0,Tmaxτ), (3.5)

    where τ is defined in Lemma 3.3.

    Proof. Integrating the first equation of (1.1) by parts and using Lemma 3.1, we find a constant C1>0 such that

    ddtΩ|u|2=2Ωuut=2ΩutΔu=2ΩΔu(dΔu+a1ua2u2a3uv)2dΩ|Δu|2+C1Ω(v+1)|Δu|for all t(0,Tmax). (3.6)

    The Gagliardo-Nirenberg inequality in Lemma 2.4, Young's inequality and Lemma 3.1 yield some constants C2,C3>0 satisfying

    Ω|u|2=u2L2(Ω)C2(ΔuL2(Ω)uL2(Ω)+u2L(Ω))d2Ω|Δu|2+C3

    and

    C1Ω(v+1)|Δu|d2Ω|Δu|2+C3Ωv2+C3for all t(0,Tmax),

    which along with (3.6) imply

    ddtΩ|u|2+Ω|u|2+dΩ|Δu|2C3Ωv2+2C3for all t(0,Tmax). (3.7)

    Then, applications of Lemma 2.2, 3.1 and 3.3 give (3.4). Finally, (3.5) can be obtained by integrating (3.7) over (t,t+τ).

    Now, the uniform-in-time boundedness of v in L2(Ω) can be established when n=2.

    Lemma 3.5. Let (u,v,w) be the solution of (1.1). If n=2, then there exists a constant C>0 such that

    Ωv2Cfor all t(0,Tmax).

    Proof. Multiplying the second equation of (1.1) by v, integrating the result by parts and using Young's inequality, we have

    ddtΩv2+2Ω|v|2=2χΩvvw+2ξΩvuv+2ρΩv22ρΩv3+2ea3Ωuv2Ω|v|2+2χ2w2L(Ω)Ωv2+2ξ2Ωv2|u|2+2ρΩv22ρΩv3+2ea3uL(Ω)Ωv2,

    which along with Lemma 3.1 and Lemma 3.2 gives some constant C1>0 such that

    ddtΩv2+Ω|v|22ξ2Ωv2|u|2+C1Ωv22ρΩv3for all t(0,Tmax). (3.8)

    Using Lemmas 3.1 and 3.3, Hölder's inequality, Lemma 2.4 and Young's inequality, we find some constants C2,C3,C4>0 such that

    2ξΩv2|u|22ξv2L4(Ω)u2L4(Ω)C2(v12L2(Ω)v12L2(Ω)+vL2(Ω))2(Δu12L2(Ω)u12L(Ω)+uL(Ω))2C3(vL2(Ω)vL2(Ω)ΔuL2(Ω)+vL2(Ω)vL2(Ω)+ΔuL2(Ω)v2L2(Ω)+v2L2(Ω))v2L2(Ω)+C4(1+Δu2L2(Ω))v2L2(Ω)for all t(0,Tmax). (3.9)

    Furthermore, Young's inequality yields some constant C5>0 such that

    C1Ωv22ρΩv3C5for all t(0,Tmax). (3.10)

    Substituting (3.9) and (3.10) into (3.8), we get

    ddtΩv2C4(1+Δu2L2(Ω))v2L2(Ω)+C5for all t(0,Tmax),

    which alongside Lemma 2.2, Lemma 3.3 and Lemma 3.4 completes the proof.

    To get the global existence of solutions in any dimensions, we derive the following functional inequality which gives an a priori estimate on u.

    Lemma 3.6. Let (u,v,w) be the solution of (1.1) and q2. If n1, then there exists a constant C>0 such that

    ddtΩ|u|2q+dqΩ|u|2(q1)|D2u|2q(n+2(q1))K22dΩ(v2+1)|u|2(q1)+Cfor all t(0,Tmax),

    where K2 is defined in (1.3).

    Proof. From the first equation of (1.1) and the fact 2uΔu=Δ|u|22|D2u|2, it follows that

    ddtΩ|u|2q=2qΩ|u|2(q1)uut=2qΩ|u|2(q1)u(dΔu+a1ua2u2a3uv)=dqΩ|u|2(q1)Δ|u|22dqΩ|u|2(q1)|D2u|2+2qΩ|u|2(q1)u(a1ua2u2a3uv)

    which implies

    ddtΩ|u|2q+2dqΩ|u|2(q1)|D2u|2=dqΩ|u|2(q1)Δ|u|2+2qΩ|u|2(q1)u(a1ua2u2a3uv)=:I1+I2for all t(0,Tmax). (3.11)

    Now, we estimate the right hand side of (3.11). Choosing s(0,12) and

    θ=12s+12nq121nq(0,1),

    we get

    12s+12n=θ(121n)+(1θ)q,

    which, along with the Gagliardo-Nirenberg inequality, Young's inequality and the embedding

    Ws+12,2(Ω)Ws,2(Ω)L2(Ω),

    gives some constants C1, C2, C3, C4>0 such that

    Ω|u|2(q1)|u|2νC1Ω|u|2q=C1|u|q2L2(Ω)C2|u|q2Ws+12,2(Ω)C3|u|q2θL2(Ω)|u|q2(1θ)L1q(Ω)+C3|u|q2L1q(Ω)2(q1)q2|u|q2L2(Ω)+C4for all t(0,Tmax).

    Therefore, it holds that

    I1=dqΩ|u|2(q1)|u|2νdqΩ|u|2(q1)|u|22d(q1)qΩ||u|q|2+C4dq4d(q1)qΩ||u|q|22d(q1)qΩ||u|q|2+C4dqfor all t(0,Tmax).

    Owning to the fact |Δu|n|D2u|, Young's inequality and Lemma 3.1, we have

    I2=2q(q1)Ω(a1ua2u2a3uv)|u|2(q2)|u|2u2qΩ(a1ua2u2a3uv)|u|2(q1)Δu2q(q1)K2Ω(v+1)|u|2(q2)||u|2||u|+2qnK2Ω(v+1)|u|2(q1)|D2u|qd(q1)2Ω|u|2(q2)||u|2|2+2q(q1)K22dΩ(v2+1)|u|2(q1)+dqΩ|u|2(q1)|D2u|2+qnK22dΩ(v2+1)|u|2(q1)=2d(q1)qΩ||u|q|2+dqΩ|u|2(q1)|D2u|2+q(n+2(q1))K22dΩ(v2+1)|u|2(q1)for all t(0,Tmax),

    where K2 is defined in (1.3). Hence, substituting the estimates I1 and I2 into (3.11), we finish the proof of the lemma.

    Now, we show the following functional inequality to derive the a priori estimate on v in any dimensions.

    Lemma 3.7. Let (u,v,w) be the solution of (1.1) and q2. If n1, we can find a constant C>0 such that

    ddtΩvq+2(q1)qΩ|vq2|2+ρqΩvq+1q(q1)ξ2Ωvq|u|2+CΩvq

    for all t(0,Tmax).

    Proof. Utilizing the second equation of (1.1) and integration by parts, we get

    ddtΩvq=qΩvq1vt=qΩvq1((v+χvwξvu)+v(ρ(1v)+ea3u))=q(q1)Ωvq2|v|2χq(q1)Ωvq1wv+ξq(q1)Ωvq1uv+ρqΩvqρqΩvq+1+ea3qΩuvq. (3.12)

    Now, we estimate the right hand side of (3.12). An application of Young's inequality and Lemma 3.2 yields some constant C1>0 such that

    χq(q1)Ωvq1wvχq(q1)supt(0,Tmax)wL(Ω)Ωvq1|v|q(q1)4Ωvq2|v|2+C1Ωvq

    and

    ξq(q1)Ωvq1uvq(q1)4Ωvq2|v|2+q(q1)ξ2Ωvq|u|2,

    which along with (3.12), Lemma 3.1 and the fact

    vq2|v|2=4q2|vq2|2

    gives a constant C2>0 such that

    ddtΩvq+2(q1)qΩ|vq2|2q(q1)ξ2Ωvq|u|2+(ρq+C1)ΩvqρqΩvq+1+ea3qΩuvqq(q1)ξ2Ωvq|u|2ρqΩvq+1+C2Ωvqfor all t(0,Tmax).

    Hence, we finish the proof of the lemma.

    Combining Lemmas 3.6 and 3.7, we have the following inequality which can help us to achieve the global existence of solutions in any dimensions.

    Lemma 3.8. Let (u,v,w) be the solution of (1.1) and p2. If n1, we can find a constant C>0 such that

    ddt(Ω|u|2p+Ωvp)+2(p1)pΩ|vp2|2+Ω|u|2p+Ωvp(K3(p)ρp2)Ωvp+1+Cfor all t(0,Tmax),

    where K3(p) is defined in (1.4).

    Proof. Combining Lemmas 3.6 and 3.7, we see for any p=q2 there exists a constant C1>0 such that for all t(0,Tmax)

    ddt(Ω|u|2p+Ωvp)+2(p1)pΩ|vp2|2+dpΩ|u|2(p1)|D2u|2+ρpΩvp+1p(n+2(p1))K22dΩv2|u|2(p1)+p(p1)ξ2Ωvp|u|2+C1Ω|u|2(p1)+C1Ωvp+C1. (3.13)

    Now, we estimate the right hand side of the above inequality. Indeed, owing to Lemma 2.3 and Young's inequality, for all t(0,Tmax), we have

    p(n+2(p1))K22dΩv2|u|2(p1)dp8(4p2+n)u2L(Ω)Ω|u|2(p+1)+2p+1(dp(p+1)8(p1)(4p2+n)u2L(Ω))p12(p(n+2(p1))K22d)p+12Ωvp+1dp4Ω|u|2(p1)|D2u|2+23p12pdp(n+2(p1)K22p+1)p+12((p1)(4p2+n)K21)p12Ωvp+1

    and

    p(p1)ξ2Ωvp|u|2dp8(4p2+n)u2L(Ω)Ω|u|2(p+1)+pp+1(dp(p+1)8(4p2+n)u2L(Ω))1p(p(p1)ξ2)p+1pΩvp+1dp4Ω|u|2(p1)|D2u|2+23pp2d1p((p1)ξ2p+1)p+1p((4p2+n)K21)1pΩvp+1,

    where K1 and K2 are defined in (1.3). Similarly, we can find a constant C2>0 such that

    C1Ω|u|2(p1)dp8(4p2+n)u2L(Ω)Ω|u|2(p+1)+C2dp4Ω|u|2(p1)|D2u|2+C2for all t(0,Tmax).

    Substituting the above estimates into (3.13), we get

    ddt(Ω|u|2p+Ωvp)+2(p1)pΩ|vp2|2+dp4Ω|u|2(p1)|D2u|2+ρpΩvp+1K3(p)Ωvp+1+C1Ωvp+C1+C2for all t(0,Tmax), (3.14)

    where K3(p) is given in (1.4). Furthermore, we can use Young's inequality and Lemma 2.3 to get a constant C3>0 such that

    (C1+1)Ωvpρp2Ωvp+1+C3,

    and

    Ω|u|2pdp8(4p2+n)u2L(Ω)Ω|u|2(p+1)+C3dp4Ω|u|2(p1)|D2u|2+C3for all t(0,Tmax),

    which together with (3.14) finishes the proof.

    Next, we shall deduce a criterion of global boundedness of solutions for the system (1.1) inspired by an idea of [33].

    Lemma 3.9. Let n1. If there exist M>0 and p0>n2 such that

    Ωvp0Mfor all t(0,Tmax), (3.15)

    then Tmax=+. Moreover, there exists C>0 such that

    u(,t)W1,(Ω)+v(,t)L(Ω)+w(,t)W1,(Ω)Cfor all t>0.

    Proof. We divide the proof into two steps.

    Step 1: We claim that there exists a constant C1>0 such that

    Ωv2p0C1for all t(0,Tmax).

    Indeed, due to Lemma 3.8, for any p=2p0, there exists a constant C2>0 such that

    ddt(Ω|u|4p0+Ωv2p0)+2p01p0Ω|vp0|2+Ω|u|4p0+Ωv2p0(K3(2p0)ρp0)Ωv2p0+1+C2for all t(0,Tmax). (3.16)

    Let

    θ=nn+22p0+22p0+1(0,1).

    Then, 2p0+12p0θ<1 due to p0>n2. By the Gagliardo-Nirenberg inequality, Young's inequality and (3.15), we can find some constants C3,C4>0 such that

    (K3(2p0)ρp0)Ωv2p0+1=(K3(2p0)ρp0)vp02p0+1p0L2p0+1p0(Ω)C3(vp02p0+1p0(1θ)L1(Ω)vp02p0+1p0θL2(Ω)+vp02p0+1p0L1(Ω))C3(M2p0+1p0(1θ)vp02p0+1p0θL2(Ω)+M2p0+1p0)2p01p0Ω|vp0|2+C4for all t(0,Tmax),

    which along with (3.16) implies

    \frac{d}{dt}\left(\int_{\Omega}|\nabla u|^{4p_0}+\int_{\Omega}v^{2p_0}\right)+\int_{\Omega}|\nabla u|^{4p_0}+\int_{\Omega} v^{2p_0} \leqslant C_2+C_4\quad\text{for all}\ t\in(0, T_{max}).

    Therefore, the claim follows from the Grönwall inequality applied to the above inequality.

    Step 2: Thanks to the regularity theorem [32,Lemma 1], we can find a constant C_5 > 0 such that \|\nabla u\|_{L^\infty(\Omega)}\leqslant C_5 due to 2p_0 > n . With (3.12) and Lemmas 3.1 and 3.2, we get a constant C_6 > 0 such that for any p\geqslant2

    \begin{equation} \frac{d}{dt}\int_\Omega v^p+p(p-1)\int_\Omega v^{p-2}|\nabla v|^2\leqslant p(p-1)(C_6\chi+C_5\xi)\int_\Omega v^{p-1}|\nabla v|+p(\rho+ea_3K_1)\int_\Omega v^p. \end{equation} (3.17)

    Thanks to Young's inequality, we find a constant C_7 > 0 such that

    \begin{equation*} p(p-1)(C_6\chi+C_5\xi)\int_\Omega v^{p-1}|\nabla v|\leqslant\frac{p(p-1)}{2}\int_\Omega v^{p-2}|\nabla v|^2+C_7p(p-1)\int_\Omega v^p, \end{equation*}

    which together with (3.17) implies

    \begin{equation} \frac{d}{dt}\int_\Omega v^p+p(p-1)\int_\Omega v^p+\frac{2(p-1)}{p}\int_\Omega|\nabla v^{\frac{p}{2}}|^2\leqslant p(p-1)C_8\int_\Omega v^p, \end{equation} (3.18)

    with C_8 = C_7+\rho+ea_3K_1+1 . Applying 1+p^n\leqslant(1+p)^n and the following inequality [34]

    \|f\|_{L^2}^2\leqslant\varepsilon\|\nabla f\|_{L^2}^2+C_9(1+\varepsilon^{-\frac{n}{2}})\|f\|_{L^1}^2,

    with f = v^{\frac{p}{2}} and \varepsilon = \frac{2}{p^2C_8} , we find a constant C_{10} > 0 such that

    \begin{equation} p(p-1)C_8\int_\Omega u^p\leqslant\frac{2(p-1)}{p}\int_\Omega \left|\nabla u^{\frac{p}{2}}\right|^2+C_{10}p(p-1)(1+p^n)\left(\int_\Omega u^{\frac{p}{2}}\right)^2. \end{equation} (3.19)

    Substituting (3.19) into (3.18), we have

    \begin{equation*} \frac{d}{dt}\int_\Omega u^p+p(p-1)\int_\Omega u^p\leqslant C_{10}p(p-1)(1+p)^n\left(\int_\Omega u^{\frac{p}{2}}\right)^2. \end{equation*}

    Then, employing the standard Moser iteration in [35] or a similar argument as in [34], we can prove that there exists a constant C_{11} > 0 such that

    \begin{equation*} \|v\|_{L^\infty(\Omega)}\leqslant C_{11}\quad\text{for all}\ t\in(0, T_{max}). \end{equation*}

    Thus, with the help of Lemma 3.2, we finish the proof.

    Now, utilizing the criterion in Lemma 3.9, we prove the global existence and boundedness of solutions for the system (1.1).

    Proof of Theorem 1.1. If n\leqslant2 , then the conclusion of the theorem can be obtained by Lemmas 3.3, 3.5 and 3.9. If n\geqslant3 and

    \rho\geqslant \frac{2K_3\left(\left[\frac{n}{2}\right]+1\right)}{\left[\frac{n}{2}\right]+1},

    then according to Lemma 3.8, by fixing p = [\frac{n}{2}]+1 we can find a constant C_1 > 0 such that

    \frac{d}{dt}\left(\int_{\Omega}|\nabla u|^{2\left[\frac{n}{2}\right]+2}+\int_{\Omega}v^{\left[\frac{n}{2}\right]+1}\right)+\int_{\Omega}|\nabla u|^{2\left[\frac{n}{2}\right]+2}+\int_{\Omega} v^{\left[\frac{n}{2}\right]+1} \leqslant C_1\quad\text{for all}\ t\in(0, T_{max}),

    which along with the Grönwall inequality gives a constant C_2 > 0 ,

    \int_{\Omega} v^{\left[\frac{n}{2}\right]+1} \leqslant C_2\quad\text{for all}\ t\in(0, T_{max}).

    Together with Lemma 3.9, we finish the proof by Lemma 2.1.

    In this section, we will employ suitable Lyapunov functionals to study the large-time behavior of u , v and w . We first improve the regularity of the solution.

    Lemma 4.1. There exist constants \theta_1, \theta_2, \theta_3\in(0, 1) and C > 0 such that

    \|u\|_{C^{2+\theta_1, 1+\frac{\theta_1}{2}}(\overline{\Omega}\times[t, t+1])}+\|v\|_{C^{2+\theta_2, 1+\frac{\theta_2}{2}}(\overline{\Omega}\times[t, t+1])}+\|w\|_{C^{2+\theta_3, 1+\frac{\theta_3}{2}}(\overline{\Omega}\times[t, t+1])}\leqslant C\quad\mathit{\text{for all}}\ t > 1.

    In particular, one can find C > 0 such that

    \|\nabla u\|_{L^\infty(\Omega)}+\|\nabla v\|_{L^\infty(\Omega)}+\|\nabla w\|_{L^\infty(\Omega)}\leqslant C\quad\mathit{\text{for all}}\ t > 1.

    Proof. The conclusion is a consequence of the regularity of parabolic equations in [36].

    We split our analysis into two cases: a_1 > a_3 and a_1\leqslant a_3 .

    We know that there are four homogeneous equilibria (0, 0, 0) , (0, 1, 0) , \left(\frac{a_1}{a_2}, 0, \frac{ra_1}{\gamma a_2}\right) and (u_*, v_*, w_*) when a_1 > a_3 , where u_*, v_* and w_* are defined in (1.5). In this case, we shall prove the coexistence steady state (u_*, v_*, w_*) is globally exponentially stable under certain conditions. Define an energy functional for (1.1) as follows:

    \mathcal{F}(t) = \varepsilon_1\int_\Omega\left(u-u_*-u_*\ln\frac{u}{u_*}\right)+\int_\Omega\left(v-v_*-v_*\ln\frac{v}{v_*}\right)+\frac{\varepsilon_2}{2}\int_\Omega\left(w-w_*\right)^2,

    where \varepsilon_1 and \varepsilon_2 are to be determined below.

    Proof of Theorem 1.2–(1). We complete the proof in four steps.

    Step 1: The parameters \varepsilon_1 and \varepsilon_2 can be chosen in the following way. First, we recall from (1.5) and (1.6) that

    \begin{equation} \Gamma = \frac{4du_*}{K_1^2v_*}, \Phi = \frac{2 a_2 \rho}{a_3^2}+e, \Psi = \frac{\gamma\eta a_3^2K_1^2}{d\rho^2 r^2u_*}. \end{equation} (4.1)

    Let

    f(y) = \frac{\Psi(\Gamma y-\xi^2)(-y^2+2\Phi y-e^2)}{y}, \quad y > 0.

    It is clear that f\in C^0((0, +\infty)) . Then, if

    \frac{\xi^2}{\Gamma} < \Phi+\sqrt{\Phi^2-e^2},

    the following holds:

    \begin{equation} \frac{\xi^2K_1^2 v_*}{4du_*} < \frac{2a_2\rho}{a_3^2}+e+\frac{2}{a_3}\sqrt{a_2\rho\left(\frac{a_2\rho}{a_3^2}+e\right)}. \end{equation} (4.2)

    Under (4.2), we let a = \max\left\{\frac{\xi^2}{\Gamma}, \Phi-\sqrt{\Phi^2-e^2}\right\} and b = \Phi+\sqrt{\Phi^2-e^2} with a < b . Then, f(y) is continuous on [a, b] with f(a) = f(b) = 0 , and consequently f(y) must attain the maximum at some point, say \varepsilon_1 , in (a, b) , namely f(\varepsilon_1) = \max \limits_{y\in [a, b]} f(y) . Then, a < \varepsilon_1 < b , or equivalently (see (4.1))

    \begin{equation} \max\left\{\frac{\xi^2u^2 v_*}{4du_*}, \frac{2a_2\rho}{a_3^2}+e-\frac{2}{a_3}\sqrt{a_2\rho\left(\frac{a_2\rho}{a_3^2}+e\right)}\right\} < \varepsilon_1 < \frac{2a_2\rho}{a_3^2}+e+\frac{2}{a_3}\sqrt{a_2\rho\left(\frac{a_2\rho}{a_3^2}+e\right)}. \end{equation} (4.3)

    Next, we assume \chi > 0 is suitably small such that

    \begin{align*} \chi^2 < f(\varepsilon_1) = &\frac{\gamma\eta a_3^2K_1^2}{d\rho r^2u_*\varepsilon_1}\left(\frac{4du_*\varepsilon_1}{v_*K_1^2}-\xi^2\right)\left(-\varepsilon_1^2+2\left(\frac{2a_2\rho}{a_3^2}+e\right)\varepsilon_1-e^2\right)\\ = &\frac{4\gamma\eta}{d\rho r^2u_*v_*\varepsilon_1}\left(4du_*\varepsilon_1-\xi^2v_*K_1^2\right)\left(a_2\rho\varepsilon_1-\frac{a_3^2(\varepsilon_1-e)^2}{4}\right), \end{align*}

    which implies

    \frac{d\chi^2u_*v_*^2\varepsilon_1}{\eta\left(4du_*v_*\varepsilon_1-\xi^2v_*^2K_1^2\right)} < \frac{4\gamma}{\rho r^2}\left(a_2\rho\varepsilon_1-\frac{a_3^2(\varepsilon_1-e)^2}{4}\right).

    Hence, there exists a constant \varepsilon_2 > 0 such that

    \begin{equation*} \label{24} \frac{d\chi^2u_*v_*^2\varepsilon_1}{\eta\left(4du_*v_*\varepsilon_1-\xi^2v_*^2K_1^2\right)} < \varepsilon_2 < \frac{4\gamma}{\rho r^2}\left(a_2\rho\varepsilon_1-\frac{a_3^2(\varepsilon_1-e)^2}{4}\right) \end{equation*}

    which along with Lemma 3.1 yields

    \begin{equation} \frac{d\chi^2u_*v_*^2\varepsilon_1}{\eta\left(4du_*v_*\varepsilon_1-\xi^2v_*^2u^2\right)} < \varepsilon_2 < \frac{4\gamma}{\rho r^2}\left(a_2\rho\varepsilon_1-\frac{a_3^2(\varepsilon_1-e)^2}{4}\right). \end{equation} (4.4)

    Step 2: We claim

    \|u-u_*\|_{L^{\infty}(\Omega)}+\|v-v_*\|_{L^{\infty}(\Omega)}+\|w-w_*\|_{L^{\infty}(\Omega)} \rightarrow0\quad\text{as}\ t\rightarrow +\infty.

    Indeed, using the equations in system (1.1) along with integration by parts, we have

    \begin{aligned} &\frac{d}{dt}\int_\Omega\left(u-u_*-u_*\ln\frac{u}{u_*}\right) = \int_\Omega \frac{u-u_*}{u}u_t\\ = &-du_*\int_\Omega\frac{|\nabla u|^2}{u^2}+\int_\Omega(u-u_*)(a_1-a_2u-a_3v)\\ = &-du_*\int_\Omega\frac{|\nabla u|^2}{u^2}-a_2\int_\Omega(u-u_*)^2-a_3\int_\Omega(u-u_*)(v-v_*). \end{aligned}

    Similarly, we obtain

    \begin{aligned} &\frac{d}{dt}\int_\Omega\left(v-v_*-v_*\ln\frac{v}{v_*}\right) = \int_\Omega \frac{v-v_*}{v}v_t\\ = &-v_*\int_\Omega\frac{|\nabla v|^2}{v^2}-\chi v_*\int_\Omega\frac{\nabla v\cdot\nabla w}{v}+\xi v_*\int_\Omega\frac{\nabla u\cdot\nabla v}{v}+\int_\Omega(v-v_*)(\rho-\rho v+ea_3u)\\ = &-v_*\int_\Omega\frac{|\nabla v|^2}{v^2}-\chi v_*\int_\Omega\frac{\nabla v\cdot\nabla w}{v}+\xi v_*\int_\Omega\frac{\nabla u\cdot\nabla v}{v}\\ &\quad-\rho\int_\Omega(v-v_*)^2+ea_3\int_\Omega(u-u_*)(v-v_*) \end{aligned}

    and

    \begin{align*} \frac{d}{dt}\int_\Omega\left(w-w_*\right)^2 = &2\int_\Omega\left(w-w_*\right)w_t = 2\int_\Omega\left(w-w_*\right)\left(\eta \Delta w+ru-\gamma w\right)\\ = &-2\eta\int_\Omega|\nabla w|^2+2r\int_\Omega(u-u_*)(w-w_*)-2\gamma\int_\Omega(w-w_*)^2\quad\text{for all}\ t > 0. \end{align*}

    Then, it follows that

    \begin{align*} \frac{d}{dt}\mathcal{F}(t) = &-du_*\varepsilon_1\int_\Omega\frac{|\nabla u|^2}{u^2}-v_*\int_\Omega\frac{|\nabla v|^2}{v^2}-\eta \varepsilon_2\int_\Omega|\nabla w|^2\\ &\quad-\chi v_*\int_\Omega\frac{\nabla v\cdot\nabla w}{v}+\xi v_*\int_\Omega\frac{\nabla u\cdot\nabla v}{v}\\ &\quad-a_2\varepsilon_1\int_\Omega(u-u_*)^2-\rho\int_\Omega(v-v_*)^2-\gamma\varepsilon_2\int_\Omega(w-w_*)^2\\ &\quad-a_3(\varepsilon_1-e)\int_\Omega(u-u_*)(v-v_*)+r\varepsilon_2\int_\Omega(u-u_*)(w-w_*)\\ = &:-X^TSX-Y^TTY, \end{align*}

    where X = (\nabla u, \nabla v, \nabla w) , Y = (u-u_*, v-v_*, w-w_*) , and

    S = \left[\begin{array}{ccc} \frac{du_*\varepsilon_1}{u^2} & -\frac{\xi v_*}{2v} & 0\\ -\frac{\xi v_*}{2v} & \frac{v_*}{v^2} & \frac{\chi v_*}{2v}\\ 0 & \frac{\chi v_*}{2v} & \eta\varepsilon_2 \end{array}\right], \quad T = \left[\begin{array}{ccc} a_2\varepsilon_1 & \frac{a_3(\varepsilon_1-e)}{2} & -\frac{r\varepsilon_2}{2}\\ \frac{a_3(\varepsilon_1-e)}{2} & \rho & 0\\ -\frac{r\varepsilon_2}{2} & 0 & \gamma\varepsilon_2 \end{array}\right].

    Note that (4.3) yields

    \frac{du_*v_*\varepsilon_1}{u^2v^2}-\frac{\xi^2v_*^2}{4v^2} > \frac{v_*^2}{4v^2}\bigg(\frac{4d u_*\varepsilon}{K_1^2}-\xi^2\bigg) > 0,

    and (4.4) gives

    \frac{\eta du_*v_*\varepsilon_1\varepsilon_2}{u^2v^2}-\frac{d\chi^2u_*v_*^2\varepsilon_1}{4u^2v^2}-\frac{\eta \xi^2v_*^2\varepsilon_2}{4v^2} > 0.

    The above results indicate that matrix S is positive definite. Using (4.3) and (4.4) again, we observe that

    a_2\rho\varepsilon_1-\frac{a_3^2(\varepsilon_1-e)^2}{4} > 0,

    and

    a_2\rho\gamma\varepsilon_1\varepsilon_2-\frac{\rho r^2\varepsilon_2^2}{4}-\frac{a_3^2\gamma(\varepsilon_1-e)^2\varepsilon_2}{4} > 0,

    which imply that matrix T is positive definite. Therefore, one can choose a constant C_1 > 0 such that

    \begin{equation} \frac{d}{dt}\mathcal{F}(t)\leqslant-C_1\left(\int_\Omega(u-u_*)^2+\int_\Omega(v-v_*)^2+\int_\Omega(w-w_*)^2\right)\quad\text{for all}\ t > 0. \end{equation} (4.5)

    Integrating the above inequality with respect to time, we get a constant C_2 > 0 satisfying

    \int_{1}^{+\infty}\int_\Omega(u-u_*)^2+\int_{1}^{+\infty}\int_\Omega(v-v_*)^2+\int_{1}^{+\infty}\int_\Omega(w-w_*)^2\leqslant C_2,

    which together with the uniform continuity of u, v and w due to Lemma 4.1 yields

    \begin{equation} \int_\Omega(u-u_*)^2+\int_\Omega(v-v_*)^2+\int_\Omega(w-w_*)^2\rightarrow0, \quad\text{as}\ t\rightarrow +\infty. \end{equation} (4.6)

    By the Gagliardo-Nirenberg inequality, we can find a constant C_3 > 0 such that

    \begin{equation} \|u-u_*\|_{L^{\infty}(\Omega)} \leqslant C_3\|u-u_*\|_{W^{1, \infty}(\Omega)}^{\frac{n}{n+2}}\|u-u_*\|_{L^{2}(\Omega)}^{\frac{2}{n+2}}, \end{equation} (4.7)
    \begin{equation} \|v-v_*\|_{L^{\infty}(\Omega)} \leqslant C_3\|v-v_*\|_{W^{1, \infty}(\Omega)}^{\frac{n}{n+2}}\|v-v_*\|_{L^{2}(\Omega)}^{\frac{2}{n+2}} \end{equation} (4.8)

    and

    \begin{equation} \|w-w_*\|_{L^{\infty}(\Omega)} \leqslant C_3\|w-w_*\|_{W^{1, \infty}(\Omega)}^{\frac{n}{n+2}}\|w-w_*\|_{L^{2}(\Omega)}^{\frac{2}{n+2}}\quad\text{for all}\ t > 1, \end{equation} (4.9)

    which along with (4.6) and Lemma 4.1 prove the claim.

    Step 3: From the L'Hôpital rule, it holds that for any s_0 > 0

    \lim\limits_{s\rightarrow s_0}\frac{s-s_0-s_0\ln\frac{s}{s_0}}{(s-s_0)^2} = \lim\limits_{s\rightarrow s_0}\frac{1-\frac{s_0}{s}}{2(s-s_0)} = \lim\limits_{s\rightarrow s_0}\frac{1}{2s} = \frac{1}{2s_0},

    which gives a constant \eta > 0 such that for all |s-s_0|\leqslant\eta

    \begin{equation} \frac{1}{4s_0}(s-s_0)^2\leqslant s-s_0-s_0\ln\frac{s}{s_0}\leqslant\frac{1}{s_0}(s-s_0)^2. \end{equation} (4.10)

    By (4.6), there exists T_1 > 1 such that

    \begin{equation*} \|u-u_*\|_{L^\infty(\Omega)}+\|v-v_*\|_{L^\infty(\Omega)}+\|w-w_*\|_{L^\infty(\Omega)}\leqslant\eta\quad\text{for all}\ t\geqslant T_1. \end{equation*}

    Therefore, by (4.10), we get

    \begin{equation} \frac{1}{4u_*}\int_\Omega(u-u_*)^2\leqslant\int_\Omega \left(u-u_*-u_*\ln\frac{u}{u_*}\right)\leqslant\frac{1}{u_*}\int_\Omega(u-u_*)^2\quad\text{for all}\ t\geqslant T_1 \end{equation} (4.11)

    and

    \begin{equation} \frac{1}{4v_*}\int_\Omega(v-v_*)^2\leqslant\int_\Omega \left(v-v_*-v_*\ln\frac{v}{v_*}\right)\leqslant\frac{1}{v_*}\int_\Omega(v-v_*)^2\quad\text{for all}\ t\geqslant T_1. \end{equation} (4.12)

    Step 4: From (4.11) and (4.12), it follows that

    \mathcal{F}(t)\leqslant\max\left\{\frac{\varepsilon_1}{u_*}, \frac{1}{v_*}, \frac{\varepsilon_2}{2}\right\}\left(\int_\Omega(u-u_*)^2+\int_\Omega(v-v_*)^2+\int_\Omega(w-w_*)^2\right),

    which alongside (4.5) yields a constant C_4 > 0 such that

    \frac d{dt}\mathcal{F}(t)\leqslant-C_4\mathcal{F}(t)\quad\text{for all}\ t\geqslant T_1.

    This immediately gives a constant C_5 > 0 such that

    \begin{equation*} \mathcal{F}(t)\leqslant C_5e^{-C_4 t}\quad\text{for all}\ t\geqslant T_1. \end{equation*}

    Hence, utilizing (4.11) and (4.12) again, one obtains a constant C_6 > 0 such that

    \begin{equation*} \int_\Omega(u-u_*)^2+\int_\Omega(v-v_*)^2+\int_\Omega(w-w_*)^2\leqslant C_6e^{-C_4 t}\quad\text{for all}\ t\geqslant T_1. \end{equation*}

    Finally, by (4.7)–(4.9) and Lemma 4.1, we get the decay rates of \|u-u_*\|_{L^\infty(\Omega)} , \|v-v_*\|_{L^\infty(\Omega)} and \|w-w_*\|_{L^\infty(\Omega)} , as claimed in Theorem 1.2–(1).

    In this case, there are three homogeneous equilibria (0, 0, 0) , (0, 1, 0) and \left(\frac{a_1}{a_2}, 0, \frac{ra_1}{\gamma a_2}\right) , and we shall show that the steady state (0, 1, 0) is global asymptotically stable, where the convergence rate is exponential if a_1 < a_3 and algebraic if a_1 = a_3 . Define an energy functional for (1.1) as follows:

    G(t) = e\int_\Omega u+\frac{\zeta_1}{2}\int_\Omega u^2+\int_\Omega\left(v-1-\ln v\right)+\frac{\zeta_2}{2}\int_\Omega w^2,

    where \zeta_1 and \zeta_2 will be determined below.

    Proof of Theorem 1.2–(2). We divide the proof into five steps.

    Step 1: We shall choose the appropriate parameters \zeta_1 and \zeta_2 . By the definitions of A and B in (1.7), since A < B , we have

    \begin{equation} \left(\frac{\xi^2}{4d}\right)^2 < \frac{\xi^2ea_2}{4da_1} < \left(\frac{ea_2}{a_1}\right)^2. \end{equation} (4.13)

    Let

    g(y) = \frac{16\eta\gamma}{dr^2}\frac{(dy-\frac{\xi^2}{4})(ea_2-a_1y)}{y}, \quad \frac{\xi^2}{4d} < y < \frac{ea_2}{a_1}.

    Then, g\in C^1\left(\left(\frac{\xi^2}{4d}, \frac{ea_2}{a_1}\right)\right) , and g(y) > 0 in \left(\frac{\xi^2}{4d}, \frac{ea_2}{a_1}\right) . We further observe that

    g\left(\frac{\xi}{2}\sqrt{\frac{ea_2}{da_1}}\right) = D\left(A+B-2\sqrt{AB}\right)

    which along with \chi^2 < D\left(A+B-2\sqrt{AB}\right) implies

    \chi^2 < g\left(\frac{\xi}{2}\sqrt{\frac{ea_2}{da_1}}\right).

    By the definition of g , one has

    g'(y_0) = \frac{16\eta\gamma}{dr^2}\left(-da_1+\frac{\xi^2ea_2}{4y_0^2}\right) = 0,

    which alongside (4.13) gives y_0 = \frac{\xi}{2}\sqrt{\frac{ea_2}{da_1}}\in\left(\frac{\xi^2}{4d}, \frac{ea_2}{a_1}\right) . Thus, g(y) is increasing in \left(\frac{\xi^2}{4d}, \frac{\xi}{2}\sqrt{\frac{ea_2}{da_1}}\right) and decreasing in \left(\frac{\xi}{2}\sqrt{\frac{ea_2}{da_1}}, \frac{ea_2}{a_1}\right) . We can find a constant \zeta_1 > 0 such that

    \begin{equation} \frac{\xi}{2}\sqrt{\frac{ea_2}{da_1}} < \zeta_1 < \frac{ea_2}{a_1} \end{equation} (4.14)

    and

    0 = g\left(\frac{ea_2}{a_1}\right) < \chi^2 < g(\zeta_1) < g\left(\frac{\xi}{2}\sqrt{\frac{ea_2}{da_1}}\right).

    With the definition of g , we get

    \frac{d\chi^2\zeta_1}{4\eta(d\zeta_1-\frac{\xi^2}{4})} < \frac{4\gamma}{r^2}(ea_2-a_1\zeta_1),

    which implies that there exists \zeta_2 > 0 such that

    \begin{equation} \frac{d\chi^2\zeta_1}{4\eta(d\zeta_1-\frac{\xi^2}{4})} < \zeta_2 < \frac{4\gamma}{r^2}(ea_2-a_1\zeta_1). \end{equation} (4.15)

    One can verify that

    \begin{equation} \eta d\zeta_1\zeta_2-\frac{d\chi^2}{4}\zeta_1-\frac{\eta \xi^2}{4}\zeta_2 > 0, \end{equation} (4.16)

    and

    \begin{equation} (ea_2-a_1\zeta_1)\rho\gamma\zeta_2-\frac{\rho r^2}{4}\zeta_2^2 > 0. \end{equation} (4.17)

    Thanks to (4.13) and (4.14), one obtains

    \begin{equation} \frac{\xi^2}{4d} < \zeta_1 < \frac{ea_2}{a_1}. \end{equation} (4.18)

    Step 2: We claim

    \begin{equation} \|u\|_{L^{\infty}(\Omega)}+\|v-1\|_{L^{\infty}(\Omega)}+\|w\|_{L^{\infty}(\Omega)} \rightarrow0\quad\text{as}\ t\rightarrow +\infty. \end{equation} (4.19)

    Indeed, if (u, v, w) is the solution of system (1.1), then we get

    \begin{gather} \frac{d}{dt}\int_\Omega u = a_1\int_\Omega u-a_2\int_\Omega u^2-a_3\int_\Omega uv, \end{gather} (4.20)
    \begin{gather} \frac{d}{dt}\int_\Omega u^2 = 2\int_\Omega uu_t = -2d\int_\Omega|\nabla u|^2+2a_1\int_\Omega u^2-2a_2\int_\Omega u^3-2a_3\int_\Omega u^2v, \end{gather} (4.21)
    \begin{gather} \begin{aligned} &\frac{d}{dt}\int_\Omega\left(v-1-\ln v\right) = \int_\Omega \frac{v-1}{v}v_t\\ = &-\int_\Omega\frac{|\nabla v|^2}{v^2}-\chi \int_\Omega\frac{\nabla v\cdot\nabla w}{v}+\xi \int_\Omega\frac{\nabla u\cdot\nabla v}{v}+\int_\Omega(v-1)(\rho-\rho v+ea_3u)\\ = &-\int_\Omega\frac{|\nabla v|^2}{v^2}-\chi \int_\Omega\frac{\nabla v\cdot\nabla w}{v}+\xi\int_\Omega\frac{\nabla u\cdot\nabla v}{v}-\rho\int_\Omega(v-1)^2+ea_3\int_\Omega uv-ea_3\int_\Omega u \end{aligned} \end{gather} (4.22)

    and

    \begin{equation} \frac{d}{dt}\int_\Omega w^2 = 2\int_\Omega ww_t = -2\eta\int_\Omega|\nabla w|^2+2r\int_\Omega uw-2\gamma\int_\Omega w^2\quad\text{for all}\ t > 0. \end{equation} (4.23)

    Then, combining (4.20), (4.21), (4.22) and (4.23), we have from the definition of G(t) that

    \begin{equation} \begin{aligned} \frac{d}{dt}G(t)\leqslant&-d\zeta_1\int_\Omega|\nabla u|^2-\int_\Omega\frac{|\nabla v|^2}{v^2}-\eta\zeta_2\int_\Omega|\nabla w|^2\\ &\quad-\chi \int_\Omega\frac{\nabla v\cdot\nabla w}{v}+\xi\int_\Omega\frac{\nabla u\cdot\nabla v}{v}+e(a_1-a_3)\int_\Omega u\\ &\quad-(ea_2-a_1\zeta_1)\int_\Omega u^2-\rho\int_\Omega(v-1)^2-\gamma\zeta_2\int_\Omega w^2+r\zeta_2\int_\Omega uw\\ = &:-X^TPX-Y^TQY+e(a_1-a_3)\int_\Omega u, \end{aligned} \end{equation} (4.24)

    where X = (\nabla u, \nabla v, \nabla w) , Y = (u, v-1, w) ,

    P = \left[\begin{array}{ccc} d\zeta_1 & -\frac{\xi}{2v} & 0\\ -\frac{\xi}{2v} & \frac{1}{v^2} & \frac{\chi}{2v}\\ 0 & \frac{\chi}{2v} & \eta\zeta_2 \end{array}\right]\quad\text{and}\quad Q = \left[\begin{array}{ccc} ea_2-a_1\zeta_1 & 0 & -\frac{r\zeta_2}{2}\\ 0 & \rho & 0\\ -\frac{r\zeta_2}{2} & 0 & \gamma\zeta_2 \end{array}\right].

    It can be checked that (4.16) and (4.18) ensure that the matrix P is positive definite while (4.17) and (4.18) guarantee that the matrix Q is positive definite. Thus, there is a constant C_1 > 0 such that if a_1 < a_3 , then

    \begin{equation} \frac{d}{dt}G(t)\leqslant-C_1\left(\int_\Omega u+\int_\Omega u^2+\int_\Omega(v-1)^2+\int_\Omega w^2\right)\quad\text{for all}\ t > 0, \end{equation} (4.25)

    and if a_1 = a_3 , then

    \begin{equation} \frac{d}{dt}G(t)\leqslant-C_1\left(\int_\Omega u^2+\int_\Omega(v-1)^2+\int_\Omega w^2\right)\quad\text{for all}\ t > 0. \end{equation} (4.26)

    Integrating the above inequalities with respect to time, we find a constant C_2 > 0 satisfying

    \int_1^{+\infty}\int_\Omega u^2+\int_1^{+\infty}\int_\Omega(v-1)^2+\int_1^{+\infty}\int_\Omega w^2\leqslant C_2,

    which together with the uniform continuity of u, v and w due to Lemma 4.1 yields

    \begin{equation} \int_\Omega u^2+\int_\Omega(v-1)^2+\int_\Omega w^2\rightarrow0, \quad\text{as}\ t\rightarrow +\infty. \end{equation} (4.27)

    Thus, (4.19) is obtained by the Gagliardo-Nirenberg inequality and Lemma 4.1.

    Step 3: By the L'Hôpital rule, we get

    \lim\limits_{s\rightarrow 1}\frac{s-1-\ln s}{(s-1)^2} = \lim\limits_{s\rightarrow 1}\frac{1-\frac{1}{s}}{2(s-1)} = \lim\limits_{s\rightarrow 1}\frac{1}{2s} = \frac{1}{2},

    which gives a constant \varepsilon > 0 such that

    \begin{equation} \frac{1}{4}(s-1)^2\leqslant s-1-\ln s\leqslant(s-1)^2 \ \ \text{for all}\ \ |s-1|\leqslant\varepsilon. \end{equation} (4.28)

    By (4.19), there exists T_1 > 0 such that

    \begin{equation} \|u\|_{L^\infty(\Omega)}+\|v-1\|_{L^\infty(\Omega)}+\|w\|_{L^\infty(\Omega)}\leqslant\varepsilon\quad\text{for all}\ t\geqslant T_1. \end{equation} (4.29)

    Therefore, it follows from (4.28) that

    \begin{equation} \frac{1}{4}\int_\Omega(v-1)^2\leqslant\int_\Omega (v-1-\ln v)\leqslant \int_\Omega(v-1)^2\quad\text{for all}\ t\geqslant T_1. \end{equation} (4.30)

    Step 4: If a_1 < a_3 , from the definition of G(t) and (4.30), one has

    G(t)\leqslant\max\left\{e, \frac{\zeta_1}{2}, \frac{\zeta_2}{2}, 1\right\}\left(\int_\Omega u+\int_\Omega u^2+\int_\Omega(v-1)^2+\int_\Omega w^2\right),

    which along with (4.25) yields a constant C_3 > 0 such that

    \frac d{dt}G(t)\leqslant-C_3G(t)\quad\text{for all}\ t\geqslant T_1.

    This gives a constant C_4 > 0 such that

    \begin{equation*} G(t)\leqslant C_4e^{-C_3 t}\quad\text{for all}\ t\geqslant T_1. \end{equation*}

    Hence, utilizing (4.30) again, we find a constant C_5 > 0 such that

    \begin{equation*} \int_\Omega u^2+\int_\Omega(v-1)^2+\int_\Omega w^2\leqslant C_5e^{-C_3 t}\quad\text{for all}\ t\geqslant T_1. \end{equation*}

    Then, by the Gagliardo-Nirenberg inequality and Lemma 4.1, we get the exponential convergence for \|u\|_{L^{\infty}(\Omega)}+\|v-1\|_{L^{\infty}(\Omega)}+\|w\|_{L^{\infty}(\Omega)} .

    Step 5: If a_1 = a_3 , we use (4.29), (4.30) and Young's inequality to find a constant C_6 > 0 :

    \begin{align*} G^2(t)\leqslant& C_6\left(\int_\Omega u+\int_\Omega u^2+\int_\Omega(v-1)^2+\int_\Omega w^2\right)^2\\ \leqslant&C_6(\varepsilon+1)^2\left(\int_\Omega u+\int_\Omega(v-1)+\int_\Omega w \right)^2\\ \leqslant&3C_6(\varepsilon+1)^2|\Omega|\left(\int_\Omega u^2+\int_\Omega(v-1)^2+\int_\Omega w^2 \right)\quad\text{for all}\ t\geqslant T_1, \end{align*}

    which alongside (4.26) implies some constant C_7 > 0

    \frac{d}{dt}G(t)\leqslant-C_7G^2(t)\quad\text{for all}\ t\geqslant T_1.

    Solving the above inequality directly yields a constant C_8 > 0 such that

    G(t)\leqslant C_8(t+1)^{-1}\quad\text{for all}\ t\geqslant T_1.

    Similar to the case a_1 < a_3 , we can use (4.30), the Gagliardo-Nirenberg inequality and Lemma 4.1 to get the convergence rate of \|u\|_{L^{\infty}(\Omega)}+\|v-1\|_{L^{\infty}(\Omega)}+\|w\|_{L^{\infty}(\Omega)} .

    The author warmly thanks the reviewers for several inspiring comments and helpful suggestions. The research of the author was supported by the National Nature Science Foundation of China (Grant No. 12101377) and the Nature Science Foundation of Shanxi Province (Grant No. 20210302124080).

    The author declares there is no conflict of interest.



    [1] W. O. Kermack, A. G. McKendrick, A contribution to the mathematical theory of epidemics, Proc. R. Soc. London, Ser. A, 115 (1927), 700-721. doi: 10.1098/rspa.1927.0118
    [2] P. A. Naik, Global dynamics of a fractional-order SIR epidemic model with memory, Int. J. Biomath., 13 (2020), 1-24.
    [3] P. A. Naik, J. Zu, M. Ghoreishi, Stability analysis and approximate solution of SIR epidemic model with Crowley-Martin type functional response and holling type-II treatment rate by using homotopy analysis method, J. Appl. Anal. Comput., 10 (2020), 1482-1515.
    [4] P. A. Naik, J. Zu, K. M. Owolabi, Global dynamics of a fractional order model for the transmission of HIV epidemic with optimal control, Chaos Solitons Fractals, 138 (2020), 109826.
    [5] P. A. Naik, J. Zu, K. M. Owolabi, Modeling the mechanics of viral kinetics under immune control during primary infection of HIV-1 with treatment in fractional order, Phys. A, 545 (2020), 123816.
    [6] J. Zu, M. Li, Y. Gu, S. Fu, Modelling the evolutionary dynamics of host resistance-related traits in a susceptible-infected community with density-dependent mortality, Discrete Contin. Dyn. Syst. Ser. B, 25 (2020), 3049-3086.
    [7] J. Zu, G. Zhuang, P. Liang, F. Cui, F. Wang, H. Zheng, et al., Estimating age-related incidence of HBsAg seroclearance in chronic hepatitis B virus infections of China by using a dynamic compartmental model, Sci. Rep., 7 (2017), 2912.
    [8] J. Zu, M. Li, G. Zhuang, P. Liang, F. Cui, F. Wang, et al., Estimating the impact of test-and-treat strategies on hepatitis B virus infection in China by using an age-structured mathematical model, Medicine, 97 (2018), e0484.
    [9] A. Chekroun, M. N. Frioui, T. Kuniya, T. M. Touaoula, Global stability of an age-structured epidemic model with general Lyapunov functional, Math. Biosci. Eng., 16 (2019), 1525-1553. doi: 10.3934/mbe.2019073
    [10] Z. Feng, W. Huang, C. Castillo-Chavez, Global behavior of a multi-group SIS epidemic model with age structure, J. Differ. Equations, 218 (2005), 292-324. doi: 10.1016/j.jde.2004.10.009
    [11] T. Kuniya, Global stability analysis with a discretization approach for an age-structured multigroup SIR epidemic model, Nonlinear Anal. Real World Appl., 12 (2011), 2640-2655. doi: 10.1016/j.nonrwa.2011.03.011
    [12] J. Wang, H. Shu, Global analysis on a class of multi-group SEIR model with latency and relapse, Math. Biosci. Eng., 13 (2016), 209-225. doi: 10.3934/mbe.2016.13.209
    [13] J. Yang, R. Xu, X. Luo, Dynamical analysis of an age-structured multi-group SIVS epidemic model, Math. Biosci. Eng., 16 (2019), 636-666. doi: 10.3934/mbe.2019031
    [14] R. S. Cantrell, C. Cosner, Spatial ecology via reaction-diffusion equations, Wiley Series in Mathematical and Computational Biology, John Wiley and Sons Ltd, Chichester, 2004.
    [15] J. Wu, Theory and applications of partial functional differential equations, Springer-Verlag, New York, 1996.
    [16] A. Chekroun, T. Kuniya, An infection age-space structured SIR epidemic model with Neumann boundary condition, Appl. Anal., 1-14.
    [17] R. Cui, Y. Lou, A spatial SIS model in advective heterogeneous environments, J. Differ. Equations, 261 (2016), 3305-3343. doi: 10.1016/j.jde.2016.05.025
    [18] A. Ducrot, P. Magal, S. Ruan, Travelling wave solutions in multigroup age-structured epidemic models, Arch. Ration. Mech. Anal., 195 (2010), 311-331. doi: 10.1007/s00205-008-0203-8
    [19] W. E. Fitzgibbon, J. J. Morgan, G. F. Webb, Y. Wu, A vector-host epidemic model with spatial structure and age of infection, Nonlinear Anal. Real World Appl., 41 (2018), 692-705. doi: 10.1016/j.nonrwa.2017.11.005
    [20] T. Kuniya, R. Oizumi, Existence result for an age-structured SIS epidemic model with spatial diffusion, Nonlinear Anal. Real World Appl., 23 (2015), 196-208. doi: 10.1016/j.nonrwa.2014.10.006
    [21] Y. Luo, S. Tang, Z. Teng, Global dynamics in a reaction-diffusion multi-group SIR epidemic model with nonlinear incidence, Nonlinear Anal. Real World Appl., 50 (2019), 365-385. doi: 10.1016/j.nonrwa.2019.05.008
    [22] X. Wang, X. Q. Zhao, J. Wang, A cholera epidemic model in a spatiotemporally heterogeneous environment, J. Math. Anal. Appl., 468 (2018), 893-912. doi: 10.1016/j.jmaa.2018.08.039
    [23] J. Yang, Z. Jin, F. Xu, Threshold dynamics of an age-space structured SIR model on heterogeneous environment, Appl. Math. Lett., 96 (2019), 69-74. doi: 10.1016/j.aml.2019.03.009
    [24] J. Yang, R. Xu, J. Li, Threshold dynamics of an age-space structured brucellosis disease model with Neumann boundary condition, Nonlinear Anal. Real World Appl., 50 (2019), 192-217. doi: 10.1016/j.nonrwa.2019.04.013
    [25] L. Zhao, Z. C. Wang, S. Ruan, Traveling wave solutions in a two-group SIR epidemic model with constant recruitment, J. Math. Biol., 77 (2018), 1871-1915. doi: 10.1007/s00285-018-1227-9
    [26] L. Liu, J. Wang, X. Liu, Global stability of an SEIR epidemic model with age-dependent latency and relapse, Nonlinear Anal. Real World Appl., 24 (2015), 18-35. doi: 10.1016/j.nonrwa.2015.01.001
    [27] L. Liu, X. Feng, A multigroup SEIR epidemic model with age-dependent latency and relapse, Math. Methods Appl. Sci., 41 (2018), 6814-6833. doi: 10.1002/mma.5193
    [28] X. Ren, Y. Tian, L. Liu, X. Liu, A reaction-diffusion within-host HIV model with cell-to-cell transmission, J. Math. Biol., 76 (2018), 1831-1872. doi: 10.1007/s00285-017-1202-x
    [29] H. L. Smith, Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems, vol. 41 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI, 1995.
    [30] K. Yosida, Functional analysis, Springer, New York, 1978.
    [31] Y. Lou, X. Q. Zhao, A reaction-diffusion malaria model with incubation period in the vector population, J. Math. Biol., 62 (2011), 543-568. doi: 10.1007/s00285-010-0346-8
    [32] J. K. Hale, Asymptotic behavior of dissipative systems, American Mathematical Society, Providence, RI, 1988.
    [33] A. Berman, R. J. Plemmons, Nonnegative matrices in mathematical sciences, Academic Press, New York, 1979.
    [34] M. A. Krasnosel'skiǐ, Positive solutions of operator equations, Noordhoff, Groningen, 1964.
    [35] M. Y. Li, Z. Shuai, Global-stability problem for coupled systems of differential equations on networks, J. Differ. Equations, 248 (2010), 1-20. doi: 10.1016/j.jde.2009.09.003
  • This article has been cited by:

    1. Luis E. Ayala-Hernández, Gabriela Rosales-Muñoz, Armando Gallegos, María L. Miranda-Beltrán, Jorge E. Macías-Díaz, On a deterministic mathematical model which efficiently predicts the protective effect of a plant extract mixture in cirrhotic rats, 2023, 21, 1551-0018, 237, 10.3934/mbe.2024011
    2. Rinaldo M. Colombo, Paola Goatin, Elena Rossi, A Hyperbolic–parabolic framework to manage traffic generated pollution, 2025, 85, 14681218, 104361, 10.1016/j.nonrwa.2025.104361
  • Reader Comments
  • © 2020 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(4353) PDF downloads(238) Cited by(7)

Figures and Tables

Figures(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog