Global dynamics for a class of reaction–diffusion multigroup SIR epidemic models with time fractional-order derivatives

. This paper investigates the global dynamics for a class of multigroup SIR epidemic model with time fractional-order derivatives and reaction–diffusion. The fractional order considered in this paper is in (0 , 1] , which the propagation speed of this process is slower than Brownian motion leading to anomalous subdiffusion. Furthermore, the generalized incidence function is considered so that the data itself can ﬂexibly determine the functional form of incidence rates in practice. Firstly, the existence, nonnegativity, and ultimate boundedness of the solution for the proposed system are studied. Moreover, the basic reproduction number R 0 is calculated and shown as a threshold: the disease-free equilibrium point of the proposed system is globally asymptotically stable when R 0 (cid:54) 1 , while when R 0 > 1 , the proposed system is uniformly persistent, and the endemic equilibrium point is globally asymptotically stable. Finally, the theoretical results are veriﬁed by numerical simulation.


Introduction
As we all know, mathematical models play an important role in researching the dynamical behavior of infectious diseases. In the classical epidemic model, it is generally considered that individuals are completely mixed, and everyone has the same possibility of infection. However, due to the differences in age, geographical distribution, and other factors, it is more realistic to divide the total population into several different populations, that is, to establish a multigroup epidemic model. Lajmanovich et al. first proposed the SIS multigroup systems and researched the stability of the endemic equilibrium point [11]. Subsequently, there are many research efforts devoted to investigating the importance of multigroup epidemic models [6,8,14]. Guo et al. were the first to successfully establish the complete global dynamics of the multigroup epidemic model based on the basic reproduction number [7]. Boosted by the work of Guo et al., many researchers discussed the stability of various multigroup systems [3,15,20,21,25].
Meanwhile, individual diffusion behavior is widespread in the actual propagation of infectious diseases. With the development of global transportation, individuals in incubation period can easily travel from one place to another, which is thought to be one of the main reasons of the global pandemic of infectious diseases. For instance, SARS first appeared in China's Guangdong Province in November 2002 and then quickly spread to other parts of China and even the world [26]. Also, COVID-19 was first detected at the end of December 2019 with successive cases occurring worldwide. Therefore, in order to better understand the impact of population mobility on the spread of infectious diseases, it is necessary to incorporate human movement into epidemic model to provide more theoretical guidance for epidemic control. Li et al. analyzed the stability and the uniform persistence of a SIRS epidemic model with diffusion [12]. Xu et al. studied the stability and the existence of traveling wave solutions of a SIS epidemic model with diffusion [28]. Recently, many diffusive epidemic models have been used to model within-group and inter-group interactions in spatially environments, for example, Wu et al. investigated a multigroup epidemic model with nonlocal diffusion and obtained the asymptotic behavior of traveling wave fronts [27].
It is worth noting in real life that the spread of infectious diseases not only depends on its current state, but also on its past state. Actually, it can be achieved that current state of fractional-order epidemic models depends on the past information since any fractional derivative contains a kernel function [30]. Furthermore, Smethurst et al. found that the patient waits for the doctor's time to follow a power law model P [J n > t] = Bt −α [24]. More importantly, Angstmann et al. proposed a infectivity SIR model with fractionalorder derivative, and they showed how fractional-order derivative arise naturally by continuous time random walk [2]. As generalized of classical integers ones, Hethcote firstly proposed a fractional-order SIR model with a constant population [8]. Then Almeida et al. considered the local stability of two equilibrium points of a fractional SEIR epidemic model [1].
Typically, the reaction term describes a birth-death reaction occurring in a habitat or reactor. The diffusion term simulates the movement of the individual in the environment in real-world applications. The diffusion is often described by a power law x 2 (t) − x(t) 2 ∼ Dt α , where D is the diffusion coefficient, and t is the elapsed time. In normal diffusion, the order α = 1. But if α > 1, particle undergoes superdiffusion, which mainly describes the process of active cell transport; if α < 1, this phenomenon is called subdiffusion, which can be the diffusion of proteins within cells or the diffusion of viruses between individuals [29]. And it results in a Caputo time-fractional reactiondiffusion system with fractional order α < 1. Meanwhile, it is pointed out in [19] that long waiting times model particle sticking, and the density of this process spreads slower than normal diffusion. Also, as shown in [19], Caputo time-fractional reactiondiffusion curve has a sharper peak and heavier tails, which can be used to describe the ability to control the transmission of the disease when only a small number of people are infected, such as COVID-19. The study of subdiffusion system has attracted widespread interest in recent years. Mahmoud et al. studied the Cauchy problem of the fractionalorder evolution equation and obtained the expression of the solution of the time fractionalorder reaction diffusion system [18]. The subdiffusive predator-prey system is discussed, and the analytical solution of the system is studied in [29]. However, few works have been devoted to studying the subdiffusion epidemic model. Motivated by this, in this work, we focus on time-fractional reaction-diffusion epidemic system, which means the spread of infectious diseases is slower than a Brown motion.
Based on the above discussion, the dynamics of the multigroup SIR epidemic model with generally incidence rates is investigated in this paper. Particularly, the susceptible individuals, infective individuals, and recovered individuals are assumed to follow Fickian diffusion.
The organization of this paper is as follows. A class of diffusive SIR epidemic model with time fractional-order derivatives is formulated and some preliminaries are introduced in Section 2. In Section 3, global dynamics of the proposed model are studied, and numerical simulations are presented to illustrate theoretical results in Section 4. Finally, a brief discussion is given in Section 5.

Model development
Before presenting a class of multigroup reaction-diffusion SIR epidemic model with time fractional-order derivatives, some necessary preliminaries are presented.

Preliminaries
This section begins with some notations, definitions, and results.
and φ i ∈ Y (i = 1, 2); X 0 be a open set of X such that X = X 0 ∪ ∂X, where ∂X is the boundary of X; X + = Y + × Y + be the positive cone of X.

Proof. Define the Lyapunov function
Calculating the fractional derivative of V (u) along the trajectories of system (2), one has and V 1 = 0 if and only if u = u * . Then according to [4], there exists a unique global asymptotic stability of constant equilibrium u * = b/µ for system (2).

System description
In [10], Korobrinikov et al. studied a multigroup SIR model as follows: https://www.journals.vu.lt/nonlinear-analysis But individual movement is not be considered in system (4) that is unrealistic, then Wu et al. considered the following SIR epidemic model with diffusion [27]: Based on the previous analysis, since fractional order has the long-term memory, which can describe the spread of infectious diseases more accurate. In addition, it is traditionally assumed that the incidence of disease transmission is bilinear with respect to the number of susceptible individuals and the number of infected individuals. But in reality, it is often difficult to obtain detailed information on the spread of infectious diseases because they may change with the surrounding environment. Therefore, the general incidence rates will be chose in this paper. Motivated by the above work, as an extension of system (5), a class of multigroup SIR epidemic model are investigated as follows: x ∈ Ω, k = 1, 2, . . . , n, where C 0 D α t implies Caputo fractional-order operator (0 < α 1); ∆ = ∂ 2 /∂x 2 denotes the Laplace operator; ∂/∂υ denotes the outward normal derivative on the smooth boundary ∂Ω; S k (x, t), I k (x, t), and R k (x, t) represent the number of the susceptible, infective, and recovered individuals in kth group at time t and spatial location x, respectively; µ ik (i = 1, 3) imply the nature death rates of S k and R k in kth, respectively; µ 2k denotes the disease-related death rates of I k in kth; b k represents the recruitment rate of the total population; r k implies the recovery rate of the infected individuals in kth group; d ik (k = 1, 2, 3) denotes the diffusion rate of S k , I k , and R k in kth group; β kj represents the infection rate of S k infected by I j . Furthermore, d ik , b k , µ ik (i = 1, 2, 3), β kk and r k are positive constants for k = 1, 2, . . . , n, and β kj (k = j) are nonnegative constants for k = 1, 2, . . . , n.
Before giving the main results, hypothesis in terms of generalized incidence rates f k (S k ) and g k (I k ) is made as follows: (H) (i) f k (S k ) and g k (I k ) satisfy the local Lipschitz condition and f k (0) = 0, g k (0) = 0 for k = 1, 2, . . . , n; Remark 1. Note that under hypothesis (H), many existing models can be regarded as a special form of system (6), such as , and other nonlinear incidence rate in [16].

Model analysis
Some dynamical behavior of system (6) are investigated in this section. Here it can be found that the susceptible class S k and the infected class I k are not effected by the recovered class R k of system (6). Hence, we will focus our attention on the following reduced system: Then some basic properties of system (7) are discussed in following parts.

Nonnegative and boundedness
It is significant to demonstrate the existence, uniqueness, and boundedness of a nonnegative solutions for system (7) before implementing its stable process. Thus, this subsection moves to the discussion of proprieties mentioned above.
Proof. Consider these two function: According to condition (i) in hypothesis (H), it is obvious that F 1k and F 2k are mixed quasimonotonous.
Consider the following auxiliary system: It is obvious that (S k , I k ) = (0, 0) is a pair of the lower solution to system (7). Then, according to Lemma 1, one has S k 0 and I k 0. Furthermore, the following auxiliary system is introduced: then the above system (8) has a solution as follows: Therefore, lim sup t→∞ S k (t) = b k /µ 1k , then there exists a constant T 0 satisfied S k (x, t) b k /µ 1k for t > T 0 . Further, consider the following auxiliary system: then the solution for the above system (9) is .
However, it is easy to see that It can be deduced from Eqs. (10) and (11) that Similar to the above analysis, it can be obtained the following equation: Based on the above analysis and Lemma 2, system (7) has a unique nonnegative global solution. Furthermore, the expression for the solution of system (7) is with ζ α (θ) represents a probability density; G 1k (t) t 0 represent generated strong continuous operator semigroups by d 1k ∆; G 2k (t) t 0 , denoting generated strong continuous operator semigroups by d 2k ∆ − µ ik − r k (k = 1, 2, . . . , n), can be rewritten by [23] G 2k (t)φ 2k = Ω T 2k (t, x, z)φ 2k (z) dz, https://www.journals.vu.lt/nonlinear-analysis where T 2k is the Green function yielded Hence, by the boundedness of the eigenfunction ϕ k j (x) one has According the upper solution S k of system (7), one has lim sup t→∞ S k = b k /µ 1k , which implies S k is ultimate bounded. Further, the ultimate bounded of I k will be analyzed. Let N k = S k + I k and P k = Ω N k dx. Adding the first two equations of system (7) and integrating it on Ω, one has Therefore, by [13] one has then there exist two constants M > 0 and T 1 > 0 satisfying P k M for t T 1 . According to [17], the operator families {T 2k α } is uniformly bounded. Hence, there exist two constants M 1 > 0 and T 2 > T 1 > 0 satisfying T 2k α M 1 for t T 2 . Finally, the uniformly boundedness of the infected group I k (x, t) can be studied as follows: where S 0 k = b k /µ 1k , thus I k is ultimate bounded. Therefore, there exists a unique positive global solution (S(x, t), I(x, t)) of system (7), and it is also ultimately bounded.

Lemma 5. The basic reproduction number
Proof. Linearizing system (7) at the disease-free equilibrium point E 0 , one has Let F = (β kj f k (S 0 k )c j ) n×n and V = diag (µ 21 + r 1 , . . . , µ 2n + r n ). Then it is easy to find M (S 0 ) = V −1 F . Obviously, we have R 0 = ρ(V −1 F ). By the definition of the basic reproduction number [5] one has R 0 = ρ(F V −1 ). Thus, by the properties of matrix eigenvalues it can be deduced that R 0 = R 0 . Therefore, R 0 is considered as a threshold parameter in place of R 0 . In the following, the uniqueness and the global stability of E 0 = (S 0 1 , 0, S 0 2 , 0, . . . , S 0 n , 0) are studied.
Calculating the time fractional derivative of L 1 along the trajectories of system (7), one has Let L 2 (I) = Ω (1 − ρ(M (S 0 ))αI dx, which is a positive definition function in Γ . Then it is concluded from [4] that E 0 is globally asymptotically stable in domain Γ .

Theorem 3. Under hypothesis (H) and
where δ > 0 is a constant.
Let Y (t, x) = (S, I) be the solution of system (7) under the initial value Y 0 = Y (0, x) = (φ 1 (x), φ 2 (x)) ∈ X + . For any t 0, it can be known that all nonnegative solutions (S(t, φ), I(t, φ)) generate a solution semiflow T (t) : , and it is obvious that T (0) = E, where E is the identity matrix. It can be deduced from Lemma 4 that Then T (t + s) = T (t)T (s). Based on the above analysis, one has T (t) is C 0 -semigroup on X + . Obviously, T (t) is compact for t 0 and point dissipative in X + . The following system is considered: It can be found from Lemma 3 that S 0 k is globally asymptotically stable. Thus, system (7) is globally asymptotically stable at the disease-free equilibrium point E 0 = (S 0 1 , 0, . . . , S 0 n , 0). It can be deduced that the disease-free equilibrium E 0 in ∂X + is a global attractor of T (t), which implies Ω( there exists a solution (S k , I k ) of system (7) with the initial value such that (S k , I k ) → (S 0 k , 0) (k = 1, 2, . . . , n) as t → ∞, then there exists a constant τ > 0 such that S k > S 0 k − ε 0 and I k > ε 0 for t τ . Since B is irreducible, M (S 0 , ε 0 ) is irreducible. Let (α 1 , α 2 , . . . , α n ) be the positive left eigenvector of M (S 0 , ε 0 ) corresponding to the spectral radius ρ(M (S 0 , ε 0 )), that is, Define the following arbitrary function: Calculating the time fractional derivative of L 1 along the trajectories of system (7), one has which leads to a contradiction with lim t→∞ I k (t) = 0. Therefore, W s (E 0 ) ∩ X 0 = ∅. Thus, it can be deduced from [25] that T (t) is uniformly persistent. It is concluded that system (7) is uniformly persistent.
The ultimate boundedness and the uniform persistence imply the existence of a positive equilibrium point of system (7). Therefore, the existence and global stability of the positive endemic equilibrium point of system (7) can be further discussed.

Remark 3.
Not considering infection between populations, that is, when β kj = 0 (k = j), the reproduction number R k 0 of group k is R k 0 = β kk c k f k (S 0 k )/(µ 2k + r k ). Furthermore, the disease-free equilibrium point E 0 k = (b k /µ 1k , 0) is globally asymptotically stable when R k 0 1, and the endemic equilibrium point E * k = (S * k , I * k ) is globally asymptotically stable when R k 0 > 1.

Numerical simulations
In order to verify theoretical results numerically, numerical simulations are presented in this section. We consider system (7) with two-group case, which is suitable for infectious diseases transmitted between two cities or communities. Furthermore, system (7) with two groups (n = 2) can be calculated by the central difference method in L 1 -type space and Alikhanov-type discretization in time [9]. Furthermore, we consider the following incidence rate as an example: f k (S k ) = S k , g k (I k ) = I k /(1 + τ I k ), which τ is a positive parameter measuring the psychological or inhibitory effect. Obviously, f k (S k ) and g k (I k ) satisfy hypothesis (H). The corresponding system can be expressed as C 0 D α t I 2 = d 22 ∆I 2 + β 21 S 2 I 1 1 + τ I 1 + β 22 S 2 I 2 1 + τ I 2 − (µ 22 + r 2 )I 2 .
Let assign the following values to the parameters of system (23) It is easy to calculate that R 0 = 0.1549 < 1. Based on by Theorem 2, the disease-free equilibrium point E 0 of system (23) is global stable which is verified by Figs. 1 and 2. Further, the following parameters is chose: System (23) has a unique equilibrium point E * = (S * 1 , I * 1 , S * 2 , I * 2 ). It can be calculated that R 0 = 1.0101 > 1, and Eq. (14) is satisfied. Based on the above analysis, the endemic equilibrium point E * of system (23) is global stable, which is verified by Figs. 3 and 4.
Further, with regard to the disease-free equilibrium point of the first group, the influence of different fractional order α on the stability of the infected are discussed. The error   images of the infected of α = 0.4 and α = 0.75, α = 0.75 and α = 0.98 are described in Fig. 5, respectively. It is easy to seen from Fig. 5 that although the infected will disappear, different order α will have a sensitive effect on the change of solution. Further, when α tends to 1, the numerical solutions of system (7) are also convergent to the solutions of the classical ones [17]. But the relationship between the change of the solution for system (7), and fractional order α is not discussed.

Discussion
In this article, incorporating the population diffusion and time fractional-order derivatives, theory analysis of a class of multigroup SIR epidemic model are investigated. Firstly, the existence and uniqueness of the nonnegative solution for system (7) are established. By using Lyapunov functions the global stability of the disease-free equilibrium point E 0 is obtained when the basic reproduction number R 0 1. Besides, when R 0 > 1, the uniform persistence and the global stability of the endemic equilibrium point E * are discussed. The proposed model, a more accurate epidemic model, can help us to understand some dynamical behaviors of infectious diseases. Moreover, theoretical results may provide some useful guidance for making effective countermeasures on infectious diseases. However, the relationship between system (7) and fractional order α is still an open question, which will be our future work.