key: cord-1001140-di0pnmw5
authors: Avila-Vales, Eric; Pérez, Ángel G. C.
title: Dynamics of a reaction–diffusion SIRS model with general incidence rate in a heterogeneous environment
date: 2021-11-17
journal: Z Angew Math Phys
DOI: 10.1007/s00033-021-01645-0
sha: c00eb6e2c445e3e3aafe88c2ce15dec5803ff97c
doc_id: 1001140
cord_uid: di0pnmw5

In this paper, we study a diffusive SIRS-type epidemic model with transfer from the infectious to the susceptible class. Our model includes a general nonlinear incidence rate and spatially heterogeneous diffusion coefficients. We compute the basic reproduction number [Formula: see text] of our model and establish the global stability of the disease-free steady state when [Formula: see text] . Furthermore, we study the uniform persistence when [Formula: see text] and perform a bifurcation analysis for a special case of our model. Some numerical simulations are presented to illustrate the dynamics of the solutions as the model parameters are varied.

The spread of infectious diseases within human populations is an important topic of research in mathematical modelling that has become widely studied in recent times. Many different models, such as the SI, SIS, SIR and SIRS types, have been proposed to better understand the underlying transmission mechanism of infections. In these models, the population under study is divided into several classes or subpopulations, such as susceptible (S), infected (I) and recovered (R).

The incidence rate of a disease, which measures how many people become infected by that disease per time unit, plays an important role in the dynamics of epidemic models. Traditionally, the incidence rate has been assumed to be bilinear with respect to the number of susceptible individuals (S) and the number of infected individuals (I). However, many other functions have been proposed to model the transmission of infectious diseases. In general, the incidence function can be written in the form F (S, I), where the function F satisfies some common properties, such as the following: • F (S, I) = βSI (bilinear) [1] ; • F (S, I) = βSI/(1 + a 1 I) (saturated with respect to infectives) [2] ; • F (S, I) = βSI/(1 + a 2 S) (saturated with respect to susceptibles) [3] ;

• F (S, I) = βSI/(1 + k 1 I + k 2 S) (Beddington-DeAngelis) [4] ; • F (S, I) = βSI/ (1 + k 1 I)(1 + k 2 S) (Crowley-Martin) [5] ;

Given the wide variety of forms that have been proposed for the incidence rate, some authors have opted to study epidemic models that include a general class of incidence functions. This has the advantage that the results proved for particular models can be generalized to a broader class of models, and the researchers can focus on other features, such as spatial heterogeneity that could yield more complicated dynamics or bifurcations.

A SIRS model with transfer from the infected to the susceptible class was studied in [9] using the system of ordinary differential equations dS dt = Λ − μS − Sf (I) + γ 1 I + δR,

where S(t), I(t) and R(t) denote the number of susceptible, infected, and recovered individuals at time t, respectively, and the parameters are interpreted as follows:

• Λ: recruitment rate of susceptible population.

• μ: natural death rate.

• α: disease-induced death rate.

• γ 1 : transfer rate from the infected class to the susceptible class.

• γ 2 : transfer rate from the infected class to the recovered class.

• δ: rate of immunity loss.

The authors in [9] considered an incidence function of the form Sf (I) and proved that the threshold dynamics for model (1.1) is completely determined by the basic reproduction number, which is denoted as R 0 . Variants of such model were later studied in [10] and [11] using more general incidence rates, which take the form f (S, I). The authors in [10, 11] showed that the generalized model retains its threshold dynamics under certain conditions of the parameters.

On the other hand, we should consider that the individuals of the population under study may move randomly in the space. Several studies on modelling of infectious diseases such as influenza [12] , cholera [13] , dengue [14] , brucellosis [15] and COVID-19 [16] have highlighted the importance of individual motility in the dynamics of an outbreak. Thus, it is appropriate to consider epidemic models based on reactiondiffusion equations, where the moving patterns for susceptible, infected and recovered individuals are modelled through the diffusion rate of each subpopulation.

Moreover, many communicable diseases occur in a heterogeneous environment due to the differences in environmental conditions such as humidity, temperature or the varying availability of medical resources. This has led some researchers to study epidemic models where some of the parameters depend on a spatial variable x. In particular, the transmission rate of the disease may be represented by a function that depends not only on the number of susceptible and infected individuals, but also on the spatial location. For example, a SIRS-type model with diffusion and heterogeneous parameters was studied in [12] with the incidence function f (x, I)S. Other reaction-diffusion models have been proposed in [17] using the general incidence function β(x)f (S, I) and in [13] using the function f (x, S, I), which take into account both the dependence on the spatial variable and the non-linearity with respect to S and I. A diffusive epidemic model based on system (1.1) was studied by Yang et al. in [18] , where the authors established the global attractiveness of the disease-free equilibrium when R 0 < 1 and the persistence of the disease when R 0 > 1. Although the authors in [18] included a general incidence rate of the form β(x)f (S, I) and spatially variable parameters, they assumed that the diffusion coefficients for the susceptible, infected and recovered individuals are all equal to a constant D. In reality, there could be differences in the motility patterns of these subpopulations due to the individuals changing their behaviour when they get the disease, and the diffusion rate may also vary with the spatial location.

Based on the above discussion, we propose here a diffusive version of the SIRS model (1.1) that includes three different diffusion coefficients with spatial heterogeneity. The model to be studied is given by the following system of partial differential equations: (1.2) subject to the initial conditions

Here, the variables S(x, t), I(x, t) and R(x, t) represent the number of susceptible, infected, and recovered individuals, respectively, at position x and time t, while ∇ is the gradient operator. We assume that the domain Ω is a connected, bounded subset of R n with smooth boundary ∂Ω. The positive functions d S (·), d I (·) and d R (·) denote the spatially heterogeneous diffusion coefficients for each subpopulation. The spatially dependent parameters Λ(·), μ(·), α(·), γ 1 (·), γ 2 (·) and δ(·) have the same meaning as for model (1.1); for biological reasons, they are assumed to be strictly positive (except γ 1 , which is nonnegative) and uniformly bounded on Ω. Furthermore, for the model to be well-posed, we assume that Λ(·), μ(·), α(·), γ 1 (·), γ 2 (·), δ(·) ∈ C 0 (Ω, R).

Based on (A1')-(A4'), we make the following assumptions on the incidence function F (x, S, I) for model (1.2):

∂S (x, S, I) > 0 and ∂F1 ∂I (x, S, I) ≤ 0 for all x ∈ Ω, S, I ≥ 0. (A4) There exists a Hölder continuous function β : Ω → R + such that F (x, S, I) ≤ β(x)SI for all x ∈ Ω, S, I ≥ 0.

The rest of this paper is organized as follows. First, in Sect. 2, we prove the existence of bounded global solutions for our model. Next, in Sect. 3, we compute the basic reproduction number and study the stability of the disease-free steady state. In Sect. 4, we study the persistence of the model and the existence of endemic equilibria. In Sect. 5, we perform a bifurcation analysis for a special case of our model. In Sect. 6, we carry out some numerical simulations. Finally, we provide some concluding remarks in Sect. 7.

Let X := C(Ω, R 3 ) be the Banach space with the supremum norm · , and define X + := C(Ω, R 3 + ). Then, (X, X + ) is a strongly ordered space. We will now prove the existence of unique global solutions for our model. For that, we introduce the notation

where ψ is any of Λ, μ, α, γ 1 , γ 2 , δ or β.

subject to the Neumann boundary condition. Then,

where T i (x, y, t) represents the Green function associated with ∇·(d i (·)∇)−π i (·) subject to the Neumann boundary condition. By [19, Corollary 7.2.3] , it follows that Γ(t) := Γ 1 (t), Γ 2 (t), Γ 3 (t) is strongly positive and compact for each t > 0. Following [20, 21] , we can see that there exist constants M i > 0 (i = 1, 2, 3) such that

where α i < 0 is the principal eigenvalue of ∇·(d i (·)∇)−π i (·) subject to the Neumann boundary condition. For every initial value function

Then, model (1.2) can be rewritten as the integral equation

where U (t) = S(t), I(t), R(t) T . It is easy to see that the subtangential condition in [22, Corollary 4] is satisfied. Thus, the model has a unique positive solution S(·, t; φ), I(·, t; φ), R(·, t; φ) on [0, τ), where 0 < τ ≤ +∞.

We will now prove that the local solution can be extended to a global one, i.e. τ = +∞. Suppose by contradiction that τ < +∞. By the theory of abstract functional differential equations (see [22, Theorem 2] ), we know that

Hence, it suffices to prove that the solution is bounded on Ω × [0, τ). To this end, we define

By the divergence theorem [23, Theorem 3.7] and the homogeneous Neumann boundary conditions, we have Thus, d dt

By the comparison principle, there exists a constant N 1 > 0 and t 1 > 0 such that

Next, denote by τ i j the eigenvalue of ∇ · (d i (·)∇) − π i (·) subject to the Neumann boundary condition corresponding to the eigenfunction ϕ i j (x), such that τ i 

Since ϕ i j is uniformly bounded, there exists ω > 0 such that

For i = 1, 2, 3, define

and −∇ · (d i (·)∇) + π i (·), respectively, subject to the Neumann boundary condition.

For all x ∈ Ω and z ∈ R, we have 

Using (2.2) together with (2.1) and (2.4), we obtain

Hence,

From the second equation of (1.2) and assumption (A4), we get

Hence, we obtain

Similarly, the third equation of (1.2) yields

and thus,

Dynamics of a reaction-diffusion SIRS model with general... Page 7 of 23 9 Hence, we conclude that S, I and R are bounded on Ω × [0, τ), which contradicts (2.3). This proves that τ = +∞.

In a similar way, we can obtain the following corollary on the boundedness of solutions on [0, ∞).

Furthermore, the solution semiflow Φ t : X + → X + is point dissipative and has a global compact attractor.

Proof. Using the same notation as in the proof of Theorem 2.1 and replacing τ by +∞, we can show that (2.5) . From this, it follows that the system is point dissipative. In addition, we know by [ 

In this section, we will study the disease-free dynamics of model (1.2) and determine the basic reproduction number R 0 , which is defined as the average number of secondary infections generated by a single infected individual introduced in a completely susceptible population.

Setting

, we obtain the following equation for the density of susceptible population in absence of the disease:

System (3.1) admits a unique positive steady state S 0 (x), which is globally asymptotically stable (see [20, Lemma 1] ). We will call E 0 := (S 0 (·), 0, 0) ∈ X + the disease-free steady state of model (1.2). Linearizing the equation for the infected population in (1.2) around E 0 , we obtain

Substituting I(x, t) = e λt ϕ(x) in the above equation yields the following eigenvalue problem:

It follows from the standard Krein-Rutman theorem that (3.3) has a principal eigenvalue

corresponding to a positive eigenfunction ϕ * (x), where s(A) denotes the spectral bound of a closed linear operator A.

It is well-known that λ * (S 0 ) is given by the variational characterization

Page 8 of 23 E. Avila-Vales and Á. G. C. Pérez ZAMP Let T 2 (t) be the semigroup on C(Ω, R) associated with ∇ · d I ∇I − (μ + α + γ 1 + γ 2 ). In order to define the basic reproduction number for our model, we assume that the population is near the diseasefree steady state E 0 , and introduce infected individuals at time t = 0, where the spatial distribution of infected population is described by φ 2 (x). Then, T 2 (t)φ 2 (x) is the distribution of infected population as time evolves. It follows that the distribution of new infections at time t is ∂F ∂I (x, S 0 , 0)T 2 (t)φ 2 (x). Thus, the distribution of total new infections becomes

where L is the next infection operator, which maps the initial distribution φ 2 of infected individuals to the distribution of total infected individuals produced during the infection period. Following [28] [29] [30] [31] , we define the basic reproduction number for model (1.2) as the spectral radius of L, that is,

Thus, model (1.2) can be rewritten as

The disease-free steady state is given by u 0 (x) = 0, S 0 (x), 0 , with the variables ordered as (I, S, R).

We can immediately verify that assumptions (A1)-(A4) in [31] hold. Moreover, if we define

then M 0 (x) and −V (x) are cooperative matrices for all x ∈ Ω, i.e. their off-diagonal elements are nonnegative. Hence, assumptions (A5) and (A6) in [31] also hold, so we can conclude the following result by [31, Theorem 3.1] .

Before proving the main result of this section, we give the following lemma.

(ii) If there exists t 0 ≥ 0 such that I(·, t 0 ; φ) ≡ 0 (respectively, R(·, t 0 ; φ) ≡ 0), then, for t > t 0 , we have I(·, t; φ) > 0 (respectively, R(·, t; φ) > 0).

Dynamics of a reaction-diffusion SIRS model with general... Page 9 of 23 9 Proof. By assumption (A4) and Corollary 2.2, we have F (x, S, I) ≤ β + SI ≤ β + M I S for all x ∈ Ω, S, I ≥ 0. Then, from the first equation of (1.2), we get

By the comparison principle, this shows that lim inf t→∞ S(x, t; φ) ≥ Λ − / (μ + + β + M I ) uniformly for x ∈ Ω.

From system (1.2), we can also get

and

Using the strong maximum principle (see [32, Proposition 13 .1]), part (ii) of the lemma is concluded.

Next, we give the main result concerning the stability of E 0 in terms of the basic reproduction number. Proof. (i) By Lemma 3.1, we have λ * (S 0 ) < 0 when R 0 < 1. Thus, there exists a sufficiently small ε such that λ * (S 0 + ε) < 0. According to Corollary 2.2, there exists a τ 1 > 0 such that

It follows from (A1)-(A3) that

for t ≥ τ 1 . Hence, by the second equation of (1.2), we have

Let ψ(x) be the eigenfunction corresponding to the principal eigenvalue λ * (S 0 + ε) < 0. Let ξ 1 > 0 be such that I(x, τ 1 ) ≤ ξ 1 ψ(x). By the comparison principle, we get

This yields that lim t→∞ I(x, t) = 0 uniformly for x ∈ Ω. Then, the third equation of (1.2) is asymptotic to

so lim t→∞ R(x, t) = 0 uniformly for x ∈ Ω. Moreover, the first equation of (1.2) is asymptotic to

which implies that lim t→∞ S(x, t) = S 0 (x) uniformly for x ∈ Ω. This completes the proof. E. Avila-Vales and Á. G. C. Pérez ZAMP (ii) Assume by contradiction that there exists a positive solution of (1.2) such that lim sup t→∞ S(·, t), I(·, t), R(·, t) − (S 0 , 0, 0) < ε 0 .

Then, there exists t 1 > 0 such that S 0 − ε 0 < S(x, t) < S 0 + ε 0 and 0 < I(x, t) < ε 0 for t ≥ t 1 . It follows from (A1)-(A3) that

For any ε ∈ 0, min x∈Ω S 0 (x) , we consider the eigenvalue problem

Define R ε as the spectral radius of the operator

As a consequence, we can fix a small ε 0 ∈ ε, min x∈Ω S 0 (x) such that λ * ε0 > 0. By assumption, I(x, t) > 0 for all x ∈ Ω and t > 0. Then, by Lemma 3.2(ii), we can choose a sufficiently small number η > 0 such that I(·, t) ≥ ηφ * ε0 (·), where φ * ε0 (·) is a strongly positive eigenfunction corresponding to λ * ε0 . Notice that u 1 (x, t)

is a solution to the linear system

It then follows from (3.11) and the comparison principle that

Since λ * ε0 > 0, this implies that I(x, t) is unbounded, which is a contradiction.

We will now study the existence and persistence of the endemic equilibrium for model (1.2). 

Moreover, (1.2) admits at least one endemic equilibrium (S * , I * , R * ).

Dynamics of a reaction-diffusion SIRS model with general... Page 11 of 23 9 Proof. Let

By Lemma 3.2(ii), it follows that for any φ ∈ W 0 , we have I(x, t; φ) > 0 and R(x, t; φ) > 0 for all x ∈ Ω,

and let ω(φ) be the omega limit set of the orbit O + (φ) := {Φ t (φ) : t ≥ 0}. We will now show that

For any given φ ∈ M ∂ , we have Φ t (φ) ∈ ∂W 0 . It then follows that for each t ≥ 0, either I(·, t; φ) ≡ 0 or R(·, t; φ) ≡ 0. In the case where I(·, t; φ) ≡ 0, we can see that R(x, t; φ) satisfies

which implies that lim t→∞ R(·, t; φ) = 0 uniformly for x ∈ Ω. Thus, for any sufficiently small ε > 0, there is a t 2 > 0 such that R(x, t; φ) < ε for t ≥ t 2 . Then, from the first equation of (1.2), we can get

Since ε is arbitrary, this shows that lim t→∞ S(x, t; φ) = S 0 (x). On the other hand, when I(·, t 3 ; φ) ≡ 0 for some t 3 ≥ 0, we obtain by Lemma 3.2(ii) that I(x, t; φ) > 0 for all x ∈ Ω, t > t 3 . Hence, R(·, t; φ) ≡ 0 for all t > t 3 , but the last equation of (1.2) implies that I(·, t; φ) ≡ 0 for t > t 3 , which is a contradiction. Thus, we have proved that ω(φ) = {(S 0 (·), 0, 0)} for all φ ∈ M ∂ . Moreover, since R 0 > 1, it follows from Theorem 3.3(ii) that E 0 = (S 0 , 0, 0) is a uniform weak repeller for W 0 . Define a continuous function p :

It follows from Lemma 3.2(ii) that p −1 (0, ∞) ⊆ W 0 , and p has the property that if p(φ) > 0 or φ ∈ W 0 with p(φ) = 0, then p Φ t (φ) > 0 for all t > 0. Thus, p is a generalized distance function for the semiflow Φ t : X + → X + (see [33] ). Note that any forward orbit of Φ t in M ∂ converges to E 0 . Furthermore, the claim above implies that E 0 is isolated in X + and W S (E 0 ) ∩ W 0 = ∅, where W S (E 0 ) is the stable set of E 0 . Moreover, it is clear that there are no cycles in M ∂ from E 0 to E 0 . By [33, Theorem 3] , it follows that there exists a constant η 1 > 0 such that min{p(ψ) : ψ ∈ ω(φ)} > η 1 for any φ ∈ W 0 . Hence,

Thus, (4.1) holds by taking ε = min{η 1 , η 2 }, which proves the uniform persistence result. By Theorem 3.7 and Remark 3.10 in [34] , it follows that Φ t : W 0 → W 0 has a global attractor. It then follows from [34, Theorem 4.7] that Φ t has a steady state (S * , I * , R * ) ∈ W 0 , which is positive by Lemma 3.2(ii). Therefore, (S * , I * , R * ) is an endemic equilibrium of model (1.2). 

In this section, we will use γ 1 (the transfer rate from the infected class to the susceptible class) as the main bifurcation parameter to perform the bifurcation analysis of model (1.2). To apply the local and global bifurcation theory by Crandall and Rabinowitz [35] , we need to assume that the diffusion coefficient d I and the parameter γ 1 are constant, where d I > 0 and γ 1 ≥ 0. However, we still allow all other parameters and the diffusion rates d S (x) and d R (x) to be spatially variable.

A steady state of model (1.2) is a solution of the elliptic problem

When d I (x) = d I and γ 1 (x) = γ 1 are constant, system (5.1) becomes

It is easy to see that (S 0 (·), 0, 0) is a semi-trivial steady state solution of (5.2). Denote by γ * 1 the principal eigenvalue of the eigenvalue problem

associated with a positive eigenfunction ψ 0 (x), which is uniquely determined by the normalization max x∈Ω ψ 0 (x) = 1. By expression (3.5), we know that R 0 is decreasing with respect to γ 1 . Furthermore, (5.3) with γ 1 = γ * 1 is equivalent to system (3.3) with γ 1 (x) = γ * 1 , d I (x) = d I and λ = 0, and by Lemma 3.1, we have that λ * = 0 if and only if R 0 = 1. This shows that the condition γ 1 = γ * 1 is equivalent to R 0 = 1.

Next, we define a function

Then, it follows that γ * 1 = H when H is a constant. We will now study the case where H is not constant and it could change sign in Ω. Consider the eigenvalue problem with indefinite weight (ii) If

We next regard γ 1 as a bifurcation parameter and investigate a local branch of positive solutions of (5.2) that bifurcates from the branch of semi-trivial solutions (S 0 (·), 0, 0, γ 1 ) : γ 1 ≥ 0 . We note that S 0 (·) is independent of the parameter γ 1 .

S = (S, I, R, γ 1 ) ∈ X × R + : (S, I, R, γ 1 ) = (S 0 (·), 0, 0, γ 1 ) , be the set of semi-trivial solutions of (5.2), where ψ 0 > 0 is the principal eigenfunction of (5.3), and φ 0 , χ 0 satisfy

x ∈ ∂Ω.

Moreover, γ 1 (0) can be calculated by

Proof. We begin by making the change of variables u =S − S 0 , v =Ĩ, w =R, so the disease-free steady state (S,Ĩ,R) = S 0 (·), 0, 0 is shifted to the origin (u, v, w) = (0, 0, 0). Let Y = L p (Ω). Define a mapping G :

Page 14 of 23 E. Avila-Vales and Á. G. C. Pérez ZAMP Letf (u, v, w, γ 1 ) = f (u + S 0 , v, w, γ 1 ),g(u, v, w, γ 1 ) = g(u + S 0 , v, w, γ 1 ) andh(u, v, w, γ 1 ) = h(u + S 0 , v, w, γ 1 ). The Fréchet derivative of G with respect to (u, v, w) at (u, v, w, γ 1 ) = (0, 0, 0, γ 1 ) is given by

In particular,

We recall that the positive eigenfunction ψ 0 of (5.3) satisfies

we can see that (φ 0 , ψ 0 , χ 0 ) satisfies G (u,v,w) (0, 0, 0, γ * 1 )[φ 0 , ψ 0 , χ 0 ] = 0. Hence, (φ 0 , ψ 0 , χ 0 ) ∈ ker G (u,v,w) (0, 0, 0, γ * 1 ) \ {0}, so ker G (u,v,w) (0, 0, 0, γ * 1 ) is non-trivial. On the other hand, it follows from the definition of γ * 1 that the principal eigenvalue λ * of (3.3) equals zero when γ 1 = γ * 1 . Since λ * = 0 is a simple eigenvalue, the eigenfunction ψ 0 is unique up to constant multiples. Moreover, χ 0 and φ 0 are uniquely determined by (5.10) and (5.11) for a fixed ψ 0 . This implies that ker G (u,v,w) (0, 0, 0, γ * 1 ) = span {(φ 0 , ψ 0 , χ 0 )}. It is easy to see that (h 1 , h 2 , h 3 ) ∈ range G (u,v,w) (0, 0, 0, γ * 1 ) if and only if there exists (φ, ψ, χ) ∈ X such that 

Then, for such ψ and χ, the first equation of (5.12) has a unique solution

Dynamics of a reaction-diffusion SIRS model with general... Page 15 of 23 9 From this, it follows that the range of G (u,v,w) (0, 0, 0, γ * 1 ) is given by

The above equation defines a constraint on the variable h 2 . Since there are no constraints on h 1 or h 3 , the solution to system (5.12) has two degrees of freedom. This shows that dim range G (u,v,w) (0, 0, 0, γ * 1 ) = 2 and thus, codim range G (u,v,w) (0, 0, 0, γ * 1 ) = 1. Furthermore, we can calculate

where G (u,v,w),γ1 denotes the second partial Fréchet derivative of G with respect to the variables (u, v, w) and γ 1 . Similarly, we denote byf u,γ1 the second partial derivative off with respect to u and γ 1 , and so on.

Since

This allows us to apply the theorem of bifurcation from a simple eigenvalue [35] to conclude that the set of positive solutions to (5.2) near (S, I, R, γ 1 ) = S 0 (·), 0, 0, γ * 1 is a curve of the form (5.6). We also note that the possibility of other bifurcation points different from γ 1 = γ * 1 is excluded by virtue of the Krein-Rutman theorem. Furthermore, according to the direction of bifurcation formula by Shi [37] , we have

where the linear functional l : Y 2 → R is defined as

while G (u,v,w), (u,v,w) and G (u,v,w),γ1 denote the second partial Fréchet derivatives of G. By direct calculation, we can see that the second component of

Page 16 of 23 E. Avila-Vales and Á. G. C. Pérez ZAMP Thus,

where L is given by (5.9) . This completes the proof.

Next, we use the global bifurcation theory to give an extension of the local bifurcation branch (5.6). We define the set of positive solutions of (5.2) as Σ = (S, I, R, γ 1 ) ∈ X × R + : (S, I, R, γ 1 ) is a solution of (5.2)with S, I, R, γ 1 > 0 .

We can thus obtain the following result. Proof. For the local bifurcation branch obtained in Theorem 5.2, let Σ 1 be any maximal extension in X × R + as a connected set of solutions of (5.2). From the remarks in [38] , we can verify that all conditions in [38, Theorem 4.4] are satisfied. Therefore, Σ 1 must satisfy one of the following:

(i) it is not compact; or (ii) it contains a point (S 0 , 0, 0,γ 1 ) withγ 1 = γ * 1 ; or (iii) it contains a point (S 0 + U, V, W, γ 1 ) where (U, V, W ) = 0 and (U, V, W ) is in the complement of span {φ 0 , ψ 0 , χ 0 }. Since the eigenvalue of (5.3) with positive eigenfunction is unique, it is clear that (ii) cannot occur. Now, by Theorem 3.3(i), system (5.2) has no positive steady state solutions for R 0 < 1, that is, when γ 1 > γ * 1 . This implies that proj γ1 Σ 1 ⊂ [0, γ * 1 ]. On the other hand, due to Theorem 4.1, there exists a positive steady state for all γ 1 < γ * 1 (which corresponds to R 0 > 1). Therefore, we can conclude that proj γ1 Σ 1 = [0, γ * 1 ]. This result, together with the boundedness of solutions, shows that Σ 1 must be compact.

The above discussion implies that case (iii) holds. Then, by a standard argument with the global bifurcation theorem, we can find a certain (Ŝ,Î,R,γ 1 ) ∈ Σ 1 such that (Ŝ,Î,R,γ 1 ) is in the boundary of the set of positive solutions Σ. Since there are no equilibrium states withŜ = 0,Î = 0 orR = 0 other than the disease-free steady state, the only possibility isγ 1 = 0 andŜ,Î,R > 0. Hence, (Ŝ,Î,R) is a positive solution of system (5.2) with γ 1 = 0, which completes the proof.

In this section, we will perform some numerical simulations for system (1.2) to investigate the dynamics of the solutions as some of the parameters are varied. For simplicity, we consider the one-dimensional domain Ω = (0, 20) ⊂ R. We will use the parameter values Λ = 10, μ = 0.001, α = 0.01, γ 1 = 0.02, γ 2 = 0.07, δ = 0.02. where β = 0.0001 and the parameter a 1 is a positive function of the spatial variable. This allows us to study the case when the incidence function saturates more rapidly in some locations, due to the different measures taken by the population to eradicate the epidemic outbreak. It is readily shown that this form of incidence function satisfies (A1)-(A4). Throughout this section, we will use the initial conditions S 0 (x) = 9700, I 0 (x) = 10 exp(−x), R 0 (x) = 0 for x ∈ [0, 20].

We will first study the effect that heterogeneity in the saturation parameter has on the model dynamics.

We define the function a 1 : Ω → R + by a 1 (x) = 0.5[1 + k sin(2πx/20)], where k ∈ [0, 1] is a number that represents the magnitude of spatial heterogeneity. The case k = 0 corresponds to a 1 (x) ≡ 0.5. Now, we fix the values of the diffusion rates as d S = d I = d R = 1 and plot in Fig. 1 the solutions to system (1.2) when k = 0, that is, when the incidence rate F (x, S, I) = βSI/(1 + a 1 I) is independent of x. In this case, the densities of susceptible, infected and recovered populations converge to positive constants.

If we increase the heterogeneity by choosing k = 0.95 while keeping all other parameter values as before, we obtain the solutions plotted in Fig. 2 . Here, the values of S(x, t), I(x, t) and R(x, t) tend to a non-homogeneous endemic equilibrium as t → ∞.

Lastly, we plot the spatial distribution of susceptible, infected and recovered populations at t = 1000 for k = 0, 0.5, 0.8, 0.95 in Fig. 3 . We can see that greater heterogeneity in the saturation parameter results in higher values of infected and recovered populations. 

We will consider now the effects of using different diffusion rates for system (1.2) while keeping all other parameter values fixed. Here, we assume that a 1 (x) = 0.5[1 + 0.95 sin(2πx/20)] and other parameters are taken as in (6.1). Suppose first that the diffusion rate of infected population is fixed as d I = 1, while the diffusion rates of susceptible and recovered populations are both equal to some constant D. Figure 4a depicts the distribution of S, I and R at t = 1000 for several values of D. The figure shows that reducing the value of the diffusion rates d S = d R = D contributes to reducing the infected population levels.

Next, we make some simulations assuming that the diffusion rates d S , d I and d R are all equal to D. In this case, the spatial distribution of solutions is shown in Fig. 4b . Our simulations show that reducing the diffusion rate in all three subpopulations results in greater heterogeneity for the distribution of the infected; however, the density of infected population actually increases at certain locations for lower values of D. Hence, restricting the motility for all individuals does not necessarily reduce the total number of infections. 

We will now simulate the bifurcation dynamics of our model when the parameter γ 1 is varied. Here, we fix the parameter values Λ = 10, μ = 0.001, α = 0.01, γ 2 = 0.07, δ = 0.02, d S = d I = d R = 1,

F (x, S, I) = βSI 1 + a 1 (x)I , a 1 (x) = 0.5 1 + 0.95 sin 2πx 20 (6.3) and let γ 1 vary from 0 to 0.8.

Using the same initial conditions as before, we obtain the distribution of susceptible, infected and recovered individuals shown in Fig. 5 for t = 4000. We can see that for γ 1 ∈ [0, 0.6], the solution tends to a positive endemic steady state, which corresponds to the case when R 0 > 1. As we increase the value of γ 1 , the infected population decreases. The basic reproduction number crosses unity at a certain value γ * 1 ∈ (0.6, 0.8), and there is no endemic equilibrium for γ 1 > γ * 1 . Figure 5 shows that the solutions of the model tend to the disease-free steady state (E 0 = (10 4 , 0, 0)) when γ 1 = 0.8.

We proposed a diffusive SIRS epidemic model that generalizes several previously studied models in the literature, such as [9] [10] [11] 18] . In particular, we extended the results obtained by Yang et al. in [18] by considering a model with three different diffusion rates, d S (x), d I (x) and d R (x), which may vary depending on the location of individuals. Moreover, our model assumes that the incidence function takes the general form F (x, S, I), which includes many different types of nonlinear or non-monotone functions and allows for spatial heterogeneity.

We established the threshold dynamics for model (1.2) with respect to the basic reproduction number: the disease-free steady state is globally asymptotically stable when R 0 < 1, while the disease persists and an endemic equilibrium appears when R 0 > 1.

Furthermore, we performed a bifurcation analysis for model (1.2) in the case when the diffusion coefficient of infected population and the rate of transfer from infected to susceptible classes are constant. We used the theory by Crandall and Rabinowitz [35, 37, 38] to determine the type of bifurcation that occurs when the basic reproduction number crosses unity. From Theorems 5.2 and 5.3, we can infer that the model always undergoes a forward bifurcation at R 0 = 1, and the existence of endemic equilibria is completely ruled out when R 0 is less than one.

We also carried out some simulations to illustrate the effects of heterogeneity on the dynamics of model (1.2). We used, as an example, the incidence rate of the form F (x, S, I) = βSI/[1 + a 1 (x)I], which does not meet the assumptions considered in [18] but is included by our general model. Our examples show that a larger variation in the incidence rate may increase the number of infections in the population. 

Population biology of infectious diseases: part I

Nilam: Stability of a time delayed SIR epidemic model along with nonlinear incidence rate and Holling type-II treatment rate

A delayed epidemic model with pulse vaccination

Dynamics of an SIR model with nonlinear incidence and treatment rate

Dynamics of an SEIR epidemic model with nonlinear incidence and treatment rates

The impact of media on the control of infectious diseases

Stability analysis of an SIR epidemic model with specific nonliner incidence rate

Bifurcations of an epidemic model with non-monotonic incidence rate of saturated mass action

Threshold dynamics of an SIRS model with nonlinear incidence rate and transfer from infectious to susceptible

Global asymptotic stability of a generalized SIRS epidemic model with transfer from infectious to susceptible

Global stability for SIRS epidemic models with general incidence rate and transfer from infectious to susceptible

Spatiotemporal transmission dynamics for influenza disease in a heterogenous environment

Analysis of a reaction-diffusion cholera epidemic model in a spatially heterogeneous environment

Spatial-temporal basic reproduction number and dynamics for a dengue disease diffusion model

Threshold dynamics of an age-space structured brucellosis disease model with Neumann boundary condition

Diffusion-reaction compartmental models formulated in a continuum mechanics framework: application to COVID-19, mathematical analysis, and numerical study

Threshold dynamics of a diffusive SIRI model with nonlinear incidence rate

Dynamical analysis of a diffusive SIRS model with general incidence rate

Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems

A reaction-diffusion within-host HIV model with cell-to-cell transmission

Analysis of a diffusive virus infection model with humoral immunity, cell-to-cell transmission and nonlinear incidence

Abstract functional-differential equations and reaction-diffusion systems

Divergence theorems and the supersphere

Partial Differential Equations of Mathematical Physics and Integral Equations

Nonlinear Elliptic Equations

Theory and Applications of Partial Functional Differential Equations

Asymptotic Behavior of Dissipative Systems

On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

A nonlocal and time-delayed reaction-diffusion model of dengue transmission

Basic reproduction numbers for reaction-diffusion epidemic models

Periodic-Parabolic Boundary Value Problems and Positivity

Robust persistence for semidynamical systems

Global attractors and steady states for uniformly persistent dynamical systems

Bifurcation from simple eigenvalues

Persistence and bifurcation of degenerate solutions

On global bifurcation for quasilinear elliptic systems on bounded domains

ZAMP Dynamics of a reaction-diffusion SIRS model with general

97119 Mérida Mexico e-mail: angel.cervantes@correo.uady.mx

Eric Avila-Vales e-mail: avila@correo.uady.mx (Received

This work was supported by Sistema Nacional de Investigadores (15284) and Conacyt-Becas [Consejo Nacional de Ciencia y Tecnología (Grant No. 777140)].Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.