Dynamics of a Diffusive Multigroup SVIR Model with Nonlinear Incidence

In this paper, a multigroup SVIR epidemic model with reaction-diﬀusion and nonlinear incidence is investigated. We ﬁrst establish the well-posedness of the model. Then, the basic reproduction number R 0 is established and shown as a threshold: the disease-free steady state is globally asymptotically stable if R 0 < 1, while the disease will be persistent when R 0 > 1. Moreover, applying the classical method of Lyapunov and a recently developed graph-theoretic approach, we established the global stability of the endemic equilibria for a special case.


Introduction
As is well-known, mathematical models have played a key role in describing the dynamical evolution of infectious diseases [1][2][3][4][5][6][7][8][9][10]. For example, during the special battle for preventing the novel coronavirus (COVID-19) spreading, not only medical and biological research but also theoretical studies based on mathematical modeling may play a crucial role in analyzing and forecasting the spreading of the disease (see, for example, [4][5][6][7][8][9][10] and references therein; also, some useful advices for controlling the disease spreading have been given by the researchers). e spatial heterogeneity has been regarded as an important role that affects the spatial spreading of disease. erefore, complex models are needed to investigate the spread of diseases. In recent years, reaction-diffusion models involving environmental heterogeneity have been proposed and studied to find reliable measures to control disease spreading [11][12][13][14][15][16][17][18][19][20][21][22]. In particular, the basic reproduction number is an important threshold value for investigation of the dynamics of epidemic models. en, Wang and Zhao [14] defined the basic reproduction number and obtained its computation formula for general reaction-diffusion epidemic models. Xu et al. [21] studied an SVIR epidemic model with reaction-diffusion and the existences of travelling waves were investigated, while Zhang and Liu [22] also studied the existence of travelling waves for an SVIR epidemic model but with nonlocal dispersal and time delay.
Notice that the essential assumption in classical compartmental epidemic models is that the individuals are homogeneously mixed, which means each individual has the same probability to get infected. However, infected probability may be different for each individual in terms of the impact of different factors, such as education levels, age, gender, and communities. To overcome this problem, multigroup models have been proposed by many researchers by dividing the individuals into different groups and most of these work focus on the global dynamics of the models (see, for example, [23][24][25][26][27][28][29][30]). Particularly, one of the earliest works of multigroup models was proposed by Lajmanovich and York [23] when studied the gonorrhea disease transmission in a nonhomogeneous population. In order to investigate the global dynamics of multigroup models, a subtle grouping method in estimating the derivatives of Lyapunov functionals guided by graph theory for multigroup models was developed in [28][29][30]. For example, Kuniya [25] considered a multigroup SVIR epidemic model with vaccination and the global stability of endemic equilibria were established by the method of Lyapunov function. However, spatially heterogeneous was not considered in the model in [25]. To the best of our knowledge, there are few results of multigroup SVIR epidemic model with reaction-diffusion. e existing results are focusing on SIS and SIR models (see [31][32][33][34]). On the contrary, incidence rate also plays an important role in modeling the epidemic models, which has been frequently used to describe the nature of certain phenomena and obtained much exact results. In addition, applying general incidence rates can obtain the unification theory by the omission of unessential detail, see, for example, [21,[35][36][37] and references therein. Hence, inspired by [21,25] and the above considerations, in this paper, we consider the following multigroup SVIR epidemic model with reactiondiffusion and general incidence rate: with the homogeneous Neumann boundary conditions (2) and the initial conditions where Ω is a domain in R n with smooth boundary Ω and ] is the outward normal vector to the boundary zΩ. Initial functions ϕ ki (x) (k � 1, 2, 3, 4, 1 ≤ i ≤ n) are nonnegative and continuously defined on � Ω.
stand for the densities of the susceptible, vaccinated, infective, and recovered individuals in the i-th group at time t and spatial location x, respectively.
, and μ R i (x) denote the natural death rates of S i , V i , I i , and R i in spatial location x, respectively; δ i (x) is the death rate induced by the disease in spatial location x; ξ i (x) is the vaccination rate of S i in spatial location x; c i (x) is the rate of recovery from infection in spatial location x; β ij (x) is the infection rate of S i infected by I j in spatial location x; β ij (x) is the infection rate of V i infected by I j in spatial location x; and d 1i (x), d 2i (x), d 3i (x), and d 4i (x) are the diffusion rate of S i , V i , I i , and R i in spatial location x, respectively. All location-dependent parameters are continuous and strictly positive defined on � Ω, and f j (I j ) denotes the force of infection. We assume the function f j (I j ) satisfies the following properties: It is natural to assume that f j (0) � 0 due to the fact that disease cannot spread if there is no infection. e disease spreads heavily with the increasing number of infected individuals; thus, it reasonable to suppose f j ′ (I j ) > 0. Since the infectivity of infected individuals cannot be unbounded, it should reach a certain level when the infected individuals are heavy. erefore, the assumption f ′ ′ j (I j ) ≤ 0 implies that there exists a peak level for the infectivity of the infected individuals at some certain time. Similar to [25,34], we also assume the n-square matrix (β ij (x)) n×n is nonnegative and irreducible.
Because the last equation of model (1) is decoupled from other equations, we indeed need to study the following subsystem of (1): 2 Complexity e rest of this paper is organized as follows. In Section 2, some preliminaries are introduced for the well-posedness of the model. In Section 3, we define the basic reproduction number R 0 . In Section 4, the threshold dynamics are established in terms of R 0 . An special case is performed as a supplementary to the theoretical results in Section 5. A brief conclusion ends the paper.
en, we give the existence of solutions of model (5).
admits a global compact attractor.
Proof. Suppose to the contrary that τ ∞ < ∞, then ‖u(t, ·, ϕ)‖ ⟶ + ∞ as t ⟶ τ ∞ by eorem 2 in [39]. It follows from the first equation of model (5) that By the comparison principle and Lemma 2 in [40], there Furthermore, similar procedure can be applied to the second equation of model (5), and then, there exists a constant Hence, from the third equation of model (5), we have Now, we consider the following comparison system: It follows from the standard Krein-Rutman theorem (see [41]) that the eigenvalue problem of system (13) admits a principle eigenvalue λ with a strongly positive eigenfunctionφ 2 � (φ 21 ,φ 22 , . . . ,φ 2n ).
us, system (13) has a solution ςe λtφ 2 (x) for t ≥ 0, where ς is a positive constant, satisfying ςφ 2 ≥ (I 1 (0, x), I 2 (0, x), . . . , I n (0, x)) for x ∈ � Ω. By using the comparison principle, we have us, there exists a constant M 3 such that which leads to a contradiction. Hence, the global existence of u(t, ·, ϕ) is derived. Now, we are in the position to show the solution is also ultimately bounded. Indeed, it follows from the comparison principle, (11), and Lemma 2 in [40] that there exist t 1 > 0 and Moreover, from the second equation of model (5) and applying similar procedures, we can find A 2 > 0 and t 2 > 0, satisfying us, there exist t 3 > 0 and A 3 > 0 such that From Chapter 5 in [42], one obtains Since φ i j (x) is uniformly bounded, there exists constant κ 3i > 0 such that Γ 3i (t, x, y) ≤ κ 3i j≥1 e τ i j t for t > 0. Moreover, assume π i j are eigenvalues of ▽ · (d � 3i Complexity subject to Neumann boundary condition, which satisfies By eorem 2.4.7 in Wang [43], one gets π i j ≥ τ i j for all j ∈ N + . Since π i j decreases like − n 2 , there exists κ 3 > 0 such that Lett � max t 1 , t 2 , t 3 . For all t ≥t, by (9) and (4), we have which yields that us, the above discussion implies that the system (5) is point dissipative. Furthermore, by eorem 2.2.6 in [44], the solution semiflow Φ(t) is compact for any t > 0. erefore, it follows from eorem 3.4.8 in [45] that Φ(t) has a global compact attractor in X + .

Basic Reproduction Number
For each 1 ≤ i ≤ n, consider the following subsystem of model (5): It follows from the Lemma 2.2 in [40] that the following system admits a unique positive steady state S 0 i (x) which satisfies the equation with zS 0 i (x)/z] � 0 for x ∈zΩ, which is globally asymptotically stable in C( � Ω, R + ). Hence, the second equation of system (22) is asymptotic to By Lemma 2.2 in [40] and Corollary 4.3 [46], there exists a globally asymptotically stable steady state V 0 i (x). erefore, concluding from the above discussion, we know that model (5) admits a unique disease-free steady state and 0 � (0, 0, . . . , 0). Furthermore, if all the parameters of model (5) are positive constants, then we have S 0 (5) at E 0 , we obtain the linearized system Substituting I i (t, x) � e λt ϕ 3i (x) into (26), we obtain From eorem 7.6.1 in [38], we have the following result.

Lemma 2. e eigenvalue problem (27) admits a principle eigenvalue λ 0 with a strictly positive eigenfunction.
Denote T 3 � (T 31 , T 32 , . . . , T 3n ). Assume that the distribution of initial infection is ϕ 3 (x) ∈ Y + . en, F(x)T 3 ϕ 3 (x) is the distribution of new infective part as time evolves. erefore, we use (28) to describe the total distribution of the new infective numbers produced during the infection period. According to [14], the basic reproduction number is defined by R 0 � r(L), where r(L) is the spectral radius of the operator L. Furthermore, following eorem 3 in [14], we have the following lemma.

Lemma 3.
e principle eigenvalue λ 0 and R 0 − 1 have the same sign, and the disease-free steady state E 0 (x) is locally asymptotically stable.

Extinction/Persistence Result
Based on the discussion above, we now investigate the extinction and persistence of the disease. Theorem 2. If R 0 < 1, then the disease-free steady state E 0 (x) is globally asymptotically stable.

Complexity
By Lemma 4 and eorem 4.7 in [49], we know that Φ(t) has a positive steady state E * (x) of model (5). e proof is complete.

Global Dynamics for a Special Case
In this section, we consider a special case with all the parameters of model (5) as constants, except for the diffusion coefficients. en, model (5) degenerates into the following model: with the homogeneous Neumann boundary conditions (2) and the initial conditions (3). We mainly focus on the global stability of the endemic steady states of model (48). e proofs of the main results applying the Lyapunov functions and a subtle grouping technique guided by graph theory were developed in [28][29][30].
It follows from the Lemma 2 in [13] that model (48) has a disease-free steady state en, it follows from [26] that the basic reproduction number is defined as the spectral radius of FV − 1 , i.e., R 0 � r(FV − 1 ). As a consequence of eorems 2 and 3, we can obtain the following results without proofs.

Theorem 4.
e disease-free steady state E 0 of model (48) is globally asymptotically stable when R 0 < 1.
uniformly for all x ∈ � Ω. Moreover, model (48) has at least one endemic steady state.
Theorem 6. If R 0 > 1, then the endemic steady state E * is globally asymptotically stable.
Proof. Define which is a Laplacian matrix whose column sums are zero [30]. en, Θ is irreducible because (β ij ) n×n is irreducible. It follows from Lemma 1 in [29] that the solution space of the linear system Θζ � 0 has dimension 1 with a base ζ � (ζ 1 , ζ 2 , . . . , ζ n ) with ζ i � c ii , where c ii > 0 is the cofactor of the i-th diagonal entry of Θ. Define where Differentiating H i (t), we obtain It follows from [50] that Complexity en, using steady state equations (52), we have